Effect of lattice inhomogeneity on collapsed phases of semi-stiff ISAW polymers

The canonical picture of polymer collapse is that long chain polymers in a poor solvent undergo a phase transition from an extended phase to a collapsed phase at a critical temperature, known as the $\theta$-point Flory (1953); Duplantier and Saleur (1987). In the interacting self-avoiding walk (ISAW) model this is induced by an interaction between non-consecutive monomers causing the polymer to energetically favour a random globule configuration at low temperatures. There are, however, a number of other lattice models of polymer collapse with different properties which may suggest the $\theta$-point transition is not universal Blote and Nienhuis (1989); Shapir and Oono (1984). Often these models involve a change to the self-avoiding condition of the lattice walk such as the vertex-interacting SAW Foster and Pinettes (2003), or the self avoiding trail Owczarek and Prellberg (1995, 2007) but can also be due to the addition of another interaction Doukas et al (2010); Bedini et al (2016). In the latter case a new collapsed phase can appear where the configurations of the walk have crystalline structure. A similar example of this kind of model is the semi-stiff ISAW Bastolla and Grassberger (1997); Krawczyk et al (2009) which allows for varying the energy cost of a bend in a walk. When bends are inhibited the typical configurations resemble $\beta$-sheet configurations and the transition between the extended phase and the crystal phase is first order, compared to the continuous $\theta$-point transition to the random globule phase in the flexible regime.

If the collapsed phase is characterised by the formation of dense configurations then it is sensible to ask how robust this phase is to some degree of inhomogeneity in the underlying medium. The effect of random disorder on otherwise well-understood statistical mechanical problems has been a topic of study for some time, going back to the Ising model Watson (1969). The study of the SAW model on randomly diluted lattices goes back almost as far and has been closely related to the problem of percolation. By modeling disorder in the underlying medium as percolation early studies showed that the fundamental scaling laws for SAWs, that of the size and number of walks, have the same universal exponents $\nu$ and $\gamma$ on inhomogeneous lattices as for homogeneous lattices, provided the disorder is above the percolation limit $p_{c}$ Kremer (1981). Change in this scaling behaviour only occurs at the percolation limit $p_{c}$ Blavatska and Janke (2008). These results have been confirmed with numerical work Lee and Nakanishi (1988); Rintoul et al (1994) and exact enumeration Lam (1990); Ordemann et al (2000); Nakanishi and Lee (1991). Renormalisation group theory also finds that Kim (1983); Blavatska et al (2005) however there is evidence of additional fixed points due to the disordered medium Meir and Harris (1989); Grassberger (1993). The addition of disorder also introduces new considerations such as how the type of averaging over disorder affects SAWs Nakanishi and Moon (1992); Birkner and Sun (2010) and when scaling laws are well-defined Janssen and Stenull (2007). Other approaches to linear polymers in disordered media include fractal and irregular lattices Chakrabarti (2005); Dhar and Singh (2005); Rammal et al (1984); Roy and Blumen (1990).

The link back to polymer collapse involves the formation of sub-clusters in a SAW on a smaller scale which can be related to cluster formation in percolation models Brak et al (1998) and different scaling if one considers SAWs on a single percolation cluster or across multiple percolation clusters Sahimi (1984). It is also known that the (effective) connective constant of SAWs on a regular lattice decreases linearly with bond dilution $1-p$ but the $\theta$-point is only weakly affected, even down to critical $p_{c}$ Barat et al (1991). In this paper we use the percolation formulation as the representation of lattice inhomogeneity, where each site has a probability $p$ to be open. From the point of view of a lattice polymer, each site has a probability $1-p$ to be an invalid step in the walk, thus serving as an impurity. We are more interested in this work on the effect of the inhomogeneity on the collapsed phase of the semi-stiff ISAW model so we work well above the percolation limit $p_{c}\approx 0.593$ Newman and Ziff (2000) on the square lattice where scaling laws are unchanged. Away from the fixed points at $1-p=0,1-p_{c}$ the disruption to the symmetry of the lattice should have a significant effect on the highly ordered crystal phase. We characterise the phases present in this system and see how they are affected by an increase in lattice inhomogeneity. In particular, we investigate whether the crystal phase is more sensitive to disruption by impurities compared to the globule phase.

