Irreversible Markov chain Monte Carlo algorithm for self-avoiding walk

Abstract

We formulate an irreversible Markov chain Monte Carlo algorithm for the self-avoiding walk (SAW), which violates the detailed balance condition and satisfies the balance condition. Its performance improves significantly compared to that of the Berretti–Sokal algorithm, which is a variant of the Metropolis–Hastings method. The gained efficiency increases with spatial dimension (D), from approximately 10 times in 2D to approximately 40 times in 5D. We simulate the SAW on a 5D hypercubic lattice with periodic boundary conditions, for a linear system with a size up to L=128, and confirm that as for the 5D Ising model, the finite-size scaling of the SAW is governed by renormalized exponents, {\nu }^{*}=2/d and \gamma /{\nu }^{*}=d/2. The critical point is determined, which is approximately 8 times more precise than the best available estimate.

Graphical abstract

Keywords

Monte Carlo algorithms / self-avoiding walk / irreversible / balance condition

1 Introduction

The self-avoiding walk (SAW) serves as a paradigmatic model in polymer physics (see, e.g., Ref. [1] and the references therein). It is equivalent to the n\to 0 limit of the O(n) model [2], and plays an important role in the study of critical phenomena. In the grand-canonical ensemble, the length of a walk can fluctuate and the SAW model is defined by the following partition sum:

\begin{array}{}\text{(1)}& \mathcal{Z}=\sum _{\omega }{x}^{|\omega |},\end{array}

where |\omega | is the length of the walk, ω, x is the weight of each unit length, and the summation is over all possible self-avoiding paths. In two and higher dimensions (D), the SAW has two distinct phases separated by a critical point, x_c. The length, |\omega |, remains finite in the dilute phase with , and becomes divergent in the dense region with x>x_c.

Markov chain Monte Carlo (MCMC) methods have been extensively used in simulation of the SAW [3]. The balance condition (BC) and ergodicity are two key factors in designing an MCMC algorithm. The BC states that the probability flow entering into a configuration equals the flow out of the configuration. Thus, it ensures a stationary distribution. Then, ergodicity ensures convergence to the distribution [4]. In practice, the BC is typically satisfied by employing the detailed balance condition (DBC), which implies that the probability flow from one configuration to another is equal to reverse flow, i.e., the dynamics is reversible.

In recent years, there have been several successful studies [5-14] that show a promising future for MCMC algorithms beyond the DBC. Geometric allocation approaches have been applied to the Potts model [5, 6]; irreversible MCMC methods have been designed for the mean-field Ising model [7, 8]; event-chain Monte Carlo (ECMC) methods have been proposed for simulation of hard-sphere systems [9, 10] and generalized to particle systems with arbitrary pairwise interactions [11], including soft-disk systems [12], the XY model [13], and the Heisenberg model [14]. Near the phase transition point, the geometric allocation method outperforms the standard Metropolis–Hastings (MH) method by 6.4 times for the q=4 square lattice Potts model, and its performance increases with q [5, 6]. For the mean-field Ising model, the irreversible MCMC method has a dynamic exponent, z\simeq 0.85, which is considerably smaller than z\simeq 1.43 for the reversible MH method [7, 8]. In comparison with the MH method, the gained efficiency of the ECMC method reaches two orders of magnitude in large systems consisting of {10}^6 hard spheres [10]. For the 3D ferromagnetic Heisenberg model, it has been reported that the ECMC method has a dynamic exponent, z\simeq 1, in contrast to z\simeq 2 for the MH method [14].

