[PYTHON] Spin map generation of Ising model by Metropolis-Hastings method (M-H algorithm)

Introduction

The spin map of the two-dimensional square lattice Ising model was created for each temperature by the Metropolis-Hastings method, so it is summarized. The Ising model is only a simplified explanation, so if you want to know more, please refer to the textbook of statistical mechanics. You can also create spin maps such as triangular lattices other than square lattices by this method, so please try it. (Just change the Hamiltonian change)

Metropolis-Hastings algorithm

The Metropolis-Hastings method is a Markov chain Monte Carlo method (MCMC) for obtaining a sequence of random samples from a probability distribution that is difficult to sample directly, and is often used in statistics and statistical physics. This sequence can be used to approximate the distribution, integrate, and calculate the expected value.

algorithm

Set an arbitrary initial state, repeat the update by the following method, and adopt the final state as a plausible state. Suppose a state $ S $ is adopted. Generate the state $ \ acute {S} $ according to the probability density function $ g (\ acute {S} \ mid S) $ as a candidate for update. $ g (\ acute {S} \ mid S) $ is a distribution function of the probability that state $ S $ will move to state $ \ acute {S} $. Here, let the likelihood of the state be $ P (S) , and the new functions $ a_1 = \ frac {P (\ acute {S})} {P (S)}, a_2 = \ frac Introduce {g (S \ mid \ acute {S})} {g (\ acute {S} \ mid S)} $ $ a = a_1a_2 $$. $ a_1 $ is the likelihood ratio (probability ratio) and $ a_2 $ is the ratio of the probability density function. The state $ S $ uses $ a $ to update the state according to the following rules. ・ When $ a \ ge 1 $ Update the state $ S $ to the state $ \ acute {S} $ and proceed to the next step. ・ When $ a \ lt 1 $ Update the state $ S $ to the state $ \ acute {S} $ with the probability $ a $ and proceed to the next step, and next without updating the state with the probability $ 1-a $ Proceed to the step of. The update is repeated until the state is not completely updated.

Ising Model

The Ising model is a model mainly used as a model of a ferromagnet. (It is not a model peculiar to ferromagnets, but works as a model for various phenomena.) This model is a theory that deals only with the interaction between the grid points of interest and their reclining grid points. Although it is such a simple theory, it is an excellent model that can describe the phase transition phenomenon. The Hamiltonian of this model is $ without external force H = -J\sum_{\}\sigma_i\sigma_j Obtained by $. Where J is the exchange interaction energy, $ \ <i, j > $ is the combination of adjacent grid points, and $ \ sigma_k $ is the spin value at the grid point k. The likelihood of the Ising model is expressed using the Hamiltonian, where the Hamiltonian of the system is $ H $ and the temperature is T, then $ \ frac {\ exp (\ frac {H} {k_BT})} {Z_ \ beta Represented as} $. In other words, if the Hamiltonian displacement due to the change of state is $ \ Delta H $, then $ a_1 = \ exp (\ frac {\ Delta H} {k_BT}) $. Further, if the update candidate is a state in which the spin of one randomly selected lattice point is inverted, the probability is constant regardless of the transition from any state to any state. In other words, since the probability density does not change depending on the direction of transition, $ a_2 = 1 $, so $ a = a_1 = \ exp (\ frac {\ Delta H} {k_BT}) $.

Two-dimensional square lattice Ising model