We consider single polymers in dilute solution modeled as self-avoiding walks (SAWs) on the square lattice. The canonical partition function for walks of length $n$ with $m$ bulk interactions and $s$ straight segments is

where $\omega$ and $\eta$ are the Boltzmann weights of a non-consecutive neighbouring pair and a straight segment, respectively, and $c_{n}(m,s)$ is the density of states. To represent the presence of impurities, each lattice site has a probability $p$ that it is a valid site for a step in the walk. Thus the inhomogeneous lattice is represented by making a fraction $1-p$ of lattice sites unavailable. An example walk is shown in Fig. 1. Details of how the lattice configuration is chosen are given below when discussing the flatPERM algorithm.

To characterise the phases of the system, particularly the collapsed phases, we calculate the internal energy-type order parameters, namely the average fraction of straight segments $\langle s\rangle/n$ and the average fraction of neighbouring steps $\langle m\rangle/n$. An additional signature of the crystal phase is the anisotropy, defined as

where $n_{x}$ and $n_{y}$ are the number of steps of the walk in each orthogonal direction of the lattice. In the case of a homogeneous lattice the collapsed phases are indicated by a non-zero value of $\langle\zeta\rangle$ with $\langle\zeta\rangle=1$ in the crystal phase. Within the collapsed regime we expect $\langle s\rangle/n$ to change value across the globule-crystal transition while $\langle m\rangle/n$ remains smooth and varies slowly. The effect of lattice inhomogeneity on the system will be seen in how these quantities vary with $1-p$.

The effect of disorder on the nature of the phase transitions determined from scaling laws. In particular, in the thermodynamic limit of long chains the specific heat near the transition temperature $T_{c}$ scales according to

controlled by the critical exponent $\alpha$. If $\alpha<1$ the transition is continuous; if $\alpha=1$ then Eq. (3) diverges, indicating a first-order transition. For a finite system we calculate the reduced specific heat from moments of the partition function, for example

and similar for $c_{n}^{(s)}(\eta)=\text{var}(s)/n$. The peak of $c_{n}^{(m)}$ is at the transition point $\omega_{c}$ and the value of the peak scales as $n^{\alpha\phi}$, where $\phi$ is the crossover exponent. By considering Eq. (4) and using the standard relation $2-\alpha=1/\phi$ we may determine the nature of the transitions. The rule of thumb is that if the peak values of $c_{n}^{(m)}/n$ coincide for a range of $n$ then $\alpha\phi=1$, thus $\alpha=1$ and the transition is first-order; else, continuous. Since we have two microcanonical parameters we generalise and calculate the Hessian covariance matrix

where $f_{n}=-\tfrac{1}{n}\log Z_{n}$ is the reduced free energy.

Walks are simulated using the flatPERM algorithm Prellberg and Krawczyk (2004), an extension of the pruned and enriched Rosenbluth method (PERM) Grassberger (1997). The algorithm works by growing a walk on a given lattice up to some maximum length $N_{\text{max}}$. At each step the number of bulk interactions $m$ and straight segments $s$ are calculated and the cumulative Rosenbluth & Rosenbluth weight Rosenbluth and Rosenbluth (1955) is compared with the current estimate of the weights of all samples $W_{n,m,s}$. If the current state has relatively low weight the walk is ‘pruned’ back to an earlier state. On the other hand, if the current state has relatively high weight, then microcanonical quantities $m$ and $s$ are measured and $W_{n,m,s}$ is updated. The state is then ‘enriched’ by branching the simulation into several possible further paths (which are explored when the current path is eventually pruned back). When all branches are pruned a new iteration is started from the origin. FlatPERM improves upon PERM by altering the prune or enrich choice such that the sample histogram is flat in the microcanonical parameters $n$, $m$ and $s$. Further improvements are made to account for the correlation between branches that are grown from the same enrichment point, which provides an estimate of the number of effectively independent samples. The main output of the simulation are the weights $W_{n,m,s}$, which are an approximation to the athermal density of states $c_{n}(m,s)$ in Eq. (1), for all $n\leq N_{\text{max}}$. In practice, thermodynamic quantities are determined by specifying $\omega$ and $\eta$ and using the weighted sum