For a Monte Carlo Markov chain, let \pi (\omega ) be the weight of a configuration, \omega, A(\omega \to {\omega }^{\prime }) be the a priori probability of proposing a transition to another configuration, {\omega }^{\prime }, and P(\omega \to {\omega }^{\prime }) be the probability of accepting the proposal. One obtains the stationary probability flow, \varphi (\omega \to {\omega }^{\prime })\equiv \pi (\omega )A(\omega \to {\omega }^{\prime })P(\omega \to {\omega }^{\prime }). The DBC states that for any pair of \omega and {\omega }^{\prime }, \varphi (\omega \to {\omega }^{\prime })=\varphi ({\omega }^{\prime }\to \omega ). Instead, the BC requires that for any ω, {\sum }_{{\omega }^{\prime }}\varphi (\omega \to {\omega }^{\prime })={\sum }_{{\omega }^{\prime }}\varphi ({\omega }^{\prime }\to \omega ), where the summation is over all possible configurations, {\omega }^{\prime }. Without the DBC, net probability flows can exist between two states, ω and {\omega }^{\prime }, i.e., \varphi (\omega \to {\omega }^{\prime })\ne \varphi ({\omega }^{\prime }\to \omega ), and the probability fluxes make circles in a phase space [7].

While local probability flow circles are introduced in the geometric allocation approach, considerably larger or even global circles can appear in other methods beyond the DBC. This is achieved through a lifting technique, which enlarges a phase space by an auxiliary variable. As a result, a system can be at different modes specified by different values of the auxiliary variable. Within a given mode, an a priori direction of Monte Carlo updates is preferred, and probability flow may exist. The BC is recovered by allowing switches between different modes, and the probability flow can satisfy a skew detailed balance condition (see, e.g., Refs. [15, 16]). For example, the phase space for the mean-field Ising model [7] is doubled and denoted as decreasing and increasing modes. In the first mode, only positive spins are flipped and total magnetization, M, decreases, whereas in the second mode, M increases only through flipping of negative spins. Updates can persist for a long time in one mode until a spin-flip proposal is rejected, after which the updates are switched to the other mode. This leads to large probability flow circles, and the diffusive feature of random updates is suppressed or even replaced by ballistic-like behavior. It is noted that while the efficiency of the mean-field Ising model improves qualitatively, the lifting technique does not help significantly in the 2D Ising model [8, 17].

In this work, we design an irreversible MCMC algorithm for the SAW by employing the lifting technique. The update direction is selected such that the length |\omega | increases in one mode and decreases in the other. The two modes are “linked” to each other through switching. For practical coding, only a few lines need be added to the widely used Berretti–Sokal (BS) algorithm, a simple variant of the MH method [18]. Nevertheless, numerical results show that the irreversible MCMC method is considerably superior to the BS algorithm. We use this new method to explore the finite-size scaling (FSS) of the 5D SAW, above the upper critical dimension, d_{\mathrm{u}}=4. The critical point is located with high precision.

The rest of this article is organized as follows: In Section 2, we review a few Monte Carlo methods for the SAW, including the conventional MH algorithm and the BS method. Section 3 describes the irreversible MCMC algorithm. We compare the performances of these algorithms in Section 4. Section 5 contains an FSS analysis for the SAW on a periodic hypercubic lattice in five dimensions. A brief conclusion and discussion are presented in Section 6.

2 Reversible MCMC algorithms

We review below the MH and BS algorithms for the SAW, which employ the DBC and are reversible.

[Metropolis–Hastings algorithm.] Given a regular d-dimensional lattice with coordination number z, the SAW, ω, with length \mathcal{N}\equiv |\omega |, is a sequence of \mathcal{N}+1 lattice sites connected through a chain of occupied edges. For simplicity, we fix an end of the walk at the origin and denote the movable end by I. The MH algorithm [19, 20] is a standard reversible MCMC algorithm, which is constructed as follows:

Metropolis–Hastings Algorithm

i) Randomly select one of the z neighboring sites of I, e.g., I^{\prime }, and propose a symmetric update of the edge in between, which flips an empty edge to be occupied, and vice versa.

ii) If the update leads to a valid SAW, {\omega }^{\prime }, accept the proposal with probability P=\min \{1,\pi ({\omega }^{\prime })/\pi (\omega )\}.

The weight \pi (\omega ) is given by Eq. (1) as \pi (\omega )=x^{|\omega |}, and the acceptance probabilities are

