Local and Global Order in Dense Packings of Semi-Flexible Polymers of Hard Spheres

The local and global order in dense packings of linear, semi-flexible polymers of tangent hard spheres are studied by employing extensive Monte Carlo simulations at increasing volume fractions. The chain stiffness is controlled by a tunable harmonic potential for the bending angle, whose intensity dictates the rigidity of the polymer backbone as a function of the bending constant and equilibrium angle. The studied angles range between acute and obtuse ones, reaching the limit of rod-like polymers. We analyze how the packing density and chain stiffness affect the chains’ ability to self-organize at the local and global levels. The former corresponds to crystallinity, as quantified by the Characteristic Crystallographic Element (CCE) norm descriptor, while the latter is computed through the scalar orientational order parameter. In all cases, we identify the critical volume fraction for the phase transition and gauge the established crystal morphologies, developing a complete phase diagram as a function of packing density and equilibrium bending angle. A plethora of structures are obtained, ranging between random hexagonal closed packed morphologies of mixed character and almost perfect face centered cubic (FCC) and hexagonal close-packed (HCP) crystals at the level of monomers, and nematic mesophases, with prolate and oblate mesogens at the level of chains. For rod-like chains, a delay is observed between the establishment of the long-range nematic order and crystallization as a function of the packing density, while for right-angle chains, both transitions are synchronized. A comparison is also provided against the analogous packings of monomeric and fully flexible chains of hard spheres.