For each case or simulation mentioned in this work the results are comprised of a composite of ten independent simulations in order to obtain some measure of statistical error.

The inhomogeneous lattice is implemented by choosing a set of lattice sites to be inaccessible to the walk. The number of impurities is drawn from the appropriate binomial distribution with $1-p$ being the probability of any particular site being a valid site for the walk. These impurities are distributed uniformly over the area of the lattice that would be accessible to a walk of length $n$ in the case of a homogeneous lattice. The set of inaccessible sites is reseeded at the beginning of each flatPERM iteration (growing the walk from the origin). The initial weight of each iteration is set to be the probability of the configuration of lattice impurities. In this way the output weights $W_{n,m,s}$ contain the sum over disorder such that any $\langle Q\rangle$ in Eq. (6) also represents a quenched-type average over disorder Nakanishi and Moon (1992). It is also possible to achieve annealed-type averages by ignoring the probability of the lattice configuration at each iteration but we found no appreciable difference between the two methods for this problem.

First we characterise the phases by looking at the order parameters for the system and the expected configurations, with and without lattice impurities. In Fig. 2 we plot, from left to right, the average fraction of straight segments $\langle s\rangle/n$, the average fraction of neighbouring steps $\langle m\rangle/n$ and the average anisotropy $\langle\zeta\rangle$. Further visualisation of the phases is given in Fig. 3 which shows typical configurations at points in the $(\omega,\eta)$ phase diagram that are indicative of each phase. These points are marked with black dots on Fig. 2(a). In both figures the top and bottom row are for $1-p=0,0.25$, respectively. In the homogenous lattice case, we see the three distinct phases of the semi-stiff ISAW system Krawczyk et al (2009). For small $\omega$ at all $\eta$ the chains are extended, with near-zero fraction of neighbouring monomers arranged isotropically and the fraction of straight segments increasing weakly with $\eta$. Configurations in this phase look like Fig. 3(a,b). For small $\eta$ and $\omega>\omega_{\text{c}}$ the chains are collapsed in a random (and thus isotropic) globule, with a small fraction of straight segments, like in Fig. 2(c). At $\eta=1$ the chains are completely flexible and there a collapse transition at $\omega=\omega_{\theta}$ so that for large $\omega>\omega_{c}$ the polymer is in an amorphous dense state where the average size scales as $\sqrt{n}$. The third phase is the crystal phase which is also collapsed since where the average size scales as $\sqrt{n}$ which is characterised by folded configurations of zero entropy. This phase occurs when $\eta$ is larger than some critical value $\eta_{c}$ such that the fraction of straight segments is high, and thus the chains have anisotropic conformations arranged as parallel lines, like in Fig. 2(d), analogous to $\beta$-sheets in three dimensions.

In the case of the inhomogeneous lattice, where a considerable fraction of the lattice is unavailable to the walks, there are several differences. The extended phase ($\omega<\omega_{c}$) is still clear in this case, marked by small but non-zero $\langle m\rangle/n$, $\langle\zeta\rangle\approx 0$ and $\langle s\rangle/n$ increasing weakly with $\eta$, as in the homogeneous case. The configurations in Fig. 3(e,f) are like those in Fig. 3(a,b). A collapsed phase still exists for $\omega>\omega_{c}$, as evidenced in Fig. 2(d) where the value of $\langle m\rangle/n$ is non-zero, albeit smaller than the homogeneous lattice case.the globule and crystal phases. This is noticeable in Fig. 2(e) where the value of $\langle s\rangle/n$ is reduced for the inhomogeneous lattice in the high stiffness and nearest neighbour interaction region ($\omega>\omega_{c}$ and $\eta>\eta_{c}$). Moreover, $\langle s\rangle/n$ increases only slowly with increasing $\eta$ and is largely independent of $\omega$ in this region. This suggests a smooth crossover as $\eta$ increases for $\omega>\omega_{c}$, rather than a phase transition as is the case for the homogeneous lattice. Similarly, $\langle\zeta\rangle$ is greatly reduced in this region, however it remains non-zero, as shown in Fig. 2(f). The corresponding configurations in Fig. 3(g,h) - for the globule, crystal phase, respectively - show a change in the global structure (compared to Fig. 3(c,d)), with the appearance of several connected smaller clusters. In particular, the straightness of the connecting strands between clusters accounts for the small but non-zero value of $\langle\zeta\rangle$.

