Monte Carlo simulation of the 2D Potts model - Part 1¶
Gen Kuroki
2017-08-15 Ver.2.2 (2017-08-13 created)
The Wolff method version
Reference
In [1]:
using PyPlot
rc("font", size=8)
rc("figure", figsize=(3,3))
using PyCall
@show PyDict(matplotlib["rcParams"])["font.size"]
@show PyDict(matplotlib["rcParams"])["figure.figsize"];
In [2]:
# ## Torus boundary condition
#
function Potts2dEnergy(sigma::AbstractArray{Int64,2})
m, n = size(sigma)
E = zero(Float64)
delta(a,b) = ifelse(a==b, one(Float64), zero(Float64))
chi(cond) = ifelse(cond, one(Float64), zero(Float64))
for x in 1:m-1
for y in 1:n-1
E += chi(sigma[x,y] != sigma[x+1,y])
E += chi(sigma[x,y] != sigma[x,y+1])
end
end
for x in 1:m-1
E += chi(sigma[x,n] != sigma[x+1,n])
E += chi(sigma[x,n] != sigma[x,1])
end
for y in 1:n-1
E += chi(sigma[m,y] != sigma[1,y])
E += chi(sigma[m,y] != sigma[m,y+1])
end
E += chi(sigma[m,n] != sigma[1,n])
E += chi(sigma[m,n] != sigma[m,1])
return E
end
function Potts2dNeachcolor(q::Integer, sigma::AbstractArray{Int64,2})
return [length(filter(x->x==c, sigma)) for c in 1:q]
end
### MCMC simulation by the Wolff method
#
function Potts2dWolff(
q::Integer,
init::Array{Int64,2},
beta::Float64,
niters::Integer,
burnin::Integer,
thin::Integer;
maxtotalms::Integer = 10^7
)
m, n = size(init)
cb = Array{Bool,2}(m,n)
cx = Array{Integer}(m*n)
cy = similar(cx)
# nbr[:,x,y] is the array of the neighbors of (x,y)
nbr = Array{Tuple{Integer,Integer},3}(4,m,n)
for x in 1:m
for y in 1:n
nbr[1,x,y] = (mod1(x-1,m), y)
nbr[2,x,y] = (mod1(x+1,m), y)
nbr[3,x,y] = (x, mod1(y-1,n))
nbr[4,x,y] = (x, mod1(y+1,n))
end
end
# Wolff Update
#
function Potts2dWolffUpdate!(q::Integer, sigma::AbstractArray{Int64,2}, beta::Float64, N::Integer)
p = exp(-2*beta)
local ms::Integer = 0 # number of modified sites
for Niter in 1:N
fill!(cb, false)
k = 1
cx[k], cy[k] = rand(1:m), rand(1:n) # location of seed
cb[cx[k],cy[k]] = true # cn[x,y] is true iff (x,y) is in cluster
ss = sigma[cx[k],cy[k]] # state of seed
cs0 = 0
cs1 = k # cluster size
while cs0 < cs1
for i in cs0+1:cs1
for d in 1:4
x, y = nbr[d,cx[i],cy[i]]
if sigma[x,y] != ss || cb[x,y] == true || rand() < p; continue; end
k += 1
cx[k], cy[k] = x, y
cb[x,y] = true
end
end
cs0 = cs1
cs1 = k
end
s = mod1(rand(1:q-1)+ss,q)
for i in 1:cs1
sigma[cx[i],cy[i]] = s
end
ms += cs1
end
return (Potts2dEnergy(sigma), maximum(Potts2dNeachcolor(q,sigma)), ms)
end
N = div(niters, thin)
@printf("N = %d\n", N)
S = Array{Int64,3}(m,n,N+1)
maxnc = maximum(Potts2dNeachcolor(q,init))
E = Potts2dEnergy(init)
totalms = 0 # number of total modified sites
@printf("initial value: maxnc = %d, E = %.0f\n", maxnc, E)
t = 0
S[:,:,t+1] = init
E, maxnc, ms = Potts2dWolffUpdate!(q, (@view S[:,:,t+1]), beta, burnin)
totalms += ms
@printf("after burnin: maxnc = %d, E = %.0f, totalms = %d\n", maxnc, E, totalms)
for t in 1:N
S[:,:,t+1] = @view S[:,:,t]
E, maxnc, ms = Potts2dWolffUpdate!(q, (@view S[:,:,t+1]), beta, thin)
totalms += ms
if mod(t, div(N,16)) == 0
@printf("t = %d, maxnc = %d, E = %.0f, totalms = %d\n", t, maxnc, E, totalms)
end
if totalms ≥ maxtotalms || maxnc == m*n
@printf("breaked: t = %d, maxnc = %d, E = %.0f, totalcs = %d\n", t, maxnc, E, totalms)
break
end
end
return S[:,:,1:t+1]
end
Out[2]:
In [3]:
### Critical inverse temperature of q-state 2D Potts model
#
# See https://mathtod.online/@ips_ips_ips_ips/495560
#
# 宮下精二・轟木義一『別冊数理科学・演習形式で学ぶ相転移・臨界現象』
#
betacritical(q) = 0.5*log(1+sqrt(q))
@show betacritical(2)
@show betacritical(3)
@show betacritical(4)
;
q=2: Ising¶
In [4]:
q = 2
m, n = 128, 128
init = rand(1:q, m, n)
niters = 2^12
burnin = 0
thin = div(niters,2^4)
beta = betacritical(q)/1.01
@time S = Potts2dWolff(q, init, beta, niters, burnin, thin)
sleep(0.2)
betastr = @sprintf("%.4f", beta)
for k in 1:2:size(S,3)
figure()
pcolormesh(S[:,:,k], cmap="bwr", vmin=1, vmax=q)
title("q = $q, beta = $betastr, t = $(k-1)")
axes()[:set_aspect]("equal")
end
q=3¶
In [5]:
q = 3 # 3-state Potts
m, n = 128, 128
init = rand(1:q, m, n)
niters = 2^9
burnin = 0
thin = div(niters,2^4)
beta = betacritical(q)/eps()
@time S = Potts2dWolff(q, init, beta, niters, burnin, thin)
sleep(0.2)
betastr = @sprintf("%.4f", beta)
for k in 1:1:size(S,3)
figure()
pcolormesh(S[:,:,k], cmap="bwr", vmin=1, vmax=q)
title("q = $q, beta = $betastr, t = $(k-1)")
axes()[:set_aspect]("equal")
end
In [6]:
q = 3
m, n = 128, 128
init = rand(1:q, m, n)
niters = 2^13
burnin = 0
thin = div(niters,2^4)
beta = betacritical(q)/0.99
@time S = Potts2dWolff(q, init, beta, niters, burnin, thin, maxtotalms=10^8)
sleep(0.2)
betastr = @sprintf("%.4f", beta)
for k in 1:2:size(S,3)
figure()
pcolormesh(S[:,:,k], cmap="bwr", vmin=1, vmax=q)
title("q = $q, beta = $betastr, t = $(k-1)")
axes()[:set_aspect]("equal")
end
