TAUP 1920-91 December 1991

Cluster Algorithm for a Solid-On-Solid Model with Constraints

Martin Hasenbusch^{1}, Gideon Lana^{2}, Mihail Marcu^{2,3} and Klaus Pinn^{4}

1 Fachbereich Physik, Universität Kaiserslautern, Postfach 3049, D-6750 Kaiserslautern, Germany

2 School of Physics and Astronomy,

Raymond and Beverly Sackler Faculty of Exact Sciences, Tel Aviv University, 69978 Tel Aviv, Israel

3 II. Institut für Theoretische Physik, Universität Hamburg, Luruper Chaussee 149, D-2000 Hamburg 50, Germany

4 Institut für Theoretische Physik I, Universität Münster, Wilhelm-Klemm-Str. 9, D-4400 Münster, Germany

Abstract

We adapt the VMR (valleys-to-mountains reflections) algorithm, originally devised by us for simulations of SOS models, to the BCSOS model. It is the first time that a cluster algorithm is used for a model with constraints. The performance of this new algorithm is studied in detail in both phases of the model, including a finite size scaling analysis of the autocorrelations.

PACS numbers: 68.35.-p, 68.35.Rh, 64.60.Fr, 64.60.Ht, 05.50.+q

Introduction

Cluster algorithms are one of the most promising new ideas to overcome the problem of critical slowing down (CSD) in computer simulations of statistical physics models [1, 2]. Each time they were sucessful, the dramatic improvement in performance came by integrating into the algorithm the relevant physical properties of the model under study. In the case of solid-on-solid (SOS) models [3], we devised a new cluster algorithm [4] that is based on reflecting “valleys” and “mountains” of the fluctuating SOS surface. We shall call this algorithm the VMR algorithm (valleys-to-mountains-reflections).

Up to now, no cluster algorithms were devised for models with constraints on the configurations. In this letter we present a first example where this is achieved: the adaptation of the VMR algorithm to the BCSOS model.

The BCSOS model is defined as follows [5, 3]. On each site i of a square lattice there is an integer-valued variable h_{i}, which can be thought of as the height of a surface above the two-dimensional lattice. Nearest neighbor variables are constrained to a height difference of ±1. Next to nearest neighbors (i.e. diagonal neighbors) can thus have a height difference of ±2 or 0, and the ±2 is assigned the smaller Boltzmann weight. The partition function is then

Z = ∑_{h}′ exp (−\frac{1}{4}K ∑_{[i,j]}(h_{i} − h_{j})^{2}) ,  (1)

where [i, j] are pairs of diagonal neighbors (∑′ is the constrained sum). The model undergoes a roughening transition at K_{r} = ln 2. At large K the SOS surface is smooth, at small K it is rough. The BCSOS model is equivalent to the F-model [3, 6, 7], which is a particular case of the 6-vertex model (the coupling K in eq. (1) is chosen such as to conform with standard 6-vertex notations).

The idea for the VMR algorithm came from the picture of the SOS configurations as landscapes with mountains and valleys. We first choose a horizontal reflection plane. The clusters are, roughly, the connected regions above the plane (mountains) and below it (valleys). Large scale changes of a configuration are performed by “flipping” valleys or mountains, i.e. by reflecting them through the plane independently and with an appropriate probability. Then a new reflection plane is chosen, the reflection procedure repeated, and so on.

In one way or another the reflection idea had turned up in other cluster algorithms [8, 9, 2]. In the case of the VMR algorithm for SOS models however, the crucial ingredient in eliminating CSD altogether turned out to be not only the reflection, but the precise procedure for choosing the reflection plane [4]. It is at this stage that the particular features of the SOS models were incorporated into the algorithm.

We shall see that the VMR algorithm for the BCSOS model eliminates CSD only half-way. It is nevertheless of considerable algorithmic significance to achieve even that much for a model with constraints. Besides, the improvement in performance was enough for the purpose of a high precision determination of the roughening transition in several SOS models [10]. In that study the renormalization group flows of the SOS models were matched to that of the critical trajectory of the BCSOS model.

The VMR algorithm for the BCSOS model

We shall describe a “single cluster flip” algorithm [9, 2]. A step consists in growing a cluster around a seed, then flipping it. The seed is a randomly chosen site i_{0}.

