# Study of a microcanonical algorithm on the ±𝐽 spin glass model in 𝑑=3.

## 1 Introduction

The Ising spin glass is the paradigm of the complex systems. It possesses two fundamental characteristics: disorder, because the couplings are random variables and frustration since the signs of the coupling are positive or negative. These two characteristics produce the main property of a complex system: a very slow dynamics.

This slow dynamics is due to the existence and competition of a large number of pure and metastable states below the critical temperature. In some cases, a large number of metastable states above the critical temperature can produce this effect even in the paramagnetic phase.

During the last two decades, the existence of a low-temperature phase in Ising spin glass in three dimensions has been investigated and it has consumed a large amount of CPU resources. At equilibrium it is possible to simulate up to $`L=16`$ lattices and the signature of the transition is very weak .

In reference  was found a behaviour for the dynamical critical exponent $`z(T)7T_c/T`$, where $`T_c1`$ means for the critical temperature and $`T`$ is the temperature, for the three dimensional Gaussian Ising spin glass using a Metropolis dynamics. This numerical fact has been corroborated in experiments using samples of CuMn (at 6%) and thiospinel . This result for the dynamical critical exponent implies Monte Carlo (MC) thermalization times proportional to $`L^7`$ near the phase transition ($`L`$ is the size of the system). One can compare this behaviour with that of the pure Ising model: the thermalization time diverges as $`L^2`$ (near its critical temperature).

From the previous discussion follows that traditional approaches using only local algorithms should fail in thermalizing a large system. Moreover, the absence of a non local update algorithm and the high value of the dynamical critical exponent for the local ones convert this problem in a very challenging computational issue. Recent works using large amounts of computational power with a standard local MC simulation  showed some evidences of a cold phase in this model. Other recent approaches  in similar models have used more sophisticated update algorithms, based on a combination of a standard local MC run and an innovative update process in the temperature. These algorithms , the simulated tempering and parallel tempering, succeed in thermalizing systems at temperatures lower than those reachable by a standard Metropolis simulation.

Moreover, in a canonical simulation, the most time-consuming task is the generation of the random numbers, and so, one possibility is to use a MC method that does not use random numbers. This calls for microcanonical methods. In particular, Creutz developed the so-called demon algorithm that does not need random number generation to work.

The aim of this paper is to investigate the behavior of this microcanonical local update algorithm on this model. Since that, we will use only this algorithm, although it is clear that other tools, as parallel tempering, have to be implemented to improve the simulation in order to try to elucidate the low temperature regime of this model. In particular, we will study numerically the ergodicity of the algorithm, efficiency (i.e. autocorrelation times) and the difference with the canonical algorithm when we work at finite volume and how these differences go to zero with the volume. It is clear, that a study of this kind is essential if we will use the microcanonical algorithm in extensive numerical simulations.

As further studies we will plan to analyze the performance of a combination of microcanonical and canonical algorithms. This combination has worked very well in the simulation of some physical systems as Quantum Chromodynamics  and it could be of great interest to check if this combination will work well in spin glasses.

Moreover, one of the authors of this paper is finalizing the construction of a dedicated machine to simulate this model , being this paper a preliminary study of the characteristics of the demon algorithm, previous to its hardware implementation. We remark that this algorithm could also be used on computers of general purpose not only in dedicated machines.

## 2 Model, observables and update algorithm

The $`\pm J`$ spin glass model is defined by the Hamiltonian

$$H\underset{<i,j>}{}\sigma _iJ_{ij}\sigma _j,$$
(1)   
where the spins $`\sigma _i`$ take values $`\pm 1`$. The nearest neighbor quenched couplings $`J_{ij}`$ take values $`\pm 1`$ with equal probability. The spins live in a cubic lattice containing $`V=L^3`$ sites. We have used helicoidal boundary conditions in two directions and periodic in the third one. The reason is because we wanted to check the special purpose computer developed for this physical model .

As usual, for every realization of the bonds or sample two independent copies of the system are studied. The main quantity to be measured is the overlap between the two copies with the same disorder, which acts like an order parameter for this model. The overlap between two spin configurations $`\sigma `$ and $`\tau `$ is given by