\begin{array}{}\text{(2)}& \begin{array}{l}P(\Delta \mathcal{N}=+1)=\text{m}\text{i}\text{n}\left\{\mathrm{1,}x\right\}\mathrm{,}\\ P(\Delta \mathcal{N}=-1)=\text{m}\text{i}\text{n}\{\mathrm{1,1}/x\}\mathrm{,}\end{array}\end{array}

where \mathrm{\Delta }\mathcal{N}=|{\omega }^{\prime }|-|\omega |. Here, the a priori probability is a constant, A(\omega \to {\omega }^{\prime })=1/z, independent of \mathrm{\Delta }\mathcal{N}=\pm 1. The DBC can be easily proven.

[Berretti–Sokal algorithm.] On a lattice with a large coordination number, z (e.g., a high-dimensional lattice), the critical point of the SAW occurs at x_c\approx 1/(z-1). This implies that, in the abovementioned MH algorithm, a typical update would attempt to flip an empty edge; however, it will most likely fail, and thus, it is ineffective. A more efficient MH algorithm for SAWs is the BS method [18], which employs a different scheme for a priori probabilities. A version of the BS algorithm, slightly different from that in Ref. [18], is described below.

Berretti–Sokal Algorithm

i) Choose Action: propose, with equal probability, the “add” action, a_{+}, for \mathrm{\Delta }\mathcal{N}=+1 and the “delete” action, a_{-}, for \mathrm{\Delta }\mathcal{N}=-1.

ii) Perform Action:

• For action a_{+}, randomly occupy, with probability P_{\mathrm{BS}}^{+}, one of the (z-1) empty edges if it leads a valid SAW.

• For action a_{-}, delete, with probability P_{\mathrm{BS}}^{-}, the last occupied edge, incident to I.

Note that an empty edge is chosen with probability 1/2(z-1) while the occupied edge is chosen with probability 1/2. The DBC condition leads to the Metropolis acceptance probabilities as

\begin{array}{}\text{(3)}& \begin{array}{l}{P}_{\text{B}S}^{+}=\text{m}\text{i}\text{n}\{\mathrm{1,}x(z-1)\}\mathrm{,}\\ P_{\text{B}S}^{-}=\text{m}\text{i}\text{n}\{\mathrm{1,1}/x(z-1)\}\mathrm{.}\end{array}\end{array}

Near the critical point, the acceptance probabilities are close to unity. In the original version [18], a priori probabilities are slightly more optimized such that the BS algorithm becomes rejection free.

It is noted that special attention is required for the case of a “null” walk with \mathcal{N}=0, in which there are z empty incident edges instead of z-1. For simplicity, we permanently delete an edge incident to the original site.

The diffusive feature of the BS algorithm is clear. Particularly for x(z-1)\approx 1, the SAW length {\mathcal{N}}^{\prime } in the next step will be {\mathcal{N}}^{\prime }=\mathcal{N}\pm 1 with approximately equal probabilities.

3 Irreversible MCMC method

A direct approach toward an irreversible method is to double the state space of the SAW by introducing an auxiliary variable with two values, (+) and (-). In the increasing mode, (+), action a_{-} is forbidden and the walk length, \mathcal{N}, increases because of action a_{+}. In contrast, \mathcal{N} decreases in the decreasing mode, (-). The balance condition is satisfied by allowing switches between the two modes. Figure 1 shows a sketch of the associated probability flows. The formulation of the irreversible algorithm is considerably similar to the abovementioned BS method, as illustrated below.

Irreversible Algorithm

i) For the increasing mode, (+), perform action a_{+} with probability P^{+}. Randomly select and occupy one of the (z-1) empty edges if this leads to a valid SAW; otherwise, switch to the decreasing mode, (-).

ii) For the decreasing mode, (-), perform action a_{-} with probability P^{-}. Delete the last occupied edge in the SAW if \mathcal{N}>0; otherwise, switch to the increasing mode, (+).

