Experimental WebGL_2 based applications work only
in Chrome Canary and Firefox Dev yet!

2D Ising model

512x512 Ising lattice is shown below. Up and down spins are
white and black pixels (squares). You can watch thermal fluctuations,
phase transition and clusters formation (or melting) for different
temperatures. Drag by mouse "thermometer" to the right
(the black bar corresponds to the phase transition
temperature T_{c}).
Antiferromagnetic spins ordering appears for negative temperature.

it
T
mag
fps
The thermostat algorithm is used. it is the iteration
per frame number. mag is "magnification".
Choose 1024x1024 lattice.

Phase transitions in magnetics

Phase transitions are observed in surprisingly simple systems, e.g. on
a lattice of interacting spins s_{i} (magnet
vectors). Interaction energy of nearest neighbours pair is
E_{ij} = -J (s_{i}s_{j }) .
Total energy E and magnetization M for a spins configuration
{s_{1}, s_{2}, ... s_{n }} is obtained
by summation throughout the lattice. In the Heisenberg model every spin
can take arbitrary direction. In the XY model spins rotate in a
plane.
In the Ising model spins have only two possible states ±1
(up or down).
As since every spin takes two values, therefore there are 2^{n}
different configurations for n spins. You see below 2^{4} =
16 spin configurations for 2x2 lattice.

For J > 0 the state of lowest energy is when all spins are aligned.
The state has macroscopic magnetizaion, i.e. it is ferromagnetic.
When temperature is low Ising spins minimize energy. Interaction
aligns all spin vectors in the same direction, giving huge total magnetic
fields. At high temperature thermal fluctuations break spins order.
The randomness of the spin configuration tends to wash out the large scale magnetism.
In the 2D Ising model there is a phase transition at T_{c} =
2.269 from disordered (non-magnetic) to ordered magnetic state (see Fig.1).

Antiferromagnetics

If the exchange interaction constant
J < 0 , then for low temperature nearest spins are anti-aligned.
In the simplest Ising antiferromagnetic on square grid they form two
ordered sub-lattices (disposed as cells on a chess-board).
Starting at random configuration you can see clusters of ordered phase
formation and movement of domain walls.
Antiferromagnetic ordering appears in the thermostat
algorithm for negative temperature.

The Monte-Carlo method

For large systems it is impossible to calculate statistical averages directly.
One could generate spins configurations at random and approximate real
thermal averages by Monte-Carlo averages. The next problem here is, most
randomly chosen states will make a negligible contribution to the sum.
Random walks are used to take into account only important spins configurations
by the introduction of a fictitious dynamics (so-called Markov chain of
configurations). This wandering generates states, which are the most probable
from energy-entropy point of view (importance sampling).

The Metropolis algorithm

In the Metropolis algorithm we try to turn over a single spin direction with
transition probability

W_{12} =
exp[(E_{1}-E_{2})/T]

if E_{1} < E_{2}

W_{12} = 1

if E_{1} > E_{2}

where E_{1}, E_{2} are energies of the old and new
configurations (see details in the Gould and Tobochnik book).
Statistical averages may then be computed as simple arithmetic means.
The Metropolis algorithm can by applied easy e.g. to the XY model simulation.
There are many possible choices for W_{12}.
To get equilibrium final spins distribution it is enough that
W_{ij} satisfy the detailed balance conditions
W_{12} exp(-E_{1}/T) =
W_{21} exp(-E_{2}/T) or
W_{12} / W_{21} =
exp[(E_{1}-E_{2})/T].

The thermostat algorithm

In the thermostat algorithm we get a single spin (all the rest are
fixed) into contact with big thermostat at temperature T. Then we
get for transition probabilities (the Glauber formula)
P_{+} = exp(-E_{+} /T) /
[exp(-E_{+} /T) + exp(-E_{-} /T)] =
exp(-(E_{+}-E_{-} )/T) /
[exp(-(E_{+}-E_{-} )/T) + 1] ,
P_{-} = exp(-E_{-} /T) /
[exp(-E_{+} /T) + exp(-E_{-} /T)] = 1 - P_{+} . The Glauber transition probabilities satisfy the detailed
balance conditions.

In this application spins are processed in parallel
on GPU (WebGL 2) in the chess-board order.
The linear congruential random number generator is used
r = 1664525u * r + 1013904223u,
where r is unsigned 32 bit integer. It is initialized
by 0xFFFFFFFF*Math.random() values.