For a given seed i_{0}, the first thing is to choose an integer reflection plane M. We tried two choices of M. One is to set M = h_{i_{0}} ± 1, the +1 or −1 being chosen randomly with probability one half. The other is to choose another random site j_{0}, and set M = h_{j_{0}}. Both choices gave the same dynamical exponent z (see definition below). The smallest autocorrelation times were obtained with the second procedure, in the case that j_{0} is taken on the sublattice that does not contain i_{0} (we have an even and an odd sublattice; if on the even sublattice all heights are even integers, then on the odd sublattice they are odd integers, and vice versa).

Starting from the seed, which is the first site in the cluster, we freeze and delete links with probabilities to be given below. We have to consider both nn links (nearest neighbor) and nnn links (next to nearest neighbor). A new site is taken into the cluster if it is connected by a frozen link to a site already in the cluster. We repeat the freeze/delete procedure until all links pointing from sites in the cluster either end in the cluster or are deleted. After the cluster is finished, we reflect it (flip it) through the reflection plane M:

h_{i} → 2M − h_{i}  (2)

for all sites i in the cluster.

The freeze/delete probabilities are defined as follows (for each link the sum of the freeze and delete probability is one). For the nn link < i, j > we have two situations. If either h_{i} = M or h_{j} = M, the link is always deleted. Otherwise (i.e. if both sites are on the same side of the reflection plane) the link is always frozen. Thus it is ensured that the constraint on the nearest neighbors is not violated by the reflection.

As a consequence of the freeze/delete procedure for nn links, the nnn link [i, j] is automatically frozen in all but three situations (remember that |h_{i} − h_{j}| is zero or two). First, if h_{i} and h_{j} are on different sides of the reflection plane, the nnn link [i, j] is always deleted. Second, if either h_{i} = M or h_{j} = M, [i, j] is also deleted with probability one. Third, if h_{i} = h_{j} = M ± 1, the nnn link [i, j] is deleted with probability exp(−K). In eq. (1) we wrote the nnn interaction in the form of the discrete Gaussian model. It can be easily seen that the delete probabilities used for the nnn links are indeed the same as those defined in [4] for the discrete Gaussian model without constraints.

Notice that for the case of the BCSOS model, the clusters are precisely valleys or mountains cut by the reflection plane from the SOS surface. In the case of the discrete Gaussian model [4], this was only approximately so.

Autocorrelations, sublattice effects

We shall call “autocorrelation time” τ the quantity often referred to in the literature [2] as the “exponential autocorrelation time”. At large “Monte Carlo times” t the autocorrelation function of some physical quantity decays exponentially as exp(−t/τ), which is related to the slowest mode in the Markov process used for the simulation.

In practice, for some physical quantities the autocorrelations immediately enter the exp(−t/τ) regime, while for other quantities this only happens after the autocorrelations have decayed much faster for a while. Of course, the theoretical expectation is that at very large t the autocorrelations decay with the same τ for all quantities (except if prevented by some conservation law).

For the BCSOS model we found that quantities defined on one sublattice tended to behave as exp(−t/τ) already for t of the order of a few to a few tens of cluster flips. On the other hand, averages over the two sublattices tended to enter this “clean exponential” regime only after the autocorrelation functions had decayed over several orders of magnitude.

As an example we discuss here the energy density, which is, up to a constant factor, < V^{−1} ∑_{[i,j]}(h_{i} − h_{j})^{2} >, with V the number of points on one sublattice, and the sum over nnn pairs [i, j] in that sublattice. The values of the energy density on the two sublattices are so strongly anticorrelated that the statistical error of the average over the two sublattices is a whole order of magnitude less than the error of the one-sublattice quantity.

These strong anticorrelations can be understood by regarding the four points at the corners of an elementary plaquette of the lattice. In the case that the height difference of the nnn variables at one pair of diagonally opposite sites is two, the constraint on the nn pairs implies that for the other pair of diagonally opposite sites the height difference is necessarily zero. The anticorrelations stem from the fact that the two nnn pairs discussed here lie on different sublattices.

In what follows, we shall present the results for the autocorrelation time of the one-sublattice energy. At each value of K and for each lattice size quoted, the statistics was between 250000 and 800000 clusters (vectorization of the cluster algorithm [11] was of great help). We used periodic boundary conditions in the directions of the coordinate axes of the sublattices (these axes are rotated by 45^{∘} with respect to the coordinate axes of the original lattice). Each sublattice has a size of L^{2}, so our lattice has volume 2L^{2}. For our study we took L to be a power of 2, and 8 ≤ L ≤ 256 (in the rough phase the largest L was 128).

As opposed to the one-sublattice energy, for the energy averaged over both sublattices the exp(−t/τ) regime is not reached before the autocorrelations become zero within errorbars, if the statistics is a few hundred thousand of clusters. In some runs with ten times as many clusters we checked that the energy averaged over the two sublattices does indeed have the same τ.

