Enumeration of Self-Avoiding Random Walks on Lattices as Model Chains in Polymer Crystals

1. Introduction

The term “soft matter” refers to a class of physical systems which includes polymers, colloids, granular media, surfactants, and gels. The common feature of all these diverse materials is that they consist of units whose size is significantly larger than the constituent atoms [1]. One of the main characteristics of soft matter is the existence of thermal fluctuations which are mainly manifested as Brownian motion of atoms, particles, and molecules. Due to this, the molecular shape and size are constantly changing even under conditions of equilibrium. Thus, to describe equilibrium local and global structure, a statistical approach is required. Additionally, soft matter systems are characterized by spontaneous self-assembly, self-organization and phase transitions. An important tool to aid in the statistical description of mesoscopic structure, but also of phenomena and processes relevant to soft matter, in general, and polymer science, in particular, is the concept of random walk (RW) and its variation in the form of self-avoiding random walk (SAW) [2].

The mathematical concept of a SAW corresponds to a trajectory of fixed step length that grows randomly on a lattice under the condition that it cannot visit the same point twice. In a strict mathematical sense, several salient aspects of SAW are still unknown or only partially resolved [2]. Still, this has not been an obstacle for its successful application in a very wide range of physical systems and topics in chemistry, physics, polymer science, material technology, mathematics, process optimization, computer science and biology [3,4,5,6,7,8,9,10,11,12]. Over the years significant advances have been made in the algorithms related to SAW identification and enumeration, further advancing its usage in diverse problems, especially ones related to synthetic and biological macromolecules. The excluded volume effect in polymer melts and solutions is in fact intimately related to the condition of self-avoidance in random walks [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].

Self-avoiding random walks and polygons (SAPs) have been studied extensively on 2-D and 3-D lattices [30] with specific examples being the honeycomb [31,32,33,34], square [35,36], triangular [37], simple cubic [35,38,39,40], body centered cubic (BCC) [41] and face centered cubic (FCC) [41,42] lattices. Due to their importance and general applicability significant algorithmic and theoretical advances have been made in the enumeration, characterization, identification and scaling behavior of SAWs, SAPs, and their variants [19,20,33,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62]. It is interesting to notice that while the simple cubic (SC), BCC and FCC crystals have been extensively studied in the literature, no such wealth of information exists for the hexagonal close packed (HCP) one. This could be attributed to a possible assumption that the HCP and FCC crystals are characterized by the same number of SAWs, independently of the step size [63], because both crystals have the same coordination number (n_coord = 12) and packing density (φ = $\pi \sqrt{2}/6$). As will be demonstrated in the continuation such claim is not valid: beyond a certain SAW length the number of SAWs compatible with the HCP crystal becomes marginally higher than the one for the FCC crystal, the difference increasing as the number of steps of the SAW grows.

In the past we employed a SAW-based analysis to enumerate the possible configurations of single-chain crystals and thus explain their thermodynamic stability in regular lattices in slits (tubes) [64] and plates [65]. These studies have been motivated by the spontaneous, entropy-driven crystallization of chains of tangent hard spheres as demonstrated in Monte Carlo (MC) simulations of dense packings under various conditions of spatial confinement [66,67,68]. Very recently, it was documented that starting from random (disordered) packings, freely-jointed chains of hard spheres in the bulk show a transition to the ordered state following Ostwald´s rule of stages [69]: initially a random hexagonal close packed (rHCP) morphology is formed of mixed HCP and FCC character [70,71]. Given enough observation time the rHCP ordered morphology is eventually succeeded by an almost defect-free FCC crystal [71]. In parallel modeling efforts, a wealth of distinct crystal polymorphs, including non-compact crystals like the body centered cubic (BCC) and simple hexagonal (HEX), has also been encountered in the (free)energy-driven self-assembly of freely-jointed chains of tangent monomers interacting with the square well potential under very dilute conditions [72].

With respect to the stability of the HCP and FCC crystals made of athermal polymer chains semi-analytical calculations were presented in [70], based on the separation of translational and conformational degrees of freedom. Resulting upper bound estimates demonstrate that the conformational entropy of an HCP crystal of hard-sphere polymers is higher than the one of the FCC crystal by a margin of 0.331 × 10^-5 k per monomer, where k is Boltzmann´s constant. However, this minute difference is significantly smaller than the free energy advantage of the FCC crystal in terms of translational entropy. The latter can be assumed to be very similar, if not identical, to the translational entropy of monomers disregarding the constraints imposed by chain connectivity. A widely accepted value for the entropic difference between FCC and HCP crystals for monomeric hard spheres corresponds to approximately 112 × 10^-5 k per monomer, has been established in the literature [73,74,75], depending further on packing density [76]. The semi-analytical predictions of [70] have been supported by extremely long molecular simulations based on Monte Carlo algorithms demonstrating FCC perfection in athermal polymer crystallization of entangled chains [66,77].