$$q(\sigma ,\tau )\frac{1}{V}\underset{i}{}\sigma _i\tau _i,$$
(2)   
which is usually denoted $`q`$. Using powers of this quantity one can compute different observables. The second and fourth power are used to build the Binder parameter

$$g\frac{1}{2}\left[3\frac{\overline{q^4}}{\overline{q^2}^2}\right],$$
(3)   
where $`(\mathrm{})`$ stands for the thermal average for a given realization of the bonds, and $`\overline{(\mathrm{})}`$ means the average over the disorder. Since this quantity is dimensionless, it obeys the finite size scaling law (near the critical point)

$$g=\stackrel{~}{g}\left(L^{1/\nu }\left(TT_c\right)\right),$$
(4)   
being independent of the volume at the critical temperature, $`T=T_c`$. This property makes it appropriate to investigate the existence of any spin glass phase transition by studying the intersections of the functions $`g(T)`$ for different lattice sizes.

In addition, one can compute the spin glass susceptibility,

$$\chi V\overline{q^2}.$$
(5)   
The algorithm we want to investigate is the demon algorithm  proposed by Creutz. For this microcanonical algorithm, the physical system is the standard lattice plus a demon, which acts like an entity able to store energy. The update algorithm keeps constant the sum of the energy of the lattice and the demon.

In order to carry out the MC simulation for a given total energy $`H`$, one can start from a spin configuration with that energy and the demon energy equal to zero. To generate new spin configurations, the spins are updated as follows: first, a spin is selected and its sign is proposed to be inverted. If the flip lowers the spin energy, the demon takes that energy and the flip is accepted. On the other hand, if the flip grows the spin energy, the change is only made if the demon has that energy to give to the spin.

At this level, the value of the temperature $`T`$ is unknown and it can be obtained from the demon energy, whose probability distribution $`p(E_d)`$ is given by the expression

$$p(E_d)e^{\beta E_d}.$$
(6)   
A fit to this function can provide the $`\beta `$ value, although a better estimate can be calculated if the mean energy of the demon on the sample, $`E_d`$, is computed. Thus, the $`\beta `$ value is obtained as

$$\beta =\frac{1}{T}=\frac{1}{4}\mathrm{log}\left(1+\frac{4}{E_d}\right).$$
(7)   
In a spin glass model, given the energy of the simulation, a different value of the temperature is obtained for every sample. The average of them all gives the temperature of the simulation corresponding to the fixed energy.

## 3 Numerical results

We have simulated some different energies per spin $`e`$ in three lattice sizes ($`L=8,12`$ and $`16`$). For every couple of parameters $`(L,e)`$ we generate a large set of samples. The initial demon energies and spin configurations are chosen to give a total energy $`H=eV`$. In order to generate them, we start with $`E_d=0`$ and all the spins equal to 1. We use the demon algorithm to change the spin configuration, but when $`Ed>4`$ we steal 4 units with a probability of 50%. In this way, we approach smoothly towards the energy desired for the simulation. When the required energy is reached, the simulation starts. The spins to be updated are sequentially chosen, having a completely updated spin configuration after every $`V`$ updates.

Table 1 shows the parameters of the different runs: lattice size, energy per spin, number of samples and total number of Monte Carlo sweeps are showed in the first columns. The next one shows the thermalization time $`t_0`$ (we will discuss in detail below how we have computed the thermalization time). The last column is the number of measures of the overlap considered in every run to compute the thermodynamical average.

The results obtained in these runs are shown in Table 2. We have computed the temperature according to Eq. (7). The second and fourth powers of the overlap have been also calculated to obtain the Binder cumulant. The work has been carried out on the RTNN computer, which holds 32 PentiumPro processors, for a total CPU time of approximately 20 days of the whole machine. The errors in the estimates of the observables have been calculated with the jack-knife method .

In order to be sure that the system is thermalized before measuring, we check the symmetry on the probability distribution of the overlap and also the MC evolution of the spin glass susceptibility. For every run reported in this paper the mean value of the overlap is zero (within the statistical error). Moreover, we have checked the symmetry around zero of the probability distributions of the overlap.

