[Calculation condition] In the simulation, a grid of $ 20 \ times 20 $ was used, and $ J = 1 $ and $ h = 0 $. The total number of Monte Carlo steps was set to 40,000, and the number of samplings for thermal near value calculation was set to about 8000 (Monte Carlo steps averaged from about 32000).

In addition, the spin-spin interaction of spins at the boundary is [Periodic boundary condition](https://ja.wikipedia.org/wiki/%E5%91%A8%E6%9C%9F%E7%9A%84% E5% A2% 83% E7% 95% 8C% E6% 9D% A1% E4% BB% B6) was used for evaluation.

** Calculation code **

Uses Python 3.5.1. Implemented and confirmed operation on Jupyter notebook 5.0.0.

"""
Monte Carlo simulation by Metropolis method of 2D Ising model

"""

from random import random, randrange
import numpy as np
import matplotlib.pyplot as plt



#Define a function to calculate the energy for a given spin arrangement alis
def Ecalc2(alis2, Nx,Ny):
    dum=0
    for i in range(-1,Nx):
        for j in range(0,Ny):
            
            l=i
            m=j
            
            if l==-1:
                l=Nx
            if l == Nx:
                l=0
            if m ==-1:
                m=Ny
            if m== Ny:
                m=0
            
            ll=i+1
            if ll == Nx:
                ll=0
            dum+=alis2[l, m]*alis2[ll, m] 
            
    for i in range(0,Nx):
        for j in range(-1,Ny):
            
            l=i
            m=j
            
            if l==-1:
                l=Nx
            if l == Nx:
                l=0
            if m ==-1:
                m=Ny
            if m== Ny:
                m=0
            
            mm=j+1
            if mm == Ny:
                mm=0
            dum+=alis2[l, m]*alis2[l, mm] 
        
    return dum 



#Random initial state generation:The number of upspins and the number of downspins are the same. M=0
def Initial_rand(s, Nx,Ny):
    NN=int(Nx*Ny/2)
    for k in range(NN):
        i=randrange(Nx-1) #A random number is chosen from 1 to N
        j=randrange(Ny-1) #A random number is chosen from 1 to N
        s[i,j]=-1*s[i,j]
    
    return s


Nx= 20  #Divided into 100 in x direction
Ny =20  #Divided into 100 in the y direction
Ntot=Nx*Ny


KBT_lis=np.linspace(0.001,6, 201) #Temperature 0.From 001 to 6(In kBT units)It fluctuates in steps of 201.
ETplot=[]
MTplot=[]
CTplot=[]
for KBT in KBT_lis:
    
    J = 1 #Spin coupling constant
    B = 0.0 #External magnetic field
    steps =40000 #MC step
    avsteps=int(steps/5)
    
        
    #Initial state generation:Random spin placement
    s= np.ones([Nx,Ny], int)  #Quantum number for N spins(Sz = +1 or -1)Set as 1
    s=Initial_rand(s, Nx,Ny)
 
    E = -J* Ecalc2(s, Nx,Ny) -B*np.sum(s) #(initial)Calculate energy.
    E2 = E**2  # E^Store 2
    
    #setup
    eplot = []  #Define a list to store energy values under each temperature kBT
    Eavplot=[] #Define a list to store the thermal mean of energy
    
    e2plot = []#Define a list to store the square of the energy value under each temperature kBT. Used in the calculation of specific heat.
    
    eplot.append(E/Ntot)
    Eavplot.append(E/Ntot)
    
    e2plot.append((E/Ntot)**2)
    
    Cvplot = [] ##Define a list to store the specific heat under each temperature kBT
    
    
    M = np.sum(s)  #Magnetization calculation
    Mplot = []
    Mavplot=[]
    Mplot.append(M/Ntot) #Thermal mean value of magnetization
    Mavplot.append(M/Ntot) #Define a list to store the thermal mean of magnetization

    #Main
    for k in range(steps):
        i=randrange(Nx-1) #0 to N-A random number is chosen up to 1
        j=randrange(Ny-1) #0 to N-A random number is chosen up to 1
    
        s_trial=s.copy()
        s_trial[i,j]= -1*s[i,j]
        delta_E=2*s_trial[i,j]*-1*J*(s[i+1,j]+s[i-1,j]+s[i,j+1]+s[i,j-1])-B*(s_trial[i,j]-s[i,j])

        E_trial =E+  delta_E

        #Status update by metropolis method
        if E_trial < E : 
            s = s_trial
            E = E_trial
        else :
            if random() < np.exp(-(delta_E)/KBT):
                s = s_trial
                E = E_trial
   
        eplot.append(E/Ntot)
        e2plot.append((E/Ntot)**2)
        M = np.sum(s)
        Mplot.append(M/Ntot)
    
    ETplot.append(np.sum(eplot[-avsteps:])/avsteps)
    MTplot.append(np.sum(Mplot[-avsteps:])/avsteps)
    CTplot.append((np.sum(e2plot[-avsteps:])/avsteps-((np.sum(eplot[-avsteps:]))/avsteps)**2)/KBT**2)  #Calculation of specific heat



plt.plot(KBT_lis,ETplot)
plt.ylabel("Totla energy per spin")
plt.xlabel("kBT")
plt.show()

plt.plot(KBT_lis,MTplot)
plt.ylabel("Total magnetic moment per spin")
plt.hlines([0], 0, KBT_lis[-1], linestyles="-")  # y=Draw a line at 0.
plt.xlabel("kBT")
plt.show()


plt.plot(KBT_lis,CTplot)

plt.ylabel("Heat capacity per spin")
plt.xlabel("kBT")
plt.show()

** Result (1) Time evolution of energy and magnetization **

An example of how the amount of substance changes with the Monte Carlo step during a Monte Carlo simulation under isothermal conditions is shown below. This cannot be obtained directly from the calculation code. To display it, write the following plot statement without looping the temperature in the calculation code (for KBT in KBT_lis :).

plt.plot(eplot)
plt.ylabel("Totla energy per spin")
plt.xlabel("kBT")
plt.show()

plt.plot(Mplot)
plt.ylabel("Total magnetic moment per spin")
#plt.hlines([0], 0, steps, linestyles="-")  # y=Draw a line at 0.
plt.xlabel("kBT")
plt.show()

$ k_BT = 1 $. It can be seen that the physics converges as the number of steps increases (the state of reaching thermal equilibrium).

Time evolution of energy per spin.

Time evolution of magnetization per spin.

** Result (2) Temperature dependence **

The figure obtained by the simulation is shown below. If you want to obtain more accurate results, you can increase the number of grids in the simulation. The theoretical value of $ T_c $, the phase transition temperature from ferromagnetism to paramagnetism, is $ T_c = 2.2691 $ [see Addendum (4-1) (12)].

Total energy per spin. As the temperature rises, pairs with opposite spins occur, increasing energy.

Temperature dependence of magnetization. Magnetization disappears with heating (phase transition from ferromagnet to paramagnetic). Although it is difficult to see in the figure, strictly speaking, it changes discontinuously at $ T = T_c $ (Secondary phase transition. % BB% A2% E7% A7% BB # .E7.AC.AC.E4.BA.8C.E7.A8.AE.E7.9B.B8.E8.BB.A2.E7.A7.BB)) is there.

Temperature dependence of specific heat. Magnetic transition is [secondary phase transition](https://ja.wikipedia.org/wiki/%E7%9B%B8%E8%BB%A2%E7%A7%BB#.E7.AC.AC.E4 .BA.8C.E7.A8.AE.E7.9B.B8.E8.BB.A2.E7.A7.BB), and the specific heat dissipates at $ T = T_c $. This tendency can be seen from this figure as well.

Addendum

** Addendum (0): Time evolution of spin state **

The change in spin state during the simulation is shown in animation below. The initial spin states were all upward (upspin, $ s_ {i, j} = 1 $), and $ k_BT = 30 $. ** At this temperature, the paramagnetic state becomes stable, and the spin state changes so that the upward and downward spins are about the same **.

Change in spin state of 100x100 grid. ** Pink represents upspin ($ s_ {i, j} = 1 ) and light blue represents downspin ( s_ {i, j} = -1 $). Downspin increases with time **.

** Addendum (1): Theory **

According to quantum statistical mechanics, the thermal mean value of the physical quantity (observable operator) $ \ hat A $ in the thermal equilibrium state of temperature T is $ \ <A > $, where the Hamilton operator of the system is $ \ hat H $. ，

\< A\> = Tr (\hat ρ \hat A) = e^{-\beta \hat H \hat A} /Z \tag{2}

Given in. Here, Tr is the diagonal sum, and $ \ beta = 1 / (k_B T) $ is the quantity called inverse temperature.

$ \ Hat ρ $ is the Density Matrix Operator And in thermal equilibrium \hat ρ = e^{-\beta \hat H} /Z \tag{3} Was used to be.

If you use the display that diagonalizes $ \ hat H $ and $ \ hat A $ at the same time, the diagonal sum is

\< A\> =\sum_{i=1}^{n} A_i \exp(\{ -\beta E_{i} \}), \ \ \ Z =\sum_{i=1}^{n} \exp(\{ -\beta E_{i} \})\tag{4}

Call. $ E_ {i} $ is the $ i $ th energy eigenvalue of $ \ hat H $, and $ Z $ is the partition function E9% 85% 8D% E9% 96% A2% E6% 95% B0).

** Addendum (2): Monte Carlo simulation by Metropolis method **

From Addendum (1), it was found that the thermal mean value of physical quantities under a finite temperature can be theoretically described. This theoretical system will be one of the monuments in theoretical physics. However, from the standpoint of ** actually calculating **, it is generally extremely difficult to calculate the thermal mean value. The reason is that the number of solutions (number of states) of the Schledinger-equation is often enormous and cannot be summed for the states of ** all **.

** Therefore, multiple appropriate states are sampled using random numbers, and $ \ <A > $ is approximately calculated as the average of them (Monte Carlo simulation). ** ** Assuming that the number of samples is N, the $ \ <A_N > _ {MC} $ obtained in this way becomes all pairs as N becomes infinite when the following conditions (1) and (2) are satisfied. It is known that the sum of angles gradually approaches $ \ <A > $, which is calculated exactly. Therefore, if a sufficiently large N is taken, it is possible to accurately obtain the thermal mean value of the physical quantity.

\lim_{N\to\infty} \< A_N \>_{MC} = \< A\> \tag{5}

conditions (1): Monte Carlo simulation is ergodic (2): ** Satisfy "Detailed Balance Conditions" **

(1) is that the state of the system changes with time evolution starting from an arbitrary initial state, but it is possible to go through all possible states (Ergodic condition. .org / wiki /% E3% 82% A8% E3% 83% AB% E3% 82% B4% E3% 83% BC% E3% 83% 89% E7% 90% 86% E8% AB% 96)).

** (2) sets the probability that the state $ i $ is realized in the thermal equilibrium state $ \ pi (i) $ and the probability of transition from the state $ i $ to the state $ j $ (transition probability) $ w (i, j). ) $ \pi(i)\ w(i,j)=\pi(j) \ w(j,i) To meet. ** **

