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)
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.
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)
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
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
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)
Metropolis–Hastings algorithm Hastings, W.K. (1970). "Monte Carlo Sampling Methods Using Markov Chains and Their Applications".