Exact Partition Functions of Interacting Self-Avoiding Walks on Lattices

Ideas and methods of statistical physics have been shown to be useful for understanding some interesting problems in physical systems, e.g. universality and scaling in critical systems. The interacting self-avoiding walk (ISAW) on a lattice is the simplest model for homopolymers and serves as the framework of simple models for biopolymers, such as DNA, RNA, and protein, which are important components in complex systems in biology. In this paper, we briefly review our recent work on exact partition functions of ISAW. Based on zeros of these exact partition functions, we have developed a novel method in which both loci of zeros and thermodynamic functions associated with them are considered. With this method, the first zeros can be identified clearly without ambiguity. The critical point of a small system can then be defined as the peak position of the heat capacity component associated with the first zeros. For the system with two phase transitions, two pairs of first zeros corresponding to two phase transitions can be identified and overlapping Cv can be well separated. ISAW on the simple cubic lattice is such a system where in addition to a standard collapse transition, there is another freezing transition occurring at a lower temperature. Our approach can give a clear scenario for the collapse and the freezing transitions.


Introduction
The purpose of statistical physics is to understand the macroscopic properties of a many body system from the interactions of the constituents of that system.It has been found that a simple model in statistical physics can often be used to understand phenomena in complex systems [1].
a cnchen@mail.ndhu.edu.twb huck@phys.sinica.edu.twIn 1952, C. N. Yang and T. D. Lee proposed that the critical behavior of the gas-liquid system can be represented by a lattice-gas model [12,13], in which every lattice site is either occupied by an atom or vacant.They showed that the lattice gas model is equivalent to the Ising model.According to their theory, the gas-liquid systems and the three dimensional (3d) Ising model should be in the same universality class and have the same set of critical exponents, a feature which was later confirmed by Monte Carlo simulations of the 3d Ising model [17] and the molecular dynamics simulation (MD) [4] of the Lennard-Jones system [3].Yang and Lee [12,13] also used the zeros of the partition function in the complex fugacity plane to study the critical behavior of the Ising model.Such zeros have been called Yang-Lee zeros or Lee-Yang zeros.
In 1965, M. E. Fisher [15] published a paper on the zeros of the canonical partition function in the complex temperature plane of the square lattice Ising model.He used the exact partition function of the Ising model on the square lattice by Kaufman [9] to show that the partition function zeros are distributed on circles in the thermodynamic limit, and that the logarithmic singularity of the twodimensional Ising model is related to the zero distribution.Such zeros have been called Fisher zeros.
Using subgraph expansion for the partition functions of the Ising model, the Potts model, and a lattice model of hydrogen-bonding in water molecules, Hu showed that the phase transitions of such models can be related to percolation transitions of the corresponding correlated percolation models [50][51][52][53].Based on the connection between the q-state Potts model and a bond-correlated percolation model, Hu and Chen developed a percolation renormalization group method [54,55] based on geometrical factors; Chen and Hu also developed a fast algorithm to calculate exact geometrical factors for the q-state Potts model [56].
To obtain geometrical factors for larger lattices so that the calculated critical point and exponents can be accurate, in early 1992 Hu proposed histogram Monte Carlo simulation methods (HMCSMs) for percolation and phase transition models [57][58][59].Instead of calculating the percolation probability P, the mean cluster size S , etc, at various discrete bond or site occupation probabilities p for percolation problems, Hu used the Monte Carlo simulation method to calculate the histograms of various important quantities from which the geometrical factors, the percolation probability P, the mean clusters S , and the existence probability E p for finite systems at any bond or site occupation probability p may be calculated [57,58].Using the percolation renormalization group method [54,55] and the data of histograms, one may obtain very accurate critical points and exponents, and the thermodynamic order parameters for percolation and Potts models on the square lattice [57,58].The HMCSM had been used to calculate finite-size scaling functions [60], the universal finite-size scaling functions [61] of bond and site percolation models on two-dimensional lattices [62][63][64][65], and three-dimensional lattices [66]; other Monte Carlo methods had been used to calculate universal finite-size scaling functions and exponents for the two-dimensional Ising model [67][68][69].
Since the method of [58] can be used to calculate the partition function of the q-state Potts model, F. Y. Wu proposed to Hu to use the method of [58] to calculate Fisher zeros [15] of the q-state Potts model.Hu considered that the partition function zeros should be highly sensitive on the value of the partition function and thus it is better to use exact partition functions of the q-state Potts model to calculate their zeros.Hu invited Chi-Ning Chen to use exact partition functions obtained by the fast algorithm of [56] to study such a problem.Chen, Hu and Wu found that such zeros are distributed on the unit circle of the complex x = (e K − 1)/ √ q plane for the self-dual square lattices when the real component of the zero is positive [70].
In the following, we will first introduce the model of interacting self-avoiding walks (ISAWs) on lattices, review our results for exact partition functions of ISAWs on the simple cubic lattice [49], then analyze physical quantities of ISAWs from zeros of the exact partition functions.

