Ising model simulation

This posting shows how to use Mathematica (version 6) to create an interactive simulation of a 2-dimensional Ising model, which is based on some work on Boltzmann machines that I did back in the mid 1980s. This nicely illustrates the use of the Manipulate function in Mathematica, which is by far the easiest way of creating interactive graphical demonstrations that I have ever used.

Define a function for randomly initialising the state of the 2-dimensional Ising model.

init[narray_] := array = RandomInteger[{0,1}, {narray,narray}];

Define a function for doing a single Monte Carlo update of the 2-dimensional Ising model:

update[{pe_,pne_,pn_,pnw_}, narray_] :=
Module[
{probpatch={{pnw,pn,pne}, {pe,0.5,pe}, {pne,pn,pnw}}, i, j, di, dj, arraypatch, freq1, freq0},
{i,j} = RandomInteger[{1,narray}, {2}];
arraypatch = Table[array[[Mod[i+di,narray,1],Mod[j+dj,narray,1]]], {di,-1,1}, {dj,-1,1}];
freq1 = Apply[Times, probpatch arraypatch+(1-probpatch)(1-arraypatch), {0,1}];
freq0 = Apply[Times, (1-probpatch)arraypatch+probpatch(1-arraypatch), {0,1}];
array[[i,j]] = If[Random[] < freq1/(freq0+freq1), 1, 0];
];

Define a function for implementing the interactive graphical demonstration.

Manipulate[
If[narray!=Length[array], init[narray], ];
If[u, Do[update[{p2en[[1]],p2nenw[[1]],p2en[[2]],p2nenw[[2]]}, narray], {nupdate}], ];
ArrayPlot[array],
{{p2en, {0.5,0.5}, "probability(\[RightArrow],\[UpArrow])"}, {e, e}, {1-e, 1-e}, ControlPlacement->Left},
{{p2nenw, {0.5,0.5}, "probability(\[UpperRightArrow],\[UpperLeftArrow])"}, {e, e}, {1-e, 1-e}, ControlPlacement->Left},
{{u,"False","update"}, {True,False}, ControlPlacement->Top},
{{narray,32,"size"}, 3, 64, 1, Appearance->"Labeled", ControlPlacement->Top},
FrameLabel->"Ising Model",
Initialization->(e=0.001; array={}; nupdate=50)
]

When Manipulate is evaluated the initial output typically looks like this:

Use the check box at the top of the output window to switch the simulation on and off. The number of Monte Carlo updates that is applied between the display of each "frame" of the simulation is hardwired to be 50.

The initial size of the 2-dimensional array of "spins" is 32 by 32, which is controlled by the 1D slider at the top of the output window. This slider can be moved at any time, even whilst a simulation is being run.

Each Monte Carlo update is implemented by selecting a spin at random, and calculating 8 probability factors due to its interaction with each of its 8 surrounding spins. Each such factor specifies a multiplicative contribution to the relative likelihood that the central spin has the same/different value compared to its neighbouring spin (this is not the most general interaction that is possible). There are 4 independent probability factors corresponding to the following directions from the central spin: east, north, north-east, and north-west. These probability factors are controlled by the two 2D sliders on the left of the output window. The top 2D slider controls the east (left/right slider movement) and north (up/down slider movement) factors, and analogously the bottom 2D slider controls the north-east and north-west factors.

Denote the probability factor for a pair of spins being the same as p, and for being different as 1-p. The central position of each 2D slider corresponds to a probability factor p=1/2 (i.e. the same contribution whether the spins are the same or different), the left and bottom edges of each 2D slider correspond to p=0 (i.e. force the spins to be different), and the right and top edges correspond to p=1 (i.e. force the spins to be the same). These sliders can be moved at any time, even whilst a simulation is being run.

So, using a 1-dimensional example, the spin configuration 0?0 gives a product of probability factors (1-p)^2 for 010 and p^2 for 000, whereas 1?0 gives p(1-p) for 110 and (1-p)p for 100 (i.e. the same in each case). A Monte Carlo update of the central spin selects ? as being 1 or 0 with a relative probability given by the ratio of the corresponding products of probability factors. This ratio is ((1-p)/p)^2 and 1 in the two examples above, respectively. The 2-dimensional case is a straightforward generalisation of the 1-dimensional case.

Each snapshot of the output below shows the typical result of running a randomly initialised simulation for a few seconds with the 2D slider settings as shown in the snapshot.

Note that extreme values are used for the probability factors in each of these simulations, and because the simulations each have a short duration they do not actually reach equilibrium. For instance, the first simulation enforces such strong positive correlations between adjacent spins that it should eventually degenerate to all the spins having the same colour.

Enjoy playing with this Ising model simulation. It is quite interesting to move the 2D sliders to vary the probability factors as the simulation is running, because the speed of the simulation is sufficiently fast that you get an almost real-time response as the Ising model dynamically adjusts its equilibrium state.