[理系大学生応援] Wolffの手法でイジングモデルを解く

概要

この記事は Wolff, "Collective Monte Carlo Updating for Spin Systems", 1989 のアルゴリズムでイジングモデルを数値的に解くことを目的とした記事です。
「いきなり複雑なアルゴリズムでイジングモデルを解くのはちょっと...」という方は、過去記事
[理系大学生応援] シングルスピンフリップでイジングモデルを解く
[理系大学生応援] Swendsen-Wangの手法でイジングモデルを解く
でイジングモデルの理解を深めてから、この記事を読んでいただければと思います。

アルゴリズム

STEP 0

適当に初期状態のスピン配置(ここでは例としてランダムとします)を用意します。

STEP 1

ランダムに座標$(i, j)$のスピンを選択します。

STEP 2

STEP 1 で選択したスピンをpocket, clusterとして登録します。

前のSTEPですでに登録されているclusterがあればこれを空にしておきます。

STEP 3

pocketの中から1つのスピンに着目し、それと隣あうスピンneighborに対して以下のループを行います。

STEP 3-1

以下の条件がすべて満たされるとき、neighborにあったスピンをpocketとclusterに登録します。

(1) 隣あうスピンの向きが同じ
(2) 隣にあったスピンがclusterに登録されていない
(3) $0<r<1$の乱数が確率$p=1-e^{-2J/(k_B T)}$より小さい

STEP 4

STEP 3で選択したスピンをpocketから取り除きます。

STEP 5

pocketが空っぽになるまでSTEP 2 ~ STEP 4を繰り返します。

STEP 6

以上のSTEPでclusterに登録されたスピンをフリップさせます。

Pythonスクリプト

以下にPythonで書かれた実装スクリプトを示します。

# !/usr/bin/env python


import matplotlib.pyplot as plt
import numpy as np
import random


class IsingWolff:
	def __init__(self):
		# set the number of spins in x-direrection, y-direction
		self.IXMAX = 200
		self.JYMAX = 200
		# set J_ij (= J, uniform)
		self.J = 1.0
		# set initial condition (-1 or 1, random)
		self.S = (np.random.rand(self.IXMAX, self.JYMAX)<0.5) * (-2) + 1
		# set 1/(k_B T)
		self.ONEOKBT = 1.0
		# set the number of iteration 
		self.ITER = 10000
		# set (i, j) neighbor
		self.NEIGHBOR = {}
		for i in range(self.IXMAX):
			for j in range(self.JYMAX):
				self.NEIGHBOR[(i, j)] = [((i+1)%self.IXMAX, j), ((i-1)%self.IXMAX, j), (i, (j+1)%self.JYMAX), (i, (j-1)%self.JYMAX)]
		
	def Compute_Probability(self):
		# calcurate probability of flip
		return 1.0 - np.exp(-2.0*self.J*self.ONEOKBT)

	def Cluster_Flip(self, p):
		# set random coordinate (i, j)
		coord = (np.random.randint(self.IXMAX), np.random.randint(self.JYMAX))
		# set pocket
		pocket = [coord]
		# set cluster
		cluster = [coord]
		# search pocket while pocket is not empty
		while pocket != []:
			# choose random tuple from pocket
			l = random.choice(pocket)
			# for loop for neightborhood
			for m in self.NEIGHBOR[l]:
				# if spins are parallel & (i, j) is not in cluster and r < p, add to pocket & cluster
				if self.S[m[0]][m[1]] == self.S[l[0]][l[1]]:
					if m not in cluster:
						if np.random.rand() < p:
							pocket.append(m)
							cluster.append(m)
			# remove choosen (i, j) from pocket list
			pocket.remove(l)
		# flip spin including cluster
		for l in cluster:
			self.S[l[0]][l[1]] *= -1

	def Plot_Spin(self, n):
		# remove previous figure
		plt.clf()
		# plot color map
		plt.imshow(self.S, cmap='PuOr', vmin=-1, vmax=1, animated=True)
		# add colorbar
		plt.colorbar(ticks=[-1, 1], orientation='vertical')
		# set pause interval time
		plt.pause(1.0e-6)
		# save figure as png format
		plt.savefig('fig'+str(n)+'.png')
		
	def run(self):
		# set cluster flip probability
		p = self.Compute_Probability()
		# main loop start
		for n in range(self.ITER):
			# plot
			if n % 1 == 0:
				self.Plot_Spin(n)
			# flip cluster
			self.Cluster_Flip(p)
		
if __name__ == '__main__':
	a = IsingWolff()
	a.run()

スクリプト詳細

__init__

self.IXMAX, self.JYMAXでそれぞれ$x, y$方向のスピンの数を指定しています。self.Jは相互作用係数を、self.ONEOKBTは$1/k_B T$の大きさ、self.ITERでモンテカルロの試行回数を決定しています。self.NEIGHBORに$(i, j)$のスピンと隣あうスピン$(i+1, j), (i-1, j), (i, j+1), (i, j-1)$を辞書型で登録しています。

Compute_Probability

確率$p$を計算しています。

Cluster_Flip

上述のアルゴリズムのSTEP 1 ~ STEP 6 を実行しています。

Plot_Spin

この関数でmatplotlibを用いてself.S($\pm 1$のスピンの並び)を2次元カラーコントアとして描画しています。また描画結果をjpegファイルとして保存しています。

run

この関数でWolffの方法によるスピン反転をself.ITERの回数だけ繰り返し行います。

結果

self.IXMAX = 200
self.JYMAX = 200
self.J = 1.0
self.ONEOKBT = 1.0
self.ITER = 10000		

上述のパラメータで計算した結果を動画で示します。

Wolffの手法でイジングモデルを解く、動画はこちら pic.twitter.com/6rVJQuLj08

— 中の人はダンブルドア (@tweet_nakasho) September 12, 2019

$k_B T = 1$のときを計算しているので、スピンが綺麗に揃う結果が示せていることがわかります。

結言

この方法もSwendsen-Wangの手法と同様、シングルスピンフリップよりも強力な方法です。ですがWolffのアルゴリズムは、Swendsen-Wangのアルゴリズムをより大きなクラスターでスピン交換が起きる確率を高くした改良版です。どんどん活用していきましょう。