The MC evolution of the spin glass susceptibility is plotted in figure 1. Every point in the plot has been computed by averaging the values for the overlaps only in its corresponding MC time. One expects to see the susceptibility rising with the Monte Carlo time until a plateau is reached. The beginning of this plateau defines the thermalization time $`t_0`$. We use this criterion for the thermalization. As we said above, the temperature corresponding to the simulation can be computed by using Eq. 7. For every sample we obtain its own temperature. Figure 2 shows the normalized probability distribution of temperatures obtained for the highest and the lowest cases in both $`L=8`$ and $`L=12`$ lattices. Note the width of the histograms.

## 4 Demon vs. Canonical: a comparison

The results exposed in Table 2 can be compared with those previously obtained  in a canonical simulation running a heat bath update algorithm. Moreover we have performed canonical numerical simulations (using a Metropolis algorithm and periodic boundary conditions) in order to run just at the temperature given by the demon algorithm (and hence, compare at the same temperatures). We report the parameters of these canonical simulations in Table 3 (the colums are: lattice size, temperature, number of samples, number of Monte Carlo sweps, thermalization time ($`t_0`$) and number of measures) and their results in Table 4. In addition to this, in Table 4 we have computed the difference between the Binder cumulant computed in the demon simulation and the canonical one.

Fig. 3 shows the second and fourth moments of the overlap. We take as reference the heat bath data from Kawashima and Young  and our own Metropolis data. As a check of our Metropolis simulation it is clear that our data match very well (one standard deviation) between those from Kawashima and Young.

Now we can confront the demon data with the canonical. While the squared overlap seems to fit perfectly in the data obtained with heat bath and Metropolis, the overlap to the fourth differs significantly, being the microcanonical data lower than the canonical. The discrepancy between the two ensembles decreases with the lattice size.

These two quantities are used to calculate the Binder parameter, which is plotted in Fig. 4. The canonical simulation gives a cut point between the two Binder parameters. In our microcanonical simulation both Binder parameters approach at low temperature. Data converge to be compatible in the error bar, but no cut point is resolved using $`L=8`$ and $`L=12`$ data (we have simulated with the demon algorithm $`L=16`$ data in the region $`TT_c=1.11`$).

To be sure of the correctness of the algorithm we have carried out some extra runs at $`e=1.650`$ and $`L=8`$. One of them was using periodic boundary conditions (to check the effect of periodic and helicoidal boundary conditions on the observables using the demon algorithm). In addition, to check the ergodicity of the algorithm, we repeated the simulation with the same realizations of the disorder but starting from spin configurations obtained from a thermalized heat bath simulation. In both cases, we obtained compatible results.

The canonical and microcanonical ensembles must agree in the thermodynamical limit. The discrepancies between them have to decrease when volume goes to infinite. To check it, we have simulated $`e=1.650`$ at $`L=16`$. In this case, we obtained $`g=0.640(4)`$, nearer to the Binder cumulant coming from the canonical simulation than the $`L=8`$ and $`L=12`$ cases. We have seen that the discrepancy of the Binder cumulant (sixth column of Table 4) for $`e=1.650`$ goes to zero following a power law (with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}1`$ where d.o.f stands for degrees of freedom): $`\mathrm{\Delta }g(e=1.65)L^{0.86(36)}`$. We can repeat this procedure for other energies. For instance, if $`e=1.700`$ then $`\mathrm{\Delta }g(e=1.70)L^{0.49(43)}`$ with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=1.56`$ with a confidence level of 21%.<sup>1</sup><sup>1</sup>1The confidence level is the probability that $`\chi ^2`$ were greater than the observed value assuming that the statistical model used is correct (in our case the power law behaviour). A very low value of this confidence level (e.g. $`<5\%`$) would imply that our statistical model is incorrect. See for instance . In any case, a more detailed study of this issue is needed.

We can study in more detail the previous issue by comparing the probability distributions of the overlap, its second power, and its fourth power obtained with the demon algorithm with $`e=1.650`$. Moreover, we have also measured the previous three probability distributions carrying out a heat bath simulation at temperature $`T=1.272`$ and $`L=8`$ with the same sets of bounds and number of iterations of the case $`e=1.650`$. We show these probability distributions in Fig. 5.(a). The different shape of the distribution is clarified in the plots of the powers of the overlap, being the demon distribution more peaked than the canonical one.

