Juliaを使って非常に単純な系でイジングモデルの数値計算をやってみた

はじめに

統計力学の講義で、3次元のイジング模型について1格子,磁場がかかってないという条件で内部エネルギーなどの物理量の温度依存性を数値計算せよという課題がありました.
実は先日永井先生のJulia言語の本
https://www.kspub.co.jp/book/detail/5282823.html
を買いまして,折角なのでこの課題をJuliaでやってみたというのがこの記事の内容です.

実行環境

一応実行環境を書いておきます.Windows11でJulia 1.7.3を動かし,VScode(1.68.1)上で編集を行いました.

何を計算したのか?

格子を構成する各粒子が$1/2$のスピンを持っているとし,$i$番目の粒子の「スピン」を$S_i=\pm1$とすると系のエネルギーは
$$\mathcal{H}=-J\mu^2\sum_{(i,j)}S_iS_j-H\mu\sum_i S_i$$と書けます。ただし$J$は結合のスピン間の相互作用の強さを表す定数,$\mu$は磁気モーメントの次元を持つ定数,$H$は磁場です.また,$1$項目の$\sum_{(i,j)}$は各隣接原子間について和をとることを表しています.この記事では磁場が$0$の場合を考えるので$2$項目は消えます.また,温度依存性が分かればそれでいいので$J\mu^2=1$とします.カノニカルアンサンブルにおいて物理量$A$の期待値は
$$＜A＞=\sum_{i}\frac{e^{-\beta E_i}}{Z}A_i$$で求まるのでこれを用いて期待値を数値計算します.$i$は各状態のラベルです.また,$Z$は分配関数で,
$$\beta=\frac{1}{k_BT}$$は逆温度です.この逆温度に対する依存性を調べます.

プログラムの説明

(冗談抜きで)初心者なのでおかしな書き方をしている部分があれば教えてください.
まず,次のコードでスピンの状態を表すベクトルを取得します.

function trans_lattice(lattice_index::Int64)::Vector{Int64} #maiking spin state ventor
	# This function returns a vector of spin state of 8 particles by bit operation
	# This function requires Int64 (0 ≤ lattice_index ≤ 255)

	statevec = zeros(Int64,8)

	for i in 0:7 

		spin = (lattice_index >> i) & 1 == 1 ? 1 : -1

		# if i th bit is 1, assign 1 to spin
		# if i th bit is 0, assign -1 to spin

		statevec[8-i]=spin
	end

	return statevec
end

$0$から$255=2^8-1$の数字を2進法で表したときには$2^0$の位から$2^7$の位まで$0$と$1$をおいたすべてのパターンが出現するわけですが,ビット演算によってその値を取得し,$1$は$S=1$に,$0$は$-1$に対応させることで8個の粒子の全スピンのパターンを取得しようというコードです.自分は最初,再帰関数を定義してこのベクトルを取得しようとしていたのですが,(少なくとも自分よりは)プログラミングが得意な友人がこの方法を教えてくれました.型をかなり厳密に指定してますが,そんなに重い計算をするわけではないので,多分要りません
次に,この関数で各状態に対応したエネルギーを求めます.

function energy(statevec::Vector{Int64})::BigFloat 
	# This function calculates energy corresponding to state
	# This function requires Vector{Int64} (state of spin)
	# This function does Not return expected value 

	energy = BigFloat
	
	Sys =Array{Int64,3}(undef,2,2,2) #arranging particles cubic
	Sys[:,:,1]=[statevec[1] statevec[2] ; statevec[3] statevec[4]]
	Sys[:,:,2]=[statevec[5] statevec[6] ; statevec[7] statevec[8]]

	energy = -1*(
		Sys[1,1,1]Sys[2,1,1] + Sys[1,2,1]Sys[2,2,1] + Sys[1,1,2]Sys[2,1,2] + Sys[1,2,2]Sys[2,2,2] +
		Sys[1,1,1]Sys[1,2,1] + Sys[2,1,1]Sys[2,2,1] + Sys[1,1,2]Sys[1,2,2] + Sys[2,1,2]Sys[2,2,2] +
		Sys[1,1,1]Sys[1,1,2] + Sys[2,1,1]Sys[2,1,2] + Sys[1,2,1]Sys[1,2,2] + Sys[2,2,1]Sys[2,2,2]
	)

	return energy
end

途中のSysで現実に即した3次元の立方格子に配置することで相互作用エネルギーを計算しやすくしています.この配置の仕方はもちろん一意ではないですが,各状態に対応したエネルギーを計算するためにこの関数のみを用いている限り問題はないはずです.Array{Int64,3}(undef,2,2,2)はInt64の型をもつ数の成分未決定の2×2×2の配列を表していて,Sys[:,:,1]でこの配列の[i,j,1]成分すべてにアクセスし,
$$\begin{align}\begin{pmatrix}{\rm statevec}[1] & {\rm statevec}[2]\\ {\rm statevec}[3] & {\rm statevec}[4]\end{pmatrix}\end{align}$$という行列(?)を代入しています.
次の関数は分配関数です.

function dist(β::BigFloat)::BigFloat 
	# This function is distribution function
	# This function requires BigFloat as inverse temperature

	sum_weight=big"0"

	for lattice_index=0:255

		Boltzmannweight=exp(-β*energy(trans_lattice(lattice_index)))
		sum_weight += Boltzmannweight

	end

	return sum_weight
end

内部で先に述べたtrans_lattice関数とenergy関数を呼び出しているので単独では存在できません.
次の関数は物理量に相当する関数と逆温度を引数にとり,その逆温度での期待値を出力してくれる関数です.この記事では内部エネルギーの期待値を求めることまでしかしませんが,熱容量や磁化,帯磁率など他の物理量を計算する際に便利です.

function Expect(f,β::BigFloat)::BigFloat # expected value
	# This function calculates expected value
	# This function requires function as physical quantity And BigFloat as inverse temperature 

	Z=dist(β) # distribution function

	sum_value = big"0"

	for lattice_index in 0:255 
		
		statevec = trans_lattice(lattice_index)

		Boltzmannweight=exp(-β*energy(statevec))
		value = f(statevec)*Boltzmannweight/Z    # exp(-βEᵢ)fᵢ/Z

		sum_value += value
	end

	return sum_value
end

以上で準備は終わりで,

Energy(β::BigFloat) = Expect(energy,β)

とすればこのEnergy関数で逆温度$\beta$での内部エネルギーの期待値を計算できます.

betas=big"0.0":big"0.01":big"2"

としておいて,(これは簡単に言えば$0$から$2$までを$0.01$区切りに得られる数字つまり$0,0.01,0.02,\ldots,1.99,2$の集合です.)

using Plots
plot(betas,Energy,label=false,xlabel=:"inverse temperature",ylabel=:"internalenergy")

とすれば次のような図が得られるはずです.
画像は先のプログラムの直後に

savefig("internalenergy.png")

などと書けば保存できます.高温でスピンの向きが乱雑になり,内部エネルギーが$0$になっていることが見て取れます.

さいごに

この記事をみてくれたあなたが少しでもJulia言語に興味を持っていただけたなら幸いです.まだコミュニティが小さく,情報が少ないのが難点ですが,個人的にどんどんと広がっていってほしいと思っています.冒頭に述べた永井先生のJuliaの本には2次元イジング模型をモンテカルロシミュレーションする方法が示されているので,気になる人は買ってみてください!

なんで？

plotするコードのところでこんな問題点がサジェストされるんですけど、なぜですかね…？

(2022-7-5追記)
VScodeの未修正のバグなのかもしれません.
https://github.com/julia-vscode/julia-vscode/issues/1576