Polymorph Stability and Free Energy of Crystallization of Freely-Jointed Polymers of Hard Spheres

The free energy of crystallization of monomeric hard spheres as well as their thermodynamically stable polymorph have been known for several decades. In this work, we present semianalytical calculations of the free energy of crystallization of freely-jointed polymers of hard spheres as well as of the free energy difference between the hexagonal closed packed (HCP) and face-centered cubic (FCC) polymorphs. The phase transition (crystallization) is driven by an increase in translational entropy that is larger than the loss of conformational entropy of chains in the crystal with respect to chains in the initial amorphous phase. The conformational entropic advantage of the HCP polymer crystal over the FCC one is found to be ΔschHCP−FCC≈0.331×10−5k per monomer (expressed in terms of Boltzmann’s constant k). This slight conformational entropic advantage of the HCP crystal of chains is by far insufficient to compensate for the larger translational entropic advantage of the FCC crystal, which is predicted to be the stable one. The calculated overall thermodynamic advantage of the FCC over the HCP polymorph is supported by a recent Monte Carlo (MC) simulation on a very large system of 54 chains of 1000 hard sphere monomers. Semianalytical calculations using results from this MC simulation yield in addition a value of the total crystallization entropy for linear, fully flexible, athermal polymers of Δs≈0.93k per monomer.

