Exact Enumeration Approach to Estimate the Theta Temperature of Interacting Self-Avoiding Walks on the Simple Cubic Lattice

We compute the exact root-mean-square end-to-end distance of the interacting self-avoiding walk (ISAW) up to 27 steps on the simple cubic lattice. These data are used to construct a fixed point equation to estimate the theta temperature of the collapse transition of the ISAW. With the Bulirsch–Stoer extrapolation method, we obtain accurate results that can be compared with large-scale long-chain simulations. The free parameter ω in extrapolation is precisely determined using a parity property of the ISAW. The systematic improvement of this approach is feasible by adopting the combination of exact enumeration and multicanonical simulations.


Introduction
With the rapid development of computer hardware and software, computational approaches [1][2][3] are increasingly important in polymer science. Among various computational methods, exact enumeration is a very primitive technique. Almost four decades passed since W.J.C. Orr [4] studied the equilibrium properties of a single polymer chain at dilute solution by exact enumeration. In this paper, we show that exact enumeration can already produce quantitative results as accurate as those of large-scale simulations. We focus on the interacting self-avoiding walk (ISAW) on the simple cubic lattice [5][6][7][8][9][10]. The ISAW model is a very basic polymer model that serves as the framework of most lattice protein models [11,12].
An ISAW is a self-avoiding walk with attraction between monomers. The energy of an ISAW chain is defined as m(-ε), where m is the number of nonconsecutive nearest-neighbor contacts, and -ε is the attractive contact energy between two monomers. The canonical partition function of a N-step ISAW is where x = e ε/k B T , c m is the number of ISAW configurations with m contacts, and M is the maximal number of contacts. Here, ε and k B can be set to 1 by adjusting the units. Another important function related to the square of the end-to-end distance is defined as follows: where R 2 is the square end-to-end distance, c m,R 2 is the number of ISAW configurations with m contacts and square end-to-end distance R 2 , and c m,R 2 satisfies ∑ N 2 R 2 =1 c m,R 2 = c m .
Both Z N (x) and R 2 N (x) are polynomials of x with positive integer coefficients. The rootmean-square end-to-end distance at a certain temperature can be expressed as follows: Two examples of the normalized end-to-end function R N (T)/ √ N are shown in Figure 1. Both of them are increasing functions, and the slope of the function with larger N is also larger around the intersection point. A polymer chain with attraction between monomers undergoes a collapse transition at theta temperature T θ . R N (T) has different scaling behaviors at different temperature regions [13]: In three dimensions, ν θ and φ can be determined by the mean-field theory to both be 1/2. ν saw can be determined with great accuracy by simulation to be 0.587597(7) [14].
The precise estimation of the theta temperature of the ISAW on the simple cubic lattice came from large-scale simulations [15][16][17][18][19][20]. P. Grassberger [17] proposed the wellknown pruned-enriched Rosenbluth method (PERM) on the basis of the Rosenbluth-Rosenbluth method and the idea of enrichment. For free chains with N = 10,000, the best estimate of the theta temperature was 3.717(3). The recursive sampling algorithm used by P. Grassberger and R. Hegger [15] was a previous version, and the estimated theta temperature was 3.721 (6) with N = 5000. H. Frauenkron and P. Grassberger [18] used the PERM to simulate polymer solutions with N = 2048 and obtained an estimated theta temperature of 3.717(2). T. Vogel et al. [20] used the new PERM with simple sampling up to N = 32,000 and performed a scaling analysis to obtain an estimate of 3.72(1). The above are all chain-growth methods. Tesi et al. [16] used two Markov chain-sampling methods, the multiple Markov chain method and umbrella sampling, to obtain an estimated theta temperature of 3.62 (8) with N = 1600. Yan et al. [19] used the expanded grand-canonical ensemble simulation for polymer solutions with N = 16,000 and obtained an estimate of 3.71(1). In general, Monte Carlo methods must be able to generate unbiased samples and overcome the trapping problems. Exact enumeration, in contrast, is much clearer and simpler. The only challenge with exact enumeration is how to count the total number of larger systems. This is more of a computational problem than a theoretical problem. The solution to this problem benefits directly from the rapid development of computer hardware and software.

Method
The computational methods used in this paper are the exact enumeration algorithm to count the total number of ISAW configurations, and the Bulirsch-Stoer algorithm to extrapolate the finite-size data.