In the present contribution we further support the semi-analytical calculations of [70] by exhaustive enumeration of the self-avoiding random walks (SAWs) compatible with the inter-site geometry as encountered in the FCC and HCP crystals. We also analyze the SAW statistics, including the local geometry (bending and torsion angles) and the global size (distribution and mean of the end-to-end distance). Two additional crystals, the HEX and BCC, both having n_coord = 8, encountered in simulations of chains whose monomers interact with the square well potential [72], are considered in the current work.

The manuscript is organized as follows: Section 2 presents the model, the reference crystals and the method employed for the SAW enumeration. Section 3 hosts the results on the local and global properties of SAWs on the HCP, FCC, HEX and BCC lattices. The discussion of the results can be found in section 4. The manuscript is concluded with section 5 summarizing the main conclusions and listing future extensions.

2. Materials and Methods

As reported in the introduction the present work forms part of an on-going modeling study of (free)energy- [72,78] and entropy-driven [70,71,79,80] self-organization of idealized systems based on polymers interacting with the hard sphere or the square well potential [81]. Four different crystals have been considered in this work: the hexagonal close packed (HCP), face centered cubic (FCC), body centered cubic (BCC or BC) and holoedric 6/mmm (simple hexagonal, HEX), whose periodic structure and salient characteristics can be found in [77,82]. The HCP and FCC ones have been selected because they compete as emerging and resulting structures in the disorder-order transition of freely-jointed chains of hard spheres at sufficiently high packing densities [70,71,79,80]. The HEX and BCC crystals, along with the HCP and FCC crystals and the quasi-crystalline Frank-Kasper [83,84] phase, are encountered as final stable morphologies in the crystallization of clusters formed from chains interacting with the square well, attractive potential [72,78] at dilute conditions. Density-based [85,86] and geometric [72] arguments can accurately explain the dominance of non-compact crystals in specific ranges of the interaction potential in two and three dimensions.

In the present work, a polymer is represented as a linear chain of strictly tangent hard spheres. This chain is grown on the sites of a crystal, so its bonded geometry is defined by the linear architecture of the chain, the tangency condition of bonded atoms and the spatial constraints imposed by the crystal sites. As in our past works [64,65] “monomer” refers to each of the hard spheres that constitute the polymer chain, while “site” refers to the each of the sites of the ideal crystal. A single chain is grown on each reference lattice, with a monomer occupying a single site, and bonded monomers lying on adjacent sites, which are thus separated by a bond/step length. With respect to bond geometry, bending (θ) and torsion (ϕ) angles are formed by successive triplets and quartets of monomers along the chain backbone. These angles must respect the connectivity and geometry of each lattice. Accordingly, bending and torsion angles along the polymer chain are different for different crystal types, as can be seen in Table 1 (bending angles) and Table 2 (torsion angles). Chains on the HCP crystal show the largest variety, with 6 and 19 distinct bending and torsion angles, respectively, while SAWs on the BCC crystal have only 3 compatible bending angles.

Following our past works on SAWs under confinement [64,65] a SAW of N steps consists of the ordered sequence of sites ω^N(0), , ω^N(1), … , ω^N(N-1), ω^N(N) with ω^N(0) being the SAW origin. In contrast to SAWs on spatially restricted lattices [64,65], the ones studied here correspond to the unconstrained bulk case. Thus, any site can be selected as the origin of the SAW, ω^N(0), due to the maximal symmetry of the bulk lattice. Self-avoidance of the random walk (or equivalently excluded volume of polymer) is fulfilled as ω^N(i) ≠ ω^N(j) for all i ≠ j. Bond tangency is satisfied as |ω^N(i+1) - ω^N(i)| = 1, with i $\in$ [0,1,….N-1] and |ω|= (ω ^. ω)^(1/2) denoting the Euclidean norm. The size of the SAW is quantified by the square end-to-end distance: |ω^N|² = (ω^N(N) - ω^N(0)) ^. (ω^N(N) - ω^N(0)).