Introduction
Over the last few decades, developments in the synthesis of novel polymers and the fabrication of polymer-based materials have turned them into key components of our daily lives. The research is ever-growing in the pursuit of polymer-based materials with enhanced properties as chain connectivity endows macromolecules with unique properties compared to monoatomic systems [1,2]. Describing polymer conformations and configurations statistically and understanding how these are connected to macroscopic properties are thus of paramount importance in technology and industry [3,4].
In parallel, many aspects of the phase behavior and self-organization of general atomic and particulate systems remain unknown or poorly understood. The analysis becomes a great deal harder when macromolecular systems are tackled. Advances in experimental, theoretical, and simulation methods continually enrich our fundamental knowledge of the phenomenon in a wide range of physical systems [5][6][7][8][9][10][11]. For example, theoretical models have been developed to predict the effect of the processing conditions on the phase transition in order to analyze the final morphologies [12,13]. Advances in the synthesis and characterization of colloidal and granular systems have provided significant insights on the phase behavior of monomeric and polymeric systems, particularly given their simplicity and large size compared to traditional polymers [14][15][16][17][18][19][20][21][22]. Simulations can further aid in the research studies of phase transition given the persistent advances at the level of hardware and software, with traditional approaches at the molecular level being based on the Molecular Dynamics (MD) or Monte Carlo (MC) algorithms [23][24][25][26][27].
The hard-sphere model is frequently used in the study of crystallization due to its simplicity and athermal nature, despite the obvious disadvantage of lacking chemical detail. Through molecular dynamics (MD) simulations, early and pioneering works have reported the entropy-driven crystallization of systems of hard spheres [12,28] once the melting point ϕ M monomers = 0.545 is reached [29] and given sufficient time for the observation of the phenomenon [30,31]. For individual hard spheres, the crystal expected to be obtained in experiments or simulations should be the face centered cubic (FCC), as it is slightly more thermodynamically stable than the hexagonal close-packed (HCP) one [32][33][34]. However, extensive experimental [35][36][37][38][39][40], analytical and simulation [30,[41][42][43][44][45] studies have clearly demonstrated the difficulty of obtaining a perfect crystal. Thus, random hexagonal close packing (RHCP) is almost always observed due to the small free energy of stacking faults [46], while the more stable FCC structure appears after a very slow transition from the RHCP morphology [47][48][49][50]. Therefore, crystal perfection in the form of pure FCC crystals is rarely encountered in computer simulations under constant volume and starting from a predominantly amorphous initial configuration [50][51][52].
The mechanism of the crystallization of hard colloidal polymers [52][53][54][55][56][57][58][59][60] differs substantially from that of traditional polymers, the latter as revealed in x-ray scattering studies of short alkane chains [61][62][63], and in MD [64][65][66][67] or MC [54,68] simulations, including simulations with highly precise, all-atom interaction potentials [69][70][71][72][73][74][75][76][77][78]. Just as for monomeric packings, off-lattice MC simulations have shown that random packings of fully flexible linear chains of tangent hard spheres are able to crystallize at high volume fractions through an entropy-driven mechanism, very similar to the one of monomeric analogs [52,[55][56][57]79]. Extremely long simulations [52] have further shown the FCC polymorph to be the stable one for fully flexible chains of tangent hard spheres. In shorter simulations, an RHCP phase with a unique stacking direction [55][56][57] appears most often. In spite of the similarities to monomeric counterparts, athermal polymer crystallization shows unique characteristics: the critical volume fraction for the phase transition (melting point) for the fully flexible (freely-jointed, FJ) chains of tangent hard spheres, ϕ M chains (FJ) (>0. 56), is higher than the melting point of monomers. Monomeric systems present diverse crystal morphologies, ranging between single FCC or HCP structures and close-stacked packings with random directions where twin defects, associated with the formation of fivefold (FIV) symmetry structures, kinetically frustrate crystallization [80,81]. However, both the FCC and RHCP crystals of tangent hard-sphere chains are usually free of twin defects due to a conformational entropic barrier [82]. It has been established that the phase transition and the ordered morphologies of hard-sphere chains are affected insignificantly by chain length, but are sensitive to factors such as gaps in bond lengths [58,83], the presence of surfaces/confinement [59,60,84] and chain stiffness [85][86][87][88].
For semi-flexible polymers [89][90][91][92], the rigidity of the chains is generally implemented in molecular simulations through bending [93][94][95][96][97][98][99][100] and/or torsional [99] potentials. A common topic of study about the phase behavior of semi-flexible polymers is the orientational (nematic) ordering of the chains. According to Onsager, thin and infinitely long rigid rods undergo a transition from the isotropic to the nematic phase, which is driven by entropy, resulting from the competition between translational and orientational ordering [101]. Building on this, Bolhuis and Frenkel mapped out the complete phase diagram of hard spherocylinders as a function of shape anisotropy [102]. The long-range, nematic order in solutions of semi-flexible polymers depends on the size ratio of the chain [103][104][105]. For tangent hard-sphere chains in the rod limit, the transition occurs at Polymers 2023, 15, 551 3 of 23 high concentrations, still lower than the solidification of the system [106,107]. Theories have been developed [108][109][110][111] and studies have been conducted [112][113][114][115][116][117] under various conditions on the isotropic-to-nematic transition of semi-flexible chains, as the effect on the nematic order and the phase separation of different types of blends through experimental studies [118][119][120] and computer simulations [117,[121][122][123][124]. Nguyen et al. studied the effect of chain stiffness and temperature on the competition between crystallization and glass formation for unentangled semi-flexible polymer melts [85,86]. Shakirov and Paul studied the crystallization in melts of short, semi-flexible hard-sphere chains employing Monte Carlo simulations [87,88]. In their works, rod-like polymer melts undergo a first-order transition at a density lower than a melt of hard-sphere monomers or flexible hard-sphere chains [88]. During crystallization, the orientational ordering (nematic) transition is accompanied by a 2-D translational ordering on the plane perpendicular to the chains. These orderings induce the formation of hexagonal crystal planes [89].
Although rod-like chains have been extensively researched in the literature, the corresponding body of work on semi-flexible chains with different bending angles is very limited. Recent MD simulations, studying the jamming and solidification of semi-flexible polymers employing the freely-rotating model varying the fixed bending angle, have revealed the different mechanisms in solidification depending on the equilibrium bending angle, θ 0 [125,126].
The present manuscript analyzes the phase behavior of linear, semi-flexible chains of tangent hard spheres, as a function of the equilibrium bending angle and packing density. This is achieved through extensive Monte Carlo (MC) simulations using the Simu-D simulator-descriptor [127], built around polymer-oriented algorithms and particularly chain-connectivity-altering moves (CCAMs) [128,129]. Motivated by previous studies on the crystallization of linear, fully flexible chains of hard spheres [55][56][57][58][59]130,131], we extend the research through the inclusion of chain stiffness via a bending potential that depends on a bending constant, k θ , and an equilibrium bending angle, θ 0 . The global order, at the level of the chains, is quantified through the orientational (nematic) order parameter, q. The local order, at the level of the monomers, is gauged through the Characteristic Crystallographic Element (CCE) norm descriptor [132,133]. Compiling the results of all the simulated systems, we present a phase diagram that reflects the combined effect of packing density and chain stiffness (equilibrium bending angle) on the local and global order of the semi-flexible systems.
This article is organized as follows: in Section 2 we describe the molecular model and the simulated systems, and we briefly explain the descriptors to gauge the long-range order and local structure; Section 3 presents the results on the phase behavior of the simulated semi-flexible systems; finally, we summarize the main conclusions and discuss the current efforts in Section 4.

Molecular Model and Systems studied
Polymers are modeled as linear chains of identical hard spherical monomers with a collision diameter σ, which is considered the characteristic length. The systems are composed of N at monomers distributed in N ch chains of average chain length of N av (in monomers). All of the systems consist of N ch = 100 chains of average chain length N av = 12, resulting in N at = 1200. The chains present dispersity in their lengths following a uniform distribution in the range N ∈ [6,18] due to the presence of chain-connectivity-altering MC moves (see below).
The non-overlapping condition of the hard monomers is adopted by employing the Hard Sphere (HS) potential to describe all of the non-bonded interactions between the monomers. According to the HS potential, the pair-wise energy, U HS r ij , is determined by: where r ij is the distance between the centers of the monomers i and j. Periodic boundary conditions are applied in all dimensions, corresponding to a bulk, unconstraint system. For non-overlapping objects' packing density (volume fraction), ϕ, is defined as the total volume occupied by the monomers over the volume of the simulation cell. This is given by: where V mon is the total volume occupied by the monomers and V cell is the volume of the cubic simulation cell. Figure 1 shows a sketch of the bond geometry encountered in linear chains, consisting of bond lengths, l, bending angles, θ, and torsion angles, φ. Bond lengths can fluctuate uniformly in the interval l ∈ [σ, σ + dl], where dl is the maximum gap between successive monomers of the chain. For all of the systems simulated here, dl = 6.5 × 10 −4 (in units of σ), practically enforcing tangency between bonded monomers [58]. Chain stiffness is introduced through a potential controlling the bending angle, θ, which is the angle formed by two successive bond vectors ( Figure 1). The bending potential is given by: where k θ is the bending constant and θ 0 is the equilibrium bending angle. For fixed bond lengths, setting k θ = 0 allows the simulation of freely-jointed chains while k θ → ∞ corresponds to freely-rotating chains. The equilibrium bending angle can vary from fully extended (θ 0 = 0 • ) to fully compact (θ 0 = 120 • ) configurations of triplets. In this work, we study five specific equilibrium bending angles: 0 • , 60 • , 90 • , 108 • , and 120 • , as they are represented, schematically, in Figure 1. Torsion angles, φ, are allowed to fluctuate freely.
monomers. According to the HS potential, the pair-wise energy, , is determined by: where rij is the distance between the centers of the monomers i and j. Periodic boundary conditions are applied in all dimensions, corresponding to a bulk, unconstraint system. For non-overlapping objects' packing density (volume fraction), , is defined as the total volume occupied by the monomers over the volume of the simulation cell. This is given by: where Vmon is the total volume occupied by the monomers and Vcell is the volume of the cubic simulation cell. Figure 1 shows a sketch of the bond geometry encountered in linear chains, consisting of bond lengths, l, bending angles, θ, and torsion angles, . Bond lengths can fluctuate uniformly in the interval ∈ , + , where dl is the maximum gap between successive monomers of the chain. For all of the systems simulated here, = 6.5 × 10 (in units of σ), practically enforcing tangency between bonded monomers [58]. Chain stiffness is introduced through a potential controlling the bending angle, θ, which is the angle formed by two successive bond vectors ( Figure 1). The bending potential is given by: where is the bending constant and is the equilibrium bending angle. For fixed bond lengths, setting = 0 allows the simulation of freely-jointed chains while → ∞ corresponds to freely-rotating chains. The equilibrium bending angle can vary from fully extended ( = 0°) to fully compact ( = 120°) configurations of triplets. In this work, we study five specific equilibrium bending angles: 0°, 60°, 90°, 108°, and 120°, as they are represented, schematically, in Figure 1. Torsion angles, , are allowed to fluctuate freely. The bending constant is set equal to ⁄ = 9 rad −2 (for simplicity in the continuation we will report it as = 9) in all of the simulated systems, where and T are the Boltzmann constant and temperature, respectively. The specific choice of the harmonic constant was made considering the computational expediency and the fact that this range of constraints leads to orientational ordering for rod-like chains, as demonstrated in [87]. This degree of stiffness limits the deviations from the equilibrium bending angle to, approximately, ±20°, as shown in Figure 2 where, for comparison, we also present the corresponding curves as obtained by increasing the bending constant to = 50. The bending constant is set equal to k θ /k B T = 9 rad −2 (for simplicity in the continuation we will report it as k θ = 9) in all of the simulated systems, where k B and T are the Boltzmann constant and temperature, respectively. The specific choice of the harmonic constant was made considering the computational expediency and the fact that this range of constraints leads to orientational ordering for rod-like chains, as demonstrated in [87]. This degree of stiffness limits the deviations from the equilibrium bending angle to, approximately, ±20 • , as shown in Figure 2 where, for comparison, we also present the corresponding curves as obtained by increasing the bending constant to k θ = 50.

Simulation Algorithm
The simulations are carried out by means of the Simu-D simulator-descriptor suite [127]. The simulations are performed under constant volume and temperature ( = 1 ⁄ ), employing Periodic Boundary Conditions (PBCs). The equilibrium simulations are performed mainly at high packing densities in the range = 0.54, 0.60 . Additionally, simulations are also carried out at significantly lower volume fractions to study the longrange orientational order, particularly for rod-like chains ( = 0°). The choice of volume fractions is driven by the fact that, in past studies, we showed that crystallization for systems of freely-jointed chains of tangent hard spheres occurs once a packing density of ≈ 0.58 is reached while at = 0.56 polymer packings remain amorphous [55][56][57]. For comparison, monomeric analogs crystallize at = 0.545. In Figure S1 (see Supplementary Materials), computer-generated configurations are visualized corresponding to dilute conditions ( = 0.10) for all of the semi-flexible systems studied here.
The systems are generated and equilibrated by employing a mix of localized algorithms and chain-connectivity-altering moves (CCAMs), as implemented in the Simu-D simulator [127]. The MC mix is composed of the following moves: reptation (10%); rotation (10%); flip (34,8%); intermolecular reptation (25%); end-segment re-arrangement (or CCB as in Refs. [134,135]) (20%); simplified end-bridging (sEB) (0.1%); and simplified intramolecular end-bridging (sIEB) (0.1%). The attempt probabilities of each MC move are reported in parentheses and remain unaltered for all of the systems. Due to the high packing densities, all of the local moves are executed according to a configurational bias (CB) pattern (see more technical explanations in Refs. [127,128]), with the number of trials per attempted move depending on the packing fraction, according to Table 1 of Ref. [128].
The initial configurations are generated at very dilute conditions ( = 0.001) as fully flexible systems. Under these conditions, the desired equilibrium bending angle, , and bending constant, are activated, followed by rapid constant-volume equilibration. The dilute systems of the semi-flexible chains are then shrunk until the desired volume fraction is reached through isotropic volume compressions attempted at regular intervals [128]. The equilibrium, constant-volume simulations are conducted between a minimum of 5 × 10 and a maximum of 8 × 10 MC steps, depending on the system. The systems' configurations ("frames" or "snapshots") are recorded every 10 MC steps.

Simulation Algorithm
The simulations are carried out by means of the Simu-D simulator-descriptor suite [127]. The simulations are performed under constant volume and temperature (T = 1/k B ), employing Periodic Boundary Conditions (PBCs). The equilibrium simulations are performed mainly at high packing densities in the range ϕ = [0.54, 0.60]. Additionally, simulations are also carried out at significantly lower volume fractions to study the long-range orientational order, particularly for rod-like chains (θ 0 = 0 • ). The choice of volume fractions is driven by the fact that, in past studies, we showed that crystallization for systems of freely-jointed chains of tangent hard spheres occurs once a packing density of ϕ ≈ 0.58 is reached while at ϕ = 0.56 polymer packings remain amorphous [55][56][57]. For comparison, monomeric analogs crystallize at ϕ M monomers = 0.545. In Figure S1 (see Supplementary Materials), computer-generated configurations are visualized corresponding to dilute conditions (ϕ = 0.10) for all of the semi-flexible systems studied here.
The systems are generated and equilibrated by employing a mix of localized algorithms and chain-connectivity-altering moves (CCAMs), as implemented in the Simu-D simulator [127]. The MC mix is composed of the following moves: reptation (10%); rotation (10%); flip (34,8%); intermolecular reptation (25%); end-segment re-arrangement (or CCB as in Refs. [134,135]) (20%); simplified end-bridging (sEB) (0.1%); and simplified intramolecular end-bridging (sIEB) (0.1%). The attempt probabilities of each MC move are reported in parentheses and remain unaltered for all of the systems. Due to the high packing densities, all of the local moves are executed according to a configurational bias (CB) pattern (see more technical explanations in Refs. [127,128]), with the number of trials per attempted move depending on the packing fraction, according to Table 1 of Ref. [128].
The initial configurations are generated at very dilute conditions (ϕ = 0.001) as fully flexible systems. Under these conditions, the desired equilibrium bending angle, θ 0 , and bending constant, k θ are activated, followed by rapid constant-volume equilibration. The dilute systems of the semi-flexible chains are then shrunk until the desired volume fraction is reached through isotropic volume compressions attempted at regular intervals [128]. The equilibrium, constant-volume simulations are conducted between a minimum of 5 × 10 11 and a maximum of 8 × 10 11 MC steps, depending on the system. The systems' configurations ("frames" or "snapshots") are recorded every 10 7 MC steps.

Global (Long-Range) Orientational Order
As we simulate semi-flexible chains, it is important to quantify the long-range orientational order of the chains [101]. In the nematic order, the chains exhibit a long-range orientation with their long axis aligned along a preferent direction, denominated nematic director n. The isotropic-to-nematic transition is influenced by parameters such as packing density [106,107] and ratio of chain size [103][104][105].
The chain orientational order is determined through averages of a second-order invariant [136]. The orientation of each molecule is defined by the unit vector u along the long axis of the molecule. The long axis of each chain is determined by the inertia tensor I, a second-order tensor calculated according to Equation (4), where N(j) is the number of monomers of the chain j, m i is the mass of the monomer (considered here as unity), x i is the coordinate vector of the monomer, x cm is the coordinate vector of the centre of mass of the chain and δ is the second order isotropic tensor. The long axis of each chain corresponds to the eigenvector, v 3 , of the lowest eigenvalue of the inertia tensor I so that the unit vector u is calculated normalizing the eigenvector v 3 .
As an average of all the molecular orientations, we define the second-order tensor Q according to: The Q tensor is: (1) symmetric, and (2) traceless (TrQ = 0). These two properties reduce the independent components of the second-order tensor from 9 to 5.
The tensor Q can be calculated for the different ideal phases (isotropic phase; prolate mesogen, nematic mesophase; and oblate mesogen, nematic mesophase; for their schemes see for example Figure 4 in [136]). For the isotropic phase, denoted here as "ISO", the molecules can present any random orientation with equal probability, so no orientation is preferred. Thus, all components of the resulting matrix Q ISO are zero. In the perfectly aligned prolate mesogen, the long axis of every molecule is aligned along a preferred orientation leading to the nematic mesophase. In a coordinate system in which the preferred orientation of the system is along the x-axis, the matrix representation of Q for the perfect prolate mesogen, nematic mesophase, denoted here as "PRO", takes the form: In the oblate mesogen, nematic mesophase, it is the short axis of the molecules that tends to be aligned with the nematic director n. In the same coordinate system, the matrix representation of Q for the oblate mesogen, nematic mesophase, denoted here as "OBL", takes the form: In order to identify the similarities with these ideal cases, we diagonalize Q and use its normalized eigenvalues as unit vectors of a new Cartesian coordinate system where Q' is diagonal and its diagonal elements are its eigenvalues (λ 1 , λ 2 , and λ 3 ), ordered by decreasing the absolute value. The eigenvector of the largest absolute eigenvalue λ 1 , in the case of a nematic mesophase (denoted here as "NEM"), would correspond to the preferred orientation of the system or nematic director n.
A scalar order parameter q is obtained by comparing the diagonal tensor Q' with the corresponding tensor Q PRO of a perfect prolate mesogen, nematic mesophase. This scalar q represents the degree of alignment between the chains. In an isotropic system, where Q is null, q → 0. For the perfectly aligned prolate mesogen, nematic mesophase system, q = 1, while q = −1/2 for the ideal oblate mesogen, nematic mesophase. The orientational order parameter, thus, takes values in the range −1/2 ≤ q ≤ 1.
Accordingly, the long-range orientational order is characterized by the scalar orientational order parameter q and the nematic director n.

Local Order: Characteristic Crystallographic Element Norm
The metric used in this work to identify the disorder-order transition at the local level and quantify the degree of crystallinity of the simulated systems is the Characteristic Crystallographic Element (CCE) norm, which is explained in detail in Refs. [132,133]. The CCE norm is integrated into the descriptor part of the Simu-D suite, which is also used for the MC simulations [127].
The CCE norm descriptor gauges the crystallinity of an atomic or particulate system in two (2D) or three (3D) dimensions by comparing the local environment around each site with a set of ideal reference crystals under the main concept that each ideal crystal is uniquely identified by a set of symmetry operations [137][138][139][140]. Therefore, for every site i of a system, the CCE norm descriptor identifies the nearest neighbors and quantifies the orientational and radial deviations of the "real" local environment with respect to the "ideal" environment of each reference crystal. This comparison provides, for the given site i, a CCE norm value with respect to each X reference crystal, ε X i . The closer the X-CCE norm is to zero, the higher the similarity of the local environment to the respective reference crystal X. Site i is identified as an X-type crystal when the calculated CCE norm is lower than a critical threshold, ε X i ≤ ε thres . In past studies, an empirical threshold of 0.245 was determined for packings of non-overlapping spheres [56,132,133] in the bulk [55][56][57][58]80,81] and under confinement [59,60]. Due to the strict concept behind the descriptor, the CCE norm is highly discriminatory, so the value of a site cannot be simultaneously very low for two different reference crystals. The current version of the CCE norm descriptor can identify similarity with respect to the following crystals: hexagonal close packed (HCP); face centered cubic (FCC); body centered cubic (BCC); and hexagonal (HEX) for 3D systems, and triangular (TRI); square (SQU); and honeycomb (HON) for 2D systems. The fivefold (FIV) and pentagonal (PEN) local symmetries can also be identified in three and two dimensions, respectively. Sites that cannot be assigned to any of the previous reference ideal environments are labeled as amorphous (AMO) or, more precisely, as "unidentified".
The process explained before is repeated over all sites of the system for each reference crystal. Once the CCE norm has been evaluated for every site and reference crystal, an order parameter for each reference crystal X, S X (S X ∈ [0, 1]), can be calculated for the snapshot as: where P ε X is the probability function of CCE-based norms of the reference crystal X. Additionally, a degree of crystallinity (or total crystallinity), τ c , can be calculated as the sum of the order parameters of all the reference crystals. For a bulk, 3D system, the degree of crystallinity is calculated as: In the present work, no appreciable population of sites with HEX or BCC similarity was detected (S HEX , S BCC → 0). Thus, in the continuation, we consider only the HCP and FCC crystals, as well as fivefold local symmetry (FIV).