The choice of $ w (i, j) $ is arbitrary. In the metropolis method, the transition probability $ w (i, j) $ from the state i to the state j is set as follows.

w(i,j) =$ 1 (E_j \le E_i)\ $, $ \ e^{-\beta \ (E_j-E_i)} \ (E_j > E_i) \tag{6}$

** This creates various states one after another, and it has been proved that after a sufficient time, the possible states of the system follow the Boltzmann distribution for the canonical ensemble in statistical mechanics [4]. The method of selecting $ w (i, j) $ by the metropolis method therefore has a high affinity with statistical mechanics and is indispensable in computational statistical physics. ** **

** Addendum (3): Procedure for Monte Carlo simulation by Metropolis method of Ising model **

1. Make an arbitrary spin array $ A_k $

2. Randomly pick up the i-th spin to generate a new spin placement $ A_ {k + 1} $

3. Create $ A_ {trial} with the spin direction of spin i inverted (spin flip) to create a trial arrangement, and calculate $ energy $ E_ {tri} $.

4. If $ E_i \ \ le E_ {trial} , accept the trial placement ( A_ {k + 1} $ = $ A_ {trial}) $

5. If $ E_i \ \ gt E_ {trial} $, accept with probability $ P_ {tr} = exp [-(\ beta \ (E_ {trial} -E_i))] $ ($ A_ {k + 1} $ = $ A_ {trial} $)