The values for τ we shall quote are in units of sweeps, i.e. we compute τ in units of clusters first, then multiply it by the average cluster size and divide it by 2L^{2}. One unit then amounts to a work proportional to the volume. As opposed to other cluster algorithms as e.g. the single cluster algorithm for the Ising model [9, 2] the average cluster size is quite stable. If we choose the reflection plane to be the height at a random site on the sublattice not containing the seed, the average cluster size is around 0.35 of the volume.

Autocorrelation study in the rough phase and at K_{r}

Since the SOS surface thickness increases only logarithmically with the lattice size, the SOS surface has in reality shallow valleys and small mountains. As the temperature increases (K decreases) the valleys become deeper and the mountains higher. This leads us to the expectation that the higher the temperature, the better our algorithm should perform.

We studied the autocorrelations at the roughening transition K_{r} = ln 2 and at the point K = K_{r}/2 which is deep in the rough phase. Since in the rough phase the correlation length is infinite, the dynamical critical exponent z for CSD is defined by the relation [2]

τ ∝ L^{z} .  (3)

The results are shown in figs. 1 and 2. Fig. 1 suggests that at K_{r} eq. (3) is well fulfilled even for small L. The best fit of our data gives

τ(K_{r}) = 1.20 ± 0.02 .  (4)

At K_{r}/2, fig. 2 suggests that the regime of eq. (3) was not yet reached, and that by increasing L further the fitted value of τ would decrease. We conclude that

τ(K_{r}/2) < 0.79 ± 0.09 ,  (5)

the number on the r.h.s. being the fitted value for 32 ≤ L ≤ 128. Thus the algorithm does indeed perform better at higher temperatures.

In [4] we analyzed the VMR algorithm for the Discrete Gaussian model. The BCSOS algorithm presented here performs similarly to the H-algorithm discussed there. Because of the constraints, we did not find a way to include the analogue of the I-algorithm that further improved the performance in the case of the Discrete Gaussian model.

Figure 1: Log-log plot of τ vs. L at K = K_{r}. The solid line is the fit with eq. (3) for 8 ≤ L ≤ 128.

Figure 2: Log-log plot of τ vs. L at K = K_{r}/2. The solid line is the fit with eq. (3) for 32 ≤ L ≤ 128.

Table 1: The exponent z for fixed L/ξ in the smooth phase.

\begin{matrix}L/ξ & 0 & 0.25 & 0.5 & 1 & 2 & 4 & 8 \\ z & 1.20(2) & 1.08(4) & 1.02(5) & 0.90(4) & 0.91(4) & 0.86(5) & 0.77(9) \end{matrix}

Finite size scaling of the autocorrelation time in the smooth phase

In [12] it is conjectured that the autocorrelation time obeys the usual finite size scaling law

τ = L^{z} f(L/ξ) ,  (6)

where f is some smooth function. We therefore took data for fixed ratios L/ξ. Fortunately, there is an exact formula for the correlation length of the BCSOS model (ref. [7], eq. 8.11.24). We made runs for L/ξ = 0.25, 0.5, 1, 2, 4 and 8, with L chosen to be a power of two, 8 ≤ L ≤ 256.

For each value of L/ξ, we fitted the values of τ, obtained from the autocorrelations of the one-sublattice energy, with eq. (3). In each case the fits were good, but, as opposed to eq. (6), the value of z was not constant. In table 1 we give the fitted values of z as a function of L/ξ. We see that z slowly increases as L/ξ decreases, and this increase seems consistent with the value z = 1.20(2) at the transition point (L/ξ = 0).

This finding is compatible with a scaling law τ ∝ ξ^{z} in the limit L → ∞, but with an exponent z smaller than 1.20(2). In all cases investigated up to now the two exponents were equal, so such a situation would be very interesting. We stress, however, that we do not have data close enough to the thermodynamic limit in order to claim that the law τ ∝ ξ^{z} was well checked.

In some cases logarithmic corrections to the finite size scaling variable L/ξ turn out to be important [13]. We tried for the same data for τ in the smooth phase the ansatz

τ = L^{z} f(\frac{L}{ξ(ln ξ)^{a}} ) ,  (7)

where we took z = 1.20. We looked for a value of a for which the values of τ/L^{z} plotted against L/(ξ(ln ξ)^{a}) collapse onto one curve. In fig. 3 we show the result for a = −0.7. The range of a for which reasonable collapse occured was −0.9 ≤ a ≤ −0.6. Thus, within our precision, this modified finite size scaling is consistent with all data and with the value z = 1.20. We caution however that a collapse of the data of roughly the same quality can be found for 1.05 ≤ z ≤ 1.25.

