GitHub

Example: 2D Ising model

The Ising model

The Ising model is a classic statistical mechanics model defined by a Hamiltonian

\[H = -J \sum_{i \sim j} S_i S_j\]

where $S_i \in (-1,1)$ is the spin of site $i$. The sum is over all adjacent pairs on a specified lattice. Here we assume that the interaction strength is in units where $J = 1$ and we ignore any external magnetic field. We will apply this model on a square grid of size $L \times L$ sites. This model is the simplest statistical mechanical system that exhibits a continuous phase transition.

Setup

To apply the Wang-Landau algorithm, we load WangLandau along with necessary definitions for the Ising model

using WangLandau, LinearAlgebra, Plots, LaTeXStrings;

include("../../examples/ising.jl");

The main definition in ising.jl is the mutable struct Ising2D, which holds an integer array representing the configuration of spins and the energy of the configuration. The rest of the code defines methods for applying the WangLandau interface to Ising2D (see WangLandauProblem). These methods define the histograms and specify how to make and accept random trial moves. For the Ising model a simple trial move is flipping a single spin at a random site on the lattice.

This example can have periodic or free boundary conditions, here we will define the problem for the case of free boundary conditions:

L = 15;

P = false;

prob = WangLandauProblem(Ising2D(L; periodic=P));
WangLandauProblem{Main.Ising2D{false, false}}(Main.Ising2D{false, false}([0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0], 15, 420))

When Ising2D is constructed it will initialise a random configuration of spins.

Simulation

Next, we run the simulation. We will use the default strategies for altering the $f$ parameter, but set a final desired value of $f < 10^{-6}$. We can also set how often to check for flatness by setting the keyword check_sweeps. A single Monte Carlo sweep is $N$ trial moves, where $N$ is the canonical size of the system, in this case $N = L^2$.

sim = solve(prob; check_sweeps = 500, final_logf = 1e-6);
WangLandauSimulation(Main.Ising2D{false, false}([0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0], 15, 420))
  log(f) = 9.5367431640625e-7 (final: 1.0e-6)
  iterations: 20 (checks: 3762)
  total steps: 423225004
  elapsed time: 1.90e+01 s

To see if the algorithm performed well, we should inspect the histogram of samples. The output sim.samples is the accumulation of the histogram at each iteration. 

scatter(sim.samples; label = :none, ylabel = "samples", msw = 0, ms = 2)
Example block output

If the histogram is flat then we can proceed to analysis.

The samples histogram is indexed by $i \in 1,\ldots,i_\mathrm{max}$, where $i_\mathrm{max} = 2 (L - P)^2 + 2 (L - P)$ is the number of accessible energy values and $P=0$ for free boundary and $P=1$ for periodic boundary. This is more efficient for internal processing and the conversion to consecutive array indexing is handled by the Ising2D code. However, for physical analysis we want a range of real energy values. The extremal energy values for an $L \times L$ lattice are $\pm E_\mathrm{max}$ where $E_\mathrm{max} = 2(L - (1 - P))^2 + 2(L - (1 - P))$, and the change in energy from flipping a single spin is $2(1+P)$. In fact, $E_\mathrm{max}$ is already stored in Ising2D so we can just define

maxE = sim.statedefn.maxE;

Es = -maxE:(2*(1 + P)):maxE;
-420:2:420

Next, we look at the density of states, normalised to its maximum value. We also apply a mask to avoid any hard-to-reach states that were not sampled.

mask = sim.samples .> 0;

logdos_shift = maximum(sim.logdos[mask]);

dos = exp.(sim.logdos[mask] .- logdos_shift);

scatter(Es[mask], dos; label = :none, xlabel = L"E", ylabel = L"g(E)", yscale = :log10, ms = 2, msw = 0)
Example block output

Critical point

The 2D Ising model has a continuous phase transition at $\beta_c J = \log (1+\sqrt{2}) / 2 \approx 0.44\ldots$, where $\beta$ is inverse temperature. The output of the simulation is the entropy $\log g(E)$, i.e. the logarithm of the density of states, from which we can calculate moments of the energy

\[\langle E^n \rangle = \frac{\sum_E E^n g(E) \exp(-\beta E)}{\sum_E g(E) \exp(-\beta E)}\]

To compute moments of the density of states from the output data we define 

function moment_E(Es, g, n, β); E_factor = @. Es^n * exp(-Es * β); return dot(g, E_factor); end;
moment_E (generic function with 1 method)

where Es is the range of energy values, g is the corresponding density of states, n is the desired moment and β is an inverse temperature value. Using this we compute the partition function, the average energy and the specific heat (energy variance) as the 0th, 1st and 2nd moments, for a range of values of $\beta$

βs = 0.1:0.01:1.0;

Z = [moment_E(Es[mask], dos, 0, β) for β in βs];

avgE = [moment_E(Es[mask], dos, 1, β) for β in βs];

avgE2 = [moment_E(Es[mask], dos, 2, β) for β in βs];

varE = @. (avgE2 / Z - (avgE / Z)^2) * βs^2;
91-element Vector{Float64}:
  4.381648733936245
  5.355894854301492
  6.443483647934456
  7.649446869440791
  8.979221631154859
 10.438887435720494
 12.035505773002306
 13.777513433743854
 15.6750966613888
 17.740467217881235
  ⋮
 19.150471712490173
 18.17991132174422
 17.26152757693942
 16.39234017549148
 15.569548124492169
 14.79051868084821
 14.052776997923909
 13.353996414142145
 12.691989329236094

Now we can plot the average energy

plot(βs, avgE ./ Z ./ maxE; label = :none, xlabel = L"\beta", ylabel = L"\langle E \rangle / E_\mathrm{max}", xlim = (0, Inf))
Example block output

and the specific heat, or variance

plot(βs, varE ./ L^2; label = :none, xlabel = L"\beta", ylabel = L"\mathrm{var}(E)", xlim = (0, Inf), ylim = (0, Inf))
Example block output

Finally, the peak of the specific heat is a signature of the location of the critical point

varE_max, varE_max_i = findmax(varE);

varE_max / L^2, βs[varE_max_i]
(1.1110660614018033, 0.49)

which compares well to the known value.