Results
The phase behavior of polymer chains under the effect of the chain stiffness is the main focus of this work, emphasizing the long-range, orientational (nematic) ordering and the local structure.
Prior to the simulation data being analyzed in a post-processing step, a preliminary visual inspection of the initial and resulting system configurations is performed. In the case of rod-like chains (θ 0 = 0 • ), the visual inspection of the initial configurations, as presented in Figure 3, suggests a transition from an isotropic (ϕ = 0.10 and 0.20) to a nematic phase (ϕ = 0.30 and 0.50) with increasing packing density, i.e., the rod-like chains align in the system along the common nematic director. This transition occurs at a packing density that is significantly lower than for solidification, in perfect qualitative agreement with past independent works on similar systems [101,106,107]. Figure S2 (see Supplementary Materials) shows snapshots at the end of the simulation for all systems at ϕ= 0.60. where ( ) is the probability function of CCE-based norms of the reference crystal X. Additionally, a degree of crystallinity (or total crystallinity), , can be calculated as the sum of the order parameters of all the reference crystals. For a bulk, 3D system, the degree of crystallinity is calculated as: In the present work, no appreciable population of sites with HEX or BCC similarity was detected (S HEX , S BCC → 0). Thus, in the continuation, we consider only the HCP and FCC crystals, as well as fivefold local symmetry (FIV).