Self-avoiding walks and interacting self-avoiding walks
On the basis of Flory theory about thermodynamics of high polymer solutions in 1942 [71], in 1947 Orr [40] first introduced the model of interacting self-avoiding walks on lattices and calculated exact partition functions for the ISAWs with N monomers on the simple cubic lattice for 4 ≤ N ≤ 7 [40].
In 1954, F. T. Wall and L. A. Hiller [41] considered a self-avoiding walk model for the macromolecule on lattices and used a Monte Carlo method to calculate mean square end-to-end separation r 2 n of self-avoiding walks of length (step) n on simple cubic (sc) and tetrahedral lattices.They found that for both lattices In 1955, Marshall N. Rosenbluth and Arianna W. Rosenbluth used a Monte Carlo method to calculate the average extension of molecular chains [42].The behavior of the chains of many molecules was investigated by solving a restricted random walk (i.e.self-avoiding walk) problem on the sc lattice and the square lattice.In the Monte Carlo simulations, a large number of chains were randomly generated subject to the restrictions of no crossing or doubling back, to give the average extension of the chain as a function of the chain length n denoting the number of links in the chain.A system of weights was used to ensure that all the possible allowed chains were counted equally.Results for the true random walk problem without weights were also obtained [42].For the sc lattice, the result is consistent with that of Eq. ( 1); for the square lattice, the exponent is 1.45 instead of 1.22 [42].

Exact partition functions for ISAWs on the SC lattice
The energy of an ISAW chain is defined as where m is the number of nonconsecutive nearest-neighbor contacts, and where c m is the number of ISAWs with m contacts, and M is the largest contact number.
Mathematical Modeling and Computational Physics 2015 The first exact enumeration of the ISAWs on the three-dimensional (3d) sc lattice with N = 7 monomers was done in 1947 [40].The exact enumeration of the ISAWs on the sc lattice increased from N = 13, 14 in the mid 1970s [43,44] to N = 24 in 2012 [48].In 2013, we extended the exact enumeration to N = 27 [49].In this paper, we will review our results [49].
To determine c m , we enumerate exhaustively all kinds of configurations of an ISAW chain on the simple cubic lattice by a direct counting algorithm.The convention is the following.We place the first monomer at the origin and place the second monomer at (1,0,0).This choice fixes the degree of freedom of six degenerate directions.Then we consider the third monomer, which has five possible places to locate, and so on.With this convention, Z 2 (x) = 1, Z 3 (x) = 5, and Z 4 (x) = 21 + 4x.Note that our Z 4 (x) is smaller than the corresponding expression in [40] by a factor of 3.This is due to the fact that we placed the second monomer at (1,0,0) while Orr placed the second monomer at (1,0,0) or (0,1,0) or (0,0,1).
We run 25 jobs with 25 different initial configurations in Z 4 (x).Using such parallel scheme and a PC cluster with 20 Intel Core i7-3770 CPUs running for one week, in 2013 we obtained the exact partition functions of an ISAW chain with 27 monomers on the simple cubic lattice [49]: The coefficient of the last term is the number of compact conformations, i.e.Hamiltonian walks, constrained in the 3 × 3 × 3 cube lattice.This number was obtained by Shakhnovich and Gutin [45] in 1990, while all other coefficients have been determined for the first time by us in 2013.
The partition functions for N = 26 and N = 25 obtained by us [49]  Z N (1) and Z N (0) for 14 ≤ N ≤ 27 are listed in Table 1.Z N (1) is the total number of self-avoiding walks (SAWs) of a chain with N monomers (length N − 1) on the sc lattice.Exact values of such numbers for the chain length from 1 (2 monomers) to 36 (37 monomers) were given in Table 1 of [72].Values of Z N (1) in Table 1 below are consistent with values of Z N in Table1 of [72]: N in the former corresponds to N − 1 in the latter, and the value of the former is six times smaller than that of the latter because we put the second monomer in a fixed direction of the sc lattice and the authors of [72] put the second monomer in 6 possible directions of the sc lattice.