For a given number of steps, N, and a reference crystal X (HCP, FCC, HEX or BCC) our algorithm is based on the direct enumeration of the total number of distinct SAWs, c_N(X), and proceeds with the calculation of all bending and torsion angles and of the square end-to-end distance of each identified SAW. Accordingly, for a SAW of N steps the mean square end-to-end distance, <|ω^N|² >(X), can be readily calculated as:

$$〈{\left|{\mathbf{w}}^N\right|}^2〉(X)=\frac{1}{c_N}{\sum }_{c_N}{\left|{\mathbf{w}}^N\right|}^2$$

(1)

The scaling of the number of distinct SAWs, c_N, and of the average SAW size, $〈{\left|{\mathbf{w}}^N\right|}^2〉$, as a function of the number of SAW steps is given by exponential-power-law asymptotic expressions [2,36,38,39,40,60,87]:

$$c_N\left(X\right)~A{\mu }^N{N}^{\gamma -1}$$

(2)

$$〈{\left|{\mathbf{w}}^N\right|}^2〉\left(X\right)~D{N}^{2v}$$

(3)

where A and D are the critical amplitudes, μ is the connective constant, and γ and ν are the critical exponents. The critical exponents are considered universal, while the values of A, D and μ are lattice dependent.

SAW enumeration and the calculation of their statistics over the four different reference crystals has been performed with the help of a home-made code. All reported numerical calculations have been executed on an Intel i9-10850K with 16 Gb of memory, running on Linux operating system. Starting from the work of Orr [22] elegant methods and efficient algorithms have been developed over the years to tackle the exponentially difficult SAW enumeration problem [6,9,36,38,39,61,62,88,89] allowing to reach various high-N SAWs in distinct lattices. We should further note that our SAW enumeration algorithm is not as efficient as the state-of-the-art methods described above, and thus our analysis is limited to SAWs of intermediate number of steps. However, this should not be considered as a potential disadvantage as extensive off-lattice simulations, under a wide variety of conditions, have clearly demonstrated that, beyond a moderate value, chain length has practically no effect on the crystallization of athermal packings of fully flexible chains. The phase behavior and the established ordered morphologies are the same independently of the average length of chains being as low as 12 or as high as 1000 monomers, the former value corresponding to short oligomers while the latter to well entangled chains, deep in the polymeric regime [67,79,80].

3. Results

Table 3 presents the total number of distinct SAWs, c_N, and the mean square end-to-end distance, $〈{\left|{\mathbf{w}}^N\right|}^2〉$, as a function of SAW steps, N, for the FCC and HCP lattices, along with the difference in the number of SAWs, Δc_N (= c_N(HCP) – c_N(FCC)). The corresponding data for the BCC and HEX crystals are presented in Table 4.

The FCC and HCP crystals show the same number of SAWs up to N = 5 steps. However, for N ≥ 6 the number of SAWs on the two lattices start to deviate with the HCP crystal showing systematically more SAWs than the FCC one. At N = 6, Δc_N = 12, the relative difference is approximately 7.6×10^-6. For the longest SAW studied here (N = 12) Δc_N = 60,617,100 which corresponds to a relative difference of 3.4×10^-5, showing increasing trends as the SAW grows in length. In parallel, as can be seen by the comparison of the related rows in Table 3 the average size of the SAWs in the HCP and FCC lattices is the same within a tolerance of 10^-5 for all values of N.

Similar trends are observed in the comparison of the BCC and HEX crystals as listed in Table 4. For N = 3 the BCC crystal has 12 more SAWs compared to the HEX crystal with the relative difference being 0.031, significantly higher than the one observed for the HCP-FCC pair at the corresponding deviation point (N = 6). For the longest SAW studied here (N = 14) the absolute and relative differences for the HEX-BCC pair increase to 89,206,013,508 and 0.20, respectively. Accordingly, even if both lattices have the same coordination number (n_coord = 8) the single chain conformations, as quantified by the number of SAWs, are significantly fewer in HEX compared to BCC. With respect to average size the polymer grown on the HEX lattice is systematically longer than the one on the BCC with the relative difference being approximately 0.053 for N = 14.