Results
The phase behavior of polymer chains under the effect of the chain stiffness is the main focus of this work, emphasizing the long-range, orientational (nematic) ordering and the local structure.
Prior to the simulation data being analyzed in a post-processing step, a preliminary visual inspection of the initial and resulting system configurations is performed. In the case of rod-like chains ( = 0°), the visual inspection of the initial configurations, as presented in Figure 3, suggests a transition from an isotropic ( = 0.10 and 0.20) to a nematic phase ( = 0.30 and 0.50) with increasing packing density, i.e., the rod-like chains align in the system along the common nematic director. This transition occurs at a packing density that is significantly lower than for solidification, in perfect qualitative agreement with past independent works on similar systems [101,106,107]. Figure S2 (see Supplementary Materials) shows snapshots at the end of the simulation for all systems at = 0.60.  [141]. Individual panels are also available as stand-alone, interactive, 3-D images (in Supplementary Materials).

Global Orientational Order
Following the preliminary visual inspection of the initial configurations for the rodlike chains ( = 0°), the long-range, orientational (nematic) ordering is analyzed through the second-order tensor Q and the scalar order parameter q, as explained previously. As an example, Figure S3 (see Supplementary Materials) shows the evolution of the exponential moving average of q as a function of the MC steps for rod-like chains at different packing densities. The dependence of the nematic order parameter q on the packing density, , is presented in Figure 4 for all of the semi-flexible systems under study ( = 0, 60, 90, 108, and 120°), averaged over all frames of the equilibrated part of the simulation trajectory. In line with Figure 3, the rod-like chains of average length Nav = 12 show a well-  [141]. Individual panels are also available as stand-alone, interactive, 3-D images (in Supplementary Materials).