The acceptance probabilities, P^{+} and P^{-}, can also be given by Eq. (3). For the case where , the SAW grows with probability P^{+}=x(z-1) in the increasing mode, (+), and it is deleted until \mathcal{N}=0 in the decreasing mode, (-). For the case where x(z-1)>1, which includes the critical point, x_c, one obtains P^{+}=1. Thus, in the increasing mode, (+), the SAW grows until it violates self-avoidance, after which it switches to the decreasing mode, (-), where the chain decreases with probability P^{-}=1/x(z-1). At the critical point, x_c, P^{-} is close to unity, for example, P^{-}\simeq 0.9823 in 5D. Thus, instead of returning to \mathcal{N}=0 and growing again, it decreases for a considerably long time, then, it switches to the increasing mode, (+), and grows again. It is clear that the diffusive feature of Monte Carlo updates in the state space is significantly suppressed.

We demonstrate below the balance condition for the case where , using the acceptance probabilities given by Eq. (3); the proof for x(z-1)\ge 1 follows the same procedure. Consider an \mathcal{N}-step SAW in the increasing mode, (+); the incoming probability flow from the (\mathcal{N}-1)-step SAW is

\begin{array}{}\text{(4)}& {\varphi }_{\mathrm{in},a_{+}}^{(+)}=[1/(z-1)]x^{\mathcal{N}-1}{P}^{+}=x^{\mathcal{N}}.\end{array}

The factor, 1/(z-1), accounts for the probability of selecting the edge leading to the current SAW. The incoming probability flow due to the switch from mode (-) is equal to zero, unless \mathcal{N}=0 and {\varphi }_{\mathrm{in},\mathrm{s}}^{(+)}=x^{\mathcal{N}}(1-P^{-})=0.

The total outgoing probability flows are clearly x^{\mathcal{N}}, because no action is allowed to keep the configuration unchanged. In the next step, there is an (\mathcal{N}+1)-step SAW in mode (+) or an \mathcal{N}-step SAW in mode (-). More specifically, suppose that occupying one of the z^{\prime }\in [0,z-1] empty edges leads to a valid SAW, the outgoing probability flow is {\varphi }_{\mathrm{out},a_{+}}^{(+)}=[z^{\prime }/(z-1)]x^{\mathcal{N}}{P}^{+}=z^{\prime }{x}^{\mathcal{N}+1} and the switch probability flow is {\varphi }_{\mathrm{out},\mathrm{s}}^{(+)}=x^{\mathcal{N}}(1-z^{\prime }x). Thus, the balance condition in the increasing mode, (+), is satisfied as

\begin{array}{}\text{(5)}& {\varphi }_{\mathrm{in},a_{+}}^{(+)}+{\varphi }_{\mathrm{in},\mathrm{s}}^{(+)}={\varphi }_{\mathrm{out},a_{+}}^{(+)}+{\varphi }_{\mathrm{out},\mathrm{s}}^{(+)}=x^{\mathcal{N}}.\end{array}

The same procedure is followed for the BC in mode (-).

4 Performance

We conducted simulations for the SAW at criticality on a d-dimensional periodic hypercubic lattice, for d=2 to 5, where the critical value, x_c, is listed in Table 1. For each linear size, L, we carried out 5\times {10}^6/2^d sweeps (L^d Monte Carlo steps) of simulations, in which one fifth were thrown for thermalization. The number of Monte Carlo steps between successive samples is L/2 for d=2,3, and L^2/4 for d=4,5.

We compare the efficiency of the algorithms according to the integrated autocorrelation time, τ, for an arbitrary observable, defined by [25]

\begin{array}{}\text{(6)}& \delta \mathcal{O}=\sqrt{\frac{1+2\tau /\mathrm{\Delta }\tau }{n-1}(\overline{{\mathcal{O}}^2}-{\overline{\mathcal{O}}}^2)},\end{array}

