Quarto demo

Test of Potts model

using WGLMakie
using LinearAlgebra
using CircularArrays
using StatsBase: Weights, sample
using Random: MersenneTwister
using Distributions
using MakiePublication
using QuartoTools

ham(J, q::Int, n1::Int, n2::Int) = J*cos(2pi / q * (n1 - n2))

ham (generic function with 1 method)

function choose_weights(_s, _q)
  # choses weights for spin flip
  weights = Weights(ones(_q)/(_q-1))
  weights[_s+1] = 0.0 # set probability to choose spin-i again to be zero
  return weights
end

function spin_flip(_si, _q)
  # randomly chooses ANOTHER spin
  # spin = round(Int,_si / 2pi * _q) # convert to integer selection

  w = choose_weights(_si, _q)
  return sample(0:q-1, w)
end

spin_flip (generic function with 1 method)

q = 4  # theta discreteness
N = 64

n_dims = (N, N)
kB = 5.0
J = -1.0
spins = CircularArray(sample(0:q-1, Weights(ones(q)/q), n_dims)) 


# neighbors
locs      = eachindex(spins)
basis_vec = eachcol(diagm(ones(Int8,length(n_dims)))) # general dimensions
neighbors = map(_n->CartesianIndex(_n...), [basis_vec; -basis_vec])

MAXITER  = 100_000
STEP     = MAXITER ÷ 1000
time     = LinRange(0, 1, MAXITER)
βmax     = 50.0
β        = @. βmax * ((1 + exp(-2*(time-1/2)*10))^-1 - 1/2)  #@. βmin  - 2βmin*time
temp     = @. 1.0 / β

spin_vec = zeros(size(spins)..., MAXITER)

for k in 1:MAXITER

  spin_vec[:,:,k] = spins.data # save state


  i = rand(locs) # choose random spin location
  si = spins[i]

  nearest = i .+ neighbors
  
  energy = mapreduce(sj -> ham(J, q, si, sj), +, spins[nearest]) # sum 

  si_new = spin_flip(si, q)

  energy_new = mapreduce(sj -> ham(J, q, si_new, sj), +, spins[nearest]) # sum 
  dE = energy_new - energy

  if dE < 0.0
    spins[i] = si_new
  elseif rand() < exp(-dE * β[k])
    spins[i] = si_new
  end
end