Global Orientational Order
Following the preliminary visual inspection of the initial configurations for the rod-like chains (θ 0 = 0 • ), the long-range, orientational (nematic) ordering is analyzed through the second-order tensor Q and the scalar order parameter q, as explained previously. As an example, Figure S3 (see Supplementary Materials) shows the evolution of the exponential moving average of q as a function of the MC steps for rod-like chains at different packing densities. The dependence of the nematic order parameter q on the packing density, ϕ, is presented in Figure 4 for all of the semi-flexible systems under study (θ 0 = 0, 60, 90, 108, and 120 • ), averaged over all frames of the equilibrated part of the simulation trajectory. In line with Figure 3, the rod-like chains of average length N av = 12 show a well-defined isotropic-nematic transition as the concentration increases. This transition takes place in the interval 0.15 ≤ ϕ ≤ 0.20 and very closely resembles the prediction of Onsager's theory. At higher volume fractions, the chains form a prolate mesogen, nematic mesophase (PRO), as seen by the monotonically increasing value of q. A perfect PRO phase (q → 1) is established at packing densities approximately equal to ϕ = 0.45. For the average length studied here (N av = 12), the transition from the isotropic to the nematic phase, ISO → PRO (i.e., nematic mesophase with prolate mesogens) for rod-like athermal chains occurs at significantly lower packing densities (ϕ nem ≈ 0.20) than the freezing (ϕ F monomers ≈ 0.495) and melting (ϕ M monomers ≈ 0.545) points of monomeric hard spheres.  In addition to this expected isotropic-to-nematic (ISO → PRO) transition exhibited by the rod-like chains, a nematic mesophase with oblate mesogens, OBL, (q < 0) appears for semi-flexible chains with = 90°. In contrast to the long-range orientational order of the rod-like chains, the isotropic-to-oblate, ISO → OBL, transition takes place at ≈ 0.57, which is higher than the melting point of HS.
For the remaining values of the equilibrium bending angles, each chain system remains in the isotropic phase in the whole concentration range. Still, a small peak is produced at very high packing densities that can be associated with chains primarily of = 60°, taking some preferred directions, corresponding to the values of the nematic order parameter of around q ≈ 0.15. This can be explained as the establishment of the FCC and HCP crystallites (see later discussion on local order) enforces specific bending angles which are compatible with these crystals. If such bending angles are not available through intra-chain arrangements, then the only other option is inter-chains ones, thus inducing partial alignment among the chains at a local level.
In the case of semi-flexible chains with = 90°, the left panel of Figure 5 hosts the evolution of the exponential running average of the orientational order parameter, q, as a function of the MC steps for the packing densities where a certain degree of nematic order In addition to this expected isotropic-to-nematic (ISO → PRO) transition exhibited by the rod-like chains, a nematic mesophase with oblate mesogens, OBL, (q < 0) appears for semi-flexible chains with θ 0 = 90 • . In contrast to the long-range orientational order of the rod-like chains, the isotropic-to-oblate, ISO → OBL, transition takes place at ϕ obl ≈ 0.57, which is higher than the melting point of HS.
For the remaining values of the equilibrium bending angles, each chain system remains in the isotropic phase in the whole concentration range. Still, a small peak is produced at very high packing densities that can be associated with chains primarily of θ 0 = 60 • , taking some preferred directions, corresponding to the values of the nematic order parameter of around q ≈ 0.15. This can be explained as the establishment of the FCC and HCP crystallites (see later discussion on local order) enforces specific bending angles which are compatible with these crystals. If such bending angles are not available through intrachain arrangements, then the only other option is inter-chains ones, thus inducing partial alignment among the chains at a local level.
In the case of semi-flexible chains with θ 0 = 90 • , the left panel of Figure 5 hosts the evolution of the exponential running average of the orientational order parameter, q, as a function of the MC steps for the packing densities where a certain degree of nematic order was observed in Figure 4. A tendency for the formation of an oblate mesogen, nematic mesophase (OBL) is observed after enough simulation time. The degree of ordering of the oblate mesogens increases with the volume fraction, as indicated by q approaching the ideal value −1/2. Compared to the crystallization at the level of monomers (see below), the ISO → OBL transition for θ 0 = 90 • takes places with an equal or even slower rate. We should note here that the "rate" used in the present manuscript has no physical meaning and it corresponds to the number of MC steps required to observe the phase transition. Thus, the nematic phase for the semi-flexible chains with θ 0 = 90 • practically coincides with crystallization, i.e., self-organization at the local and global levels appears synchronized. the ISO → OBL transition for = 90° takes places with an equal or even slower rate. We should note here that the "rate" used in the present manuscript has no physical meaning and it corresponds to the number of MC steps required to observe the phase transition. Thus, the nematic phase for the semi-flexible chains with = 90° practically coincides with crystallization, i.e., self-organization at the local and global levels appears synchronized. The right panels of Figure 5 show the semi-flexible chains with = 90° at = 0.59 at the end of the MC simulation, once the system has reached the OBL state. The top snapshot presents the unwrapped representation of the semi-flexible chains of the complete system. As in the previous snapshots, the monomers are color-coded according to the parent chain. For clarity, the bottom snapshot contains only 10 randomly selected chains of the same system. The semi-flexible chains tend to form flat layers, interrupted occasionally by right-angle jumps between planes, consistent with the employed constraint of = 90°. These flat chain configurations have, in general, a common behavior that further explains the OBL phase in Figure 4 and the left panel of Figure 5. Images created with VMD visualization software [141]. Individual panels are also available as stand-alone, interactive, 3-D images (in Supplementary Materials). The right panels of Figure 5 show the semi-flexible chains with θ 0 = 90 • at ϕ = 0.59 at the end of the MC simulation, once the system has reached the OBL state. The top snapshot presents the unwrapped representation of the semi-flexible chains of the complete system. As in the previous snapshots, the monomers are color-coded according to the parent chain. For clarity, the bottom snapshot contains only 10 randomly selected chains of the same system. The semi-flexible chains tend to form flat layers, interrupted occasionally by right-angle jumps between planes, consistent with the employed constraint of θ 0 = 90 • . These flat chain configurations have, in general, a common behavior that further explains the OBL phase in Figure 4 and the left panel of Figure 5.