where \delta \mathcal{O} is the standard deviation, \mathrm{\Delta }\tau denotes the number of sweeps between successive samples, and n is the number of samples. For \mathrm{\Delta }\tau \gg \tau, i.e., when the successive samples are effectively independent, Eq. (6) is simplified as \delta \mathcal{O}=\sqrt{(\overline{{\mathcal{O}}^2}-{\overline{\mathcal{O}}}^2)/(n-1)}.

In the simulations, we sampled the walk length, \mathcal{N}, and the observable, {\mathcal{D}}_0, that describes the event of a null SAW; {\mathcal{D}}_0=1 for \mathcal{N}=0 and {\mathcal{D}}_0=0 otherwise. The statistical average, D_0=⟨{\mathcal{D}}_0⟩, accounts for the probability that the walk end, I, returns to the original site.

The autocorrelation time, τ, is measured for \mathcal{N} and {\mathcal{D}}_0. Figures 2 and 3 compare the autocorrelation times, \tau (\mathcal{N}) and \tau ({\mathcal{D}}_0), for the MH, BS, and irreversible algorithms. It can be seen that the performance of the irreversible algorithm is considerably superior to that of the BS and MH algorithms. The gained efficiency becomes more pronounced as d increases. For d=2, it outperforms the MH and BS algorithms by approximately 12 and 8 times, respectively, while they become approximately 200 and 40 for d=5. As an illustration, Table 2 shows \tau (\mathcal{N}) for various algorithms, for the maximum linear size, L_{\max }, in this comparative study. We also compare the efficiency in terms of CPU time using a desktop computer with 3.8 GB memory and four Intel i5 cores. As shown in the last column of Table 2, the results are similar to those for the walk length, \mathcal{N}.

5 Finite-size scaling above the upper critical dimension

The method of FSS, which is derived from the renormalization group theory, plays a fundamental role in the numerical study of critical phenomena. It predicts that near x_c, the energy-like quantity, specific-heat-like quantity, and magnetic susceptibility, i.e., E, C, and χ, respectively, scale as

\begin{array}{}\text{(7)}& \begin{array}{l}E(x\mathrm{,}L)=L^d{E}_{\text{r}}(x)+L^{1/\nu }{E}_{\text{s}}[L^{1/\nu }(x-x_c)]\mathrm{,}\\ C(x\mathrm{,}L)=C_{\text{r}}(x)+L^{\alpha /\nu }{C}_{\text{s}}[L^{1/\nu }(x-x_c)]\mathrm{,}\\ \chi (x\mathrm{,}L)={\chi }_{\text{r}}(x)+L^{\gamma /\nu }{\chi }_{\text{s}}[L^{1/\nu }(x-x_c)]\mathrm{,}\end{array}\end{array}

where the critical exponents, α and γ, are for the thermodynamic (i.e., infinite system size) quantities, C and χ, respectively. In Eq. (7), the first term accounts for the regular functions that are size-independent (except that the regular part of energy is proportional to L^d), while E_{\mathrm{s}}, C_{\mathrm{s}}, and {\chi }_{\mathrm{s}} are universal functions that account for singular behavior.

The FSS formula (7) is correct in dimensions lower than the upper critical dimensionality, d_{\mathrm{u}}. For d>d_{\mathrm{u}}, the thermodynamic critical exponents take mean-field values; however, the FSS behavior is considerably more complicated. Different FSS behaviors occur for different boundary conditions, and for the \mathit{k}=0 and the \mathit{k}\ne 0 fluctuations. Extensive studies have been carried out for the 5D Ising model, i.e., the n=1 case of the O(n) model, for which the renormalization group theory gives d_{\mathrm{u}}=4, and mean-field exponents \nu =1/2, \gamma =1 and \alpha =0. For periodic boundary conditions, the FSS formula (7) holds if the critical exponent, \nu =1/2, is replaced by a renormalized exponent, {\nu }^{*}=2/d, even though there are debates on the involved physical scenarios [26-29].