We do not claim that eq. (7) is the “correct” finite size scaling law, but it is certainly closer to the truth than eq. (6).

Figure 3: Finite size scaling plot with z = 1.20 and a = −0.7. Notice the logarithmic axes. Different marker symbols are used for different values of L/ξ.

Conclusions and outlook

We have presented for the first time a successful cluster algorithm for a model with constraints. We carefully measured and analyzed the performance of our algorithm. We also performed a finite size scaling analysis for the autocorrelation time.

We would like to point out that a completely different algorithm for vertex models was also developed recently [14]. Another research direction we pursue is the generalization of the VMR algorithm to more complicated constraints, which may be of relevance for comparing to interface experiments.

While further algorithmic studies are interesting by themselves, we point out that the algorithm described in this letter was used in the (to date) most precise determination of the roughening (Kosterlitz-Thouless) transition for a variety of SOS models (including the dual of the XY model) [10].

Finally, we point out that our cluster algorithm can be used, in slightly modified form, for accelerating quantum Monte Carlo simulations for XXZ chains [15].

Acknowledgements

This work was supported in part by the German-Israeli Foundation for Research and Development (GIF), by the Basic Research Foundation of The Israel Academy of Sciences and Humanities, and by the Deutsche Forschungsgemeinschaft. We would like to express our gratitude to the HLRZ in Jülich, where most of our computer runs were performed.

References

[1] R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58 (1987) 86.

[2] For reviews on cluster algorithms see e.g.:

U. Wolff, in Lattice ’89 , Capri 1989, N. Cabbibo et al. editors, Nucl. Phys. B (Proc. Suppl.) 17 (1990) 93; D. Kandel and E. Domany, Phys. Rev. B43 (1991) 8539; A. Sokal, in Lattice ’90 , Tallahassee 1990, U. M. Heller, A. D. Kennedy and S. Sanielevici editors, Nucl. Phys. B (Proc. Suppl.) 20 (1991).

[3] For reviews on SOS type of models see e.g.:

D. B. Abraham, Surface Structures and Phase Transitions – Exact Results, in Phase Transitions and Critical Phenomena Vol. 10, C. Domb and J. L. Leibowitz editors (Academic, 1986) p. 1; H. van Beijeren and I. Nolden, The Roughening Transition, in Topics in Current Physics, Vol. 43: Structure and Dynamics of Surfaces II , W. Schommers and P. von Blanckenhagen editors (Springer, 1987) p. 259.

[4] H. G. Evertz, M. Hasenbusch, M. Marcu, K. Pinn and S. Solomon, Phys. Lett. 254B (1991) 185, in Lattice 90 (see above) p. 80, and in Workshop on Fermion Algorithms, Jülich 1991, H. J. Herrmann and F. Karsch editors, Int. J. Mod. Phys. C3 (1992) 235.

[5] H. van Beijeren, Phys. Rev. Lett. 38 (1977) 993.

[6] E. H. Lieb, Phys. Rev. Lett. 18 (1967) 1046; E. H. Lieb and F. Y. Wu, Two-dimensional Ferroelectric Models, in Phase Transitions and Critical Phenomena Vol. 1, C. Domb and M. S. Green editors (Academic, 1972) p. 331.

[7] R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic, 1982).

[8] R. C. Brower and P. Tamayo, Phys. Rev. Lett. 62 (1989) 1087.

[9] U. Wolff, Phys. Rev. Lett. 62 (1989) 361, Nucl. Phys. B322 (1989) 759, and Phys. Lett. 228B (1989) 379.

[10] M. Hasenbusch, M. Marcu and K. Pinn, High precision verification of the Kosterlitz-Thouless scenario, preprint DESY 92-010, to appear in the proceedings of Lattice 91 Tsukuba (November 1991).

[11] H. G. Evertz, preprint FSU-SCRI-91-183.

[12] R. G. Edwards, J. Goodman and A. D. Sokal, Nucl. Phys. B354 (1991) 289.

[13] E. Brézin and J. Zinn-Justin, Nucl. Phys. B257 (1985) 867.

[14] H. G. Evertz, G. Lana and M. Marcu, in preparation.

[15] M. Suzuki (editor), Quantum Monte Carlo methods in equilibrium and nonequilibrium systems, Taniguchi symposium, Springer (1987).