Local Order: Crystal Nucleation and Growth
The local structure is gauged through the CCE norm descriptor [132,133]. For consistency with our past studies on the flexible chains of hard spheres, we employ a threshold of ε thres = 0.245 to label a site as X-type, where X is the reference crystal or local symmetry. In the present work, given that no appreciable population of BCC or HEX sites is detected in any of the simulated systems, X corresponds to HCP, FCC, and FIV. In the continuation, and throughout the manuscript, the corresponding colors to be used for the representation of the HCP, FCC, and FIV sites (in snapshots) and curves (in figures) are blue, red, and green, respectively. The results from the semi-flexible systems are also compared with the ones of fully flexible (freely-jointed) chains of tangent hard spheres simulated and analyzed through the Simu-D software under the same conditions of volume fraction and average chain length [55][56][57]130,131]. Figure 6 presents the CCE-based order parameters for HCP (S HCP ), FCC (S FCC ), and FIV (S FIV ), along with the degree of crystallinity, τ c , as a function of the MC steps at a packing density ϕ = 0.59 for all of the semi-flexible systems, including the fully flexible one as reference. This volume fraction is higher than the melting transition of the freely-jointed hard-sphere chains, as established in [52,[55][56][57]79]. An inspection of all of the panels shows clearly that the mechanism of the phase transition, here in the form of crystallization, is very similar and rather independent of the equilibrium bending angle. The initial packings are amorphous (AMO) as the fraction of sites with non-ordered local structure is vastly dominant: the percentage of close packed sites does not exceed 5% and, in most cases, the FIV population is commensurate or even exceeds the combined HCP + FCC one. This initial state of the AMO athermal chain packings agrees with the past works on dense monomeric and polymer assemblies [57,60,[80][81][82]. Out of all of the systems, the rod-like one is characterized by the lowest population of FIV sites at the beginning, as FIV appears to be incompatible with the perfect nematic ordering exhibited by the chains at this range of volume fractions.
Once crystallites of the HCP and FCC characters start growing, at almost identical rates, the fraction of the FIV-like sites decreases gradually until it practically disappears. The structural competition, observed here for all of the semi-flexible chain systems, between the FIV local symmetry and crystallization in the form of the HCP and FCC sites is in perfect match with identical observations in simulations of fully flexible [52,[55][56][57][58][59][60] and monomeric HS systems [80,81], including bead-spring chains under quenching [85,142]. Once the population of FIV sites is eliminated, a trend manifestly valid for all of the simulated packings, the relative fractions of the HCP and FCC crystals undergo a sharp variation, which is followed by the establishment of the final, stable ordered structures of the crystalline (CRY) character. The average degree of crystallinity in the final CRY phase ranges between 0.65 and 0.85, with the lowest and highest values corresponding to the fully flexible and rod-like chains, respectively. At ϕ = 0.59, as seen in Figure 6, all of the crystals contain appreciable fractions of FCC and HCP sites, and no perfection is registered. Accordingly, one expects that the formed morphologies correspond to fivefold-free but defect-ridden RHCP crystals. Figure 7 shows, for all of the systems whose phase transition is presented in Figure 6, the snapshots at the beginning (top panels) and the end (bottom panels) of the corresponding MC simulations. As already quantified by the data in Figure 6, all of the initial configurations show amorphous structures with a remarkable population of FIV-like sites. At the end of the MC simulations, the crystalline phase of every system shows a stable configuration of mixed HCP/FCC structures, being defect-ridden and fivefold-free, exactly as expected by the fractions of sites in Figure 6. In their majority, the semi-flexible chain packings form RHCP structures of alternating layers of unique FCC and HCP character with a single stacking direction, as in the freely-jointed systems. Additionally, the semiflexible systems can also form crystal structures with multiple, random stacking directions, as it is observed for the system with θ 0 = 90 • .  Figure 7 shows, for all of the systems whose phase transition is presented in Figure  6, the snapshots at the beginning (top panels) and the end (bottom panels) of the corresponding MC simulations. As already quantified by the data in Figure 6, all of the initial configurations show amorphous structures with a remarkable population of FIV-like sites. chain packings form RHCP structures of alternating layers of unique FCC and HCP character with a single stacking direction, as in the freely-jointed systems. Additionally, the semi-flexible systems can also form crystal structures with multiple, random stacking directions, as it is observed for the system with = 90°.
The formation of RHCP crystals of mixed HCP/FCC structures, with unique or multiple stacking directions, is also observed for all of the equilibrium bending angles and packing densities studied where crystallization takes place. Although almost all of the semi-flexible systems crystallize in RHCP crystals of mixed HCP/FCC layers, two important exceptions exist for the rod-like chains ( = 0°) at = 0.58 and 0.60. Figure 8 shows the CCE-based local order parameter, , and the total crystallinity, , as a function of the MC steps, while the corresponding snapshots at various simulation instances can be found in Figure 9. Both systems show very similar trends: after a very short initial period, characterized by the rapid increase in the HCP and FCC sites and the parallel reduction of the FIV population, the growth rates stop being the same and one type grows in favor of the other. At = 0.58, the resulting morphology is an almost perfect FCC crystal ( = 0.83, = 0.83), while the opposite occurs at = 0.60, where a (less) perfect HCP crystal emanates ( = 0.77, = 0.74). In the case of the HCP crystal at = 0.60, the resulting structure is not as stable as that of the FCC crystal, alternating between the almost perfect HCP crystal ( → 0) and an HCP crystal with a small population of FCC-like sites ( ≈ 0.1). For the latter case, when FCC impurities appear in the HCP crystal, these are produced on the border of the crystal with a regime of defects in the form of amorphous (AMO) sites. Figure S4 (see Supplementary Materials) hosts the CCE-based snapshots of the final configuration for the corresponding systems. While perfection in the form of an FCC crystal has been reported very recently from extremely long simulations on freely-jointed chains of hard spheres [52], it is the first time that an almost perfect HCP crystal, made of hard-sphere chains, in extended rod-like conformations, is observed. The easiness with which the crystal perfection is observed could Figure 7. Snapshots of the initial (top) and final (bottom) configuration of the MC simulations for the 100-chains N av = 12 systems at ϕ = 0.59 of (from left to right) freely-jointed chains and semiflexible chains (k θ = 9 ) with equilibrium bending angle θ 0 = 0, 60, 90, 108, and 120 • . Monomers are color-coded according to CCE norm. Blue, red, cyan, purple, and green correspond to HCP-, FCC-, BCC, HEX-, and FIV-like sites, respectively. Amorphous (AMO) sites are colored in yellow and are shown with reduced dimensions for visual clarity. Individual panels are also available as stand-alone, interactive, 3-D images (in Supplementary Materials).
The formation of RHCP crystals of mixed HCP/FCC structures, with unique or multiple stacking directions, is also observed for all of the equilibrium bending angles and packing densities studied where crystallization takes place.
Although almost all of the semi-flexible systems crystallize in RHCP crystals of mixed HCP/FCC layers, two important exceptions exist for the rod-like chains (θ 0 = 0 • ) at ϕ = 0.58 and 0.60. Figure 8 shows the CCE-based local order parameter, S X , and the total crystallinity, τ c , as a function of the MC steps, while the corresponding snapshots at various simulation instances can be found in Figure 9. Both systems show very similar trends: after a very short initial period, characterized by the rapid increase in the HCP and FCC sites and the parallel reduction of the FIV population, the growth rates stop being the same and one type grows in favor of the other. At ϕ = 0.58, the resulting morphology is an almost perfect FCC crystal (τ c = 0.83, S FCC = 0.83), while the opposite occurs at ϕ = 0.60, where a (less) perfect HCP crystal emanates (τ c = 0.77, S HCP = 0.74). In the case of the HCP crystal at ϕ = 0.60, the resulting structure is not as stable as that of the FCC crystal, alternating between the almost perfect HCP crystal ( S FCC → 0 ) and an HCP crystal with a small population of FCC-like sites (S FCC ≈ 0.1). For the latter case, when FCC impurities appear in the HCP crystal, these are produced on the border of the crystal with a regime of defects in the form of amorphous (AMO) sites. Figure S4 (see Supplementary Materials) hosts the CCE-based snapshots of the final configuration for the corresponding systems.
While perfection in the form of an FCC crystal has been reported very recently from extremely long simulations on freely-jointed chains of hard spheres [52], it is the first time that an almost perfect HCP crystal, made of hard-sphere chains, in extended rod-like conformations, is observed. The easiness with which the crystal perfection is observed could be related to the nematic ordering exhibited by the rod-like chains. This trend, particularly compared to the fully flexible model, will be explored in more detail in future studies. be related to the nematic ordering exhibited by the rod-like chains. This trend, particularly compared to the fully flexible model, will be explored in more detail in future studies.   be related to the nematic ordering exhibited by the rod-like chains. This trend, particularly compared to the fully flexible model, will be explored in more detail in future studies.

Local Order: Total Crystallinity
From the local order parameters S X , the total degree of crystallinity, τ c , can be gauged for all the simulated semi-flexible systems. As in the previous analysis, the present results are compared against the ones of fully flexible (freely-jointed) analogs. Figure S5 (see Supplementary Materials) presents the evolution of the degree of crystallinity as a function of the MC steps at different packing densities, ϕ.
Based on these results and utilizing the value of the total crystallinity as established in the final, stable part of the MC simulation, we can further extract the one-dimensional phase diagram of local order, quantified by the degree of crystallinity, as a function of the packing density and equilibrium bending angle ( Figure 10). It is worth keeping in mind here that the melting point for monomers is ϕ M monomers = 0.545, while freely-jointed chains of tangent hard spheres crystallize at higher volume fractions, ϕ M chains (>0.56). Some important conclusions can be drawn from the diagram. First, all of the semi-flexible chain packings eventually crystallize, independently of the equilibrium angle. Second, the equilibrium angle profoundly affects the onset of the phase transition of the systems. For example, the rod-like chains present semi-crystalline and crystalline phases at volume fractions significantly lower than the melting point for freely-jointed chains, even lower than the one for monomeric hard spheres, due to the effect of the nematic (global) ordering that precedes the local order (as will be demonstrated in the continuation). On the other hand, chains with θ 0 = 60 • crystallize later than all of the other systems, as at ϕ = 0.57 no crystal nucleation and growth is observed, even after 7 × 10 11 MC steps. Finally, based on the above, the melting transition shows the trend: From the local order parameters , the total degree of crystallinity, , can be gauged for all the simulated semi-flexible systems. As in the previous analysis, the present results are compared against the ones of fully flexible (freely-jointed) analogs. Figure S5 (see Supplementary Materials) presents the evolution of the degree of crystallinity as a function of the MC steps at different packing densities, .
Based on these results and utilizing the value of the total crystallinity as established in the final, stable part of the MC simulation, we can further extract the one-dimensional phase diagram of local order, quantified by the degree of crystallinity, as a function of the packing density and equilibrium bending angle ( Figure 10). It is worth keeping in mind here that the melting point for monomers is = 0.545, while freely-jointed chains of tangent hard spheres crystallize at higher volume fractions, (> 0.56). Some important conclusions can be drawn from the diagram. First, all of the semi-flexible chain packings eventually crystallize, independently of the equilibrium angle. Second, the equilibrium angle profoundly affects the onset of the phase transition of the systems. For example, the rod-like chains present semi-crystalline and crystalline phases at volume fractions significantly lower than the melting point for freely-jointed chains, even lower than the one for monomeric hard spheres, due to the effect of the nematic (global) ordering that precedes the local order (as will be demonstrated in the continuation). On the other hand, chains with = 60° crystallize later than all of the other systems, as at = 0.57 no crystal nucleation and growth is observed, even after 7 × 10 MC steps. Finally, based on the above, the melting transition shows the trend:  [29].
In general, with the sole exception of the rod-like chains, the presence of constraints, related to bond geometry in the form of bonds or bending angles, increases the melting transition with respect to the monomeric analogs. Nevertheless, it is quite surprising that the systems with obtuse bending angles ( = 120 and 108°) crystallize earlier, and the right angle ( = 90°) at a very similar volume fraction as the fully flexible chains, which In general, with the sole exception of the rod-like chains, the presence of constraints, related to bond geometry in the form of bonds or bending angles, increases the melting transition with respect to the monomeric analogs. Nevertheless, it is quite surprising that the systems with obtuse bending angles (θ 0 = 120 and 108 • ) crystallize earlier, and the right angle (θ 0 = 90 • ) at a very similar volume fraction as the fully flexible chains, which are completely free of bending constraints. Obtuse angles are favored at high volume fractions as they minimize the local volume compared to acute ones. However, one should further consider the fact that only specific geometric arrangements of polymer chains are compatible with the sites of ideal FCC and HCP crystals. For example, the FCC and HCP share the bending angles of 0, 60, 90, and 120 • , while the angles at 33.5 and 70.5 • exist exclusively on the HCP crystal [52]. The lack of bending angles compatible with these crystals would require the adjacent sites to be covered by monomers belonging to other chains, imposing specific arrangements at the intermolecular level.
The phase behavior at the local (CCE-based crystallinity, τ c ) and global (nematic orientational order parameter, q) levels as a function of the density for rod-like and rightangle chains is presented in Figure 11. As was explained previously, and based on the data of Figure 4, the rod-like chains present an isotropic-nematic transition at ϕ nem ≈ 0.20, a volume fraction significantly lower than the freezing point for monomeric hard spheres (ϕ F monomers = 0.494), in accordance with the trends presented in [106,107]. At volume fractions close to the freezing point, the rod-like systems reach a practically perfect PRO phase. According to Shakirov and Paul [87], the rod-like chains self-arrange into a 2D hexagonal crystal structure in the plane perpendicular to the nematic director at a volume fraction ϕ = 0.47, driven by 2D translational entropy. Thus, those 2D hexagonal crystal structures evolve into a 3D semi-crystalline phase of RHCP structures with the increase in concentration in the range between the freezing and melting point for monomeric hard spheres. After reaching the melting point for monomeric hard spheres, the nematic systems start to crystallize in the HCP, FCC, or mixed HCP/FCC structures, maintaining the nematic order. Thus, rod-like systems pack at high densities into Nematic Close Packed (NCP) structures, as was observed by independent researchers through MD simulations [85,86].

Discussion and Conclusions
We present the results from extensive Monte Carlo simulations on the phase behavior of semi-flexible chains of tangent hard spheres as a function of the packing density and equilibrium bending angle while fixing the temperature, spring bending constant, and average chain length. The local structure is quantified through the Characteristic Crystallographic Element (CCE) norm, while the global structure is gauged through the nematic order parameter. A rich one-dimensional phase diagram as a function of packing density is identified where chains crystallize in close-packed morphologies, including random hexagonal close (RHCP) ones of single or multiple stacking directions, or in almost perfect HCP and FCC crystals in the case of rod-like chains. The analysis of the long-range orientational tensor reveals the formation of prolate mesogen, nematic mesophase (PRO) for rod-like chains at rather low volume fractions and of oblate mesogen, nematic mesophase (OBL) at high packing densities. Although all of the systems of semi-flexible chains crystallize, the equilibrium bending angle significantly affects the melting point. While equi- Figure 11. Orientational order parameter, q ((left) y-axis, blue color), and crystallinity, τ c ((right) y-axis, red color), as a function of packing density, ϕ, for 100-chains N av = 12 systems of semi-flexible chains with bending constant k θ = 9 and equilibrium bending angle θ 0 = 0 • (left) and θ 0 = 90 • (right). Dashed lines connecting points are used as visual support. Vertical dashed lines correspond to the freezing (ϕ F monomers ≈ 0.494 ) and melting (ϕ M monomers = 0.545 ) transition for monomeric hard spheres [29].
Based on the local and global phase behavior, the one-dimensional diagram as a function of the packing density can be split into three distinct regions, marked by XX-YY, where the first index XX corresponds to the local structure (XX = AMO or CRY) and the second index YY to the global structure (YY = ISO, PRO or OBL). Accordingly, for the rod-like chains studied here (N av = 12, k θ = 9), we have the following approximate domains of phase behavior (left panel of Figure 11): (i) AMO-ISO (ϕ ≤ 0.15), where the system is amorphous at the local level and isotropic at the global; (ii) AMO-PRO (0.20 ≤ ϕ ≤ 0.45), where the packing is amorphous locally and nematic globally; and (iii) CRY-PRO (0.50 ≤ ϕ), where the system shows crystallinity in the form of HCP and FCC close packed morphologies of varied levels of perfection, and perfect nematic ordering at the global orientational level. Comparing the local and global trends, a delay is evident as the local order is established at volume fractions significantly higher than the ones for the ISO → PRO transition.
In the case of the right-angle chains (θ 0 = 90 • ), the phase behavior is quite different compared to the rod-like chains, as seen in the right panel of Figure 11. First, we should note that the small decrease in the long-range orientational order can be attributed to higher statistical uncertainty as the higher the packing density, and the closer to jamming, the more difficult the sampling, even when such advanced MC algorithms are employed. Accordingly, significantly longer simulations are required and the statistical uncertainty is higher. The right-angle systems do not crystallize until after they reach volume fractions higher than the melting point for monomeric hard spheres (ϕ M monomers = 0.545). Then, the ISO → OBL transition occurs simultaneously to the crystal nucleation and growth. The formation of RHCP structures requires a semi-nematic OBL phase for the right-angle systems as the right-angle chain exists but is not dominant in the FCC and HCP crystals. Accordingly, for their creation, specific chain alignments are required. As a result, the local and global orders are simultaneously established, and no delay is observed. The phase diagram consists of only two regimes: (i) AMO-ISO (ϕ ≤ 0.56), where the system is amorphous at the local level and isotropic at the global level; and (ii) CRY-OBL (ϕ > 0.56), where the system shows RHCP morphologies and a nematic order of oblate mesogens.

Discussion and Conclusions
We present the results from extensive Monte Carlo simulations on the phase behavior of semi-flexible chains of tangent hard spheres as a function of the packing density and equilibrium bending angle while fixing the temperature, spring bending constant, and average chain length. The local structure is quantified through the Characteristic Crystallographic Element (CCE) norm, while the global structure is gauged through the nematic order parameter. A rich one-dimensional phase diagram as a function of packing density is identified where chains crystallize in close-packed morphologies, including random hexagonal close (RHCP) ones of single or multiple stacking directions, or in almost perfect HCP and FCC crystals in the case of rod-like chains. The analysis of the long-range orientational tensor reveals the formation of prolate mesogen, nematic mesophase (PRO) for rod-like chains at rather low volume fractions and of oblate mesogen, nematic mesophase (OBL) at high packing densities. Although all of the systems of semi-flexible chains crystallize, the equilibrium bending angle significantly affects the melting point. While equilibrium angles of 108 • and 120 • degrees favor crystallization compared to the freely-jointed model, chains with 90 • show a behavior that almost coincides with the fully flexible chains, and the acute angle of 60 • hinders crystallization, enforcing nucleation and growth to take place at higher concentrations.
In particular, for rod-like chains, three distinct regimes can be identified in the onedimensional phase diagram, where a delay is observed between local and global selforganization: (i) AMO-ISO (ϕ ≤ 0.15), where the system is amorphous at the local level and isotropic at the global; (ii) AMO-PRO 0.20 ≤ ϕ ≤ 0.45), where the packing is amorphous locally and nematic globally; and (iii) CRY-PRO (0.50 ≤ ϕ), where the system shows crystallinity in the form of HCP and FCC close packed morphologies of varied levels of perfection and perfect nematic ordering at the global orientational level. Right-angle systems show a synchronous establishment of long-range nematic orientational order and formation of RHCP crystallites, thus splitting the behavior into two distinct regimes: (i) AMO-ISO (ϕ ≤ 0.56), where the system is amorphous at the local level and isotropic at the global level; and (ii) CRY-OBL (ϕ > 0.56), where the system is crystalline at the local level and shows a nematic mesophase of oblate mesogens at the long-range.
The present simulations are currently expanded to treat semi-flexible chains of tangent hard spheres in composites with nanofillers, under confinement, and in mixtures with different species in the form of linear chains and monomeric counterparts.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/polym15030551/s1, Figure S1: Snapshots of systems of semi-flexible chains of tangent hard spheres (N ch = 100, N = 12) at a packing density ϕ = 0.10 with a bending constant k θ = 9 and different equilibrium bending angles, θ 0 ; Figure S2: Snapshots of the final configurations from MC simulations on semi-flexible chains of hard spheres at ϕ = 0.60 and different equilibrium bending angles, θ 0 ; Figure S3: Exponential moving average of the orientational order parameter, q, as a function of MC steps for 100-chains N av = 12 semi-flexible system with bending constant and equilibrium bending angle θ 0 = 0 • at different packing densities, ϕ; Figure S4: Final configuration for the 100-chains N av = 12 semi-flexible system chains with bending constant k θ = 9 and equilibrium bending angle θ 0 = 0 • at ϕ = 0.58 and 0.60; Figure S5: Degree of crystallinity, τ c , as a function of MC steps at different packing densities, ϕ. Funding: This research was funded by MICINN/FEDER (Ministerio de Ciencia, Innovación y Universidades, Fondo Europeo de Desarrollo Regional), grant numbers "RTI2018-097338-B-I00" and "PID2021-127533NB-I00" and by UPM and Santander Bank, "Programa Propio UPM Santander".

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used throughout the manuscript: 2D Two