Similar to the Ising model, the SAW is a special case of the O(n) model in the n\to 0 limit. Using the irreversible algorithm, we performed extensive simulations for the SAW on a 5D periodic hypercubic lattice, up to L=128. This provides an independent and accurate study of the d>d_{\mathrm{u}} FSS behavior for the O(n) universality. In comparison with the Ising model, such a study of the 5D SAW has a few advantages. The regular terms, C_{\mathrm{r}}(x) and {\chi }_{\mathrm{r}}(x), in Eq. (7) vanish for x\le x_c, as the average walk length per site, N/L^d=⟨\mathcal{N}⟩/L^d, approaches zero as L increases (N is the energy-like quantity with a zero regular part). Further, the simulation can reach a relatively large linear size, whereas it is mostly restricted to for the 5D Ising model.

We sampled the average walk length, N, the specific-heat-like quantity, C=L^{-d}(⟨{\mathcal{N}}^2⟩-{⟨\mathcal{N}⟩}^2), and the returning probability, D_0=⟨{\mathcal{D}}_0⟩, which is the inverse of magnetic susceptibility, D_0=1/\chi [30]. The Taylor expansion of Eq. (7) leads to

\begin{array}{}\text{(8)}& \begin{array}{l}N=L^{1/{\nu }^{*}}\left[n_0+{\displaystyle \sum _{k=1}^2}{n}_k{(x-x_c)}^k{L}^{k/{\nu }^{*}}+b_n{L}^{y_1}\right]\mathrm{,}\\ C=c_0+{\displaystyle \sum _{k=1}^2}{c}_k{\left(x-x_c\right)}^k{L}^{k/{\nu }^{*}}+b_c{L}^{y_1}\mathrm{,}\end{array}\end{array}

where the term with exponent y₁ accounts for finite-size corrections and the coefficients of each term are non-universal. Besides the renormalized exponent, 1/{\nu }^{*}=d/2, the leading correction exponent is predicted as y_1=(d_{\mathrm{u}}-d)/2=-1/2 [29, 31]. The FSS behavior of C in Eq. (8) is based on the prediction \alpha /{\nu }^{*}=0. The data of N and C are shown in Fig. 4, in which x_c is represented by the approximate intersection point for different sizes, and the rescaled plots in the right panels imply the correctness of Eq. (8). The MC data were fitted to Eq. (8) according to the least-squared criterion, and the results are shown in Table 3. The correction exponent, y₁, is fixed at a predicted value of -1/2 in the fits. As a precaution against correction-to-scaling terms not included in the fit ansatz, we imposed a lower cut off, L\ge L_{\min }, on the data points, and observed the change in {\chi }^2/DF, where “DF” represents the number of degrees of freedom. The estimates of 1/{\nu }^{*} agree well with the predicted value, 5/2. We take the final determination, x_c=\mathrm{0.113\hspace{0.33em}140\hspace{0.33em}843}(5), which is significantly improved compared to the best available result, x_c=\mathrm{0.113\hspace{0.33em}140\hspace{0.33em}81}(4) [24]. The reliability of the estimate of x_c is further demonstrated in Fig. 5, which clearly shows that x=\mathrm{0.113\hspace{0.33em}140\hspace{0.33em}823} and \mathrm{0.113\hspace{0.33em}140\hspace{0.33em}863} are below and above the critical point, respectively.

At the critical point, x_c, the FSS behavior of the returning probability, D₀, should be

\begin{array}{}\text{(9)}& D_0=L^{-\gamma /{\nu }^{*}}(d_0+d_1{L}^{y_1}),\end{array}

with a mean-field value of \gamma =1. This is confirmed by Fig. 6.

6 Conclusion and discussion