Other interesting issue is to compare the thermalization times needed in the Metropolis simulation and in the demon algorithm. All these thermalization times have been reported in the fifth column of Tables 1 and 3. We remark that these times have been obtained by monitoring the growth of the non linear susceptibility as a function of the Monte Carlo time. We take as $`t_0`$ (in the canonical as well as in the microcanonical simulations) the Monte Carlo time in which the numerical data achieve the equilibrium plateau (see figures 1 and 6 ). From the values of $`t_0`$ for the lowest temperatures follows that demon algorithm thermalizes slower (with a factor between 2 and 3) than Metropolis.

Obviously the cost of introducing a random number generator (i.e. to run a canonical simulation instead of a microcanonical one) is less than the factor two and three found in the autocorrelation times, and so one of the conclusions of this paper is that (in an general purpose computer) the efficiency of canonical algorithm are bigger than that of the microcanonical one (bigger thermalization times).

It is clear that a dedicated machine with programmable logic (running the microcanonical algorithm) with a speed 16 times bigger than a supercomputer (running a Metropolis algorithm) clearly compensates (for lattices of order $`16`$) the excess of thermalization time of the demon algorithm (which is between 2 and 3). This is the situation if we compare in a tower of APE100 supercomputer (which has a peak performance of 25 GigaFlops) , with a real performance of 5000 ps per spin. and a machine with programmable logic (312 ps per spin). Of course, that special purpose computer running a Metropolis algorithm (at the same speed) would be even more efficient because the smaller thermalization times. As a matter of fact, the special purpose computer referred in this work  is able to run both algorithms at the same speed thanks to a fast random number generator implemented in hardware. The study of the efficiency of a combination of the two algorithms is left to a later work.

We finally report our last check of the demon algorithm by checking the Guerra relations which seem to be fulfilled within a 0.5% precision in a canonical simulation .

One of the Guerra’s relations  reads

$$\overline{q^2^2}=\frac{1}{3}\overline{q^4}+\frac{2}{3}\overline{q^2}^2.$$
(8)   
This relation has been shown exact for the Gaussian model . This relation can be rigourously demonstrated in the infinite volume limit. However, one would expect finite corrections in the Gaussian case. Even though there is no proof for the $`\pm J`$ mode, the difference between the two sides of the equation has to decrease with the volume. Table 5 shows our results for the left hand side (lhs) and the right hand side (rhs) of Eq. 8. The errors are calculated by a jack-knife analysis.

In the sixth column of Table 5 we report the difference between the lhs and the rhs of Eq. 8. The maximum deviation is 2.5 standard deviations ($`L=8`$ and $`e=1.650`$). The rest of the differences of the Table 5 have fluctuations less than two standard deviations. Hence, we can conclude that the Guerra relation is satisfied in the demon algorithm.

## 5 Conclusions

We have studied a microcanonical algorithm running on the three dimensional Ising spin glass in three dimensions.

We have obtained compatible (within the statistical errors) values among the results of a canonical numerical simulation and the demon algorithm for the second and fourth moments of the overlap whereas the values of the Binder cumulant are different (but with the discrepancy going to zero following a power law). We remark that microcanonical and canonical algorithms should provide the same numerical results only in the thermodynamic limit. Moreover the microcanonical algorithm satisfies one of the Guerra relations.

Finally we have shown that the thermalization times needed for the demon algorithm are two or three times larger than those for the Metropolis ones (for the larger simulated lattice $`L=16`$).

We remark that we have checked numerically the ergodicity and the efficiency of the algorithm.

From the point of view of the efficiency we have shown that the cost of introducing random numbers is less than the excess of thermalization which need the microcanonical simulation. Obviously, if we can design a dedicated machine where only it is possible to implementate (via hardware) a microcanonical algorithm, and if this dedicated machine runs to a speed which is bigger than 5 times the speed of a canonical code in a supercomputer the use of the microcanonical algorithm will be welcome. Obviously this work shows that if we can implementate random numbers with the cost of a factor two in time we should use the canonical algorithm instead of the microcanonical one.

We wish to thank J. M. Carmona, L. A. Fernández, D. Iñiguez, G. Parisi and A. Tarancón for useful discussions. CLU is a DGA fellow. We also wish to thank P. Young for providing us the numerical data of his reference .