Another quantity of interest for polymers is the mean-squared end-to-end distance $\langle R_{e}^{2}\rangle$, where for a single chain $R_{e}^{2}=\mathbf{r}_{n}.\mathbf{r}_{n}$ assuming the chain starts at the origin $\mathbf{r}_{0}=(0,0)$. In Fig. 4 we plot $\langle R_{e}^{2}\rangle$ for two regions of the phase diagram, $\omega=1$ and $\omega=4$ corresponding to the extended phase and the collapsed phases, respectively, as well as with and without lattice inhomogeneity, $1-p=0$ and $1-p=0.25$, respectively. The top two series (blue and red) are for the extended phase and show only a small difference between the homogeneous and inhomogeneous lattice. This small difference suggests a slight swelling of the chain for small $\eta$ when inhomogeneity reduces the number of available sites for the chain. Interestingly, this difference reverses at larger $\eta$ where samples of very long straight chains are inhibited by the inhomogeneous lattice. There is a larger distinction between the homogeneous and inhomogeneous lattice cases in the collapsed phases (green and orange). On the inhomogeneous lattice the chains are significantly larger than on on homogeneous lattice though presumably still obeying $\sqrt{\langle R_{e}^{2}\rangle}\sim\sqrt{n}$; our data is not sufficient to verify this scaling. However, apart from a small effect due to the homogeneous lattice globule-crystal transition near $\eta\approx 2$, the size of the chains are not significantly affected by stiffness in the collapsed phases for both the homogeneous and inhomogeneous lattices.

The effect of lattice inhomogeneity on the phase boundaries can be seen more clearly in Fig. 5 which plots the largest eigenvalue of the covariance matrix $H_{n}$ of Eq. (5). For $p=0$ (left) the phase boundaries are clear, especially between the extended and crystal phases. The value of $\omega_{c}$ that marks this boundary decreases slightly as $\eta$ increases, but corresponds to the $\theta$-point at $\omega_{\theta}\approx 1.94$ for SAWs on the square lattice Meirovitch et al (1989). In contrast, lattice impurities ($1-p=0.25$, right) soften the boundaries. While the variance that marks the boundary between the extended and collapsed phases is smaller it is still clear. To determine the nature of these transitions we look at slices of Fig. 5 at fixed values of $\omega$ or $\eta$. Away from the multicritical point, off-diagonal terms of the Hessian matrix are negligible and so the largest eigenvalue is equivalent to the largest specific heat. Thus, in Fig. 6 we plot the specific heats scaled by $n$ to determine the scaling exponent $\alpha$. The (scaled) variance of $m$ across the extended-crystal transition is shown in Fig. 6(a) and (c) at fixed $\eta=5$ for $1-p=0$ and $1-p=0.25$, respectively. On the homogeneous lattice (a), the peaks of $\text{var(m)}/n^{2}$ coincide, indicating $\alpha\approx 1$, consistent with a first-order transition. In the presence of disorder (c), the scaling indicates $\alpha<1$, and thus the transition is now continuous. This is a clear example of quenched disorder causing the rounding of a first-order transition in 2D Hui and Berker (1989); Aizenman and Wehr (1989). For the globule-crystal transition at fixed $\omega=3$ we show the variance of $s$ scaled by $n$, again for $1-p=0$ and $1-p=0.25$, in Fig. 6(b) and (d), respectively. For $1-p=0$ (b) the variance of $s$ indicating a globule-crystal transition is broader and smaller than the extended-crystal transition, but nevertheless it is clear enough to indicate a continuous transition, with $\alpha<1$. However, in the presence of disorder (d), there is no clear signal of a transition at all. Along with the effect on the order parameters and $\langle R_{e}^{2}\rangle$, we can confidently argue that that the distinction between the globule and crystal phase regions no longer exists for SAWs on an inhomogeneous lattice. We conclude that the crystal phase only truly exists for $p=1$.