The left panel of Figure 1 shows the logarithm of the total number of distinct SAWs, c_N, as a function of the logarithm of the total number of SAW steps, for all lattices studied here. Also shown are the differences of the two pairs (c_N(HCP) – c_N(FCC)) and c_N(BCC) – c_N(HEX)). On the right panel of Figure 1 we can observe the dependence of the logarithm of the mean square end-to-end distance, <|ω^N|² >(X), as calculated from Eq. 1, on the logarithm of the number of SAW steps, N. Furthermore, we have fitted all available data with best linear fits corresponding to the scaling formula of Eq. 3. Figure 1 shows one such best linear fit corresponding to the data for the HEX crystal, with the reliability fitting coefficient being practically equal to 1. The complete set of the parameter values, according to Eqs. 2 and 3, as obtained from best linear fits on current SAW data are reported in Table 5. Very little variation is observed for both the critical amplitude and exponent between the different lattices, especially when the FCC and HCP crystals are compared. A minimal trend suggests that the exponent increases slightly for the crystals of the lower coordination number, while the opposite behavior is observed for the amplitude. Given the short- or intermediate length of the studied SAWs the universal exponent of 0.588 is nicely matched for the HCP and FCC crystals. Between the HCP and the FCC crystal there is no appreciable difference suggesting that the corresponding SAW sizes will be very similar, but strictly not identical, for the limit of infinitely long chains (N → ∞).

The distribution of the discrete bending angles which are compatible with each lattice studied here is given in Figure 2 as a function of the number of SAW steps. For all crystal types as N grows the population of obtuse angles experiences small decreases in favor of the acute ones because of the self-avoidance condition. The fraction of bending angles with 90^o, where available (HCP, FCC and HEX), rapidly reaches a stable plateau. Extrapolating the current trends to longer N the most probable bending angle is 60º, 60º, 70.53º and 90^o for the FCC, HCP, BCC and HEX crystals, respectively.

Figure 3 hosts the corresponding results for the distribution of discrete torsion angles which are compatible with each lattice crystal. For FCC as N increases the populations at 54.7, 125.3 and 180^o increase, while the one at 90^o shows a significant reduction. Compared to other crystals the SAWs on the HCP lattice show a richer behavior with the primary angles resting at 54.7, 70.5, 109.5 and 180^o. The HEX crystal is the one that shows the smallest variation with increasing number of SAW steps.

The distribution of SAW size, as quantified by the square end-to-end distance, |ω^N|² is shown in Figure 4 and Figure 5 for the HCP – FCC and BCC – HEX pairs, respectively. The probability distribution function (PDF) is presented in the main figure while the cumulative distribution function (CDF) is shown in the inset. Given the major differences in the bonded geometry between HCP and FCC lattices it is not surprising that the distribution of the SAW size shows significant deviations. However, it is interesting to notice that in spite of these variations the global SAW size, on average, is indistinguishable between the close packed FCC and HCP crystals as shown in Table 1 and further confirmed by the CDF trends in the inset of Figure 4. On the other hand, SAWs on the BCC lattice have smaller size than the ones on the HEX as indicated by the concentrated maxima of the BCC distribution at small values of the square end-to-end distance.

4. Discussion

We enumerate and describe the self-avoiding random walks, SAWs, grown on different crystals, characterized by the same coordination number (n_coord = 12 or 8). We observe that for the HCP and FCC pair once a critical number of steps is reached (N = 6) the distinct number of SAWs becomes different. In relative numbers this difference is minute but definitely non-zero with the HCP crystal showing a larger number of SAWs compared to FCC (c_N(HCP) > c_N(FCC)). The present finding is clear evidence that the SAW behavior in the HCP and FCC crystals should not be considered as identical even though both are characterized by the same coordination number and packing density. In parallel, the average size of SAWs on these two crystals is indistinguishable for the whole range of studied SAW steps even if the internal chain/SAW bonded geometry, as quantified by the bending and torsion angles, is significantly different with the HCP crystal showing a much richer variety and same is true for the size distribution as quantified by the end-to-end distance.

The main objective of the present work is to provide a quantitative basis for the study of the thermodynamic stability of the HCP and FCC crystals made of fully flexible, athermal polymers. As mentioned in the introduction by assuming independence of the translational and conformational degrees of freedom, an argument which is analyzed in detail in [70], the total entropy of the crystal can be considered as the summation of two distinct contributions: the conformational one, ΔS_conf, dictated by chain connectivity respecting the geometric constraints of the reference crystal, and the translational one, ΔS_tran, which should be very similar to the one of monomeric systems of hard spheres, free of any constraints imposed by chain connectivity. Accordingly, the entropy difference of the HCP and FCC crystals can be written as [70]:

