- Email: [email protected]
Finite-size calculations on the three-dimensional Ising model
Physica A 196 (1993) 455-460 North-Holland SDI:
0378-4371(93)E0035-D
Finite-size calculations Ising model
on the three-dimensional
H.W.J. Bl6te and G. Kamieniarz’ Laboratorium voor Technische Natuurkunde, 2600 GA Delft, The Netherlands
Technische
Universiteit Delft, P. 0.
Box 5046,
Received 9 February 1993
Recent numerical results for the critical behaviour of the three-dimensional Ising model, based on finite-size calculations, show significant discrepancies. We review possible sources of systematic deviations and strategies to avoid these or to track them down. We present some new finite-size results for the Binder cumulant, obtained by means of a perturbation expansion and by cluster Monte Carlo methods.
The three-dimensional
Ising model continues to be a subject of investigations, in particular its critical behaviour. Whereas the most accurate results traditionally came from series expansions, (see e.g. refs. [l-3] and references therein) and from the c-expansion [4,5], the error margins quoted in Monte Carlo based approaches are shrinking considerably [6-181. Due to the availability of fast computers, and the invention of ingenious cluster algorithms by Swendsen and Wang [19] and Wolff [20], such statistical accuracy is now becoming available that the results are comparable with series and &-expansions. However, the numerical results of such simulations are not always in a satisfactory agreement with one another. For instance, the critical points of the simple cubic Ising model as found by Baillie et al. [18] (K, = 0.221652 (3)) and by Livet [17] (K, = 0.2216544 (10)) are not inside the error margins of a recent result by Ferrenberg and Landau [16] (K, = 0.2216595 (26)). Furthermore, the critical exponent Y obtained by Baillie et al. does not agree satisfactorily with others. Although these differences are not large, they invite the consideration of possible sources of systematic errors. In our own analysis, we made efforts to avoid, where possible, the following ones: 1. Effects due to a poor random generator; ’ Permanent address: Institute of Physics, A. Mickiewicz University, Poznan, Poland; present address: Institute of Theoretical Physics, University of Leuven, Celestijnenlaan 2OOD, Leuven, Belgium. 0378-4371193 /$06.00 0 1993 - Elsevier Science Publishers B .V. All rights reserved
456
H. W. J. Bhte,
G. Kamieniarz
I Finite-size
calculations
on 30
Ising model
2. Effects introduced by histogramming methods; 3. The use of a single row of pseudo-random numbers for a number of parallel simulations; 4. Neglecting corrections to scaling in fitting procedures of simulation results. A significant source of errors is undoubtedly associated with the randomnumber generators used during the simulations (see e.g. refs. [21, S-111). It is well known that the linear congruential method with a period of about 10’ (using 32-bit integers) is inadequate for long simulations. But also a number of algorithms based on 127-bit feedback shift registers with a period in the order of 1O38 have been rejected on the basis of long tests [22]. Recent tests by Ferrenberg et al. [23] have shown that cluster algorithms are very susceptible to systematic deviations due to shift-register based random generators. A practical difficulty to determine such systematic errors is obvious: once the statistical accuracies involved become better than those of all other existing data, there is nothing left to compare the results with, at least in the absence of exact solutions. This would suggest to repeat at least a part of the simulations with different realizations of the random-number generator or the spinupdating algorithm. Another possibility to test Monte Carlo algorithms is offered by numerical exact results for small systems, so that one has strict criteria that the Monte Carlo results should satisfy. Problems resulting from the use of histogram methods were mentioned by M.E. Fisher in a recent address to the statistical physics community [24]. Whereas histograms can be valuable tools in the analysis of Monte Carlo results, it is quite easy to underestimate the statistical errors or to introduce deviations as a consequence of the finite resolution of the histogram. Since the generation of pseudo-random numbers consumes a considerable fraction of the total time needed for a simulation, it is tempting to use a single random number for a number of spin updates in parallel simulations [25,15]. However, the results of the parallel simulations will then be correlated and may not be considered independent. The same effect may occur if different runs use the same seed for the random-number generator. Neglecting corrections to scaling during a fitting procedure may lead to a bias in the results, and to a gross underestimation of the statistical errors in the fitted parameters. It is obvious that assumptions concerning the form of the fitted expression are necessary. These assumptions are only verifyable up to the limits imposed by the statistical precision. Even if the fit satisfies the chi-square criterion, there is no certainty that additional corrections to scaling are absent; the error in the fitted parameters may thus be larger than that predicted by the fit. The need to run additional tests on the random-number generator used
H.W. .I. B&e,
G. Kamieniarz
I Finite-size calculations on 30
Ising model
457
during our simulations of the 3D Ising model became apparent when we noticed systematic differences between our results for the Binder cumulant [26] and those of Livet 1171. These results, which apply to simple cubic lattices with periodic boundary conditions, are shown in table I in terms of the amplitude ratio g, = (m4)L/(m2)t, where m is the magnetization and L the linear system size. Both sets of data were generated using the Swendsen-Wang cluster method. However, the coupling strength K = 0.221656 used in ref. [17] is slightly different from K = 0.221653 used in our simulations. Therefore we have used finite-size scaling to correct for the linear dependence of g, on the coupling strength in the immediate neighbourhood of the critical point. Writing g, in terms of field derivatives of the free energy, a scale transformation by a factor L yields
gl’(K,) = gZ’(K,)+ q(K, - K,)Ly’ , Table I Finite-size results for the universal ratio g, of the three-dimensional Ising model at a coupling K = 0.221656. The first column shows the finite size L, and the second one our Monte Carlo data with standard errors in the last decimal place between parentheses. The next column shows the length of the simulations. This length expresses the total number of cluster configurations. Most of these were generated by the Swendsen-Wang algorithm, but also included are 5.2, 1.6, 2.0, 1.2, 1.2, 0.8, 0.2, 0.2, and 0.8 X 10’ configurations for system sizes L = 3, 5, 12, 13, 14, 18, 22, 24, and 32 respectively, that were obtained by means of the largest-cluster method. The fourth column lists our exact numerical results, and the last one contains, for comparison, part of the numerical results of ref. [17]. L
g, (present)
# conf.
gL (exact)
g, (ref. P71)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 20 22 24 28 32
1.49598 (8) 1.51570 (8) 1.52956 (14) 1.54032 (19) 1.5483 (2) 1.5545 (2) 1.5597 (2) 1.5636 (2) 1.5666 (2) 1.5696 (2) 1.5722 (3) 1.5740 (3) 1.5756 (4) 1.5766 (4) 1.5790 (4) 1.5825 (4) 1.5837 (4) 1.5857 (7) 1.5868 (7) 1.5875 (7)
1.0 x los 1.6 x lo8 6.4 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.0 x 10’ 4.0 x 10’ 2.0 x 10’ 2.0 x 10’ 2.0 x 10’ 2.0 x 10’ 2.0 x 10’ 1.0 x 10’ 1.0 x 10’ 1.0 x 10’
1.496016 1.515605
1.5167 (6) 1.5431(6) 1.5572 (9)
1.5709 (14)
1.5791(13)
1.5820 (25) 1.5943 (28)
458
H.W.J.
Bliite, G. Kamieniarz
I Finite-size calculations on 30
Ising model
with (Ye= 0.86 and yt = 1.59, as determined from fits to data taken farther away from Kc. Our results appear to be systematically lower than those of ref. [17]. The random-number generator used in our algorithm has been selected with some care, in view of past experiences in our group [21]. Such problems can be avoided by using very long feedback shift registers [27] but here we have adopted a different method. It works as follows, The row aj = ai_92,9 @ 4i_9689, where the ai are 32-bit integers and @ stands for bitwise modulo-2 addition, is in fact composed of 32 independent feedback shift registers [28]. The row b, = b,_, X 32781, truncated to 32 bits, is a simple multiplicative random generator. Both rows are known to contain undesired correlations, but of different types. Thus we have used the row ri = a, @ bj in which we expect that these correlation effects are suppressed. Indeed, we did not find any significant differences when we replaced the shift register by a shorter one [28] of length 4423 and feedback position 1393 (runs of 8 million Swendsen-Wang steps for system size L = 20, and 10 million steps for L = 22). Neither did we find deviations when we applied instead the largest-cluster [29] variant of the Swendsen-Wang algorithm (see table I). Ref. [17] does not specify the random generator used. As a further test, we have computed g, independently for 33 and 43 systems with periodic boundaries. Exact numerical calculation of g, for the smaller system hardly poses any problems, but the 43 system is already too timeconsuming unless treated carefully. A way to tackle this problem was shown by Pearson [30], but we found that a perturbation expansion as described by Saleur and Derrida [31] is even more efficient. The idea is to write the partition sum as Z = Tr T4, where T is the 65536 x 65536 transfer matrix, to expand T in powers of the field h, and to keep track of powers of h during the matrix multiplications. Thus one obtains the second and fourth derivatives of Z to h, in which g, can be expressed. It is fortunate that the multiplications can be performed using a sparse-matrix decomposition [32], and that out of the 65536 terms in the trace only 433 are independent. The rest follows by the 16-fold translation symmetry, the 4-fold rotation symmetry, the spatial inversion symmetry and the spin-inversion symmetry of a periodic 4 X 4 square of Ising spins. For more details of the calculations, see also ref. [33]. The results for L = 3 and L = 4 arc included in table I, and are in a satisfactory agreement with our Monte Carlo results. For L = 4, the difference with ref. [17] is 1.8 standard deviations as listed there, and 2.4 standard deviations according to our determination of the statistical error for a run of the quoted number of Swendsen-Wang steps (5 x 106). For this determination we divided runs of approximately that length in 1000 partial results, and performed an analysis on the partial averages. The new data in table I are restricted to rather small system sizes, but have a
H.W.J.
B&e,
G. Kamieniarz
I Finite-size calculations on 30
Ising model
459
greater statistical accuracy than other results known to us for these systems. We have analyzed our data according to g,(K)
= g + bLY’ + a,(K - KJLY’ + a,(K -
KJ2L2Y+ * **
)
(2)
with yt = 1.59 [l-17], in order to determine KC, yi, and g. A least-squares fit, which included some additional data points for KZ0.221653 in order to find the coefficients ai in eq. (2), showed that the data point for I!. = 4 deviates significantly, and the fits thus apply to L 3 5. We found that g = 1.609 + 0.004 and yi = -0.85 +- 0.04, in a good agreement with a recent series result of Nickel and Rehr [3]: yi = -0.83 + 0.05. If we would assume a fixed value yi = -0.8, as was done by Livet [17], we would obtain KC = 0.2216445 (14) for the critical point. Including yi as a variable, our final estimate of the critical point is KC = 0.221648 (4). If the critical point is indeed within these error bounds, Rosengren’s conjecture [34] tanhK, = (fi - 2) COS(IT/~) or Kc = 0.2216586 is incorrect. Our result for KC is more than two standard errors smaller than that of Ferrenberg and Landau [16], and three or more standard errors larger than values listed by Liu and Fisher [l]. Therefore we briefly review the uncertainties in our result for Kc with regard to the sources of error mentioned above. The fact that our random-number generator has passed a number of tests does not prove that our data in table I are free of systematic errors for L 3 5. However, if there are such errors, they will be clearly recognizable as such in the future. Furthermore, additional corrections to scaling may very well play a role. We have also made fits to our data with a term b, L2” included in eq. (2), describing a nonlinearity in the irrelevant field. The result was that the irrelevant exponent yi became poorly determined. Using the value found by Nickel and Rehr [3] as input leads again to KC = 0.221648 (4), where the error margin includes the uncertainty in yi. It is, of course, possible that another correction with a different exponent is significant. In that case, it will undoubtedly reveal itself when accurate data over a wider range of systems sizes become available. In order to gain insight in the discrepancies as mentioned above, it is essential that verifiable numerical data are published, as was done by Livet [17]. Likewise, it is important to specify the details of the random-number generator. This holds especially in the case of the three-dimensional Ising model, where we are confronted with an increasing body of mutually conflicting results. It is a pleasure to thank A.M. Ferrenberg, J.R. Heringa, A. Kooiman and M. Novotny for valuable discussions. We are much indebted to H.J.M. van
460
H. W. J. Bltite, G. Kamieniarz 1 Finite-size calculations on 30 Ising model
Grol, Y.T.J.C. Fonk and M.J. Neefs for teaching us how to operate workstations. This work is part of the research program of the “Stichting voor Fundamenteel Onderzoek der Materie (FOM)“, which is financially supported by the “Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)“.
References [l] [2] [3] [4] [5] [6] [7] [8] [9] [lo] [ll] [12]
A.J Liu and M.E. Fisher, Physica A 156 (1989) 35. J. Adler, J. Phys. A 16 (1983) 3585. B.G. Nickel and J.J. Rehr, J. Stat. Phys. 61 (1990) 1. G.A. Baker, B.G. Nickel, M.S. Green and D.I. Meiron, Phys Rev. B 17 (1978) 1365. J.C. Le Guillou and J. Zinn-Justin, Phys. Rev B 21 (1980) 3976. H.W.J. Blijte and R.H. Swendsen, Phys. Rev. B 20 (1979) 2077. G.S. Pawley, R.H. Swendsen, D.J. Wallace and K.G. Wilson, Phys. Rev. B 29 (1984) 4030. M.N. Barber, R.B. Pearson, D. Toussaint and J.L. Richardson, Phys. Rev. B 32 (1985) 1720. G. Parisi and F. Rapuano, Phys. Lett. B 157 (1985) 301. A. Hoogland, A. Compagner and H.W.J. Bliite, Physica A 132 (1985) 593. G. Bhanot, D. Duke and R. Salvador, Phys. Rev. B 33 (1986) 7841. H.W.J. Blote, J.A. de Bruin, A. Compagner, J.H. Croockewit, Y.T.J.C. Fonk, J.R. Heringa, A. Hoogland and A.L. van Willigen, Europhys. Lett. 10 (1989) 105. (131 H.W.J. Bliite, A. Compagner, J.H. Croockewit, Y.T.J.C. Fonk, J.R. Heringa, A. Hoogland, T.S. Smit and A.L. van Willigen, Phyica A 161 (1989) 1. [14] N.A. Alves, B.A. Berg and R. Villanova, Phys. Rev. B 41 (1990) 383. [15] N. Ito and M. Suzuki, J. Phys. Sot. Jpn. 60 (1991) 1978. [16] A.M. Ferrenberg and D.P. Landau, Phys. Rev. B 44 (1991) 5081. [17] F. Livet, Europhys. Lett. 16 (1991) 139. (181 C.F. Baillie, R. Gupta, K.A. Hawick and G.S. Pawley, Phys. Rev. B 45 (1992) 10438. [19] R.H. Swendsen and J.S. Wang, Phys. Rev. Lett. 58 (1987) 86. [20] U. Wolff, Phys. Rev. Lett. 60 (1988) 1461. [21] A. Hoogland, J. Spaa, B. Selman and A. Compagner, J. Comput. Phys. 51 (1983) 250. [22] A. Hoogland, A. Compagner and H.W.J. Blote, in: Architecture and Performance of Specialized Computers, B. Alder, ed., Computational Techniques (Academic Press, New York, 1988) Ch. 7, The Delft Ising System Processor. [23] A.M. Ferrenberg, D.P. Landau and Y.J. Wang, Phys. Rev. Lett. 69 (1992) 3382. [24] M.E. Fisher, comment on a lecture by R.H. Swendsen at the Statphys 18 Conference (Berlin, 1992). [25] N. Ito and Y. Kanada, Supercomputer 5 (3) (1988) 31. [26] K. Binder, Z. Phys. B 43 (1981) 119. [27] J.R. Heringa, H.W.J. Blote and A. Compagner, Int. J. Mod. Phys. C 3 (1992) 561. [28] N. Zierler, Inform. Control 15 (1969) 67. [29] C.F. Baillie and P.D. Coddington, Phys. Rev. B 43 (1991) 10617. [30] R.B. Pearson, Phys. Rev. B 26 (1982) 6285. [31] H. Saleur and B. Derrida, J. Phys. (Paris) 46 (1985) 1043. [32] M.P. Nightingale, Proc. K. Ned. Akad. Wet. B 82 (1979) 23.5. [33] G. Kamieniarz and H.W.J. Bliite, J. Phys. A 26 (1993) 201. [34] A. Rosengren, J. Phys. A 19 (1986) 1709.
0378-4371(93)E0035-D
Finite-size calculations Ising model
on the three-dimensional
H.W.J. Bl6te and G. Kamieniarz’ Laboratorium voor Technische Natuurkunde, 2600 GA Delft, The Netherlands
Technische
Universiteit Delft, P. 0.
Box 5046,
Received 9 February 1993
Recent numerical results for the critical behaviour of the three-dimensional Ising model, based on finite-size calculations, show significant discrepancies. We review possible sources of systematic deviations and strategies to avoid these or to track them down. We present some new finite-size results for the Binder cumulant, obtained by means of a perturbation expansion and by cluster Monte Carlo methods.
The three-dimensional
Ising model continues to be a subject of investigations, in particular its critical behaviour. Whereas the most accurate results traditionally came from series expansions, (see e.g. refs. [l-3] and references therein) and from the c-expansion [4,5], the error margins quoted in Monte Carlo based approaches are shrinking considerably [6-181. Due to the availability of fast computers, and the invention of ingenious cluster algorithms by Swendsen and Wang [19] and Wolff [20], such statistical accuracy is now becoming available that the results are comparable with series and &-expansions. However, the numerical results of such simulations are not always in a satisfactory agreement with one another. For instance, the critical points of the simple cubic Ising model as found by Baillie et al. [18] (K, = 0.221652 (3)) and by Livet [17] (K, = 0.2216544 (10)) are not inside the error margins of a recent result by Ferrenberg and Landau [16] (K, = 0.2216595 (26)). Furthermore, the critical exponent Y obtained by Baillie et al. does not agree satisfactorily with others. Although these differences are not large, they invite the consideration of possible sources of systematic errors. In our own analysis, we made efforts to avoid, where possible, the following ones: 1. Effects due to a poor random generator; ’ Permanent address: Institute of Physics, A. Mickiewicz University, Poznan, Poland; present address: Institute of Theoretical Physics, University of Leuven, Celestijnenlaan 2OOD, Leuven, Belgium. 0378-4371193 /$06.00 0 1993 - Elsevier Science Publishers B .V. All rights reserved
456
H. W. J. Bhte,
G. Kamieniarz
I Finite-size
calculations
on 30
Ising model
2. Effects introduced by histogramming methods; 3. The use of a single row of pseudo-random numbers for a number of parallel simulations; 4. Neglecting corrections to scaling in fitting procedures of simulation results. A significant source of errors is undoubtedly associated with the randomnumber generators used during the simulations (see e.g. refs. [21, S-111). It is well known that the linear congruential method with a period of about 10’ (using 32-bit integers) is inadequate for long simulations. But also a number of algorithms based on 127-bit feedback shift registers with a period in the order of 1O38 have been rejected on the basis of long tests [22]. Recent tests by Ferrenberg et al. [23] have shown that cluster algorithms are very susceptible to systematic deviations due to shift-register based random generators. A practical difficulty to determine such systematic errors is obvious: once the statistical accuracies involved become better than those of all other existing data, there is nothing left to compare the results with, at least in the absence of exact solutions. This would suggest to repeat at least a part of the simulations with different realizations of the random-number generator or the spinupdating algorithm. Another possibility to test Monte Carlo algorithms is offered by numerical exact results for small systems, so that one has strict criteria that the Monte Carlo results should satisfy. Problems resulting from the use of histogram methods were mentioned by M.E. Fisher in a recent address to the statistical physics community [24]. Whereas histograms can be valuable tools in the analysis of Monte Carlo results, it is quite easy to underestimate the statistical errors or to introduce deviations as a consequence of the finite resolution of the histogram. Since the generation of pseudo-random numbers consumes a considerable fraction of the total time needed for a simulation, it is tempting to use a single random number for a number of spin updates in parallel simulations [25,15]. However, the results of the parallel simulations will then be correlated and may not be considered independent. The same effect may occur if different runs use the same seed for the random-number generator. Neglecting corrections to scaling during a fitting procedure may lead to a bias in the results, and to a gross underestimation of the statistical errors in the fitted parameters. It is obvious that assumptions concerning the form of the fitted expression are necessary. These assumptions are only verifyable up to the limits imposed by the statistical precision. Even if the fit satisfies the chi-square criterion, there is no certainty that additional corrections to scaling are absent; the error in the fitted parameters may thus be larger than that predicted by the fit. The need to run additional tests on the random-number generator used
H.W. .I. B&e,
G. Kamieniarz
I Finite-size calculations on 30
Ising model
457
during our simulations of the 3D Ising model became apparent when we noticed systematic differences between our results for the Binder cumulant [26] and those of Livet 1171. These results, which apply to simple cubic lattices with periodic boundary conditions, are shown in table I in terms of the amplitude ratio g, = (m4)L/(m2)t, where m is the magnetization and L the linear system size. Both sets of data were generated using the Swendsen-Wang cluster method. However, the coupling strength K = 0.221656 used in ref. [17] is slightly different from K = 0.221653 used in our simulations. Therefore we have used finite-size scaling to correct for the linear dependence of g, on the coupling strength in the immediate neighbourhood of the critical point. Writing g, in terms of field derivatives of the free energy, a scale transformation by a factor L yields
gl’(K,) = gZ’(K,)+ q(K, - K,)Ly’ , Table I Finite-size results for the universal ratio g, of the three-dimensional Ising model at a coupling K = 0.221656. The first column shows the finite size L, and the second one our Monte Carlo data with standard errors in the last decimal place between parentheses. The next column shows the length of the simulations. This length expresses the total number of cluster configurations. Most of these were generated by the Swendsen-Wang algorithm, but also included are 5.2, 1.6, 2.0, 1.2, 1.2, 0.8, 0.2, 0.2, and 0.8 X 10’ configurations for system sizes L = 3, 5, 12, 13, 14, 18, 22, 24, and 32 respectively, that were obtained by means of the largest-cluster method. The fourth column lists our exact numerical results, and the last one contains, for comparison, part of the numerical results of ref. [17]. L
g, (present)
# conf.
gL (exact)
g, (ref. P71)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 20 22 24 28 32
1.49598 (8) 1.51570 (8) 1.52956 (14) 1.54032 (19) 1.5483 (2) 1.5545 (2) 1.5597 (2) 1.5636 (2) 1.5666 (2) 1.5696 (2) 1.5722 (3) 1.5740 (3) 1.5756 (4) 1.5766 (4) 1.5790 (4) 1.5825 (4) 1.5837 (4) 1.5857 (7) 1.5868 (7) 1.5875 (7)
1.0 x los 1.6 x lo8 6.4 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.8 x 10’ 4.0 x 10’ 4.0 x 10’ 2.0 x 10’ 2.0 x 10’ 2.0 x 10’ 2.0 x 10’ 2.0 x 10’ 1.0 x 10’ 1.0 x 10’ 1.0 x 10’
1.496016 1.515605
1.5167 (6) 1.5431(6) 1.5572 (9)
1.5709 (14)
1.5791(13)
1.5820 (25) 1.5943 (28)
458
H.W.J.
Bliite, G. Kamieniarz
I Finite-size calculations on 30
Ising model
with (Ye= 0.86 and yt = 1.59, as determined from fits to data taken farther away from Kc. Our results appear to be systematically lower than those of ref. [17]. The random-number generator used in our algorithm has been selected with some care, in view of past experiences in our group [21]. Such problems can be avoided by using very long feedback shift registers [27] but here we have adopted a different method. It works as follows, The row aj = ai_92,9 @ 4i_9689, where the ai are 32-bit integers and @ stands for bitwise modulo-2 addition, is in fact composed of 32 independent feedback shift registers [28]. The row b, = b,_, X 32781, truncated to 32 bits, is a simple multiplicative random generator. Both rows are known to contain undesired correlations, but of different types. Thus we have used the row ri = a, @ bj in which we expect that these correlation effects are suppressed. Indeed, we did not find any significant differences when we replaced the shift register by a shorter one [28] of length 4423 and feedback position 1393 (runs of 8 million Swendsen-Wang steps for system size L = 20, and 10 million steps for L = 22). Neither did we find deviations when we applied instead the largest-cluster [29] variant of the Swendsen-Wang algorithm (see table I). Ref. [17] does not specify the random generator used. As a further test, we have computed g, independently for 33 and 43 systems with periodic boundaries. Exact numerical calculation of g, for the smaller system hardly poses any problems, but the 43 system is already too timeconsuming unless treated carefully. A way to tackle this problem was shown by Pearson [30], but we found that a perturbation expansion as described by Saleur and Derrida [31] is even more efficient. The idea is to write the partition sum as Z = Tr T4, where T is the 65536 x 65536 transfer matrix, to expand T in powers of the field h, and to keep track of powers of h during the matrix multiplications. Thus one obtains the second and fourth derivatives of Z to h, in which g, can be expressed. It is fortunate that the multiplications can be performed using a sparse-matrix decomposition [32], and that out of the 65536 terms in the trace only 433 are independent. The rest follows by the 16-fold translation symmetry, the 4-fold rotation symmetry, the spatial inversion symmetry and the spin-inversion symmetry of a periodic 4 X 4 square of Ising spins. For more details of the calculations, see also ref. [33]. The results for L = 3 and L = 4 arc included in table I, and are in a satisfactory agreement with our Monte Carlo results. For L = 4, the difference with ref. [17] is 1.8 standard deviations as listed there, and 2.4 standard deviations according to our determination of the statistical error for a run of the quoted number of Swendsen-Wang steps (5 x 106). For this determination we divided runs of approximately that length in 1000 partial results, and performed an analysis on the partial averages. The new data in table I are restricted to rather small system sizes, but have a
H.W.J.
B&e,
G. Kamieniarz
I Finite-size calculations on 30
Ising model
459
greater statistical accuracy than other results known to us for these systems. We have analyzed our data according to g,(K)
= g + bLY’ + a,(K - KJLY’ + a,(K -
KJ2L2Y+ * **
)
(2)
with yt = 1.59 [l-17], in order to determine KC, yi, and g. A least-squares fit, which included some additional data points for KZ0.221653 in order to find the coefficients ai in eq. (2), showed that the data point for I!. = 4 deviates significantly, and the fits thus apply to L 3 5. We found that g = 1.609 + 0.004 and yi = -0.85 +- 0.04, in a good agreement with a recent series result of Nickel and Rehr [3]: yi = -0.83 + 0.05. If we would assume a fixed value yi = -0.8, as was done by Livet [17], we would obtain KC = 0.2216445 (14) for the critical point. Including yi as a variable, our final estimate of the critical point is KC = 0.221648 (4). If the critical point is indeed within these error bounds, Rosengren’s conjecture [34] tanhK, = (fi - 2) COS(IT/~) or Kc = 0.2216586 is incorrect. Our result for KC is more than two standard errors smaller than that of Ferrenberg and Landau [16], and three or more standard errors larger than values listed by Liu and Fisher [l]. Therefore we briefly review the uncertainties in our result for Kc with regard to the sources of error mentioned above. The fact that our random-number generator has passed a number of tests does not prove that our data in table I are free of systematic errors for L 3 5. However, if there are such errors, they will be clearly recognizable as such in the future. Furthermore, additional corrections to scaling may very well play a role. We have also made fits to our data with a term b, L2” included in eq. (2), describing a nonlinearity in the irrelevant field. The result was that the irrelevant exponent yi became poorly determined. Using the value found by Nickel and Rehr [3] as input leads again to KC = 0.221648 (4), where the error margin includes the uncertainty in yi. It is, of course, possible that another correction with a different exponent is significant. In that case, it will undoubtedly reveal itself when accurate data over a wider range of systems sizes become available. In order to gain insight in the discrepancies as mentioned above, it is essential that verifiable numerical data are published, as was done by Livet [17]. Likewise, it is important to specify the details of the random-number generator. This holds especially in the case of the three-dimensional Ising model, where we are confronted with an increasing body of mutually conflicting results. It is a pleasure to thank A.M. Ferrenberg, J.R. Heringa, A. Kooiman and M. Novotny for valuable discussions. We are much indebted to H.J.M. van
460
H. W. J. Bltite, G. Kamieniarz 1 Finite-size calculations on 30 Ising model
Grol, Y.T.J.C. Fonk and M.J. Neefs for teaching us how to operate workstations. This work is part of the research program of the “Stichting voor Fundamenteel Onderzoek der Materie (FOM)“, which is financially supported by the “Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)“.
References [l] [2] [3] [4] [5] [6] [7] [8] [9] [lo] [ll] [12]
A.J Liu and M.E. Fisher, Physica A 156 (1989) 35. J. Adler, J. Phys. A 16 (1983) 3585. B.G. Nickel and J.J. Rehr, J. Stat. Phys. 61 (1990) 1. G.A. Baker, B.G. Nickel, M.S. Green and D.I. Meiron, Phys Rev. B 17 (1978) 1365. J.C. Le Guillou and J. Zinn-Justin, Phys. Rev B 21 (1980) 3976. H.W.J. Blijte and R.H. Swendsen, Phys. Rev. B 20 (1979) 2077. G.S. Pawley, R.H. Swendsen, D.J. Wallace and K.G. Wilson, Phys. Rev. B 29 (1984) 4030. M.N. Barber, R.B. Pearson, D. Toussaint and J.L. Richardson, Phys. Rev. B 32 (1985) 1720. G. Parisi and F. Rapuano, Phys. Lett. B 157 (1985) 301. A. Hoogland, A. Compagner and H.W.J. Bliite, Physica A 132 (1985) 593. G. Bhanot, D. Duke and R. Salvador, Phys. Rev. B 33 (1986) 7841. H.W.J. Blote, J.A. de Bruin, A. Compagner, J.H. Croockewit, Y.T.J.C. Fonk, J.R. Heringa, A. Hoogland and A.L. van Willigen, Europhys. Lett. 10 (1989) 105. (131 H.W.J. Bliite, A. Compagner, J.H. Croockewit, Y.T.J.C. Fonk, J.R. Heringa, A. Hoogland, T.S. Smit and A.L. van Willigen, Phyica A 161 (1989) 1. [14] N.A. Alves, B.A. Berg and R. Villanova, Phys. Rev. B 41 (1990) 383. [15] N. Ito and M. Suzuki, J. Phys. Sot. Jpn. 60 (1991) 1978. [16] A.M. Ferrenberg and D.P. Landau, Phys. Rev. B 44 (1991) 5081. [17] F. Livet, Europhys. Lett. 16 (1991) 139. (181 C.F. Baillie, R. Gupta, K.A. Hawick and G.S. Pawley, Phys. Rev. B 45 (1992) 10438. [19] R.H. Swendsen and J.S. Wang, Phys. Rev. Lett. 58 (1987) 86. [20] U. Wolff, Phys. Rev. Lett. 60 (1988) 1461. [21] A. Hoogland, J. Spaa, B. Selman and A. Compagner, J. Comput. Phys. 51 (1983) 250. [22] A. Hoogland, A. Compagner and H.W.J. Blote, in: Architecture and Performance of Specialized Computers, B. Alder, ed., Computational Techniques (Academic Press, New York, 1988) Ch. 7, The Delft Ising System Processor. [23] A.M. Ferrenberg, D.P. Landau and Y.J. Wang, Phys. Rev. Lett. 69 (1992) 3382. [24] M.E. Fisher, comment on a lecture by R.H. Swendsen at the Statphys 18 Conference (Berlin, 1992). [25] N. Ito and Y. Kanada, Supercomputer 5 (3) (1988) 31. [26] K. Binder, Z. Phys. B 43 (1981) 119. [27] J.R. Heringa, H.W.J. Blote and A. Compagner, Int. J. Mod. Phys. C 3 (1992) 561. [28] N. Zierler, Inform. Control 15 (1969) 67. [29] C.F. Baillie and P.D. Coddington, Phys. Rev. B 43 (1991) 10617. [30] R.B. Pearson, Phys. Rev. B 26 (1982) 6285. [31] H. Saleur and B. Derrida, J. Phys. (Paris) 46 (1985) 1043. [32] M.P. Nightingale, Proc. K. Ned. Akad. Wet. B 82 (1979) 23.5. [33] G. Kamieniarz and H.W.J. Bliite, J. Phys. A 26 (1993) 201. [34] A. Rosengren, J. Phys. A 19 (1986) 1709.