2次元イジングモデル

ホップフィールドネットワークのプログラムを書いたので、2次元イジングモデルのプログラムも書いてみました。イジングモデルは磁性体を表したモデルです。磁性体は磁石だと考えてもらって結構だと思います。磁石は、熱を与えるとその磁性を失います（つまり、鉄がくっつかなくなります）。そんな磁石を理論的に表したモデルがイジングモデルですね。私の記憶では2次元までイジングモデルの厳密解が得られていたと思います。つまり、磁性を失う温度がシミュレーションせずにわかると言うことです。あくまでもモデルなので、おもちゃの一種だと考えてください。

そもそも2Dイジングモデルとは

例えば、10x10のイジングモデルを考えると、その要素は-1,1になります。2次元の配列上のある要素i(x,y)を選ぶと、隣り合う要素は(x+1,y),(x-1,y),(x,y+1),(x,y-1)の4つになります。イジングモデルでは隣り合う要素のみの相互作用をエネルギーに加算します。

エネルギー関数

\( H=-J\sum_{ij} S_iS_j - h\sum_i S_i \)

ijは隣り合う要素とします。Jは相互作用定数で、hは外部地場になります。今回のシミュレーションでは、h=0としています。

シミュレーションの流れ

ランダムに一点を択び符号を逆にします。その後、エネルギーを計算し低ければ状態を遷移させ、高ければメトロポリス法に基づく状態遷移を行います。
例えば、詳細釣合の原理（遷移1->2と遷移2->1の確率は等しいとする原理）を導入すると、

\(\frac{w_{S1S2}}{w_{S2S1}} = exp(-\beta\Delta E_{S1S2}) \)

となるので、メトロポリス法による遷移確率は

\(w_{S1S2}=min[1,exp(-\beta\Delta E_{S1S2})]\)

となります。

• 境界条件
  • 境界条件は周期境界条件となっております。

コード

#!/usr/bin/env python
import copy

import numpy as np
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 = 1

    def main(self):
        self.MC(500)
        
    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 MC(self,loop):
        for i in xrange(loop):
            mp = np.random.randint(0,self.dimx,2)
            txy = copy.deepcopy(self.xy)
            txy[mp[0],mp[1]] = -txy[mp[0],mp[1]]
            Epre = self.calcE(self.xy)
            Eafter = self.calcE(txy)
            if Eafter < Epre:
                self.xy = txy
            else:
                pro = np.random.random()
                if pro < np.exp(-(Eafter-Epre)/self.T):
                    self.xy = txy
            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()
    tmp.main()

20x20の2次元イジングの結果

start->End

2015/11/12:追記

self.neighbors = [[-1,-1],[1,-1],[-1,1],[1,1]]

上記の部分を以下に訂正

self.neighbors = [[-1,0],[1,0],[0,-1],[0,1]]