Exact Enumeration
To determine c m and c m,R 2 in Equations (1) and (2), we used a direct counting algorithm that we developed [10] to exhaustively enumerate all configurations of an ISAW chain on the simple cubic lattice. The original goal of this algorithm was to generate enough typical sequences for the 27-mers [21] to study the relationship between protein sequences and structures. Using this algorithm to count all configurations (not just the ground states) of a protein sequence with 27 monomers on the simple cubic lattice now only takes the order of days. The algorithm includes not only direct counting but also reduction in the degrees of freedom in the beginning and final stages. These procedures lead to a significant reduction in computation time.
In the following, we explain some details of the counting process. Figure 2 shows all 22 representative configurations of a four-step ISAW on the simple cubic lattice. The convention is the first monomer being placed at the origin and the second monomer at (1, 0, 0), which fixes the first direction (the darker bond in every configuration). The third monomer has four or five possible directions to choose from, and so on. The three numbers above each configuration in Figure 2 are the number of contacts m, the square end-to-end distance R 2 , and degeneracy coming from the symmetry of the next direction taken. If the numbers are 0, 10, and 4, variable c 0,10 would be incremented by 4 in the program. Let us consider the six configurations in the last row of Figure 2 as an example. They are configurations with one contact, and they contribute to the coefficients of the linear terms in Z 4 (x) and R 2 4 (x). c 1 = 4 + 4 + 8 + 8 + 4 + 4 = 32. c 1,2 = 8 + 8 + 4 + 4 = 24. c 1,4 = 4 + 4 = 8. The coefficient ∑ R 2 R 2 c 1,R 2 = 2 × 24 + 4 × 8 = 80. The final results are Z 4 (x) = 32x + 89 and R 2 4 (x) = 80x + 592. The root-mean-square end-to-end distance is thus The direct counting algorithm has the advantage of easily integrating different ideas and techniques. One straightforward parallel implementation for this direct counting algorithm runs 22 jobs with 22 initial configurations shown in Figure 2. The number of jobs can be flexibly adjusted by choosing a different number of initial configurations depending on how many computer cores are available. Lastly, all results are collected and summed up to be the exact coefficients of Z N (x) and R 2 N (x). All counting jobs can be completed within a few weeks using a small PC cluster.

Bulirsch-Stoer Extrapolation
The Bulirsch-Stoer algorithm [22][23][24][25] may be the most powerful extrapolation method in existence. Its idea was borrowed from recursive algorithms such as the Richardson and Neville algorithms, but the result is much more general. In this subsection, we briefly introduce the Bulirsch-Stoer algorithm and explain how to use its main formula. The detailed derivation and proof can be found in [25].
In this paper, the data of the finite-size theta temperature of the ISAW needs to be extrapolated. Suppose its finite-size scaling behavior can be described by power-law functions: where N is the number of steps of the ISAW chain, φ is the leading exponent, 0 < φ < ω 2 < ω 3 < . . . , and T(N) is the finite-size scaling function. We may first consider the . becomes a polynomial function, to which the Neville algorithm can be applied. Below, we use five data points, {(h i , T(h i ))}, i = 0, . . . , 4, as an example to illustrate explicitly the basic procedures of the Neville and Bulirsch-Stoer algorithms, while general formulas are also provided. First, a triangular 5 × 5 matrix T i,j is prepared, and its first column is filled with {T(h i )}, i = 0, . . . , 4: Equation (9) is the main formula of the Neville algorithm. The final output for this five-point example is T 0,4 , the fourth-degree Lagrange polynomial passing through all five data points. It is an extrapolation function that can be used to approximate the finite-size scaling function T(h). Neville showed that the Lagrangian polynomials can be generated in such an iterated way. The Neville algorithm is very efficient in evaluating function values for interpolation or extrapolation.
Bulirsch and Stoer inserted an additional term: ), in the denominator in Equation (9): The appearance of T i,j−2 requires {T i,−1 }, i = 0, . . . , 4 to be defined first. They can be set to zero in the beginning. The effect of the additional term is that T 0,4 now becomes a rational function P(h)/Q(h), where both P(h) and Q(h) are polynomials. This rational function also passes through all the data points and is another extrapolation function that can be used to approximate the finite-size scaling function T(h). Since T(h) is actually not a polynomial, it is expected that the more general rational functions are more suitable to model T(h) than polynomials. Our extrapolation results also show that the Bulirsch-Stoer algorithm is always more accurate than the Neville algorithm. More details on the properties of such rational functions appearing in Bulirsch-Stoer extrapolation can be found in [25].
In Equation (6), T θ is reached as N → ∞, or h → 0. Substituting h = 0 and h i = N −ω i into Equation (10), Equation (11) is the main formula used in this paper to extrapolate the finite-size theta temperatures. Now T 0,4 is the extrapolation estimation of the theta temperature with five data points. T 0,3 and T 1,3 are also the estimations of the theta temperature but with only four data points used. They are less accurate than T 0,4 , while their difference can serve as an estimate of the error of T 0,4 . A simple argument is as follows. Suppose T 0,3 = T θ ± σ 0,3 , T 1,3 = T θ ± σ 1,3 , and T 0,4 = T θ ± σ 0,4 , where σ 0,3 , σ 1,3 , and σ 0,4 are the statistical errors of Its range is larger than ±σ 0,4 , since both σ 0,3 and σ 1,3 are larger than σ 0,4 . Thus, |T 0,3 − T 1,3 | can roughly play the role of σ 0,4 . Although there are actually no statistics for T i,j , this error estimation is appropriate for most testing examples where the answers are known. The premise is that T 0,3 , T 1,3 , and T 0,4 should be close to each other to indicate that the extrapolation values have entered a stable region. In this paper, we adopt the following definition of the error of T 0,j : In Equation (11), ω is a free parameter. ω does not need to be equal to φ, but it won't be very different either. Each ω is associated with an extrapolation function that passes through all the data points. A different ω results in a different extrapolation value. It is, thus, very important to choose ω carefully. A reasonable choice of ω is the one that minimizes the error defined in Equation (12). In this paper, we use the parity property of the ISAW on the simple cubic lattice to more precisely determine ω.