So far we have compared the homogeneous lattice case to a case with a significant amount of inhomogeneity. While the value $1-p=0.25$ is not close to the percolation limit $1-p_{c}$ it still makes a significant fraction of the lattice unavailable to the walk. We now turn our attention to the scaling region around $1-p\to 0$ and $n$ large. For a better look at this scaling region we also ran a set of additional simulations with a fixed interaction weight $\omega=3$ corresponding to a slice in the phase diagram that avoids the extended phase and looks for and change in the collapsed phase as $\eta$ is increased. This value of $\omega$ is large enough to avoid the collapse transition without being overwhelmingly biased towards collapsed configurations. This is achieved by altering the weights of the samples by a factor $\omega^{m}$ so that the sum over $m$ in Eq. (6) is performed within the simulation. This method produces histograms over $n$ and $s$ only and can sample longer walks up to $n=1024$. Figure 7 shows (a) average fraction of straight segments $\langle s\rangle/n$, (b) average fraction of neighbouring steps $\langle m\rangle/n$ and (c) average anisotropy $\langle\zeta\rangle$. As expected, increasing $\eta$ is marked by a change in the behaviour of $\langle s\rangle/n$ and $\langle\zeta\rangle$, but not $\langle m\rangle/n$. This is broadly true for all values of $1-p$ but other effects due to increasing $1-p$ are also evident. For $1-p$ near or equal to zero $\langle s\rangle/n$ has an approximately linear $\eta$ dependence at small values of $\eta$ and $\eta\approx 2.5$ $\langle s\rangle/n$ tends to a constant value. This is the continuous globule-crystal transition seen earlier and we see that it softens considerably for larger $1-p$ where there is no crystal phase. At small $1-p$ the limiting value at large $\eta$ is close to unity and it decreases slightly as $1-p$ increases. This indicates that although the available fraction of the lattice may diminish it is still possible for configurations to have large straight segments and $\langle s\rangle/n$ does not necessarily say anything about the mean size of the chains; compare the configurations in Fig. 3(d,h), which have similar values of $s/n$. In Fig. 7(b) we also see that $\langle m\rangle/n$ decreases smoothly with increasing $1-p$ at all values of $\eta$. However, the effect is stronger compared to $\langle s\rangle/n$ since each impurity prevents multiple near-neighbour bonds, presumably as a function of the coordination number of the lattice. This corresponds partly to a swelling in the size of the chain but also to the formation of smaller clusters, as seen in Fig. 3(g). Although we know that the crystal-globule transition has disappeared, it is not clear from the plots of $\langle s\rangle/n$ and $\langle m\rangle/n$ whether there is a significant change in the properties of the collapsed phase in the large $\eta$ region as the disorder increases.

The more interesting behaviour is in the anisotropy $\langle\zeta\rangle$ in Fig. 7(c). For $1-p\approx 0$ there is a steep increase from the isotropic globule state $\langle\zeta\rangle\approx 0$ to the crystal phase $\langle\zeta\rangle\approx 1$. As $1-p$ increases the sharp change near the transition disappears and the limiting value for large $\eta$ also reduces significantly. In order to focus on the effect of impurity density on the anisotropy in the collapsed phase we look more closely at the low-temperature (high $\eta$) limiting values of $\langle\zeta\rangle$. Because the anisotropy Eq. (2) is a function of the global structure of each walk we should also account for the finite size of the walks by comparing the relevant length-scales, namely the size of the chains and the mean separation of impurities. As such we define a scaling parameter $\chi=\sqrt{n\,(1-p)}$, which follows from the mean-squared radius quantities that scale as $\langle R_{e}^{2}\rangle\sim n^{2\nu}$, where $\nu=1/2$ for the collapsed phases in two dimensions. Fig. 8 shows the values of $\langle\zeta\rangle$ at a large value of $\eta$ representing the low-temperature limit, namely $\eta=5$, as a function of $\chi$. The values of $\chi$ are achieved by taking data from a several values of $n$ for each value of $1-p$. As seen before, for small $\chi$ the limiting value of $\langle\zeta\rangle$ is close to unity indicating the crystal phase. At large values of $\chi$ the limiting value of $\langle\zeta\rangle$ is reduced, but this appears to be more a function of $1-p$ only rather than $n$ via $\chi$. This is the same behaviour as $\langle s\rangle/n$ and $\langle m\rangle/n$. However, there is an intermediate regime $5\lesssim\chi\lesssim 8$ where the limiting value of $\langle\zeta\rangle$ decreases sharply. This indicates that at some critical amount of inhomogeneity the anisotropic configurations which typify the crystal phase are inhibited. The only remaining distinction from the globule phase is that $\langle\zeta\rangle>0$ for large $\chi$, or at least it tends to zero slowly, reflecting that the smaller sub-clusters retain some anisotropic order, but this is uncorrelated with other sub-clusters on the scale of the entire walk, as per Fig. 3(h).

