# Weird results of Monte Carlo simulation

I'm simulating the 3D Ising Model using the Wolff update algorithm. I am using the Mersenne Twister RNG.

When the lattice size is $$L = 50$$, the specific heat curve looks very weird!!

I want to know is it possible that this error is because of RNG (random number generator) effects? or my run configuration such as thermalization and sampling and time steps ... is incorrect?

run details:

define 1 time step = N(=L^3) spin updated(at least)

T steps = 0.001

thermalization time = 10000 time step

number of samples in each T = 10^5

total number of runs = 50 --> total samples in each T = 5 000 000 .

time interval between 2 samples = 1 time step.

RNG = SIMD-oriented Fast Mersenne Twister (SFMT)

update:

this problem is more noticeable when time interval between 2 samples is L^3(like this run) but in my previous run with time interval between 2 samples = L^2 , its not significant!!

• MT has a huge period, compared to simple old-fashioned linear congruential PRNGs. But you can still get artifacts from it if you do silly things with it, like reseeding it with system time in a fast loop. ;) – PM 2Ring Nov 15 '18 at 13:17
• I seed RNG only once by intel RDRAND. – mehrdad Nov 15 '18 at 15:22
• The RNG seems unlikely to be the problem. Do you get more sensible results for other L? Do any of your 50 runs look particularly suspicious (e.g., slow to thermalize)? What do you get with simple spin-flip updates? That last approach won't be so efficient, but it's easier to debug and can provide some reference benchmarks once debugged. – David Schaich Nov 15 '18 at 18:54
• @ David Schaich , YES all of my 50 runs have this behavior, you can see that the specific heat curve is very smooth(this is because each 50 run have this behavior ) – mehrdad Nov 15 '18 at 19:47
• @ David Schaich , concerning other L: right now I am running for L=40 and L=60 with the same settings. – mehrdad Nov 15 '18 at 19:51

There is a celebrated paper by Ferrenberg et al, Phys Rev Lett, 69, 3382 (1992) who reported errors using the Wolff method for the Ising model, which they attributed to poor random number generators. However this pre-dated the development of the Mersenne twister, and it's not completely clear what recommendations came out of this, beyond "try several different RNGs": see, for instance, this preprint.

One of the most recent studies of the Ising model by the same group uses the Wolff method with the Mersenne twister RNG: see Ferrenberg et al Phys Rev E, 97, 043301 (2018) also available as a preprint. So I would guess that they are happy with this RNG.

You might consult Deng et al Comp Phys Commun, 178, 401 (2008) who extol the virtues of "twisting and combining" random numbers, in their paper which is also targetted at the Wolff method. Unfortunately I can't find an open access version of this paper.

I think it's more likely that there is an error in your program. Or possibly you just need to simulate for much longer: you are, after all, simulating at around the critical temperature $$k_BT_c/J\approx 4.5$$. This is the wrong site to ask people for advice on debugging your program; for information on "how long to run" you should look at the literature (for example, the 2018 Ferrenberg paper cited above).

[EDIT following OP comment]

It is hard to judge whether your run parameters are satisfactory or not. The 2018 paper I mentioned is clearly "state of the art", using considerable computer power to get very precise results. But reasonable results could be obtained for this model, more than 25 years ago, using conventional updating schemes, see for example Ferrenberg and Landau Phys Rev B, 44, 5081 (1991) and, very roughly, your runs look to be of comparable length (measured in "complete lattice updates"). Your equilibration periods look somewhat short, but you seem to be systematically sweeping through the temperatures in very small steps, so maybe the equilibration period does not need to be too long. Even in that 2018 paper, they mention equilibration periods of $$\sim 10^5$$ Wolff steps, and my reading of their paper is that the actual equilibration was happening over a shorter period.

If you are not doing this already, you should at least be storing the histogram of energies at each temperature, and using histogram reweighting to compute results (average energy and heat capacity) at nearby temperatures, which will give you additional possibilities to check your results for self-consistency.

I cannot add any more to my answer: it would be guesswork, and nobody but you can check your program and your output in detail.

• I updated question, added run details. – mehrdad Nov 15 '18 at 15:02
• @ LonelyProf , thank you for introducing references! – mehrdad Nov 15 '18 at 16:02
• I have edited my answer following your comment. I hope this is helpful. Can I encourage you to re-edit the run details in your question so as to use MathJax for the variables, equations, and symbols (such as the rightarrow)? I shall be deleting my original comment on your question, and later today this comment as well, to keep things tidy. – user197851 Nov 15 '18 at 16:07
• yes the 2018 paper you mentioned is really the state of art, I have read it, and also used some of there idea such as using quad-precision floating point in my program. – mehrdad Nov 15 '18 at 16:18
• concerning about thermalization , the thermalization is reached very soon! after just 100 of time steps the energy and magnetization seems fluctuate between a constant value – mehrdad Nov 15 '18 at 16:22