Result
On the basis of Equation (4), the solution T * of the following equation is an estimation of theta temperature T θ : The choice of N and N should be both odd or both even. The random walks on the square or simple cubic lattice are naturally divided into two groups, namely, random walks with the even number of steps and the odd number of steps. On the simple cubic lattice, even-step walks only stop at the point (i, j, k) with even (i + j + k), and odd-step walks only stop at the point with odd (i + j + k). For this reason, random walks with the same parity are more similar to each other. When N and N are closer, the two walks would also be more similar, and Equation (13) would give more accurate estimations. Our numerical results confirmed this expectation. Therefore, we only discuss the case of N − N = 2 in this paper. Equation (13) becomes We generated square end-to-end functions R 2 N (x) for 10 ≤ N ≤ 27 (listed in Appendix A), and used them and Equation (14) to calculate T * (N). The results were divided into two groups according to parity and are listed in Table 1. Figure 3 also shows that the data points were clearly divided into two groups. The two lines passing through the data points are seventh-degree Lagrangian polynomials. They will merge as N → ∞. The estimation of the theta temperature could be reasonably set as (T I + T I I )/2 ≡ (T * (N odd → ∞) + T * (N even → ∞))/2. Its error is defined as 2 I + 2 I I /2, where I and I I are errors of T I and T I I (Equation (12)) in the Bulirsch-Stoer extrapolation method. Two extrapolations need to be performed here, and the chosen ω needs to make both I and I I smaller. Since there are two constraints, the extrapolation value would be less biased.
To determine an optimal ω, a wide range of ω is scanned first to find the region with small errors. The adjacent area is then zoomed in until the results do not change. Table 2 lists the estimations and errors of the theta temperature in the range of 0.814 ≤ ω ≤ 0.825. The error is minimal as ω = 0.820, so the best estimation was T θ = 3.709 (2). This result is consistent with the long-chain results of large-scale simulations [15][16][17][18][19][20]. We also used the Neville algorithm (Equation (9)) instead of the Bulirsch-Stoer algorithm to extrapolate the same data. The best estimation was T θ = 3.501(4) with ω = 1.404. For all data that we examined, the results of the Bulirsch-Stoer algorithm were always better than the results of the Neville algorithm. Extrapolation with rational functions is expected to be better than extrapolation with polynomials unless the finite-size scaling function is inherently a polynomial. Table 1. Estimation of the theta temperature T * (N) determined by Equation (14).   (14) for odd and even N. The two lines that pass through eight odd data points and eight even data points are seventh-degree Lagrangian polynomial graphs. We also considered a correction term predicted by the field-theoretic renormalization group calculation [26]: This correction term is small and may not agree with the simulations [15,17]. We followed the same procedure as above to calculate T * (N) with Equation (15), and list the results in Table 3 for comparison. Table 3. The estimation of theta temperature T * (N) determined by Equation (15). The best estimation from the data in Table 3 was T θ = 3.713 (2) with ω = 0.8617. Since T * (N) listed in Table 3 was closer to the value of the "real" T θ than those in Table 1, the extrapolation result seemed to be also improved.
We may use the same procedure to calculate other quantities. For example, we can verify the value of the crossing exponent φ in Equation (4) in the following way.
Equation (17) is a ratio method. We took 3.713 as the value of T θ . The results are listed in Table 4. The best estimation of φ was 0.531 (9) with ω = 1.76, while the exact value is 1/2. This result is not as accurate as that of theta temperature estimation. The correction term seems to be important when computing critical exponents. We have tried to include a correction term and obtained much more accurate φ results. However, since this paper mainly regards the efficiency of direct computation, we do not discuss the effect of correction further.  (17).

