Simulating 2D Ising model.

ising2d.jl

#!/usr/bin/env julia

using Printf: @printf;

function bc(x, L)
  mod(x-1, L) + 1;
end

function totalEnergy(grid ::Matrix{Int64}) ::Float64
  (x, y) = size(grid);

  local En = 0
  for i in 1:x
    for j in 1:y
      S = grid[i, j];
      # neighbors
      WF = (grid[bc(i+1, x), j] + 
            grid[bc(i-1, x), j] + 
            grid[i, bc(j+1, y)] + 
            grid[i, bc(j-1, y)])
      En += -WF * S
    end
  end

  return En / 2
end

function exponent_cache(T::Float64) ::Vector{Float64}
  cache = zeros(Float64, 9)
  cache[5+4] = exp(-4.0*2/T)
  cache[5+2] = exp(-2.0*2/T)
  cache[5+0] = 1.0
  cache[5-2] = exp( 2.0*2/T)
  cache[5-4] = exp( 4.0*2/T)
  return cache
end

function sample(nrun::Int64, grid::Matrix{Int64}, T::Float64, warm::Int64, measure::Int64)
  En = totalEnergy(grid);
  Mn = sum(grid);
  (L, _) = size(grid);

  N1 = 0;
  E1, M1, E2, M2 = 0.0, 0.0, 0.0, 0.0;
  L2 = L*L;

  expcache = exponent_cache(T);

  for irun in 1:nrun
    i = rand(1:L)
    j = rand(1:L)

    S = grid[i, j];
    WF = (grid[bc(i+1, L), j] + 
          grid[bc(i-1, L), j] + 
          grid[i, bc(j+1, L)] + 
          grid[i, bc(j-1, L)])
    P = expcache[5 + S*WF]

    if P > rand()
      grid[i, j] *= -1
      En += 2*S*WF
      Mn -= 2*S
    end

    if irun > warm && irun%measure == 0
      N1 += 1
      E1 += En
      E2 += En * En
      M1 += Mn
      M2 += Mn * Mn
    end
  end

  E = E1 / N1
  M = M1 / N1

  cv = (E2/N1 - E^2) / T^2;
  chi = (M2/N1 - M^2) / T;

  return (M/L2, E/L2, cv/L2, chi/L2)
end

function main(gsize=0)
  TRange = LinRange(0.015, 4.5, 300);
  if gsize == 0
    print("Input lattice grid size: ")
    N = parse(Int64, (readline() |> strip))
  else 
    N = gsize
  end
  println("Start simulating with grid size $N x $N")
  grid = rand([-1, 1], (N, N));

  f = open("F$N.txt", "w");
  @printf f "#%6s %20s %20s %20s %20s\n" "T" "M" "E" "cv" "chi"
  @printf "#%6s %20s %20s %20s %20s\n" "T" "M" "E" "cv" "chi"

  for T in TRange
    M, E, cv, chi = sample(4*10^5, grid, T, 10^5, 10);
    @printf " %6.3f %20.10f %20.10f %20.10f %20.10f\n" T M E cv chi
    @printf f " %6.3f %20.10f %20.10f %20.10f %20.10f\n" T M E cv chi
  end

  close(f);
end

main()