Classical Monte Carlo Simulations for the 2D Ising Model and Wolff Algorithm

1. Classical Monte Carlo Simulations for the 2D Ising Model

Model
The Hamiltonian of Ising model is written as
$\begin{matrix} H = - J \sum_{⟨ i j ⟩} s_{i} s_{j} \end{matrix}$
where $⟨ i, j ⟩$ sums over all nearest-neighbor pairs on the lattice and each site $i$ has a spin variable $σ_{i}$ that can take two values: +1 (representing an up spin $| ↑ ⟩$ ) or -1 (representing a down spin $| ↓$ ). The behavior of the model is determined by the competition between the interaction energy and the thermal fluctuations. At low temperatures, the system will exhibit ferromagnetic order, where all spins align in the same direction. At high temperatures, the system will exhibit disordered behavior, where the spins are randomly aligned. Therefore, it is expected that there may exist a finite-temperature phase transition in this model. This classical model despite its simpleness provides us with profound insights into the behavior of magnets through its concise mathematical structure and important physical properties.

Markov chain Monte Carlo simulation

Given a discrete probability distribution ${p (x)}$ , the expectation value of the observable $O (x)$ is given by
$⟨ O (x) ⟩ = \sum_{x} p (x) O (x)$
However, this may be expensive or even non-available when the value range of $x$ is huge.

One alternative way to address it is sampling, which gives an approximate estimation of $O (x)$ . To be specific, if we can generate different samples (configurations) of $x$ according to their corresponding probabilities, and sample $x$ for $N$ times, then the weak law of large numbers guarantees that
$⟨ O (x) ⟩ = lim_{N \to \infty} \frac{1}{N} \sum_{j} O (x_{j})$
where $x_{j}$ is our $j$ th sample.

The Monte Carlo (MC) simulation is basically to generate different samples on computer with some pseudo random number generator and calculate the observable. One of the popular ways to generate samples is to use the Markov chain, and this kind of MC is also called the Markov chain Monte Carlo (MCMC).

Mathematically, the Markov chain means the ( $n + 1$ )th sample $x_{n + 1}$ only depends on the $n$ th sample $x_{n}$ and some random term $ξ_{n}$ for generating samples, i.e.

$x_{n + 1} = F (x_{n}, ξ_{n})$
the master equation for this stochastic process is thus

