続*3：2次元イジングモデル - おいも貴婦人ブログ

2015/11/13修正

def calcDeltaEの部分

return -2*dE

def MCのエネルギー差の足し算部分

E+=dE

修正済みのコードは文中。

本文

今までのプログラムは、エネルギーを計算する際に全ての要素に対してエネルギーの再計算を行っていました。しかしイジングの場合、あるスピンの符号が反転するとそのスピンに関するエネルギーの符号も反転します。そうすると、スピンをひっくり返えす前の状態iとひっくり返した後の状態jのエネルギー差は、あるスピン周りのエネルギー差の２倍、つまりdE=Ej-Ei=2Ejとなります。この結果により、Pythonでも1000x1000のイジングモデルのシミュレーションが一瞬で終わってしまいます。以下、pythonのコードです。気付くのが遅くなり恥ずかしい限りです!!!

#!/usr/bin/env python                                                                                                             
import copy

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

class Ising2D:
    def __init__(self,x=20,y=20,seed=0):
	np.random.seed(seed)
	self.dimx = x
	self.dimy = y
	self.xy = np.random.choice([-1,1],x*y)
	self.xy = self.xy.reshape(x,y)
	self.E = 0
	self.h = 0
	self.J = -1
	self.neighbors = [[-1,0],[1,0],[0,1],[0,-1]]
	self.T = 0.1

    def main(self):
	self.MC(10000)

    def calcE(self,xy):
	E = 0
	for x in xrange(self.dimx):
            for y in xrange(self.dimx):
		for neighbor in self.neighbors:
                    nx = (x + neighbor[0]) % self.dimx
                    if(nx < 0):
			nx = self.dimx
                    ny = (y + neighbor[1]) % self.dimy
                    if(ny < 0):
			ny = self.dimy
                    E += -self.J*xy[x][y]*xy[nx][ny]
	return E

    def calcDeltaE(self,point):
	dE = 0
	x = point[0]
	y = point[1]
	for neighbor in self.neighbors:
            nx = (x + neighbor[0]) % self.dimx
            if(nx < 0):
		nx = self.dimx
            ny = (y + neighbor[1]) % self.dimy
            if(ny < 0):
		ny = self.dimy
            dE += -self.J*self.xy[x][y]*self.xy[nx][ny]
	preE = dE
	afterE = -dE
	dE = afterE-preE
	return -2*dE

    def MC(self,loop):
	E = self.calcE(self.xy)
	for i in xrange(loop):
            mp = np.random.randint(0,self.dimx,2)
            dE = self.calcDeltaE(mp)

            if dE<=0:
		self.xy[mp[0]][mp[1]] = -self.xy[mp[0]][mp[1]]
		E+=dE
            else:
                pro = np.random.random()
		if pro < np.exp(-(dE)/self.T):
                    self.xy[mp[0]][mp[1]] = -self.xy[mp[0]][mp[1]]
                    E+=dE
            print E
        self.show(i)

    def show(self,index):
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.matshow(self.xy)
        filename = "ising%05d.png" % (index)
        plt.savefig(filename)


if __name__ == "__main__":
    tmp = Ising2D(x=1000,y=1000)
    tmp.main()

pythonでのプログラム

1024x1024のイジングモデルを1000ステップ計算してみました。最初の総エネルギーの計算がなければもっと早く済みます。

結果

環境(1024*1024サイズ)	python
Time	6.950s