Decomposition of physical quantities by zeros of the exact partition function
We used Mathematica to calculate M zeros of Z N (x).
where c M is the ground state degeneracy, and x k is the k-th zero, i.e. root, of Z N (x).Both and k B can be set to unity by changing the units of energy and temperature so that x = x(T ) = e 1/T , and the ground state energy is −M.Following the recipe of the canonical ensemble theory, all thermodynamical functions can be expressed by a summation of contributions of partition function zeros [49].
), ( 6) Since Z N (x) is a polynomial of x with positive coefficients, its roots appear as complex conjugate pairs, except very few roots located on the negative real axis [12].For every complex or real x k , 1 2 and the heat capacity is given by EPJ Web of Conferences Thus the heat capacity C V (x) is decomposed into components C V,k (x) which are real functions of the real argument x and contributed by zeros x k .Or simply speaking, each zero x k carries a heat capacity component C V,k (x).Although all thermodynamical functions are singular at zeros, the response function C V (x)/N will diverge as zeros approach the positive real axis.All these thermodynamical functions may also be suitable for our analysis, but in this paper we only focus on using partition function zeros and C V (x) to analyze the critical phenomena of the 3d ISAW model.We will show that {C V,k (x)} k=1,M provide additional useful information in the analysis of critical phenomena.
Figure 2 shows the partition function zeros with two different variables.Figure 2(a) is the figure with the variable x, and Figure 2(b) are produced from Figure 2(a) by direct variable transformations T = 1/ln(x).The purpose of putting these two zero diagrams together is to show that the patterns of partition function zeros in small systems are quite different when different variables are used.Because the loci of zeros change from one variable to another, the statement based only on the loci of zeros may be misleading.For example, the first zero is usually defined as the zero closest to the real axis.This definition is ambiguous in small systems since the zero closest to the real axis may become the farthest after the transformation of variables.In Figures 2(a) and 2(b), the pairs of the black dots are the same zeros because they can transform to each other.The black dots in Figure 2(a) look like the first zeros, while the black dots in Figure 2(b) don't look like the first zeros at all.
The ambiguity of identifying the first zeros can be easily solved with the heat capacity component C V,k .Figures 3(a) and 3(b) show the zeros represented by dots in the first quadrant in Figures 2(a) and 2(b), respectively, with their corresponding C V,k .In both Figures 3(a) and 3(b), the darker curves correspond to the darker dots, the lighter curves with C V,k peak height 0.798 correspond to the middle lighter dots, and the lighter curves with the peak height 0.642 correspond to the leftmost lighter dots.To avoid using too many symbols, we adopt the shorthand notation: ).Because the darker curve has the largest peak, i.e. its corresponding C V,k contributes significantly to the whole C V , one may define the darker zero to be the first zero.By naive consideration, this largest peak is due to the smaller (x − x k ) 2 term in the denominator of C V,k when x k is getting closer to the real axis.In small systems, nevertheless, it's not always the case.In Figure 2(b), the black dot is not the closest point to the real axis, but its corresponding C V,k has the largest peak, therefore it can still be identified as the first zero.Obviously this definition based on C V,k for the first zero is independent of the variable used.The positions of darker diamonds now can be denoted by x C max V,1 and T C max V,1 .