$P_{n + 1} (x^{'}) = \sum_{x_{n}} ω (x^{'} | x) P_{n} (x)$
where $P_{n} (x_{n})$ denote the probability that we obtain $x_{n}$ as the $n$ th sample and
$ω (x^{'} | x)$ is the conditional probability to transfer from $x$ to $x^{'}$ . If we can expect a equilibrium on the Markov chain, i.e. $lim_{n \to \infty} P = P_{eq}$ then one brutal restriction is to require

$ω (x^{'} | x) P_{eq} (x) = ω (x | x^{'}) P_{eq} (x^{'})$
This constraint is called the detailed balance condition. If this can be satisfied, then the equilibrium must exists.

There are many ways to satisfy it, and we consider the Metropolis-Hasting algorithm here. First we rewrite
$ω (x^{'} | x) = T (x^{'} | x) A (x^{'} | x)$
where $T (x^{'} | x)$ is a defined trial probability. $A (x^{'} | x)$ is called the acceptance probability and is defined as
$A (x^{'} | x) = min {1, \frac{P_{eq} (x^{'}) T (x | x^{'})}{P_{eq} (x) T (x^{'} | x)}}$
If $T (x^{'} | x) = T (x | x^{'})$ , one can check that the detailed balance condition is naturally satisfied. Otherwise, we should further determine the functional form of $T(x’|x)$. As we will see later, in our simulation for the Ising model, we can make $T (x^{'} | x) = T (x | x^{'})$ throughout the simulations, therefore
$A (x^{'} | x) = min {1, \frac{P_{eq} (x^{'})}{P_{eq} (x)}}$
Another important thing for a MCMC scheme is the ergodicity, i.e. we have to make sure that for any possible $x$ , we can sample it from the Markov chain in a finte (discrete) time $n$ . However, this may not be easy to prove in some cases.

MCMC on the Ising model
With MCMC, we can now calculate the expectation values of observables such as energy of the Ising model by sampling different configurations (classical states). At the equilibrium, the probability distribution for different configurations $α_{j}$ obeys the Boltzmann distribution

$P_{eq} (α_{j}) = \frac{e^{- E_{j} / k_{b} T}}{\sum_{j} e^{- E_{j} / k_{b} T}}$
where $k_{b}$ is the Boltzmann constant and $E_{j}$ is the energy corresponding to $α_{j}$ . Suppose we have a configuration $x_{n} = α_{i}$ at time $n$ , then to generate a new configuration $α_{j}$ , we randomly choose one spin and flip it. Therefore, $T (α_{j} | α_{i}) = T (α_{i} | α_{j})$ is satisfied and for $α_{i}$ and $α_{j}$ , the acceptance probability reduces to
$A (α_{j} | α_{i}) = min {1, e^{- Δ E / k T}}$
where $Δ E = E_{j} - E_{i}$ . It means that if the new configuration has a lower energy ( $Δ E < 0$ ), then we accept it with 100% probability to assign $x_{n + 1} = α_{j}$ . Otherwise, we accept it with probability equal to $e^{- Δ E / k T}$ . It is straightforward that the single-flip strategy enable us to reach every possible configuration and the ergodicity can be saitisfied.

Practically, we have to assign an initial configuration, which could be $| ↑↑ \dots ↑ ⟩$ or any random one. In low-temperature case, the initial state configuration may be some high-level excited state, which has small probability to appear statistically. Consequently, if we do measurements at the beginning, it will introduce errors since the initial states are out of equilibrium. Therefore, we first take $N_{Thm}$ steps for thermalizing the system, then do measurements for $N_{Stat}$ times.

Observables
To investigate the property of the Ising model, several observables are measured. In the case of the ferromagnetic Ising model,the most important quantity to calculate is the magnetization. It should be computed using all the spins to take advantage of self-averaging to improve the statistics，which is
$m = \frac{1}{N} \sum_{i = 1}^{N} σ_{i}$
At low temperature, two spin configurations $| ↑ \dots ↑ ⟩$ and $| ↓ \dots ↓ ⟩$ lead to $m = + 1$ and $m = - 1$ respectively. Therefore we can take $⟨ | m | ⟩$ to detect the phase transition.

Another observable is magnetic susceptibility, which describes the response of a material to an external magnetic field, and it measures the degree of magnetization of a material under the action of a magnetic field. The magnetic susceptibility is usually indicated by the symbol $χ$ , which is defined as

$χ = \frac{N}{T} (⟨ m^{2} ⟩ - {⟨ | m | ⟩}^{2})$

Magnetic susceptibility plays an important role in material research and application. By measuring the magnetic susceptibility, we can understand the magnetic properties, phase transition behavior and magnetic permeability of the material. Magnetic susceptibility is also closely related to physical phenomena such as magnetic interaction, spin structure and magnetic field effect of materials.

Binder ratio is a dimensionless physical quantity used to characterize scale invariance in phase transitions. Near the critical point of the phase transition, it can be used to measure the statistical characteristics of the magnetization of systems of different sizes. It can be defined as

$R_{2} = \frac{⟨ m^{4} ⟩}{⟨ m^{2} ⟩^{2}}$

where $⟨ m^{2} ⟩$ is the second moment representing the magnetization and $⟨ m^{4} ⟩$ is the fourth moment of magnetization. At $T_{C}$ , the power-laws cancel out, and the ratio is L-independent (and also universal), up to sub-leading finite-size corrections.

The Binder ratio can be used to determine whether the system has undergone a phase transition near the phase transition point. When the system is in an ordered phase, the Binder ratio tends to 1, and when the system is in a disordered phase, the Binder ratio tends to some constant (usually much less than 1). By calculating the Binder ratio of a system of different sizes and observing whether it converges, it is possible to determine whether the system has a phase transition.

Binder Cumulant is a statistical physical quantity used to study the behavior of phase transitions, often used to describe isotropic magnetic systems or other phase transitions. This physical quantity is proposed in order to determine whether a phase transition occurs in a system and to estimate the position of the phase transition point in numerical simulation.In the case of a scalar order parameter (e.g., for the Ising model), the Binder cumulant is defined as

$U_{2} = \frac{2}{3} (1 - \frac{1}{3} R_{2})$

The main characteristic of the Binder cumulant is that it tends to a constant near the point of the phase transition and to a different value in the non-phase transition region. By observing the change of Binder cumulants with temperature or other external parameters, the phase transition points can be well determined.

2. Program

The following is the program of the classical monte carlo simulations for the 2D Ising Model in fortran 90.

!--------------------------------------------------!
Module configuration
!--------------------------------------------------!
save
integer :: L ! system size
integer :: N ! total number of sites
real(8) :: beta ! inverse temperature
real(8) :: Temperature ! temperature 
real(8),allocatable :: conf(:,:) ! store Ising configuration

end Module configuration

!--------------------------------------------------!
Module measurement
!--------------------------------------------------!
save
real(8) :: E ! energy
real(8) :: sz ! magnetization
real(8) :: sz2 ! magnetization squared 
real(8) :: sz4 ! magnetization to the fourth power
real(8) :: Sus ! susceptibility
real(8) :: Binder ! binder cumulant
real(8),allocatable :: den(:) ! density on each site

end Module measurement

!--------------------------------------------------!
Module control
!--------------------------------------------------!
save
integer :: nbin ! number of bins
integer :: nsweep ! number of sweeps in one bin
integer :: nequil ! number of bins to reach equilbrium

end Module control

!==================================================!
Program IsingMC
!==================================================!
use configuration
use measurement
implicit none

call Ising()
end program

!==================================================!
Subroutine Ising()
!==================================================!
use configuration
use measurement
use control
implicit none
integer :: nb,ns,ne

call initran(1)
call initial()

do while(Temperature<5.0)

    call tcollect()
    beta=1/Temperature

    call initialconf()

    do ne=1,nequil
        do ns=1,nsweep
            call update()
        enddo
    end do

    do nb=1,nbin
        call initialobser()
        do ns=1,nsweep
            call measure()
            call update()
        enddo
        call measureave()
        call output()
    enddo

    Temperature=Temperature+0.05
end do

call deallocateall()

end Subroutine Ising

!--------------------------------------------------!
Function energy()
!--------------------------------------------------!
use configuration
implicit none
real(8) :: energy
integer :: x,y
integer :: i1,i2,j3,j4

energy=0.d0
do x=1,L ; do y=1,L 

    i1=modulo(x,L)+1
    i2=modulo(x-2,L)+1
    j3=modulo(y,L)+1
    j4=modulo(y-2,L)+1

    energy=energy-conf(x,y)*(conf(i1,y)+conf(i2,y)+conf(x,j3)+conf(x,j4))
enddo ; enddo
energy=energy/N 

return
end Function energy 

!--------------------------------------------------!
Function energy_diff(i,j)
!--------------------------------------------------!
use configuration
implicit none
real(8) :: energy_diff
integer :: i,j
integer :: i1,i2,j3,j4

i1=modulo(i,L)+1
i2=modulo(i-2,L)+1
j3=modulo(j,L)+1
j4=modulo(j-2,L)+1

energy_diff=2*conf(i,j)*(conf(i1,j)+conf(i2,j)+conf(i,j3)+conf(i,j4))

return
end Function energy_diff

!--------------------------------------------------!
Subroutine initialconf()
!--------------------------------------------------!
use configuration
implicit none
real(8),external :: ran
integer :: i,j

do i=1,L ; do j=1,L
    if(ran()>0.5d0) then
        conf(i,j)=1
    else
        conf(i,j)=-1
    endif
enddo ; enddo 

end Subroutine initialconf

!--------------------------------------------------!
Subroutine initial()
!--------------------------------------------------!
use configuration
use measurement
use control
implicit none

open(11,file='para.in',status='old')
read(11,*) L
read(11,*) Temperature
read(11,*) nbin,nsweep,nequil
close(11)

N=L*L  

allocate(conf(L,L)) 
allocate(den(N))

end Subroutine initial

!--------------------------------------------------!
Subroutine initialobser()
!--------------------------------------------------!
use measurement

E=0.d0
sz=0.d0
sz2=0.d0
sz4=0.d0
Sus=0.d0
Binder=0.d0
den=0.d0

end Subroutine initialobser

!--------------------------------------------------!
Subroutine update()
!--------------------------------------------------!
use configuration
implicit none
real(8),external :: ran
real(8),external :: energy_diff
integer :: i,j

do i=1,L ; do j=1,L

    if(ran()<min(1.d0,exp(-beta*energy_diff(i,j)))) then
        conf(i,j)= - 1.d0 * conf(i,j)
    end if 

enddo ; enddo

end Subroutine update

!--------------------------------------------------!
Subroutine measure()
!--------------------------------------------------!
use configuration
use measurement
implicit none
integer :: i,j
real(8),external :: energy
real(8) :: sz_tmp

E=E+energy()

sz_tmp=0
do i=1,L ; do j=1,L
    sz_tmp=sz_tmp+conf(i,j)
enddo ; enddo
sz=sz+abs(sz_tmp/N)
sz2=sz2+(sz_tmp/N)**2
sz4=sz4+(sz_tmp/N)**4

Sus=Sus+(sz2-sz**2)/Temperature
Binder=Binder+(1-sz4/3/(sz2**2))

end Subroutine measure

!--------------------------------------------------!
Subroutine measureave()
!--------------------------------------------------!
use configuration
use measurement
use control
implicit none

E=E/nsweep

sz=sz/nsweep
sz2=sz2/nsweep
sz4=sz4/nsweep

Sus=N*(sz2-sz**2)/Temperature
Binder=(1-sz4/3/(sz2**2))*3/2

end Subroutine measureave

!--------------------------------------------------!
Subroutine output()
!--------------------------------------------------!
use configuration
use measurement
use control
implicit none

open(12,file='energy.dat',position='append')
write(12,'((e16.8,2x))') dble(E)
close(12)

open(13,file='M.dat',position='append')
write(13,'((e16.8,2x))') dble(sz)
close(13)

open(14,file='M2.dat',position='append')
write(14,'((e16.8,2x))') dble(sz2)
close(14)

open(15,file='M4.dat',position='append')
write(15,'((e16.8,2x))') dble(sz4)
close(15)

open(16,file='Sus.dat',position='append')
write(16,'((e16.8,2x))') dble(Sus)
close(16)

open(17,file='Binder.dat',position='append')
write(17,'((e16.8,2x))') dble(Binder)
close(17)

end Subroutine output

!--------------------------------------------------!
Subroutine tcollect()
!--------------------------------------------------!
use configuration
implicit none

open(18,file='T.dat',position='append')
write(18,'((e16.8,2x))') dble(Temperature)
close(18)

end Subroutine tcollect

!--------------------------------------------------!
Subroutine deallocateall()
!--------------------------------------------------!
use configuration
use measurement

deallocate(conf)
deallocate(den)

end Subroutine deallocateall

!--------------------------------------------------!
real(8) function ran()
!--------------------------------------------------!
! 64-bit congruental generator                     !
! iran64=oran64*2862933555777941757+1013904243     !
!--------------------------------------------------!
implicit none

real(8)    :: dmu64
integer(8) :: ran64,mul64,add64
common/bran64/dmu64,ran64,mul64,add64

ran64=ran64*mul64+add64
ran=0.5d0+dmu64*dble(ran64)

end function ran
!--------------------------------------------------!

!--------------------------------------------------!
subroutine initran(w)
!--------------------------------------------------!
implicit none

integer(8) :: irmax
integer(4) :: w,nb,b

real(8)    :: dmu64
integer(8) :: ran64,mul64,add64
common/bran64/dmu64,ran64,mul64,add64

irmax=2_8**31
irmax=2*(irmax**2-1)+1
mul64=2862933555777941757_8
add64=1013904243
dmu64=0.5d0/dble(irmax)

open(10,file='seed.in',status='old')
read(10,*)ran64
close(10)

if (w.ne.0) then
open(10,file='seed.in',status='unknown')
write(10,*)abs((ran64*mul64)/5+5265361)
close(10)
endif

end subroutine initran
!--------------------------------------------------!

3. Wolff Algorithm

The Metropolis algorithm with single spin inversion will suffer from critical slowing down when processing near the critical point. The critical region will form a large number of clusters or magnetic domains and the acceptance probability is very small. Effective inversion will only occur when picking to the boundary, so the correlation length between the two sampled configurations will be very large.

Therefore we use the classic cluster algorithm, wolff algorithm, to solve these problems. It works in the following way.

Choose a random spin as the initial cluster.
If a neighboring spin is parallel to the initial spin it will be added to the cluster with probability $1 - e^{- 2 β J}$ .
Repeat this step for all points newly added to the cluster and repeat this procedure until no new points can be added.
Perform measurements using improved estimators.
Flip all spins in the cluster.

For the program, the only difference between the wolff algorithm and the single spin flip algorithm is the subroutine update. The subroutine update is modified as follows.

!--------------------------------------------------!
Subroutine update()
!--------------------------------------------------!
use configuration
implicit none
real(8),external :: ran
integer :: i,j,k,ii,jj
integer :: i1,i2,j3,j4
integer :: stack(L*L,2),visited(L,L)
integer :: wi,wj,top
real(8) :: p

p=1-exp(-2*beta)

wi = min(int(ran()*L)+1,L)
wj = min(int(ran()*L)+1,L)

stack = 0
visited = 1

top = 1
stack(top,1) = wi
stack(top,2) = wj

visited(wi,wj) = -1

do while(top>0) 

    i = stack(top,1)
    j = stack(top,2)
    top = top-1

    i1=modulo(i,L)+1
    if ((conf(i1,j) == conf(i,j)) .and. (visited(i1,j) == 1) .and. ran()<p) then
        visited(i1,j) = -1
        top = top+1
        stack(top,1) = i1
        stack(top,2) = j
    endif

    i2=modulo(i-2,L)+1
    if ((conf(i2,j) == conf(i,j)) .and. (visited(i2,j) == 1) .and. ran()<p) then
        visited(i2,j) = -1
        top = top+1
        stack(top,1) = i2
        stack(top,2) = j
    endif

    j3=modulo(j,L)+1
    if ((conf(i,j3) == conf(i,j)) .and. (visited(i,j3) == 1) .and. ran()<p) then
        visited(i,j3) = -1
        top = top+1
        stack(top,1) = i
        stack(top,2) = j3
    endif

    j4=modulo(j-2,L)+1
    if ((conf(i,j4) == conf(i,j)) .and. (visited(i,j4) == 1) .and. ran()<p) then
        visited(i,j4) = -1
        top = top+1
        stack(top,1) = i
        stack(top,2) = j4
    endif

enddo

do i = 1,L; do j = 1,L
    if (visited(i,j) == -1) conf(i,j) = -conf(i,j)
enddo; enddo

end Subroutine update