さて、重みつきモンテカルロ法について色々調べたので、次は実際にイジング模型でシミュレーションを実行してみよう。
プログラムを順番に作っていきながら、磁化のヒストグラムやスピン分布のアニメーションを作る。
まず、イジング模型をおさらいする。
まず、古典スピン系であるイジング模型のハミルトニアンは
$$
H = -J \sum_{\langle i j \rangle} \sigma_i \sigma_j - h \sum_i \sigma_i
$$
である。第二項は磁場の効果である。ここで、$\langle i j \rangle$は、最隣接格子点のみで和を取ることを意味しており、一次元系であれば、$j=i+1$などである。
$\sigma_i$は$i$番目の格子点のスピンを表し、$+1$か$-1$を取る。
統計力学において、物理量$A$の期待値は
$$
\langle A \rangle = \frac{1}{Z} \sum_{\cal C} \left[ \exp \left( -\frac{H({\cal C})}{k_{\rm B} T} \right) A({\cal C}) \right]
$$
と書ける。ここで、$H({\cal C})$は、あるスピン配置${\cal C}$でのハミルトニアン、$A({\cal C})$はその時の物理量$A$の値である。
$k_{\rm B}$はボルツマン定数、$T$は温度である。
$Z$は分配関数であり、
$$
Z = \sum_{\cal C}\exp \left( -\frac{H({\cal C})}{k_{\rm B} T} \right)
$$
と定義されている。
つまり、すべての可能なスピン配置に関して和を取れば、物理量が計算できる。
すべての可能なスピン配置の数${\cal N}$は、$N$個の格子点を持つ系の場合、各サイトで$-1$か$1$の二通りの状態を取れるので、
$$
{\cal N} = 2^N
$$
である。
$L_x \times L_y$の正方格子の二次元イジング模型を考えよう。この時、あるサイトを${\bf i} = (i_x,i_y)$とすると、 その最近接格子は、 $$ {\bf d}_1 = (i_x+1,i_y),{\bf d}_2 = (i_x-1,i_y),{\bf d}_3 = (i_x,i_y+1), {\bf d}_4 = (i_x,i_y-1) $$ の4点である。この時、イジング模型は $$ H = -\frac{J}{2} \sum_{\bf i}^{L_x L_y} \sum_{l=1}^4 \sigma_{\bf i} \sigma_{{\bf i} + {\bf d}_l}- h \sum_{\bf i} \sigma_{\bf i} $$ となる。ここで、本来一つしかない$\sigma_1 \sigma_2$を$\sigma_1 \sigma_2 = (\sigma_1 \sigma_2 + \sigma_2 \sigma_1)/2$と 分離したため、$1/2$の因子が先頭についた。そして これは、 $$ H = -\frac{J}{2} \sum_{\bf i}^{L_x L_y} \sigma_{\bf i} \sum_{l=1}^4 \sigma_{{\bf i} + {\bf d}_l}- h \sum_{\bf i} \sigma_{\bf i} $$ $$ = -\frac{J}{2} \sum_{\bf i}^{L_x L_y} \sigma_{\bf i} S_i - h \sum_{\bf i} \sigma_{\bf i} $$ $$ S_i = \sum_{l=1}^4 \sigma_{{\bf i} + {\bf d}_l} $$ と書くことができる。よって、あるサイト${\bf i}$の隣接格子点におけるスピンの和$S_i$がそれぞれわかれば、全エネルギー$H$を計算できる。
あるサイト${\bf i}$の一つのスピンをフリップ$(\sigma_{\bf i} \rightarrow -\sigma_{\bf i})$した時、そのエネルギー差$\Delta E$は $$ \Delta E = 2J \sigma_{\bf i} S_i+ 2h \sigma_{\bf i} $$ となる。
周期境界条件を持つ4格子からなる1次元のイジング模型を考える。簡単のため磁場はゼロ($h=0$)とする。この時、ハミルトニアンは
$$
H = -J (\sigma_1 \sigma_2 + \sigma_1 \sigma_4 + \sigma_2 \sigma_3 + \sigma_3 \sigma_4)
$$
である。これは、
$$
H = -\frac{J}{2} (\sigma_1 \sigma_2 + \sigma_2 \sigma_1 + \sigma_1 \sigma_4 + \sigma_4 \sigma_1 + \sigma_2 \sigma_3 + \sigma_3 \sigma_2 + \sigma_3 \sigma_4+ \sigma_4 \sigma_3)
$$
$$
= -\frac{J}{2} \left( \sigma_1 (\sigma_2 + \sigma_4) + \sigma_2 (\sigma_1 + \sigma_3) + \sigma_3 (\sigma_4 + \sigma_2) + \sigma_4 (\sigma_1 + \sigma_3)\right)
$$
のように、$i$を固定して隣接格子の和をとったものに変更できる。
また、サイト$2$をフリップした場合は、エネルギーは
$$
H' = -J (-\sigma_1 \sigma_2 + \sigma_1 \sigma_4 - \sigma_2 \sigma_3 + \sigma_3 \sigma_4)
$$
となり、その差は
$$
H'-H = -J (-2\sigma_1 \sigma_2 - 2\sigma_2 \sigma_3 ) = 2J \sigma_2 (\sigma_1 + \sigma_3)
$$
となり、ちゃんと上で計算した$\Delta E$となっている。
まず、隣接格子点を求める関数を
function 隣接格子点(ix,iy,Lx,Ly,周期境界)
const lmax = 4 #最近接格子点の数
ls = Array{Array{Int64,1}}(lmax)
ds = ((+1,0),(-1,0),(0,+1),(0,-1))
count = 0
for d in ds
jx = ix + d[1]
jy = iy + d[2]
if 周期境界
jx += ifelse(jx>Lx,-Lx,0)
jx += ifelse(jx<1,Lx,0)
jy += ifelse(jy>Ly,-Ly,0)
jy += ifelse(jy<1,Ly,0)
end
if 1 <= jx <= Lx && 1 <= jy <= Ly
count += 1
ls[count] = [jx,jy]
end
end
return ls[1:count]
end
隣接格子点 (generic function with 1 method)
と定義しよう。これは、サイト${\bf i}$が与えられた時、${\bf i} + {\bf d}_l$の四つを返す。
そして、あるサイト${\bf i}$の周辺スピンの和を
function 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
#ls = 隣接格子点(ix,iy,Lx,Ly,周期境界)
S = 0
if ix == 1
left = スピン配置[Lx,iy]
else
left = スピン配置[ix-1,iy]
end
if ix == Lx
right = スピン配置[1,iy]
else
right = スピン配置[ix+1,iy]
end
if iy == 1
down = スピン配置[ix,Ly]
else
down = スピン配置[ix,iy-1]
end
if iy == Ly
up = スピン配置[ix,1]
else
up = スピン配置[ix,iy+1]
end
S = (left+right+down+up)
return S
end
周辺スピン和 (generic function with 1 method)
で計算してみよう。その結果、あるスピン配置における全エネルギーは
function 全エネルギー(スピン配置,周期境界,J,h)
Lx = size(スピン配置,1)
Ly = size(スピン配置,2)
E = 0.0
for ix in 1:Lx
for iy in 1:Ly
σi = スピン配置[ix,iy]
S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
E += -(J/2)*σi*S-h*σi
end
end
return E
end
全エネルギー (generic function with 1 method)
とかける。
書いた当初は隣接格子点という関数を読んで周辺スピン和の計算を行っていたが、そのまま書いたほうが3倍速かったので、隣接格子点という関数は使わないことにした。
あるサイトをランダムに選び、そのスピンをフリップさせることでスピン配置を更新するとする。
この時、全サイト数を$N = L_x \times L_y$とすると、確率$1/N$でサイト数を選ぶ。そして、このようなスピンフリップであればプロセスは対称であるので、メトロポリス法でのある配置$C_i$から$C_j$への採択率は
これは
$$
A(C_i \rightarrow C_j) = {\rm min} \left(1, \frac{P(C_j)}{P(C_i)}\right)
$$
となる。ここで、$P(C_i)$をボルツマン重み
$$
P(C_i) = \exp \left( -\frac{H({\cal C}_i)}{k_{\rm B} T} \right)
$$
と選べば、物理量を重みつきモンテカルロ法で計算できる。
スピン配置${\cal C}_k$のあるサイト${\bf i}$のスピン$\sigma_{\bf i}$をフリップさせスピン配置${\cal C}_k'$とする時、
採択率に現れる重みの比は
$$
\frac{P(C_k')}{P(C_k)} = \exp \left( -\frac{(H({\cal C}_k)-H({\cal C}_k'))}{k_{\rm B} T} \right)
$$
$$
= \exp \left( -\frac{\Delta E({\cal C}_k,{\bf i})}{k_{\rm B} T}\right)
$$
となる。
よって、エネルギー差$\Delta E$が与えられた時、メトロポリス法は
function メトロポリス法!(σi,ΔE,T,E,mz)
σ_new = ifelse(rand() <= exp(-ΔE/T),-σi,σi)
E,mz,accept = ifelse(σ_new == -σi,(E+ΔE,mz-2*σi,1),(E,mz,0))
return σ_new,E,accept,mz
end
メトロポリス法! (generic function with 1 method)
となる。ここで、フリップしたスピンが受け入れられた時、そのエネルギー$E$、全体磁化$m_z$を計算している。
熱浴法では、遷移確率は条件付き確率で計算される。
ランダムにサイトが選ばれる場合の条件付き確率を考えよう。
ランダムに選んだサイト${\bf i}$でのスピン状態を$\sigma_{\bf i}$、それ以外のサイトのスピン状態を${\bf \sigma}({\cal C}_{k,{\bf i/}})$とする。
このスピン配置${\cal C}_k$は${\cal C}_k = (\sigma_{\bf i},{\bf \sigma}({\cal C}_{k,{\bf i/}}))$とおける。
サイト${\bf i}$以外のスピン配置が${\bf \sigma}({\cal C}_{k,{\bf i/}})$の時、サイト${\bf i}$でスピン状態$\sigma_{\bf i}$が選ばれる条件付き確率を$P(\sigma_{\bf i}|{\bf \sigma}({\cal C}_{k,{\bf i/}}))$すれば、
スピン配置${\cal C}_k = (\sigma_{\bf i},{\bf \sigma}({\cal C}_{k,{\bf i/}}))$から
${\cal C}_k' = (+1,{\bf \sigma}({\cal C}_{k,{\bf i/}}))$に遷移する確率は
$$
T_{{\cal C}_k \rightarrow {\cal C}_k'} = P(+1|{\bf \sigma}({\cal C}_{k,{\bf i/}}))
$$
である。ここで、条件付き確率は、
$$
P(+1|{\bf \sigma}({\cal C}_{k,{\bf i/}})) = \frac{P((+1,{\bf \sigma}({\cal C}_{k,{\bf i/}})))}{P({\bf \sigma}({\cal C}_{k,{\bf i/}}))}
$$
とかける。ここで、$P({\bf \sigma}({\cal C}_{k,{\bf i/}}))$はスピン配置${\bf \sigma}({\cal C}_{k,{\bf i/}})$が実現する確率で、スピン配置${\bf \sigma}({\cal C}_{k,{\bf i/}})$を持ちサイト${\bf i}$の可能な状態の確率の和となるため、
$$
P({\bf \sigma}({\cal C}_{k,{\bf i/}})) = \sum_{\sigma_{\bf i}=\pm 1} P(\sigma_{\bf i},{\bf \sigma}({\cal C}_{k,{\bf i/}}))
$$
となる。よって、
$$
P(+1|{\bf \sigma}({\cal C}_{k,{\bf i/}}) = \frac{P((+1,{\bf \sigma}({\cal C}_{k,{\bf i/}})))}
{P((+1,{\bf \sigma}({\cal C}_{k,{\bf i/}}))) + P((-1,{\bf \sigma}({\cal C}_{k,{\bf i/}})))
}
$$
$$
= \frac{1}
{1 + P((-1,{\bf \sigma}({\cal C}_{k,{\bf i/}}))/P((+1,{\bf \sigma}({\cal C}_{k,{\bf i/}}))))}
$$
となる。そして、
$$
P((-1,{\bf \sigma}({\cal C}_{k,{\bf i/}}))/P((+1,{\bf \sigma}({\cal C}_{k,{\bf i/}})))) =
\exp \left( -\frac{\Delta E({\cal C}_k,{\bf i};+1 \rightarrow -1)}{k_{\rm B} T}\right)
$$
となる。ここで、$\Delta E({\cal C}_k,{\bf i};+1 \rightarrow -1)$はサイト${\bf i}$のスピンが$+1$から$-1$となった時のエネルギー差で、
$$
\Delta E({\cal C}_k,{\bf i};+1 \rightarrow -1) = 2J S_i+ 2h = \Delta E({\cal C}_k,{\bf i}) \sigma_{\bf i}
$$
となる。
少し式を整理すると、熱浴法におけるサイト${\bf i}$のスピンが$+1$になる確率は
$$
T_{{\cal C}_k \rightarrow {\cal C}_k'} = \frac{1}{2}\left( 1 + \tanh \left( -\frac{\Delta E({\cal C}_k,{\bf i};+1 \rightarrow -1)}{2k_{\rm B} T}\right) \right)
$$
となる。つまり、一様乱数$r$がこの値以下であればスピンを$+1$に、そうでなけば$-1$にすればよい。
以上より、熱浴法を実行する関数は
function 熱浴法!(σi,ΔE,T,E,mz)
α = σi*ΔE
ratio = (1 + tanh(α/2T))/2 #アップスピンの受け入れ確率
σ_new = ifelse(rand() <= ratio,+1,-1)
E,mz,accept = ifelse(σ_new == -σi,(ΔE,-2*σi,1),(E,mz,0))
return σ_new,E,accept,mz
end
熱浴法! (generic function with 1 method)
そして、ランダムにサイト${\bf i}$を選び、配置のアップデートを試みる関数は
function 配置アップデート!(スピン配置,周期境界,flipmethod,T,J,h,E,mz)
Lx = size(スピン配置,1)
Ly = size(スピン配置,2)
ix = rand(1:Lx)
iy = rand(1:Ly)
σi = スピン配置[ix,iy]
S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
ΔE = 2J*σi*S+2h*σi
スピン配置[ix,iy],E,accept,mz = flipmethod(σi,ΔE,T,E,mz)
#println(accept)
return スピン配置,E,mz,accept
end
配置アップデート! (generic function with 1 method)
using ProgressMeter
function モンテカルロ法!(スピン配置,周期境界,flipmethod,T,J,h,最大ステップ,熱化ステップ)
prog = Progress(最大ステップ,1)
E = 全エネルギー(スピン配置,周期境界,J,h)
mz = sum(スピン配置)
accept = 0
受けいれ総数 = 0
totalcount = 最大ステップ-熱化ステップ
Mz = Array{Float64}(totalcount)
Energy = Array{Float64}(totalcount)
比熱 = Array{Float64}(totalcount)
受けいれ確率 = Array{Float64}(totalcount)
Mz2 = Array{Float64}(totalcount) #磁化の二乗
mzヒストグラム = zeros(Int64,length(スピン配置)*2+1)
mzsum = 0.0
mz2sum = 0.0
energy = 0.0
energy2 = 0.0
count = 0
for i=1:最大ステップ
for j=1:length(スピン配置)
スピン配置,E,mz,accept = 配置アップデート!(スピン配置,周期境界,flipmethod,T,J,h,E,mz)
end
if rand() < 0.01
スピン配置 = -スピン配置
mz = -mz
accept = 1
end
if i > 熱化ステップ
mzヒストグラム[mz+length(スピン配置)+1] += 1
count += 1
mzsum += mz/length(スピン配置)
mz2sum += (mz/length(スピン配置))^2
energy += E
energy2 += E^2
受けいれ総数 += accept
Mz[count] = mzsum/count
Mz2[count] = mz2sum/count
Energy[count] = energy/count
比熱[count] = (energy2/count-Energy[count]^2)/T^2
受けいれ確率[count] = 受けいれ総数/count
end
next!(prog)
end
return Mz[1:count],Mz2[1:count],Energy[1:count],比熱[1:count],受けいれ確率[1:count],mzヒストグラム
end
using Plots
となる。 関数の中では、時々全反転をするようにしている。また、格子点の数だけフリップを試みることを「1モンテカルロステップ」と呼び、スピンフリップの回数は「最大ステップ」$\times$LxLyとなる。
さて、実際に計算してみよう。
まず、二次元イジング模型には$h=0$での厳密解が存在しており、強磁性転移の温度$T_c$は
$$
T_c = \frac{2J}{\ln (1 + \sqrt{2})}
$$
Tc = 2/log(1+sqrt(2))
2.269185314213022
である。この温度より低いと、強磁性相となる。
まず、磁化のヒストグラムを見てみよう。毎回毎回に磁化がどのような値になっているかを見る。転移温度以下の場合、格子点あたりの磁化は$+1$と$-1$に偏るはずである。
転移温度以下の場合、
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = 1.0
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
となり、-1と+1に集中していることがわかる。つまり、強磁性転移が起きている。
さて、この計算は手元のパソコンで1分ほどで終わった。しかし、どうやらこれよりも高速化できるようである。Twitterで@kikumacoさんから教わった方法はもっと速いようである。
まず、乱数の中で無数にif文を呼ぶのは、どんな言語でも明らかに遅そうである。よって、ここを改良する。
function 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
S = スピン配置[(ix-2+Lx)%Lx+1,iy]+スピン配置[ix%Lx+1,iy]+スピン配置[ix,(iy-2+Ly)%Ly+1]+スピン配置[ix,iy%Ly+1]
return S
end
周辺スピン和 (generic function with 1 method)
もう一度計算してみよう。
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = 1.0
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 99%|█████████████████████████████████████████| ETA: 0:00:00
28.486910 seconds (600.35 M allocations: 11.938 GiB, 4.58% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:28
速くなった。 次に、周辺スピン和に関して、ifelse文を使ってみる。余りを計算するのとどちらが速いか
function 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
S = 0
jx = ifelse(ix==Lx,1,ix+1)
jy = iy
S += スピン配置[jx,jy]
jx = ifelse(ix==1,Lx,ix-1)
jy = iy
S += スピン配置[jx,jy]
jy = ifelse(iy==Ly,1,iy+1)
jx = ix
S += スピン配置[jx,jy]
jy = ifelse(iy==1,Ly,iy-1)
jx = ix
S += スピン配置[jx,jy]
# S = スピン配置[(ix-2+Lx)%Lx+1,iy]+スピン配置[ix%Lx+1,iy]+スピン配置[ix,(iy-2+Ly)%Ly+1]+スピン配置[ix,iy%Ly+1]
return S
end
周辺スピン和 (generic function with 1 method)
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = 1.0
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 99%|████████████████████████████████████████ | ETA: 0:00:00
25.459910 seconds (600.35 M allocations: 11.938 GiB, 5.49% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:25
こちらの方が少し速いようだ。
しかし、妙にメモリをたくさん使っている。そこで、1スイープを関数にまとめてみよう。
function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,J,h,E,mz)
accept = 0
acc = 0
for j=1:Lx*Ly
スピン配置,E,mz,acc = 配置アップデート!(スピン配置,周期境界,flipmethod,T,J,h,E,mz)
accept += acc
end
return スピン配置,E,mz,accept
end
using ProgressMeter
function モンテカルロ法!(スピン配置,周期境界,flipmethod,T,J,h,最大ステップ,熱化ステップ)
prog = Progress(最大ステップ,1)
E = 全エネルギー(スピン配置,周期境界,J,h)
mz = sum(スピン配置)
accept = 0
受けいれ総数 = 0
totalcount = 最大ステップ-熱化ステップ
Mz = Array{Float64}(totalcount)
Energy = Array{Float64}(totalcount)
比熱 = Array{Float64}(totalcount)
受けいれ確率 = Array{Float64}(totalcount)
Mz2 = Array{Float64}(totalcount) #磁化の二乗
mzヒストグラム = zeros(Int64,length(スピン配置)*2+1)
mzsum = 0.0
mz2sum = 0.0
energy = 0.0
energy2 = 0.0
count = 0
Lx = size(スピン配置,1)
Ly = size(スピン配置,2)
for i=1:最大ステップ
スピン配置,E,mz,accept = スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,J,h,E,mz)
if rand() < 0.01
スピン配置 = -スピン配置
mz = -mz
accept = 1
end
if i > 熱化ステップ
mzヒストグラム[mz+length(スピン配置)+1] += 1
count += 1
mzsum += mz/length(スピン配置)
mz2sum += (mz/length(スピン配置))^2
energy += E
energy2 += E^2
受けいれ総数 += accept
Mz[count] = mzsum/count
Mz2[count] = mz2sum/count
Energy[count] = energy/count
比熱[count] = (energy2/count-Energy[count]^2)/T^2
受けいれ確率[count] = 受けいれ総数/count
end
next!(prog)
end
return Mz[1:count],Mz2[1:count],Energy[1:count],比熱[1:count],受けいれ確率[1:count],mzヒストグラム
end
モンテカルロ法! (generic function with 1 method)
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = 1.0
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 99%|█████████████████████████████████████████| ETA: 0:00:00
27.350637 seconds (600.54 M allocations: 11.941 GiB, 5.58% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:27
まだメモリを使っている。もう少しいじってみる。
function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,J,h,E,mz)
accept = 0
acc = 0
#for j=1:Lx*Ly
for jx in 1:Lx
for jy in 1:Ly
ix = rand(1:Lx)
iy = rand(1:Ly)
σi = スピン配置[ix,iy]
S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
ΔE = 2J*σi*S+2h*σi
σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz)
スピン配置[ix,iy] = σ_new
accept += acc
end
end
return スピン配置,E,mz,accept
end
スイープ! (generic function with 1 method)
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = 1.0
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 96%|███████████████████████████████████████ | ETA: 0:00:00
11.506572 seconds (437.58 k allocations: 18.585 MiB, 0.05% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:11
これで、メモリ使用量が激減した。
さらに、磁場がゼロのときは、ΔEが整数なので、あらかじめボルツマン因子を計算しておくことで高速化が期待できる。ここでJでエネルギーの次元を測ることにして、$J=1$とする。
function boltzmann(T)
w = zeros(Float64,8*2+1)
for ΔE in -8:8
w[ΔE+9] = exp(-ΔE/T)
end
return w
end
boltzmann (generic function with 1 method)
Juliaでは、同じ関数名で引数の型や数などが違うものを定義できる(多重ディスパッチ)。したがって、磁場がゼロのときに関数を以下に定義する。
function メトロポリス法!(σi,ΔE,T,E,mz,w)
σ_new = ifelse(rand() <= w[ΔE+9],-σi,σi)
E,mz,accept = ifelse(σ_new == -σi,(E+ΔE,mz-2*σi,1),(E,mz,0))
return σ_new,E,accept,mz
end
function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w) #磁場hとカップリングJを除去。wを追加
accept = 0
acc = 0
#for j=1:Lx*Ly
for jx in 1:Lx
for jy in 1:Ly
ix = rand(1:Lx)
iy = rand(1:Ly)
σi = スピン配置[ix,iy]
S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
ΔE = 2σi*S
σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz,w)#wを追加
スピン配置[ix,iy] = σ_new
accept += acc
end
end
return スピン配置,E,mz,accept
end
スイープ! (generic function with 2 methods)
これを使って、「モンテカルロ法!」を定義する。
function モンテカルロ法!(スピン配置,周期境界,flipmethod,T,最大ステップ,熱化ステップ)
prog = Progress(最大ステップ,1)
w = boltzmann(T) #追加した部分
E = 全エネルギー(スピン配置,周期境界,J,h)
mz = sum(スピン配置)
accept = 0
受けいれ総数 = 0
totalcount = 最大ステップ-熱化ステップ
Mz = Array{Float64}(totalcount)
Energy = Array{Float64}(totalcount)
比熱 = Array{Float64}(totalcount)
受けいれ確率 = Array{Float64}(totalcount)
Mz2 = Array{Float64}(totalcount) #磁化の二乗
mzヒストグラム = zeros(Int64,length(スピン配置)*2+1)
mzsum = 0.0
mz2sum = 0.0
energy = 0.0
energy2 = 0.0
count = 0
Lx = size(スピン配置,1)
Ly = size(スピン配置,2)
for i=1:最大ステップ
スピン配置,E,mz,accept = スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w)
if rand() < 0.01
スピン配置 = -スピン配置
mz = -mz
accept = 1
end
if i > 熱化ステップ
mzヒストグラム[mz+length(スピン配置)+1] += 1
count += 1
mzsum += mz/length(スピン配置)
mz2sum += (mz/length(スピン配置))^2
energy += E
energy2 += E^2
受けいれ総数 += accept
Mz[count] = mzsum/count
Mz2[count] = mz2sum/count
Energy[count] = energy/count
比熱[count] = (energy2/count-Energy[count]^2)/T^2
受けいれ確率[count] = 受けいれ総数/count
end
next!(prog)
end
return Mz[1:count],Mz2[1:count],Energy[1:count],比熱[1:count],受けいれ確率[1:count],mzヒストグラム
end
モンテカルロ法! (generic function with 2 methods)
この関数を使って、改めて計算すると、
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = 1.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 34%|██████████████ | ETA: 0:00:06
8.818178 seconds (420.74 k allocations: 18.504 MiB, 0.08% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:09
20パーセントほど速くなった。
さて、次にTc以上の温度を見てみよう。
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = 3.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 99%|█████████████████████████████████████████| ETA: 0:00:00
9.115108 seconds (389.03 k allocations: 16.669 MiB, 0.07% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:09
磁化の分布が0に集まっており、明らかに有限の磁化はない。
転移温度ではどうなるだろうか?
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200
T = Tc
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 94%|███████████████████████████████████████ | ETA: 0:00:00
8.468671 seconds (396.89 k allocations: 16.759 MiB, 0.03% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:08
ヒストグラムがギザギザしていて、サンプル数が足りない感じになっている。したがって、サンプル数を10倍増やしてみよう。
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000
熱化ステップ = 200
T = Tc
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 100%|█████████████████████████████████████████| ETA: 0:00:00
82.377186 seconds (4.06 M allocations: 163.545 MiB, 0.03% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:01:22
少し分布が綺麗になった。しかし、計算時間がかかるので、もう少し速くしたい。
乱数の発生が遅い可能性を考え、ランダムに格子点を選ぶのではなく、順番に選んでみよう。また、乱数をif文の中で生成しないようにした。
function メトロポリス法!(σi,ΔE,T,E,mz,w,r)
σ_new = ifelse(r <= w[ΔE+9],-σi,σi)
E,mz,accept = ifelse(σ_new == -σi,(E+ΔE,mz-2*σi,1),(E,mz,0))
return σ_new,E,accept,mz
end
function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w) #磁場hとカップリングJを除去。wを追加
r = rand(Lx,Ly)
accept = 0
acc = 0
#for j=1:Lx*Ly
for jx in 1:Lx
for jy in 1:Ly
ix = jx#rand(1:Lx)
iy = jy#rand(1:Ly)
σi = スピン配置[ix,iy]
S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
ΔE = 2σi*S
#σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz,w)#wを追加
σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz,w,r[ix,iy])#wを追加
スピン配置[ix,iy] = σ_new
accept += acc
end
end
return スピン配置,E,mz,accept
end
スイープ! (generic function with 2 methods)
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000
熱化ステップ = 200
T = Tc
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 94%|███████████████████████████████████████ | ETA: 0:00:01
17.058670 seconds (4.24 M allocations: 7.616 GiB, 4.99% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:17
さらに高速化できた。
最後に、100万ステップで計算してみる。
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 1000000
熱化ステップ = 200
T = Tc
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
push!(time,i-Lx*Ly-1)
end
plot(time./(Lx*Ly),mzヒストグラム)
Progress: 100%|█████████████████████████████████████████| ETA: 0:00:01
170.671098 seconds (42.15 M allocations: 76.154 GiB, 4.69% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:02:51
100万ステップでも5分以内に計算することができるようになった。
臨界点において、分布にピークが存在している。これは、磁化の絶対値の期待値をとった場合に、有限の値が残ることを意味している。
では、このピークのサイズ依存性はどうなっているだろうか?
サイズを変化させてプロットしてみる。なお、システムサイズが変わるとヒストグラムのビンのサイズが変わるため(サイズが大きいと一つあたりのビンに入る数が少なくなる)、数がサイズに比例するようにした。
Lx = 50
Ly = 50
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
T = Tc
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム1 = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time1 = Int64[]
mzヒストグラム1 = length(mzヒストグラム1)*mzヒストグラム1/sum(mzヒストグラム1)
for i in 1:length(mzヒストグラム1)
push!(time1,i-Lx*Ly-1)
end
time1 = time1./(Lx*Ly)
plot(time1,mzヒストグラム1)
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム2 = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
mzヒストグラム2 = length(mzヒストグラム2)*mzヒストグラム2/sum(mzヒストグラム2)
time2 = Int64[]
for i in 1:length(mzヒストグラム2)
push!(time2,i-Lx*Ly-1)
end
time2 = time2./(Lx*Ly)
plot!(time2,mzヒストグラム2)
Lx = 150
Ly = 150
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム3 = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
mzヒストグラム3 = length(mzヒストグラム3)*mzヒストグラム3/sum(mzヒストグラム3)
time3 = Int64[]
for i in 1:length(mzヒストグラム3)
push!(time3,i-Lx*Ly-1)
end
time3 = time3./(Lx*Ly)
plot!(time3,mzヒストグラム3)
Progress: 98%|████████████████████████████████████████ | ETA: 0:00:01
44.810197 seconds (42.12 M allocations: 19.740 GiB, 5.45% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:45 Progress: 100%|█████████████████████████████████████████| ETA: 0:00:01
171.863945 seconds (42.15 M allocations: 76.154 GiB, 4.57% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:02:52 Progress: 100%|█████████████████████████████████████████| ETA: 0:00:01
458.226056 seconds (42.39 M allocations: 170.231 GiB, 6.10% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:07:38
システムサイズが大きくなるに従って、ピークの位置は0に近づいていく。
それ以外に特徴はあるだろうか?
実は、臨界温度では、この分布は横軸が$L_{x}^{(-1/8)}$でスケールされるらしい。これを確認してみよう。
Lx = 50
Ly = Lx
plot((time1)*(Lx^(1/8)),mzヒストグラム1,title="1000000 steps",label=string(Lx)*"x"*string(Ly))
Lx = 100
Ly = Lx
plot!((time2)*(Lx^(1/8)),mzヒストグラム2,title="1000000 steps",label=string(Lx)*"x"*string(Ly))
Lx = 150
Ly = Lx
plot!((time3)*(Lx^(1/8)),mzヒストグラム3,title="1000000 steps",label=string(Lx)*"x"*string(Ly))
綺麗に重なった!
Juliaでは、@animateで簡単にGIFアニメーションファイルを作ることができる。そこで、
function モンテカルロ法可視化!(スピン配置,周期境界,flipmethod,T,最大ステップ,filename)
E = 全エネルギー(スピン配置,周期境界,J,h)
w = boltzmann(T) #追加した部分
mz = sum(スピン配置)
accept = 0
prog = Progress(最大ステップ,1)
maps = @animate for i=1:最大ステップ
スピン配置,E,accept,mz = スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w)
heatmap(1:Lx, 1:Ly, スピン配置,aspect_ratio=:equal)
next!(prog)
end every 100
gif(maps, "./"*filename, fps = 15)
end
モンテカルロ法可視化! (generic function with 1 method)
として、スピン配置の発展を見てみよう。
低温では
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
T = 1.0
J = 1.0
filename = "ising_fps15T1.gif"
@time モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:21
21.945034 seconds (58.81 M allocations: 3.786 GiB, 3.78% gc time)
INFO: Saved animation to /Users/nagai/Isingmodel/ising_fps15T1.gif
臨界温度では
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
T = Tc
J = 1.0
filename = "ising_fps15Tc.gif"
@time モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:22
22.960519 seconds (58.83 M allocations: 3.786 GiB, 3.67% gc time)
INFO: Saved animation to /Users/nagai/Isingmodel/ising_fps15Tc.gif
高温では
Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
T = 3.0
J = 1.0
filename = "ising_fps15T3.gif"
@time モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:22
23.541917 seconds (58.84 M allocations: 3.787 GiB, 3.67% gc time)
INFO: Saved animation to /Users/nagai/Isingmodel/ising_fps15T3.gif
このように、簡単にシミュレーションが実行できる。 大きさサイズでの振る舞いを見てみよう。低温で、先ほどより長い時間シミュレーションしている。
Lx = 200
Ly = 200
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000
T = 1.0
J = 1.0
filename = "ising_fps15T1_200.gif"
@time モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)
Progress: 100%|█████████████████████████████████████████| Time: 0:06:09
378.791396 seconds (823.06 M allocations: 88.898 GiB, 5.27% gc time)
INFO: Saved animation to /Users/nagai/Isingmodel/ising_fps15T1_200.gif
転移温度より低いので、全体が同じスピンの向きに揃っている。
そして、最後に物理量の温度依存性を見てみよう。量としては、磁化の絶対値を見てみる。
srand(123)
Lx = 50
Ly = 50
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000
熱化ステップ = 200
Mztdep = Float64[]
比熱tdep = Float64[]
tdep = Float64[]
@time for i in 1:20
T = (21-i)*0.2
Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
push!(Mztdep,sqrt(Mz2[end]))
push!(比熱tdep,比熱[end])
push!(tdep,T)
end
plot(tdep,Mztdep)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:05 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 100%|█████████████████████████████████████████| Time: 0:00:04 Progress: 49%|████████████████████ | ETA: 0:00:02
87.623169 seconds (82.50 M allocations: 39.453 GiB, 6.65% gc time)
Progress: 100%|█████████████████████████████████████████| Time: 0:00:04
磁化の絶対値はちゃんと転移温度以下で立ち上がっている。比熱は、
plot(tdep,比熱tdep)
となり、綺麗なλ転移をしている。
Juliaによる数値計算の速度は、Fortranとどのくらい違うのか?
http://bb.phys.se.tmu.ac.jp/~bb/pukiwiki/index.php?MCsim
に二次元イジング模型のFortranコードがあったので、使わせていただくことにする。これをifortでコンパイルして、速度を見てみよう。計算状況は上と同じにする。つまり、熱化ステップを200、ステップ数を100000とする。正確にベンチマークをとったわけではないので、参考ということでみておいてほしい。
Fortranコードは手動でパラメータを入力する部分があり、そこに数秒ほどかかっているとして、time a.outで測定してみた。
filename = "./dc.dat"
fd = open( filename, "r" )
温度f = zeros(Float64,20)
比熱f = zeros(Float64,20)
磁化f = zeros(Float64,20)
cnt = 0
while !eof(fd)
cnt += 1
line = readline(fd)
if cnt >3
u = split(line)
温度f[cnt-3] = parse(Float64,u[1])
比熱f[cnt-3] = parse(Float64,u[3])
磁化f[cnt-3] = parse(Float64,u[7])
end
end
plot(温度f,磁化f,marker=:circle,label="Fortran")
plot!(tdep,Mztdep,marker=:circle,label="Julia")
結果はほとんど一致しているので問題はないだろう。少しずれているのはきになるが、Fortranコードの詳細を読んでいないので、ここでは差異については目をつぶるとする。
そして、Fortranの計算時間は、インプットする時間を10秒と見積もると、ちょうど100秒ほどであった。
物理量を測定するタイミングなどによっても速度は変化するので、どちらが速いかははっきり結論付けられないが、JuliaがFortranと同程度のスピードを出していることがわかった。