Discussion
The longest ISAW chain used in the calculation of this paper is the 27-step walk. The total number of its configurations is huge: 431,645,810,810,533,429 ∼ 4 × 10 17 . From such a large number and others, we obtained an accurate estimate of the theta temperature of ISAW on the simple cubic lattice. Nevertheless, a 27-step ISAW is very short in three dimensions. The finite-size effect should be obvious. There are two reasons why short chains can still give accurate results. The first reason is the careful determination of the free parameter ω in the Bulirsch-Stoer algorithm. For ISAWs on the simple cubic lattice, we could just take advantage of their natural division into odd and even groups. There are two constraints on the choice of ω, resulting in more stable and accurate extrapolation results. This experience may be extended to other problems where data points are not smooth. The way of classifying the data points and combining the extrapolation results could work. If only one series of data is extrapolated, unless the data points are very smooth, the ω corresponding to the minimal error does not necessarily produce extrapolation values very close to the true answer.
Another reason is that Equation (13) may be regarded as a fixed-point equation of the real-space renormalization group transformation for the coarse-graining process, where each N-step segment of the ISAW is replaced by a N -step segment (N < N). In this view, N and N are the block sizes of the transformation and not the length of the whole ISAW chain. The finite-size effect could, therefore, be less obvious. Equation (13) can be rewritten as the equation to transform T to T : R N (T )/ √ N = R N (T)/ √ N, which will drive the temperature away from the intersectional temperature T * . Figure 1 shows this feature: if T > T * , T is larger than T, and T < T * , T is smaller than T. A simple explanation is that, so T needs to be larger to balance the transformation equation. A similar logic could be used for the case of T < T * . Only in the critical region around T * , T ≈ T, and the normalized end-to-end distances of the two ISAW segments with length N and N are approximately the same. The larger and closer N and N are, the more similar these two segments are, and the transformation equation is more accurate, so that T * is closer to T θ . Traditional real-space renormalization approaches for polymer models usually focus on the variable fugacity, e.g., [27]. The above scenario of the real-space renormalization along the ISAW chain is only for the variable temperature and can be developed in more detail. It is also similar to phenomenological renormalization, e.g., [28], in which the correlation length plays the same role as R N (T) and can be calculated by the transfer matrix method.
The number of configurations of an N-step ISAW grows exponentially with N, e.g., µ N with µ ∼ 4.7 for the simple cubic lattice. From 27 to 30 steps, the number of configurations is increased by about 100 times. To go further, the degree of freedom of the problem must be reduced in some way. The transfer matrix technique is efficient for the counting problem of the two-dimensional ISAW [29,30], but difficult to extend to the three-dimensional case. Matrix methods have been used in polymer science for many years, e.g., [31,32], and deserve further development. The length-doubling method [33,34] reduces the degree of freedom by counting and saving the information of two short walks and then joining them to form a longer walk. It is very feasible for us to use this approach because we had developed a counting algorithm with a similar idea [35]. We counted the number of graphs in the percolation and Potts models by dividing the graphs into two parts (or several layers) and then combining the two parts of the data to obtain the complete result. This algorithm was used to study the partition function zeros of the Potts model on the self-dual lattices [36].
We can push to the limit the size of the system that can be exactly enumerated to obtain more accurate results. However, this direction may become very technical. A more practical approach is to perform the multicanonical simulation [37,38] to count the configurations of longer chains. The multicanonical simulation accumulates the density of states during Monte Carlo sampling and is very different from the Metropolis algorithm, in which sampling is just for averaging. This approach can be called Monte Carlo counting, corresponding to exact counting. Among multicanonical methods, the Wang-Landau algorithm [39] is commonly used because of its concise steps and high efficiency. It has been applied to polymer simulations of various systems, including the ISAW [40]. We have tried to combine the data of medium-length chains by Wang-Landau sampling with the data of short chains by exact enumeration. Although the numbers of configurations from the Wang-Landau method are approximate, they carry information about the longer chains and would stabilize the extrapolation curve around the medium-length region. Therefore, the data from the Wang-Landau method would certainly improve the accuracy of exact enumeration. If there is no need to reach the limit, the total computation time spent on exact enumeration and Wang-Landau sampling together could still be less than that of large-scale simulations. This approach of combining exact counting and Monte Carlo counting might become a practical and general technique in computational science.
In summary, we used exact enumeration and Bulirsch-Stoer extrapolation to obtain an accurate estimate of the theta temperature of the ISAW on the simple cubic lattice. The systematic improvement of this approach by increasing the chain length is achievable. Research in this direction is in progress.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Square end-to-end distance function R 2 N (x) was used in the calculation of this paper. Z N (x) can be found in [8][9][10] and references therein. The convention adopted in these references was Z n (x), where n is the number of monomers, i.e., n = N + 1.