$$∆S^{FCC-HCP}=∆S_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{f}}^{FCC-HCP}+∆S_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}^{FCC-HCP}=\mathrm{l}\mathrm{n}\frac{\left|{\Xi }^{FCC}\right|}{\left|{\Xi }^{HCP}\right|}+∆S_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}}^{FCC-HCP}$$

(4)

where Ξ is the complete set (partition function) of all multichain configurations compatible with the polymer model (here freely-jointed chains of tangent hard spheres) and the lattice geometry (here HCP or FCC), and || denote the cardinality of a set [70]. The translational contribution is known from past studies on monomeric analogs [73,74,75]. Thus, the evaluation of Ξ for the multichain configurations is required to arrive at an accurate estimate for $∆S^{FCC-HCP}$ but such information is not available for any of the crystals. Accordingly, we need to establish an upper-bound, single-chain estimation. The concept of random walks (RWs) would not shed any light given that both crystals have the same coordination number and as such the same number of RWs over the whole range of steps. Toward this end, as a more refined criterion we resort to self-avoiding random walks (SAWs), providing a tighter and more discriminating upper bound estimate for the difference in the configurational entropy as [70]:

$$\frac{\left|{\Xi }^{HCP}\right|}{\left|{\Xi }^{FCC}\right|}<{\left(\frac{c_N\left(HCP\right)}{c_N\left(FCC\right)}\right)}^N\underset{}{\iff }\mathrm{l}\mathrm{n}\left(\frac{\left|{\Xi }^{HCP}\right|}{\left|{\Xi }^{FCC}\right|}\right)=\frac{1}{N}ln\left(\frac{c_N\left(HCP\right)}{c_N\left(FCC\right)}\right)$$

(5)

The data presented here in Table 3 and Table 4, and Figure 1 for the HCP-FCC and HEX-BCC pairs clearly identify a difference in the number of SAWs once a critical number of steps is reached. The ratio of the SAWs in Eq. 5 can be fitted using an asymptotic formula:

$$\frac{\frac{c_N\left(X\right)}{c_N\left(Y\right)}-1}{N}=A-B{e}^{-dN}$$

(6)

where X and Y are the two crystals to be compared and A, B and d are the fitting parameters. Fittings using the exponential formula of Eq. 6 on available SAW enumeration data are reported in Figure 6.

Both pairs suggest an asymptotic behavior of the entropy difference per monomer. For the HCP and FCC crystals, the established behavior is far from the asymptotic regime and significantly longer SAWs are required pointedly increasing the computational complexity of the problem. However, application of the exponential formula (Eq. 6) allows the prediction of the ratio for very long SAWs from short and intermediate values as the ones presented here. The fitting parameters for the HCP/FCC and BCC/HEX are summarized in Table 6. Although the maximum chain length for which we could exhaustively enumerate SAWs is moderate, the quality of the fit in Figure 6 suggests that the values of the parameter A in Table 6 are a quite accurate approximation to the asymptotic value, i.e. for infinite chain length, which is the physically relevant quantity for long polymers.

Concentrating on the HCP-FCC pair, which is the main focus of the present study, placing the asymptotic formula in Eq. 5 and considering the limit of very long chains provides an upper bound estimate, $∆S_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{f}}^{FCC-HCP}\approx -0.331\times {10}^{-5}k$, as the HCP polymer crystal has a higher conformational entropy than the FCC one. However, this value is significantly lower than the translational entropy $∆S_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{FCC-HCP}\approx 112\times {10}^{-5}k$. Accordingly, the FCC is the most stable crystal among the polymorphs for freely-jointed polymers of tangent hard spheres.

5. Conclusions

Prediction of the thermodynamic stability of crystals made of athermal polymers is significantly more complicated than for monomeric analogs due to the presence of constraints imposed by chain connectivity. Here, we demonstrate that once a critical number of steps is reached the HCP crystal has more self-avoiding random walks than the FCC one. Accordingly, a tight upper bound estimate of the conformational entropic advantage of the HCP crystal suggests that this lead is not sufficient to overcome the significantly larger translational advantage of the FCC crystal [70]. Accordingly, the FCC crystal is the thermodynamically most stable crystal and should prevail among competing polymorphs, as confirmed by recent off-lattice simulations [71]. The present work is currently being extended to tackle freely-rotating chains on crystals in two and three dimensions.