The structure of crystals of uniform, monomeric hard colloidal spheres, as investigated experimentally via light/X-ray scattering and confocal microscopy, is often found to be a random stacking of 2D hexagonal compact layers (rHCP) [12,[33][34][35][36][37]. In some cases, depending strongly on experimental conditions, including factors like size polydispersity, shear, and gravity [34,35], quite perfect face-centered cubic (FCC) crystals are obtained [38], typically in samples grown over weeks or months, to allow for a slow annealing or aging for the transition rHCP→FCC to take place [12,13,39,40].
A widely accepted value of the entropy difference between the FCC and HCP polymorphs, supported by simulations [41][42][43][44][45], is ≈ 112(±4) × 10 −5 k per particle, where k is Boltzmann's constant. Variations from this value depend on the conditions under which the estimation is made, and especially on packing density (volume fraction) [41,42,44,46,47]. This very small value is qualitatively consistent with the experimentally observed sluggishness of the rHCP→FCC transformation. Accordingly, under a constant volume, most studies identify the rHCP polymorph as the final crystal structure [21,28,48,49] and perfection in the form of the FCC crystal is only rarely encountered [50,51]. Interestingly, neither perfect nor defective HCP crystals seem to have been observed in experiments and/or numerical simulations starting from densely packed amorphous samples.
Fully flexible polymers of hard spheres lacking any other type of interactions (i.e., without bending, torsional or bond length energetic contributions to the Hamiltonian) are athermal so that their equilibrium phase behavior is driven solely by entropy. Our previous work [75][76][77][78][79][80] showed that, rather counter-intuitively, starting from amorphous packings, freely-jointed chains of tangent hard spheres do indeed undergo spontaneous entropy-driven crystallization under a variety of conditions. In constant-volume simulations starting from an amorphous packing crystallization sets in when the increase in translational entropy of the monomers compensates for the loss of conformational chain entropy. The increase in monomer translational entropy is caused by an increase of the space available to monomers and of its isotropy [75,76], much as it happens in the transition from the isotropic fluid to the nematic liquid crystal mesophase in Onsager's theory [24,25]. Due to the high complexity and computational demands of the algorithms needed to simulate long, athermal polymers at very high densities, the polymer crystals obtained in previous Monte Carlo (MC) simulations were made of short chains and tended to display a fivefold-free, but defect-ridden rHCP structure, so that it was not possible to establish which of the two competing crystal forms, HCP or FCC, was the thermodynamically stable one for chains of hard spheres.
Very recently we extended the studies by conducting unprecedentedly long MC simulations to study the entropy-driven crystallization of a large system of 54 chains of average length 1000 (in a number of spheres). In these isochoric simulations the initially amorphous configuration, after an early dominance of the HCP polymorph, passes to a transitory rHCP morphology and eventually reaches a stable FCC crystal of very high perfection [78]. In this work, which can be considered as a companion to [78], we present semi-analytical calculations of the free energy of crystallization of the stable FCC crystal and of its free energy advantage with respect to the HCP polymorph to support the observed simulation trends.

Free Energy Difference between FCC and HCP Polymorphs
Athermal polymers are represented here as linear, freely-jointed chains of hard spheres of uniform size, σ, which is further the characteristic (unit) length of the system. According to the hard core model, spheres i and j interact with a pair-wise energy u HS (r ij ), given by: The spontaneous crystalization of polymers is conveniently studied in the (trivially) athermal version of the isochoric semigrand canonical ensemble [VTN sites µ * i ] in which total volume V, total number of sites N sites and a spectrum µ * i of chemical potentials are specified [82,83]. This ensemble naturally allows for polydispersity and for any desired distribution of polymer chain lengths spanning an arbitrary interval l ∈ [l min , l max ], so that results are more generally valid than by assuming a specific distribution (Flory, uniform, etc). We consider a system of N such polymer chains comprising a total number of monomers (also "sites") N sites . The chain length distribution is given by the number of chains N l of length l ∈ [l min , l max ]. In the following and for compactness of notation we also use N l to refer to the entire distribution. Although both N sites and N are fixed, the numbers N l of chains of each length l are fluctuating variables. The desired polymer length distribution results from imposing a suitable spectrum of chemical potentials, so that the constraints hold, with l the counter over chain lengths. The partition function in the [VTN sites µ * i ] ensemble is [84]: where q l refers to the translational and internal contributions to the partition function for chain length l, and Z(V, T, N l min , . . . , N l max ) = d 3N r exp(−βU(r)) is the classical configurational integral. The athermal ensemble is trivially obtained by setting U(r) → ∞ whenever at least one monomer overlap exists, and 0 otherwise (Equation (1)), so that the classical configurational integral reduces to a summation over equally probable microstates. The analytical evaluation of the general partition function (Equation (3)) is not feasible. It is however the starting point both for an analytical estimate of an upper bound of the entropy difference between crystal polymorphs, and for the development of MC algorithms. Previous work [75,76] on the crystallization of chains of hard spheres has unequivocally established that the positions of chain monomers in the crystal fluctuate about the most probable sites of a well-defined polymorph (either HCP or FCC; other, non-compact crystal polymorphs do not appear in experiments or in simulations), while chains adopt random conformations compatible with the monomers fluctuating about these sites of the perfect crystal. Crystallization takes place because the loss of conformational entropy of the chains is more than compensated for by the increase in the positional/translational entropy of their monomers, even if these are forced to remain close to the sites of the crystal. The increase in monomer translational entropy is due to the larger and more isotropic volume translationally available to the monomers, in exact parallelism to what happens for monomeric hard spheres and in lyotropic phase transitions in liquid crystals as pioneered in [24] and further extended in [85][86][87].
If we ignore spatial fluctuations around their average positions, a crystal of hard sphere monomers consists of a single microstate (i.e., either the ideal FCC or HCP crystal, as specified by their lattices and bases). However, an interesting feature of the polymer crystal is that, again ignoring spatial monomer fluctuations, it possesses a large number of equally probable microstates: all possible multi-chain configurations in which the chains connect adjacent crystal sites without overlap and without leaving empty sites, while respecting monomer tangency along the chain backbone, are valid microstates. Thus, polymer crystal microstates are obtained to a very high degree of approximation as the product set of monomeric positional microstates (about the sites of the perfect crystal) and a highly multiply degenerate set of chain conformational microstates (all possible chain conformations that join the monomers in chains of the specified length).
The above description of the polymer crystal makes it possible to develop an accurate approximation to the partition function (Equation (3)) which is amenable to analytical calculations for the present case of crystals of freely jointed chains of hard spherical monomers by factorizing (Equation (3)) into a translational part for the individual monomers, as in a crystal of monomers, and a configurational part which accounts for chain connectivity and conformational variability. To this end we denote by R X i the coordinates of the center of the i-th spherical monomer in a perfect crystal of monomers of polymorph X (where X = FCC or HCP) and by R X ≡ N sites i=1 R X i the set of all monomer positions in the perfect crystal (i.e., R X is a list of 3N sites Cartesian coordinates). The lowercase versions, namely r X i and r X , denote the coordinates of the i-th spherical monomer in a specific configuration of the real crystal (i.e., subjected to fluctuations), and the set of all r X i , respectively. Given R X for a finite sample consisting of N sites of polymorph X there exists a finite set of polymer chain configurations that are obtained by tracing all sets of N nonoverlapping (i.e., simultaneously self-avoiding and mutually-avoiding) paths of the prescribed chain length distribution that connect all the N sites points of coordinates R X i . We now denote by I X jk the j-th multichain configuration for the given R X and for the k-th chain size (discrete) distribution N l , constrained by (Equation (2)). The integer counter k enumerates all possible chain distributions compatible with these constraints. Notice that it is not necessary to actually have the explicit complete finite set of possible chain distributions, nor its cardinality, because this set is independent of, and therefore the same for, all crystal polymorphs. For any given, fixed numbering scheme of the N sites monomers and of the N chains, each I X jk is (up to permutations of the site labels) a list of N sites + N − 1 integers which specify which and in which order crystal sites R X are occupied by monomers belonging to which chain (N sites integers), and which the chain lengths are (N − 1 integers, because of the first constraint above). The union set ξ X k ≡ j I X jk is then the set of all possible multichain configurations for a given chain length distribution, and the double union Ξ X ≡ k j I X jk is the complete, finite set of all possible multichain configurations for all possible chain length distributions compatible with N l , l ∈ [l min , l max ], and constrained by (Equation (2)).
Accepting the separation of translational and conformational degrees of freedom in (Equation (3)) (an assumption whose plausibility will be quantitatively discussed below), and in terms of the previous definitions, the free energy/entropy difference ∆S X−Y ≡ S X − S Y between polymorphs X and Y is, in units of k: where Z X m (V, r X i ) denotes the classical partition function for the monomeric crystal X and | | denotes the cardinality of a set.
In physical terms, (Equation (4)) splits the evaluation of the entropy of the crystal of chains in two independent, additive contributions: the first one, ∆S X−Y m (corresponding to the first fraction in (Equation (4)) due to translational degrees of freedom of the monomers, as if they were not connected to form chains; the second one ∆S X−Y ch due to the number of ways the individual monomers can be connected into N chains of the specified length distribution under the condition that the monomers occupy all sites of the perfect crystal, i.e., the chains are both self-avoiding and mutually-avoiding, simultaneous random walks on the sites of polymorph X, so that no crystal sites are left unoccupied.
One advantage of (Equation (4)) is that ∆S X−Y m is the entropy difference between the FCC and HCP crystals of monomeric hard spheres, which is precisely known [41,43,44]. The problem is then reduced to the calculation of ∆S X−Y ch , i.e., the chain configurational and conformational contribution to entropy. If ∆S HCP−FCC ch > |∆S FCC−HCP m | then the HCP crystal would be the stable polymorph, while FCC would be the thermodynamically preferred crystal form in the opposite situation.
Neither the cardinality ξ X k nor consequently Ξ X are known for any crystal type; furthermore, their evaluation by exhaustive enumeration is conjectured to be an NP-complete problem [88] so that (Equation (4)) and ∆S FCC−HCP cannot be evaluated exactly in the general case. It is however possible to establish an upper bound for ∆S X−Y ch , which, if tight enough, would be sufficient to prove the stability of the FCC polymorph. An upper bound for ξ X k is given by the product of the number of N neither self-avoiding nor mutually avoiding random walks of the given length distribution on the crystal sites of the polymorph, given by ∏ l max l min (C X ) N l , with C X the coordination number of the crystal. However, since C FCC = C HCP = 12 for both FCC and HCP, this bound predicts ∆S FCC−HCP = 0 and does not allow to discriminate between the two crystal types.
A tighter and discriminating bound can be obtained by considering the finite set c X (l) of single self-avoiding random walks (SAWs) of length l on the sites of a given crystal X. The asymptotic dependence of c X (l) on l is of the exponential-power law type [89]: where A is the critical amplitude, µ the connective constant, and γ the critical exponent. For freely jointed chains, the asymptotic regime of (Equation (5)) is already attained for chains of very moderate length (l ≈ O(10)), so that (Equation (5)) is valid with excellent accuracy in the polymeric regime for which l O (10) Although it was conjectured that c HCP (l) = c FCC (l) , direct enumeration of SAWs demonstrates that this is exact only for l ≤ 6 and first-order accurate for l > 6; above that value of l, and hence in the polymeric regime of interest here, c HCP (l) > c FCC (l) by a very small amount. We have calculated the cardinalities c FCC (l) and c HCP (l) as a function of SAW length l by exhaustive enumeration and also established their asymptotic behavior, which is accurately given by (5) with A = 1.19, µ = 10.07 and γ = 1.134 for c FCC (l) . For the HCP crystal, the values of these parameters are so similar to the FCC parameters, that it is numerically much preferable to express the minute difference between both in the expected decaying exponential form: with a 1 = 3.31 × 10 −6 , a 2 = 8.63 × 10 −6 , and a 3 = 0.24, for l > 5 and a 1 = a 2 = a 3 = 0 for l ≤ 5.
Since chains in the polymer crystal are random walks that are simultaneously mutuallyavoiding and self-avoiding, and display ideal (non-SAW) statistics [90], their conformational ensemble is guaranteed to be a proper subset of the Cartesian product of the sets of all possible (self-avoiding and non-self-avoiding) individual chain conformations. Thus, for any given chain length distribution, N l , the monotonicity of (Equation (5)) guarantees that the cardinality ξ X k is strictly bounded from above by: If we now denote by N dist the number of possible chain length distributions N l compatible with N l , l ∈ [l min , l max ], and constrained by (Equation (2)), the following bound for Ξ X results: where the sum is carried out over all possible chain length distributions N l , with l ∈ [l min , l max ], and constrained by (2). Note that N dist is independent, i.e., the same, for all crystal polymorphs and will cancel in any ratio of the type |Ξ X | |Ξ Y | , such as Equation (9) below.
Because c HCP (l) is strictly larger than c FCC (l) , because Ξ HCP ≈ Ξ FCC to first order, and because chain distributions are independent of polymorph type: An upper bound for the entropy difference per monomer of the two polymorphs X and Y, ∆s X−Y ch , is then given by: where l is the number average chain length. The smallness of a 1 and a 2 , and the polymeric regime (large l) allow us to only retain the leading term in the expansion of the logarithm, so that the last inequality follows from the asymptotic behavior of (Equation (6)). Then, to the second order: which implies that the chains in the HCP polymorph have higher conformational entropy than in the FCC crystal. However, this difference in conformational entropy is insufficient by more than two orders of magnitude to offset the translational entropic advantage ≈ 112(±4) × 10 −5 k of the FCC polymorph with respect to the HCP crystal. Thus, the upper bound (Equation (9)), although conservative, is tight enough to unequivocally show that the FCC is the thermodynamically stable polymorph for crystals of linear, freely-jointed polymers of tangent hard spheres.

Monte Carlo Simulations
In order to provide computational support for the above results, based on semianalytical calculations, in a companion paper [78] we carried out unprecedentedly extensive MC calculations of very large systems deep in the polymeric regime (N = 54 chains comprising N sites = 54, 000 monomers, with a flat chain length distribution in the interval [l min = 600, l max = 1400]), at a volume fraction ϕ = 0.56, and starting from a totally amorphous (random) packing. Observing spontaneous crystallization in a system of such size at that high volume fraction requires very efficient and proper configurational sampling, which is a very challenging task due both to the high density of the system and to the great length of the polymer chains.
One of the main motivations for simulating such a large system is to minimize potential finite size effects of small cells under periodic boundary conditions, which might conceivably lead to a crystallization advantage for the FCC polymorph due to the incommensurability of the cubic cell with the HCP crystal. Cell incommensurability is an irrelevant factor for such a large cubic cell, which is for all practical purposes just as compatible with the HCP crystals that appear in simulations as the cell usually known as "rhombic" (strictly speaking, hexagonal). Our simulation cell is significantly larger than individual chains and incipient HCP crystallites have ample possibility to freely nucleate and grow in the bulk of the system, irrespective of cell shape. For example, the size of the cubic cell is approximately 37 (measured in units of σ), more than twice the average radius of gyration R g = 16.5, where denote the average over all chains and system configurations.
In fact, as clearly seen in the following snapshots ( Figure 1) HCP crystalline domains do actually form, and even temporarily surpass in abundance the FCC crystalline regions during part of the simulation. The disappearance of HCP domains and the final formation of a pure FCC crystal is not a consequence of using a cubic simulation cell. If anything, the cubic simulation cell would make a presumptively stable HCP crystal only very slightly more defect prone than a hexagonal one, but this effect would be completely obscured by the many imperfections that inevitably appear in all computer simulations of spontaneous crystallization. As a matter of fact, relatively perfect HCP crystals appear in simulations of crystallization in confined cubic cells where spatial restrictions appear in the form of flat parallel and impenetrable walls [79,80]. The cubic shape of the box can definitely be discarded as the cause of the final dominance and greater stability of the FCC polymorph.
System configurations are first generated and equilibrated through an MC scheme based on algorithms specially designed for the efficient sampling of polymer-based systems [83,91,92] and then characterized through the Characteristic Crystallographic Element (CCE) norm descriptor [93], both modules as implemented in the Simu-D software suite [94]. The MC trajectory is generated in the [VTN sites µ * ] ensemble, where V is the volume of the cell, T is the temperature (inactive here due to the athermal nature of the model), and µ * is the spectrum of chemical potentials used to control the distribution of chain lengths. The MC scheme consists of the following moves: (i) reptation (10%), (ii) rotation (10%), (iii) flip (34.8%), (iv) intermolecular reptation (25%), (v) configurational bias (20%), (vi) simplified end-bridging (0.1%) and (vii) simplified intermolecular end-bridging (0.1%), where numbers in parenthesis denote attempt percentages. Due to the very high volume fraction all local moves (i-v) are executed in a configurational bias pattern with the number of candidate configurations set equal to 50. Trial MC moves are accepted or rejected according to the modified Metropolis criteria as explained in Ref. [83]. The total duration of the MC simulation is 1.4 × 10 12 steps with a record frequency of snapshots (frames) set at 1 × 10 8 leading to a final trajectory being composed of 14,000 frames. The radial and orientational similarity of the local environment around each site with respect to reference crystals, as quantified through the CCE norm [93], classifies them as HCP, FCC, FIV (fivefold) or AMO (amorphous, or more precisely "not identified") character. Throughout the present manuscript, blue, red, green and yellow will be used to represent the HCP, FCC, FIV and AMO sites and curves, respectively. The exact details on the methodology followed for the MC simulations and the successive structural analysis of the computer-generated system configurations can be found in the companion publication [78].  [93] with blue, red and green corresponding to sites with HCP, FCC and FIV character, respectively. Amorphous (AMO) sites are shown in yellow and with reduced dimensions for visual clarity. The stable FCC crystal (fourth, rightmost snapshot) is obtained in the steady state (up to fluctuations) MC production phase, after approximately 9 × 10 12 MC steps; (B) Sites are colored according to their parent chain and are shown with wrapped coordinates, subjected to periodic boundary conditions; (C) Sites are colored according to their parent chain and are shown with coordinates fully unwrapped in space; (D) Two randomly selected chains are shown in red and blue with sphere coordinates fully unwrapped in space. Image panels created with the VMD software [95]. Details on the MC simulation and the corresponding trajectory can be found in [78].
The evolution of the fraction of sites with HCP, FCC and FIV characters as a function of MC steps can be found in Figures 1 and 2 of [78]. Crystallinity is simply the summation of HCP and FCC fractions, while the degree of disorder of the system can be directly mapped into the fraction of AMO sites. Based on the observed trends the phase transition can be naturally split into 4 regions: (I) the rapid nucleation and growth of crystals with HCP and FCC character and the parallel decrease in the population of FIV and AMO sites, (II) the induction period where crystallization slows down, the population of HCP sites remains constant, the one of FCC increases very slowly while surviving FIV sites form characteristic linear assemblies corresponding to cyclic twin structures, (III) the FCC growth period which is accompanied by the elimination of sites with HCP similarity leading eventually to the formation of a single FCC crystal of very high perfection, (IV) the final, steady-state region, where within fluctuations, the established FCC crystal remains unaltered. The thresholds between the regions are marked by the significant slowdown in the crystal growth (I → II), the disappearance of fivefold sites (II → III), and the formation of the perfect FCC crystal (III → IV). A video showing the evolution of crystallization and the transition between the different crystal polymorphs can be found in the Supplementary Material, while system snapshots, corresponding to the four distinct regions, are presented in Figure 1. As can be seen in the supplementary video material and in the snapshots, starting from a purely amorphous, statistically homogeneous and isotropic configuration, the MC simulation is able to evolve the system through intermediate states of increasing crystallinity until a stable FCC polymorph of remarkable perfection is formed. In this final steady state the percentage of sites with FCC similarity exceeds 90% of the total population. A detailed explanation of the determination of the crystallographic type, of fivefold-symmetric, and of amorphous sites as implemented in the CCE norm descriptor is given in the companion paper [78] together with an analysis of the evolution including the entropic origins of crystal perfection. Based on the observed trends it is clear that driven by the minute FCC entropic advantage, the system shows a transition from the original disordered solid to the ordered crystal and spontaneously generates microstates of increasingly FCC character, until after approximately 9 × 10 12 MC steps when all HCP crystalline sites disappear and crystallization reaches completion.
The initial predominance of the HCP polymorph (Regime I) and its eventual complete disappearance during the evolution towards a stable, perfect FCC polymer crystal is distinctly different from the evolution during crystallization of single spheres, where crystals of mixed FCC/HCP character are invariably obtained. This phenomenon is related to the different relative (meta)stability of intermediate states for crystals of single hard spheres, and of polymers of hard spheres, and is fully explained in the accompanying work [78].

Decorrelation of Translational and Conformational Degrees of Freedom
We now examine the postulated decorrelation of translational and conformational (torsion and bending angles) degrees of freedom (d.o.f.'s), which is the basis of the calculations leading to the result of (Equation (11)). In order to test this hypothesis, we extract from the MC results the correlation between translational (monomer displacement) and conformational (bending and torsion angles of the chains) d.o.f.'s. To this end, we plot (Figure 2), for all monomers, and in Region IV the value of |R| against the torsion angles in which each monomer is involved, where |R| is the distance between a monomer and the centroid of its Voronoi cell: 12) where N V is the number of vertices of the given Voronoi cell, r i are the position vectors of its N V vertices, and r m is the position vector of the monomer. Figure 3 is a similar plot of the value of |R| versus the bending angle whose vertex is located at the monomer.  In both figures the size of the symbols is proportional to the probability density of the particular value of the torsion (Figure 2) or bending (Figure 3) angle, and information on the fluctuation amplitude is included as error bars. For reference, the probability distributions of both φ and θ are included in the same plots. For both types of confor-mational degrees of freedom, the flatness of the curves attests to the independence of translational and torsional and bending conformational d.o.f.'s. For the torsional angle φ, the value of |R| is constant to better than 2%, and better than 5% for the bending angle. The values which deviate the most from the average (e.g., for θ = 30 o ) correspond to angles that occur with extremely low frequency: p bend (θ = 30 o ) ≈ 0, for which the symbol is invisibly small in Figure 3. Only very few instances of bending angles around this value exist, as the local chain geometry they produce is incompatible with the FCC cell, so they are very strongly suppressed in the crystal (flat portion around θ = 30 • in Figure 3) and are for all practical purposes irrelevant.
From these two figures it can be concluded that the decorrelation between translational (|R| or its individual components, not shown) and conformational (φ, θ) d.o.f.'s is fulfilled to a sufficient degree to warrant their separation in the calculation of ∆s HCP−FCC , and the result (Equation (11)), between the two competing polymorphs HCP and FCC in the crystal.
In summary, although the free energy (entropy) advantage of the FCC over the HCP crystals of monomeric hard spheres is small (≈ 112(±4) × 10 −5 k), the higher conformational entropy of chains in the HCP crystal with respect to FCC is even tinier ≈ 0.331 × 10 −5 k and hence insufficient to make HCP the preferred crystal for chains. The results of the MC simulations strongly support the analytical result that the FCC polymorph is the thermodynamically stable phase for crystals of fully flexible chains of hard spheres.

Free Energy of Crystallization
Crystallization of fully flexible chains of hard spheres entails a reduction (relative to chains in an amorphous solid) in the number of conformations available to, and hence in the entropy of, the polymeric chains. However, just as in the case of monomeric hard spheres, the positional entropy of the individual monomers increases upon crystallization by a larger amount, so the overall result is still favorable for crystallization.
Thanks to the separation of translational and conformational d.o.f.'s, the increase in (monomeric) translational entropy upon crystallization of the chains can be taken as identical to the translational entropy increase of crystallization of monomeric hard spheres, which has been known for a long time: ∆s trans m = 1.17k per monomer [32,96,97]. In order to compute the free energy (entropy) of crystallization of fully flexible polymers of hard spheres, it is now necessary to subtract from ∆s trans m the value of the conformational entropy lost by chains when they crystallize.
While an analytically exact calculation of the entropy of crystallization is not feasible, the results from the MC simulation suggest a feasible approximation method. First, the Kuhn length b 0 of the chains does not significantly change upon crystallization, as we have shown in [78]. The Kuhn length b 0 = 1.52 ± 0.05 turns out to be only 50% longer than the bond length, implying that the chains behave as identical freely jointed ideal (non-self-avoiding) chains when observed at length scales beyond a few bonds, in both the disordered solid and the ordered crystal.
In Figure 4, we present the characteristic ratio C n (main panel), and the ratio , as a function of chain length, l. R 2 ee is the mean square end-to-end distance, and R 2 g is the mean square radius of gyration. Characteristic ratio is defined as C n = Ree 2 (l−1) b len , where b len is the average bond length, as explained in Ref. [78]. Both are, within statistical uncertainty, independent of chain length l ∈ [l min , l max ] and the same for both Regions, indicating that chains in the initial disordered solid (early Region I) and in the almost perfect FCC crystal (Region IV) are indistinguishable in the degree of coiling/flexibility and also equally ideal, since the ratio R 2 ee R 2 g adopts the value of 6 [98,99], as expected for ideal unperturbed polymers in the long-chain limit, as demonstrated originally by Debye [100] and then by Flory [90]. Furthermore, the distribution of the end-to-end vector of the chains, P(|R ee |) whose lengths lie within the small interval l ∈ [970, 1030] (chosen as a representative, narrow interval of chain lengths) is also the same in the disordered solid and in the ordered crystal ( Figure 5); the polymers neither shrink nor swell upon crystallization. These results unambiguously demonstrate that the large-scale features of the chains remain unchanged upon crystallization. In the two distinct phases (disordered solid and ordered crystal) chains differ only in their very small-scale features (a few bonds).
This statement may seem paradoxical at first sight but it is actually natural: the way chains crystallize (by having their monomers occupying on average the most probable sites of a crystal) uniformly (at all chain lengths) selects the chains that fulfill this condition from the other members of the ensemble. In slightly more precise terms, for a given crystal polymorph X and for a given chain length distribution, denoted by k (Section 2.1), the union set ξ X k is (up to fluctuations) an equivalence class of the ensemble [VTN sites µ * i ] under the relation "monomer coordinates in a configuration of the ensemble [VTN sites µ * i ] coincide with site coordinates of crystal X", a class for which any of the I X jk can be taken as a representative. As a consequence of the underlying uniform (average) spatial density of crystal sites, the equivalence class retains all large-scale features of the full ensemble.  Figure 1 snapshots of the whole system with chains being subjected to periodic boundary conditions (second row), being fully unwrapped in space (third row) for all four Regions exhibited during crystallization. In addition, only two, randomly selected, chains are visualized for each system for clarity in the bottom (fourth row) of Figure 1. In all cases, the chain features, down to a few bond lengths, are indistinguishable so that it is impossible to tell apart polymers in the original amorphous solid and in the various crystal polymorphs including the perfect FCC crystal.
Since conformational differences between chains in the disordered solid and the FCC crystal are overwhelmingly local, a calculation of the loss of conformational entropy upon crystallization ∆s con f m can be made, inspired by Flory's Rotational Isomeric State (RIS) theory [90] by considering continuously varying bending and torsional angles. Although the absolute value of the conformational entropy cannot be computed exactly either in the amorphous or in the crystal phases, the previous arguments make it possible to estimate the difference in conformational entropy between the amorphous solid and the crystal by considering up to next-next-next-nearest neighbors, i.e., a four-monomer portion of chain, beyond which chains remain essentially unaltered upon crystallization.
We denote by f am bt (θ 1 , φ, θ 2 ) and f cr bt (θ 1 , φ, θ 2 ) the joint orientational distribution functions of bending and torsion angles for a four-monomer chain segment (see Figure 6) in the amorphous solid and in the crystal, respectively. The angles θ 1 , θ 2 are two consecutive bending angles defined by monomers 1-2-3, and 2-3-4 (see Figure 6), and φ is the torsional angle defined by monomers 1-2-3-4. The configurational entropy difference upon crystallization (in units of Boltzmann's constant k) can be expressed in terms of the integrals of the joint orientational functions [101] by: where the integral operation . . . du is equivalent to . . . sin θdθdφ with integration ranges 0 ≤ θ ≤ π and 0 ≤ φ < 2π. The integration over the unit vectors u (1) , u (2) and the torsion angle φ are tantamount to carrying out the integration over all possible states of the four-monomer segment.
The orientational functions f cr bt (θ 1 , φ, θ 2 ) and f am bt (θ 1 , φ, θ 2 ) are unknown a priori but can be evaluated numerically as averages over Region IV (final FCC crystal) and over the initial amorphous state of Region I, respectively, and discretized on regular integration meshes of increasing resolution ranging from 20 × 20 × 20 up to 60 × 60 × 60 for all θ 1 , φ, θ 2 in order to ensure numerical convergence of the integrals.
Figures 7 (as 3D isosurfaces), 8 and 9 (as sections of these isosurfaces at several values of φ) illustrate how the preferred bending and torsional angles differ between the amorphous and the crystalline states. As can be seen in Figure 7, specific triads of θ 1 , φ, θ 2 appear in the stable crystalline polymorph with significant frequency. The highly probable combinations of θ 1 , φ, θ 2 (disconnected, high probability regions in the right panel of Figure 7) are responsible for the instantaneous positions of the monomers fluctuating about the sites of the ideal FCC crystal.
A numerical evaluation of Equation (13) yields an entropy loss (i.e., the entropy of the chains in the crystal state is lowered by the loss of conformational freedom) ∆s con f = −0.24k ± 0.04 per monomer. This figure is significantly lower (in absolute value) than the translational entropy increase due to the formation of the crystal, which is ∆s trans m = 1.17k per monomer, as mentioned above. The free energy of crystallization of fully flexible chains of hard spheres (measured by the net increase in entropy per monomer upon crystallization of the athermal system) is then ∆s m = ∆s trans m + ∆s con f ≈ 0.93k, which is the result we sought. This value is still more than sufficient to drive the phase transition, as in monomeric hard spheres, in spite of the loss of chain conformational entropy.

Conclusions
We have presented semianalytical calculations of the free energy of crystallization of linear, freely jointed chains of tangent hard spheres, as well as of the free energy difference between the FCC and HCP polymorphs. The calculations are based on a separation of degrees of freedom, i.e., a decoupling of chain conformational and monomer translational degrees of freedom. This postulate is confirmed to hold within narrow error margins. The calculations predict a small advantage of the FCC crystal over the HCP, which makes FCC the thermodynamically stable one. This prediction is consistent with the results of very long MC simulations on a large simulation cell comprising 54,000 monomeric sites assembled in 54 chains, as presented in the companion publication [78].
In addition, chain conformations in the initial amorphous and the final crystal phases are found to be identical in all large-scale features and differ only locally, at quite small length scales of very few bonds. This fact allows a calculation of the loss of conformational entropy per monomer upon crystallization, ∆s con f = −0.24k ± 0.04 . This decrease is less than the increase in translational entropy of the monomers, (∆s trans m = 1.17k per monomer), so crystallization still increases the overall free energy (entropy) in agreement with the observed spontaneous formation of crystals of polymeric chains. Thus, the calculated entropy of crystallization of freely jointed chains of hard spheres turns out to be ∆s m ≈ 0.93k.

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