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