We develop an irreversible MCMC algorithm for the SAW. It violates the DBC, and satisfies the weaker BC. In comparison with the standard MH algorithm and one of its variants, i.e., the BS algorithm, the irreversible method is considerably more efficient. While the BS algorithm is approximately d times more efficient than the standard MH algorithm, the irreversible algorithm is considerably superior to the BS algorithm. The higher is the spatial dimension, the more is the gain in efficiency. This is because the critical SAW is more similar to the ordinary random walk in higher dimensions, and thus, the diffusive feature is more suppressed in the irreversible algorithm.

Using the irreversible MCMC algorithm for the SAW, we perform an independent and accurate test of the renormalized exponents, 1/{\nu }^{*}=d/2 and \gamma /{\nu }^{*}=d/2, in the FSS behavior of systems in the O(n) universality class with periodic boundary conditions. Additionally, we calculate an estimate, x_c=\mathrm{0.113\hspace{0.33em}140\hspace{0.33em}843}(5), for the 5D SAW, which is 8 times more precise than the best available estimate.

We believe that the irreversible MCMC algorithm for the SAW will make an important contribution toward a deeper understanding of the FSS behavior above the upper critical dimension, particularly for systems with free boundary conditions. Moreover, this algorithm will be considerably valuable in the study of interacting SAWs [3], polymeric systems, and other soft matters. Similar irreversible techniques can be introduced into other algorithms such as the worm algorithm [30, 32, 33]. These are ongoing research activities.

Several efficient Monte Carlo methods exist for simulating SAWs [3, 34-36]. In particular, the pivot algorithm [35] and the pruned-enriched Rosenbluth method (PERM) [36] are known to be considerably efficient in the canonical ensemble where the SAW chain has a fixed length. For example, the critical point, x_c, for d=3,4,5 in Table 1 was obtained using the PERM [23, 24]. It would be interesting to study whether it is possible to implement the irreversible technique in these state-of-the-art algorithms for SAWs.

References

[1] P. G. de Gennes, Scaling Concepts in Polymer Physics, Cornell University, Ithaca, 1979

[2] P. G. de Gennes, Exponents for the excluded volume problem as derived by the Wilson method, Phys. Lett. A 38(5), 339 (1972)

[3] E. J. Janse van Rensburg, Monte Carlo methods for the self-avoiding walk, J. Phys. A Math. Theor. 42(32), 323001 (2009)

[4] J. R. Norris, Markov Chains, Cambridge University Press, 1998

[5] H. Suwa and S. Todo, Markov chain Monte Carlo method without detailed balance, Phys. Rev. Lett. 105(12), 120603 (2010)

[6] S. Todo and H. Suwa, Geometric allocation approaches in Markov chain Monte Carlo, J. Phys. Conf. Ser. 473, 012013 (2013)

[7] K. S. Turitsyn, M. Chertkov, and M. Vucelja, Irreversible Monte Carlo algorithms for efficient sampling, Physica D 240(4–5), 410 (2011)

[8] H. C. M. Fernandes and M. Weigel, Non-reversible Monte Carlo simulations of spin models, Comput. Phys. Commun. 182(9), 1856 (2011)

[9] E. Bernard, W. Krauth, and D. Wilson, Event-chain Monte Carlo algorithms for hard-sphere systems, Phys. Rev. E 80(5), 056704 (2009)

[10] M. Engel, J. A. Anderson, S. C. Glotzer, M. Isobe, E. P. Bernard, and W. Krauth, Hard-disk equation of state: First-order liquid-hexatic transition in two dimensions with three simulation methods, Phys. Rev. E 87(4), 042134 (2013)

[11] M. Michel, S. C. Kapfer, and W. Krauth, Generalized event-chain Monte Carlo: Constructing rejection-free global balance algorithms from infinitesimal steps, J. Chem. Phys. 140(5), 054116 (2014)

[12] S. C. Kapfer and W. Krauth, Soft-disk melting: From liquid-hexatic coexistence to continuous transitions, Phys. Rev. Lett. 114(3), 035702 (2015)

[13] M. Michel, J. Mayer, and W. Krauth, Event-chain Monte Carlo for classical continuous spin models, EPL 112(2), 20003 (2015)