Where 5 is 5-1. Generate a uniform random number $ r $$ (0 \ le r \ le 1) $ 5-2. Compare the magnitude relationship between $ r $ and $ P_ {tr} $, and if $ P_ {tr} \ ge r $, then $ A_ {k + 1} $ = $ A_ {trial} , otherwise If so, reject the state change (( A_ {k + 1} $ = $ A_k $))

And. 6. Repeat steps 1-5 above until the physical quantity is not so different (until the thermal equilibrium is reached). Energy E, magnetization M, and (equal volume) specific heat Cv are, respectively.

E= -J \ \sum_{i=1}^{n} s_i s_{i+1}\tag{7} M= -J \ \sum_{i=1}^{n} s_i \tag{8} $C_v =\frac{E^2-<E>^2}{k_B T^2}\tag{9} $ Can be given. $ n $ is the number of particles. 7. After reaching thermal equilibrium, calculate the average physical quantity in multiple steps. The number of steps depends on the temperature and the number of particles, but the higher the temperature, the more steps may be required (because the energy fluctuation becomes larger).

** Addendum (4): Other **

** Addendum (4-1) Exact solution of 2D Ising model **

There is an exact solution in the no magnetic field (h = 0) 2D Ising model. For the two-dimensional Hamiltonian given by Eq. (1), the asymptotic solutions of energy E, specific heat Cv, and magnetization M when the number of lattice points is increased are as follows [3]. \frac{E}{N}=-J\ coth(j)\ [1+\frac{2}{\pi}\kappa'\ K(\kappa)] \tag{10}

\frac{Cv}{Nk_B}=\frac{1}{\pi} \frac{j}{sinh(j)}[2K(\kappa)-\frac{2}{\pi}-cosh^2(j) \ E(\kappa)+\kappa' sinh^2(j)\ K(\kappa)]\tag{11}

Magnetization M is $ T \ lt T_c (= \ frac {J} {0.4407 \ k_B}) $

M=[\frac{cosh^2(j)}{sinh^4(j)}\ (sinh^2j-1)]^{\frac{1}{8}} \tag{12} And in $ T \ gt T_c $, $ M = 0 $. ** This $ T_c $ is the phase transition temperature from ferromagnetism to paramagnetism, which simply shows that the two-dimensional Ising model can describe the magnetic phase transition. ** In one dimension, $ M $ is a continuous function and there is no magnetic transition (same thing, but it is sometimes written as $ T_c $ = 0K $).

Where j, $ \ kappa $, $ \ kappa'$ are respectively

j=\frac{2J}{k_BT}, \kappa=2\frac{sinh(j)}{cosh^2(j)}, \kappa'=2\frac{sinh^2(j)-1}{cosh^2(j)} \tag{13} Given in. Also, $ K (\ kappa) $ and $ E (\ kappa') $ are first-class and second-class perfect elliptic integrals, respectively.

{K(\kappa)=\int_0^{\pi/2}{\frac{1}{\sqrt{1-k^2\sin^2\theta}}}d\theta}\tag{13}

{E(k)=\int_0^{\pi/2}{\sqrt{1-k^2\sin^2\theta}}d\theta}\tag{14}

(See [5] for calculations in Python).

** Addendum (4-2) Magnetism of matter **

It is in the behavior of "spin", which is one of the quantum mechanical degrees of freedom of electrons, and is the result of spin-spin competition (spin-spin interaction). Therefore, quantum mechanics is necessary to understand the essence of magnets [1]. The governing equation is the Schrodinger equation (in the non-relativistic limit). Therefore, (A) ** In order to simulate the magnetism of matter, it is necessary to solve the Schrodinger equation considering the spin-spin interaction. ** **

Statistical mechanics is a discipline that deals with the thermodynamic properties of matter at finite temperatures. According to it, quantitative information about magnetism at finite temperature is (B) ** knowing all the solutions of the (many-body) Schrodinger equation with many spins, and averaging with some weight ** You can get it.

Unfortunately, it is extremely difficult to do (A) and (B). The following three points are the main reasons.

** (C) The substance of spin-spin interaction is called exchange interaction derived from the statistical nature of electrons, and the functional form of this interaction is unknown **,

** (D) It is difficult to solve the many-body Schrodinger equation even if the functional form of the interaction is known **

** (E) Even if the Schrodinger equation can be solved, it is very difficult to calculate the thermal statistical mean (calculation of partition function, etc.) **

Because.

Therefore, in the study of magnetism, ** model interaction (model Hamiltonian) is set for (C) and (D), and the Schrodinger equation considering it is solved by numerical calculation. For (E), this paper It has often been done for a long time to reduce the amount of calculation by performing ** efficient sampling ** as dealt with in.

Phase by many-body interaction (electronic correlation) Descripting the transition is still a very challenging task. The Ising model and the Heisenberg model (two-dimensional XY model) that handles it quantum mechanically have played an important role in discussing the thermodynamics (phase transition, etc.) of magnetic materials together with the (quantum) Monte Carlo method [2]. ..

References

1. Kei Yoshida, ["Magnetic"](https://www.amazon.co.jp/%E7%A3%81%E6%80%A7-%E8%8A%B3%E7%94%B0-% E5% A5% 8E / dp / 4000054422), Iwanami Publishing, 1991.

2. Seiji Miyashita, ["Thermal and Statistical Dynamics"](https://www.amazon.co.jp/%E7%86%B1-%E7%B5%B1%E8%A8%88%E5%8A % 9B% E5% AD% A6-% E7% 89% A9% E7% 90% 86% E5% AD% A6% E5% 9F% BA% E7% A4% 8E% E3% 82% B7% E3% 83% AA% E3% 83% BC% E3% 82% BA-% E5% AE% AE% E4% B8% 8B-% E7% B2% BE% E4% BA% 8C / dp / 4563023841 / ref = sr_1_6? S = books & ie = UTF8 & qid = 1502402477 & sr = 1-6 & keywords =% E5% AE% AE% E4% B8% 8B% E7% B2% BE% E4% BA% 8C), Bakufukan, 1993.

3. Akira Morita and Yutaka Okabe, ["Simulation Statistical Physics"](https://www.amazon.co.jp/Windows%E3%83%A6%E3%83%BC%E3%82%B6%E3% 83% BC% E3% 81% AE% E3% 81% 9F% E3% 82% 81% E3% 81% AE-% E3% 82% B7% E3% 83% 9F% E3% 83% A5% E3% 83 % AC% E3% 83% BC% E3% 82% B7% E3% 83% A7% E3% 83% B3% E7% B5% B1% E8% A8% 88% E7% 89% A9% E7% 90% 86 -% E6% A3% AE% E7% 94% B0-% E7% AB% A0 / dp / 4621042033 / ref = sr_1_1? S = books & ie = UTF8 & qid = 1502402499 & sr = 1-1 & keywords =% E3% 82% B7% E3% 83% 9F% E3% 83% A5% E3% 83% AC% E3% 83% BC% E3% 82% B7% E3% 83% A7% E3% 83% B3% E7% B5% B1% E8% A8% 88% E7% 89% A9% E7% 90% 86), Maruzen Co., Ltd., 1996.

4. Qiita's article by kaityo256 "Detailed balance conditions in Monte Carlo method and how to determine transition probability" is very easy to understand and helpful. It was. recommendation.

5. Qiita article "[Scientific / technical calculation by Python] List of usage of (special) functions used in physics by using scipy" by Python I am writing about the calculation of perfect elliptic integrals.

Postscript

August 12, 2017 We received the following valuable comments from Mr. kaityo256. It will be helpful. Thank you very much.

*** "As you may know, spin-based Monte Carlo simulations generally use cluster update methods such as the Swendsen-Wang and Wolff algorithms. It's a super-powerful method compared to single-spin flips. It's a method. I wrote a simple article about Swendsen-Wang, so I hope you find it helpful. ***

http://qiita.com/kaityo256/items/6539261993e282edc5aa "