A model in which electrons are present at the grid points of a space stretched in a grid pattern is called a two-dimensional square lattice Ising model. The Hamiltonian of this model is $$, assuming $ J = 1 $ H = -\sum_{<i,j>}\sigma_i\sigma_j $$ Can be expressed as. (This does not change in any grid.) Using this, the change in Hamiltonian when the spin value of one point $ (k, l) $ is inverted is $. \Delta H = -\sum_{i=0}^{1}(\sigma_{k,l}\sigma_{k+2i-1,l}+\sigma_{k,l}\sigma_{k,l+2i-1}) $ Will be. (This depends on the grid.) As this periodic boundary condition, $ \ sigma_ {N, j} = \ sigma_ {0, j}, \ sigma_ {-1, j} = \ sigma_ {N-1, j}, \ sigma_ {i, N} = \ sigma_ {i, 0} and \ sigma_ {i, -1} = \ sigma_ {i, N-1} $ are provided. 二次元正方格子例として左図のような系を考え、格子点(i,j)を反転させたときハミルトニアンはどれほど変化するかを求める。先述の通りイジングモデルでは最隣接格子点のみ考えるので、格子点(i,j)が反転することでスピンの商が変わる格子点の組み合わせは(i,j)と(i+1,j)、(i,j)と((i-1,j)、(i,j)と((i,j+1)、(i,j)と((i,j-1)それぞれ4つのみである。幸いスピンの値は1、-1のみを取るので、ハミルトニアン変化は$\Delta{H}=2(\sigma_{i,j}\sigma_{i+1,j}+\sigma_{i,j}\sigma_{i-1,j}+\sigma_{i,j}\sigma_{i,j+1}+\sigma_{i,j}\sigma_{i,j-1})$の簡単な形式で表せる。

Implementation

The implementation was done in Python. The whole code looks like this:

GenerateSpinState.py


import numpy as np

class GSS:
    def __init__(self, N, T):
        # lattice size
        self.N = N
        # temperature list
        self.T = T
        # generate N*N spin state by random
        self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1

    # calc hamiltonian displacement　which is square lattice, after flipping one site(x,y)
    def calc_Hamiltonian_square(self, x, y):
        delta_H = 2 * self.state[x, y] * (self.state[(x+1)%self.N, y] + self.state[x, (y+1)%self.N] + self.state[(x-1)%self.N, y] + self.state[x, (y-1)%self.N])
        return delta_H

    # flip random spin site
    def flip_spin(self):
        new_state = self.state.copy()
        [x, y] = np.random.randint(0, self.N, (2))
        new_state[x, y] *= -1
        return new_state, x, y

    # calc specious spin state by metropolis method
    def calc_spin_state(self, t):
        n = 10000
        for i in range(n):
            # get flip site
            new_state, x, y = self.flip_spin()
            # calc hamiltonian displacement
            deltaH = self.calc_Hamiltonian_square(x, y)
            # decide whether to adopt
            if np.random.rand() <= np.exp(-deltaH/t):
                self.state = new_state

    def calc_each_temperature(self):
        # save spin state
        X = []
        for t in self.T:
            # init spin state
            self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1
            # generate spin state which temperature t
            self.calc_spin_state(t)
            # append generate spin state which temperature t
            X.append(self.state)
        return np.array(X)

main.py


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

import GenerateSpinState

# save path
# train test data
npy_path = "./result//"
# spin map
spin_map_path = "./result/"


# lattice size
N = 32

# number of data (square lattice)
count = 100

# temperaturelist/Significant digits is 3
T_list = np.linspace(0.1, 5.5, count).round(3)

# make gif which write spin state transition
def create_spinmap_gif(data, T, save_path):
    save_path = save_path + "squarelattice.gif"
    fig = plt.figure()

    def Plot(num):
        plt.cla()
        plt.imshow(data[num])
        t_len = len(str(T[num]))
        title = "temperature:" + str(T[num]) + "0"*(5-t_len)
        plt.title(title)

    anim = FuncAnimation(fig, Plot, frames=100)
    anim.save(save_path, writer='imagemagick')


# call module GSS
GSS = GenerateSpinState.GSS(N, T_list)

# create spin map
X_train = GSS.calc_each_temperature()

# create gif
create_spinmap_gif(X_train, T_list, spin_map_path)

First, I will explain about GenerateSpinState.py. Let the size and temperature of the grid be instance variables. It also generates the initial state of the spin map.

def __init__(self, N, T):
    # lattice size
    self.N = N
    # temperature list
    self.T = T
    # generate N*N spin state by random
    self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1

Find the Hamiltonian change $ \ Delta H $ when the spin value of any grid point $ (x, y) $ is inverted. As a feature of the Ising model, only the re-adjacent grid points are reflected in the calculation, so it is sufficient to extract and calculate only the grid points around the inverted grid points.

# calc hamiltonian displacement　which is square lattice, after flipping one site(x,y)
def calc_Hamiltonian_square(self, x, y):
    delta_H = 2 * self.state[x, y] * (self.state[(x+1)%self.N, y] + self.state[x, (y+1)%self.N] + self.state[(x-1)%self.N, y] + self.state[x, (y-1)%self.N])
    return delta_H

One grid point is randomly selected, and the spin map with the one point inverted is saved as new_state until it is decided whether or not to update it. Don't forget .copy when moving self.state to new_state as a spinmap before inversion. Without .copy, updating new_state will also update self.state and the spinmap will continue to update.

# flip random spin site
def flip_spin(self):
    new_state = self.state.copy()
    [x, y] = np.random.randint(0, self.N, (2))
    new_state[x, y] *= -1
    return new_state, x, y

After calling the above two functions, the state is updated by adopting new_state with a probability of $ \ exp {\ frac {\ Delta {H}} {k_B \ beta}} $. If it is not adopted, the status will not be updated. The number of repetitions is set to 10000, but there is no physical or mathematical basis, and it was set because it was a good value. As a result, I got the most likely result by chance, so I adopted it. (It was not completely updated after 1000 times.)

# calc specious spin state by metropolis method
def calc_spin_state(self, t):
    n = 10000
    for i in range(n):
        # get flip site
        new_state, x, y = self.flip_spin()
        # calc hamiltonian displacement
        deltaH = self.calc_Hamiltonian_square(x, y)
        # decide whether to adopt
        if np.random.rand() <= np.exp(-deltaH/t):
            self.state = new_state

List that stores the temperature you want to check first self.T I made spin maps at various temperatures.

def calc_each_temperature(self):
    # save spin state
    X = []
    for t in self.T:
        # init spin state
        self.state = 2 * np.random.randint(0, 2, (self.N, self.N)) - 1
        # generate spin state which temperature t
        self.calc_spin_state(t)
        # append generate spin state which temperature t
        X.append(self.state)
    return np.array(X)

This concludes the explanation of GenerateSpinState.py. Next, for main.py, the grid size and temperature array are set, and GenerateSpinState.py is called to obtain the spin map after temperature. In this case, the size of the grid is 36, and the temperature array is 0.1-5.5 divided into 100 equal parts. The temperature sequence was created considering that the theoretical phase transition temperature is 2.26. I made the following gif image with create_spinmap_gif. (Spin map adopted at each temperature)![Squarattice.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/591069/bba500b1-f3b6-3561-113b -d7fa35054fbb.gif)

References

Metropolis–Hastings algorithm Hastings, W.K. (1970). "Monte Carlo Sampling Methods Using Markov Chains and Their Applications".

Recommended Posts

Spin map generation of Ising model by Metropolis-Hastings method (M-H algorithm)
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method