[14] Y. Nishikawa, M. Michel, W. Krauth, and K. Hukushima, Event-chain algorithm for the Heisenberg model: Evidence for z~1 dynamic scaling, Phys. Rev. E 92(6), 063306 (2015)

[15] K. Hukushima and Y. Sakai, An irreversible Markov-chain Monte Carlo method with skew detailed balance conditions, J. Phys. Conf. Ser. 473, 012012 (2013)

[16] Y. Sakai and K. Hukushima, Dynamics of one-dimensional Ising model without detailed balance condition, J. Phys. Soc. Jpn. 82(6), 064003 (2013)

[17] R. D. Schram and G. T. Barkema, Monte Carlo methods beyond detailed balance, Physica A 418, 88 (2015)

[18] A. Berretti and A. D. Sokal, New Monte Carlo method for the self-avoiding walk, J. Stat. Phys. 40(3–4), 483 (1985)

[19] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equations of state calculations by fast computing machines, J. Chem. Phys. 21(6), 1087 (1953)

[20] W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 57(1), 97 (1970)

[21] I. Jensen, A parallel algorithm for the enumeration of self-avoiding polygons on the square lattice, J. Phys. Math. Gen. 36(21), 5731 (2003)

[22] J. L. Jacobsen, C. R. Scullard, and A. J. Guttmann, On the growth constant for square-lattice self-avoiding walks, arXiv: 1607.02984 (2016)

[23] H. P. Hsu and P. Grassberger, Polymers confined between two parallel plane walls, J. Chem. Phys. 120(4), 2034 (2004)

[24] A. L. Owczarek and T. Prellberg, Scaling of self-avoiding walks in high dimensions, J. Phys. Math. Gen. 34(29), 5773 (2001)

[25] H. Müller-Krumbhaar and K. Binder, Dynamic properties of the Monte Carlo method in statistical mechanics, J. Stat. Phys. 8, 1 (1973)

[26] K. Binder, M. Nauenberg, V. Privman, and A. P. Young, Finite-size tests of hyperscaling, Phys. Rev. B 31(3), 3 (1985)

[27] B. Berche, R. Kenna, and J. C. Walter, Hyperscaling above the upper critical dimension, Nucl. Phys. B 865(1), 115 (2012)

[28] F. Flores-Sola, B. Berche, R. Kenna, and M. Weigel, Role of Fourier modes in finite-size scaling above the upper critical dimension, Phys. Rev. Lett. 116(11), 115701 (2016)

[29] M. Wittmann and A. P. Young, Finite-size scaling above the upper critical dimension, Phys. Rev. E 90(6), 062137 (2014)

[30] Y. Deng, T. M. Garoni, and A. D. Sokal, Dynamic critical behavior of the worm algorithm for the Ising model, Phys. Rev. Lett. 99(11), 110601 (2007)

[31] E. Brézin and J. Zinn-Justin, Finite size effects in phase transitions, Nucl. Phys. B 257, 867 (1985)

[32] N. Prokofev, B. Svistunov, and I. Tupitsyn, Worm algorithm in quantum Monte Carlo simulations, Phys. Lett. A 238(4–5), 253 (1998)

[33] N. Prokofev and B. Svistunov, Worm algorithms for classical statistical models, Phys. Rev. Lett. 87(16), 160601 (2001)

[34] P. P. Nidras, Grand canonical simulations of the interacting self-avoiding walk model, J. Phys. Math. Gen. 29(24), 7929 (1996)

[35] N. Madras and A. D. Sokal, The Pivot algorithm: A highly efficient Monte Carlo method for the self-avoiding walk, J. Stat. Phys. 50(1–2), 109 (1988)

[36] P. Grassberger, Pruned-enriched Rosenbluth method: Simulation of q-polymers of chain length up to 1000 000, Phys. Rev. E 56(3), 3682 (1997)

Rights & permissions

Higher Education Press and Springer-Verlag Berlin Heidelberg 2017