Juliaで学ぶ古典モンテカルロシミュレーション

イジング模型のモンテカルロシミュレーション

さて、重みつきモンテカルロ法について色々調べたので、次は実際にイジング模型でシミュレーションを実行してみよう。
プログラムを順番に作っていきながら、磁化のヒストグラムやスピン分布のアニメーションを作る。

まず、イジング模型をおさらいする。

ハミルトニアンと分配関数

まず、古典スピン系であるイジング模型のハミルトニアンは $$ 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} $$ となる。

1次元の例

周期境界条件を持つ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$となっている。

隣接格子点の情報

まず、隣接格子点を求める関数を

In [1]:
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
Out[1]:
隣接格子点 (generic function with 1 method)

と定義しよう。これは、サイト${\bf i}$が与えられた時、${\bf i} + {\bf d}_l$の四つを返す。
そして、あるサイト${\bf i}$の周辺スピンの和を

In [2]:
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
Out[2]:
周辺スピン和 (generic function with 1 method)

で計算してみよう。その結果、あるスピン配置における全エネルギーは

In [3]:
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
Out[3]:
全エネルギー (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$が与えられた時、メトロポリス法は

In [4]:
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
Out[4]:
メトロポリス法! (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$にすればよい。
以上より、熱浴法を実行する関数は

In [5]:
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
Out[5]:
熱浴法! (generic function with 1 method)

そして、ランダムにサイト${\bf i}$を選び、配置のアップデートを試みる関数は

In [6]:
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
Out[6]:
配置アップデート! (generic function with 1 method)

となる。ここで、メトロポリス法と熱浴法の両方を選んで用いたいので、flipmethodという関数を引数とした。

モンテカルロ法

以上より、モンテカルロ法の関数は

In [7]:
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})} $$

In [8]:
Tc = 2/log(1+sqrt(2))
Out[8]:
2.269185314213022

である。この温度より低いと、強磁性相となる。
まず、磁化のヒストグラムを見てみよう。毎回毎回に磁化がどのような値になっているかを見る。転移温度以下の場合、格子点あたりの磁化は$+1$と$-1$に偏るはずである。
転移温度以下の場合、

In [9]:
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ヒストグラム)
Out[9]:
-1.0 -0.5 0.0 0.5 1.0 0 250 500 750 1000 y1

となり、-1と+1に集中していることがわかる。つまり、強磁性転移が起きている。

計算の高速化

さて、この計算は手元のパソコンで1分ほどで終わった。しかし、どうやらこれよりも高速化できるようである。Twitterで@kikumacoさんから教わった方法はもっと速いようである。
まず、乱数の中で無数にif文を呼ぶのは、どんな言語でも明らかに遅そうである。よって、ここを改良する。

In [10]:
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
Out[10]:
周辺スピン和 (generic function with 1 method)

もう一度計算してみよう。

In [11]:
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
Out[11]:
-1.0 -0.5 0.0 0.5 1.0 0 250 500 750 1000 y1

速くなった。 次に、周辺スピン和に関して、ifelse文を使ってみる。余りを計算するのとどちらが速いか

In [12]:
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
Out[12]:
周辺スピン和 (generic function with 1 method)
In [13]:
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
Out[13]:
-1.0 -0.5 0.0 0.5 1.0 0 250 500 750 1000