Collapse and freezing transitions
A 3d ISAW chain undergoes a collapse transition and another freezing transition as temperature changes [47].Here the "transition" point for finite systems is defined as the point where the heat capacity has a peak.According to the scenario proposed in [12,13,15], the partition function zeros of a 3d ISAW chain would approach the real axis at two different places, while the heat capacity has two peaks.However, in small systems, these two peaks sometimes overlap with each other.In Figure 4, the heat capacity C V (T ) of an ISAW chain with 27 monomers is illustrated as the black curve, which had been plotted in [47] and [73] by the improved pruned-enriched-Rosenbluth method (PERM) [74] and the Wang-Landau method [75].This C V curve is a typical one observed frequently in small systems [47,76], where a peak and a shoulder indicate that there are more than one transition.In this situation, C V decomposition shows its value.In Figure 4  zeros associated with the collapse transition.The gray curve is contributed by all remaining zeros.Dotted and dashed curves can be slightly negative for small T , but they are positive in the peak region and the shoulder region, respectively.Thus slight negative values of these curves do not influence the identification of freezing and collapse transitions.Figure 4 demonstrates the necessity of C V decomposition.The overlapping of the freezing transition and the collapse transition can only be separated by the peaks of dotted C V,1 (T ) for the freezing transition and dashed C V,1 (T ) for the collapse transition.The peak positions, T C max V,1 , may be simply defined as the critical points of small systems.For the freezing transition, the peak position of C V (T ), T C max V , can also be defined as the critical point.Since C V (T ) contains contributions from all other zeros, its peak will deviate from the peak of the dotted C V,1 (T ).These two peaks will coincide when N goes to infinity.In finite systems, T C max V,1 from the dotted curve should be closer to the critical point of the infinite system than T C max V .This argument is true for the two-dimensional Potts model on self-dual lattices, where partition function zeros near the exact critical point are located in an unit circle [70].In this case, T C max V,1 (L) is always closer to the the exact critical point than T C max V (L) for any lattice size L. As to the collapse transition, since C V (T ) doesn't have a peak at higher temperature, the only candidate here for the collapse transition temperature is T C max V,1 .
We have calculated the average m and the end-to-end distance d ee (defined in Figure 1) as functions of the temperature T for various N and found that the average m is a decreasing function of T for any N, and the average d ee is an increasing function of T for most N, but is not an increasing function of T for some magic numbers N=9, 13, and 19.We predict that the next magic number is N = 28.The details of such results will be presented in another paper.

Discussion
In summary, we have introduced a C V decomposition scheme by partition function zeros and showed that it can give a clear separation of the collapse and freezing transitions of the ISAW on a simple cubic lattice.This method can also be applied to lattice models for structures, folding, and aggregation of proteins [73,[76][77][78], in which more "phases" are possible for some specific sequences to cause the pattern of zeros scattered.Our approach will be still valuable to the study of such systems.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 201

Figure 1 .
Figure 1.Two examples of interacting self-avoiding walk on the square lattice with N = 17 monomers, and m = 2 (left) and m = 6 (right) nearest-neighbor interacting bonds (contacts).The end-to-end distances d ee of the left and right figures are √ 2a and 2a, respectively, where a is the lattice constant of the square lattice.
is the attractive contact energy.Two examples of the ISAW on the square lattice with N = 17 monomers, and m = 2 (left) and m = 6 (right) nearest-neighbor interacting bonds are shown in Fig. 1.The partition function of an ISAW chain with N monomers is thus a polynomial of the variable x ≡ e k B T ∈ [1, ∞).

Figure 3 .
Figure 3. Zeros represented by dots in the first quadrant in Figs.2(a) and 2(b) are plotted with corresponding heat capacity components C V,k (x) and C V,k (T ).The darker dots indicate the first zeros.The darker diamonds locate at x C max V,1 = 2.131 in 3(a) and T C max V,1 = 1.322 in 3(b), andx C max V,1 = 1/ln(T C max V,1).Three peak heights are 1.44, 0.798, and 0.642 in both (a) and (b).From Fig.2of[49].
, C V,k (T ) is shown together with C V (T ).The dotted curve indicates the heat capacity component contributed by the first zeros associated with the freezing transition.The dashed curve indicates the heat capacity component contributed by the first Mathematical Modeling and Computational Physics 2015 01005-p.7

Figure 4 .
Figure 4.The heat capacity C V (T ) (black curve) of an ISAW chain with 27 monomers.The dotted, dashed and gray curves indicate the heat capacity component C V,k (T ) contributed by the first zeros associated with the freezing transition, the first zeros associated with the collapse transition and all remaining zeros, respectively.

Table 1 .
Numbers of zeros M, real zeros n r , Z N (1), and Z N (0) for the ISAW with N monomers on the sc lattice.