Taken together, the behaviour of the order parameters indicate that as disorder, measured by $\chi$, increases the crystal phase disappears. Within the remaining collapsed phase there is still a dependence on $\eta$, i.e. stiffness. At one end, $\eta=1$, the system is most like the globule phase from the homogeneous lattice case where only $\langle m\rangle/n$ decreases as $1-p$ increases, representing the formation of clusters at a smaller scale. As $\eta$ is increased to where the crystal phase was, the overall anisotropy remains small but $\langle s\rangle/n$ increases, albeit to a smaller value than the homogeneous lattice case. This indicates that there is an increase of stiffness within each smaller cluster but with no correlation between the smaller clusters since the impurities are distributed uniformly. Furthermore, since there is no sharp change in any of the order parameters at some intermediate $\eta_{c}$ then the transition between these regimes is smooth and not critical. Thus, the globule-crystal transition has disappeared and any remnant indications of a transition that appear in Fig. 5(b) are therefore argued to be due to finite-size effects, such as conformational rearrangements, rather than due to the finite-size rounding of a thermodynamic transition.

We have simulated semi-stiff ISAWs on the square lattice with variable amount of lattice inhomogeneity represented as impurities, or sites that are unavailable to the walks. We have characterised the three expected phases, extended, globule and crystal, and looked at how the latter two collapsed phases are affected by increasing the inhomogeneity. When a large number of impurities are present the walks prefer a configuration of several smaller clusters as opposed to a global swelling. This is consistent with other studies of long chain polymers in disordered media Goldschmidt and Shiferaw (2003). For finite length walks the effect on the crystal phase is strongest, where we observe changes to the order parameters. In particular, there is a significant decrease in the average anisotropy as well as a softening of the crystal-globule transition. The latter is seen as a rounding of the singular behaviour in $\langle s\rangle/n$ to the point where it is no longer clearly distinguishable as a critical transition. We conclude that the crystal phase and the critical transition between the crystal phase and globule phase disappear when the mean distance between impurities is comparable to the size of the chains, as judged by the scaling of mean-squared radius $\langle R_{e}^{2}\rangle$. We identify that the ratio of these length scales in the system $\chi$, or rather $\chi^{2}=n\,(1-p)$, serves as the appropriate scaling variable near the onset of disorder in the lattice, i.e. $1-p>0$. Hence in the thermodynamic limit there is no crystal phase for any amount of inhomogeneity. We also note that the scaling variable depends on the length-scale exponent $\nu=1/2$ and thus the softening of the transition is consistent with the Harris criterion for the effect of quenched disorder on critical phenomena Harris (1974).

Beyond the exploration of the phase diagram the natural question is to look at the scaling behaviour of the order parameters, fluctuations of these and size measures in each of the phases and also at the critical transitions between them. While the machinery of finite-size scaling has proved to be powerful for other SAW problems, in this case the data gathered is not good enough to be used for such analysis. In particular we are unable to locate critical points with sufficient precision. The addition of disorder will always present a challenge to computational methods but in this case the stronger effect is due to the challenge of simulating dense configurations with enough fidelity to distinguish between different behaviours in the collapsed phase.