Try   HackMD

Ferromagnetic Ising model on a square lattice

No external field

Ising model on square lattice are

L×L spin that can take 2 values, either
Si=+1
or
Si=1
. Each spin interact with neighbors and with external field. The Hamiltonian of the system is
H=<ij>JijSiSjiHiSi

A particular solution for no external field are obtained by Lars Onsager where the second-order phase transition is happened at
JkBTc=ln(1+2)2
. With this in mind, we can use Metropolis algorithm to simulate Ising model with good accuracy. Following four figures show different MC steps, the spin of a
100×100
system with blue represent spin up and red represent spin down. Since the temperature happened around the critical temperature, the magnetization taken very long time to come to equilibrium.
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

First of all we can measure the mean energy at different temperature after reaching equilibrium by using following formula.

<E>=1N<<ij>JijSiSjiHiSi>
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

After equilibrium we can measure some physical quantity by averaging over time. For calculating the magnetization we use following formula.

m=1NiSi
No surprise to see that the magnetization become zero when temperature are greater than
Tc
. The simulation result are agreed with Onsager solution of magnetization.
m(T)=(1[sinh2βJ]4)1/8,T<Tc

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Heat capacity and susceptibility are best represent a phase transition at

Tc, since it increase around the critical temperature. Important to note that both heat capacity needs to taken time average since it measure the fluctuation of magnetization and energy.
Cv=kbβ2[(EE)2]

and
χ=βB[(mm)2]

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Then I plot the distribution of mean magnetization at different temperature with different size of system, the larger the system the smother and accurate the distribution be. The magnetization are agreed with the mean magnetization we fond earlier.

T=1 [J/kB]

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

T=2 [J/kB]
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

T=3 [J/kB]

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Finally, I plot the energy distribution which are indeed Boltzmann distribution and similarly the larger the system is the smoother and accurate it assemble to Boltzmann distribution which tells us that the greater the temperature the wider the spread and higher the mean.

P(E)=eβEEeβE
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

With external field

Then for the external field aren't zero, we can observe the first-order transition for

βJ>12ln(1+2) when H changes from negative to positive. Mean energy show a transition at
H=0
and has the minimum when
H=±1

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Similarly mean magnetization show first-order transition at

H=0 and the magnetization are aligned with the external field.

Heat capacity also has a transition at

H=0, and compare to no external field, the greater the external field the smaller the heat capacity.

However susceptibility are smaller compare to no external field.

The magnetization probability are totally expected that was aligned with the external field.

H=1

Unknown environment 'figure'

Unknown environment 'center'
Unknown environment 'figure'

Unknown environment 'figure'

Energy distribution didn't show Boltzmann distribution like the one
without external field, it has uniform distribution of energy
value with wider spread the distribution when the external field
is greater.

Unknown environment 'figure'
Unknown environment 'figure'

​​​​\subsection*{Ising.f95}

\begin{minted}[frame=lines, linenos, fontsize=\large]{fortran}
program Ising2d
implicit none
double precision, external :: mag, chi, E, dE, C
! spin up = +1, spin down = -1
integer, parameter :: L=64
double precision :: ss(L, L), ii, jj, z, beta, H, kT, k, T
integer :: S(L, L), i, j, n, nmax, sample, m, point, o, p
double precision :: start, finish
double precision :: E1, E2, M1, M2, avgM, avgE

​​​​call random_seed()
​​​​! open(10, file='data2/M_8.dat', status='unknown')
​​​​! open(20, file='data2/S_8.dat', status='unknown')
​​​​! open(30, file='data2/T2/Pm8.dat', status='unknown')
​​​​! open(40, file='data2/E_8.dat', status='unknown')
​​​​! open(50, file='data2/C_8.dat', status='unknown')
​​​​open(60, file='data2/PE/5.dat', status='unknown')
​​​​! open(70, file='data2/config.dat', status='unknown')
​​​​! parameter
​​​​! nmax = 25000000
​​​​nmax = 500000
​​​​k = 1.d0
​​​​H = 0.d0
​​​​! sample = 2000000
​​​​sample = 200000
​​​​start = -1.d0
​​​​finish = 1.d0
​​​​point = 30

​​​​! do m = 0, point
​​​​    ! T = start + (finish - start) / dble(point) * dble(m)
​​​​    ! H = start + (finish - start) / dble(point) * dble(m)
​​​​    H = 1.0d0
​​​​    T = 2.d0
​​​​    write(6, '(A, f5.2)') 'T=', T
​​​​    write(6, '(A, f5.2)') 'H=', H

​​​​    kT = k * T
​​​​    beta = 1.d0 / kT
​​​​    E1 = 0.d0
​​​​    E2 = 0.d0
​​​​    M1 = 0.d0
​​​​    M2 = 0.d0
​​​​    avgM = 0.d0
​​​​    avgE = 0.d0

​​​​    !!!!! step1: random start config !!!!!
​​​​    call random_number(ss)
​​​​    ! S = nint(ss)
​​​​    ! ground state start (cold start)
​​​​    S = 1

​​​​    do n = 1, nmax
​​​​        ! write config
​​​​        if (n .eq. nmax) then
​​​​            do o = 1, L
​​​​                do p = 1, L
​​​​                    write(70, *) S(o, p)
​​​​                end do
​​​​            end do
​​​​        end if

​​​​        !!!!! step2: randomly choose (i, j) one and try to flip !!!!!
​​​​        call random_number(ii)
​​​​        call random_number(jj)
​​​​        ii = ii * dble(L-1)
​​​​        jj = jj * dble(L-1)
​​​​        i = nint(ii) + 1
​​​​        j = nint(jj) + 1

​​​​        !!!!! step3: determine wether to flip !!!!!
​​​​        if (dE(S, i, j, L, H) .le. 0) then
​​​​            S(i, j) = - 1 * S(i, j)
​​​​        else
​​​​            call random_number(z) 
​​​​            if (dexp(- beta * dE(S, i, j, L, H)) .gt. z) then
​​​​                S(i, j) = - 1 * S(i, j)
​​​​            end if
​​​​        end if

​​​​        !!!!! step4: measure quantity !!!!!
​​​​        if (n .ge. (nmax-sample)) then
​​​​            avgM = avgM + mag(S, L)
​​​​            avgE = avgE + E(S, H, L)
​​​​            E1 = E1 + E(S, H, L) * dble(L) * dble(L)
​​​​            E2 = E2 + E(S, H, L) * E(S, H, L) * dble(L) * dble(L) * dble(L) * dble(L)
​​​​            M1 = M1 + mag(S, L)
​​​​            M2 = M2 + mag(S, L) * mag(S, L)
​​​​            ! P(m)
​​​​            ! write(30, *) mag(S, L)
​​​​            ! P(E)
​​​​            write(60, *) E(S, H, L)
​​​​        end if
​​​​    end do
​​​​    ! specific heat
​​​​    E1 = E1 / sample
​​​​    E2 = E2 / sample
​​​​    ! susceptibility
​​​​    M1 = M1 / sample
​​​​    M2 = M2 / sample
​​​​    ! mean magnetization 
​​​​    avgM = avgM / sample
​​​​    ! mean energy 
​​​​    avgE = avgE / sample

​​​​    ! write: <M>, <X>, <E>, <Cv>
​​​​    ! write(10, *) H, avgM
​​​​    ! write(20, *) H, beta * dble(L) * dble(L) * (M2 - M1 * M1)
​​​​    ! write(40, *) H, avgE
​​​​    ! write(50, *) H, k * beta * beta * (E2 - E1 * E1) / dble(L) / dble(L)
​​​​! end do

end program Ising2d

! calculate magnetization per spin
function mag(S, L)
implicit none
double precision :: mag
integer :: L, S(L, L), i, j
mag = 0.d0
do i = 1, L
do j = 1, L
mag = mag + dble(S(i, j))
end do
end do
mag = mag / dble(L) / dble(L)
end function

! calculate susceptibility
function chi(S, L, beta)
implicit none
double precision, external :: mag
double precision :: chi, beta
integer :: L, S(L, L), i, j
chi = 0.d0
do i = 1, L
do j = 1, L
chi = chi + dble(S(i, j)) * dble(S(i, j))
end do
end do
chi = chi / dble(L) / dble(L)
chi = chi - mag(S, L) * mag(S, L)
chi = beta * dble(L) * dble(L) * chi
end function chi

! calculate Hamiltonian
function E(S, H, L)
implicit none
double precision :: E, H
double precision :: J=1.d0
integer :: L, S(L, L), i, k
E = 0.d0
do i = 1, L-1
do k = 1, L-1
E = E - J * S(i, k) * S(i, k+1)
E = E - J * S(i, k) * S(i+1, k)
E = E - H * S(i, k)
end do
end do
do i = 1, L
E = E - J * S(i, L) * S(i, 1)
E = E - H * S(i, k)
end do
do k = 1, L
E = E - J * S(L, k) * S(1, k)
E = E - H * S(i, k)
end do
E = E / dble(L) / dble(L)
end function E

! calculate delta E (only five possible value)
function dE(S, i, k, L, H)
implicit none
double precision, external :: E
integer :: L, S(L, L), i, k
double precision :: J=1.d0, dE, H
dE = 0.d0
if (i.ne.1) then
dE = dE + S(i-1, k)
dE = dE + S(i+1, k)
else
dE = dE + S(L, k)
dE = dE + S(2, k)
end if
if (k.ne.1) then
dE = dE + S(i, k-1)
dE = dE + S(i, k+1)
else
dE = dE + S(i, L)
dE = dE + S(i, 2)
end if
dE = dE * 2.d0 * J * S(i, k) + H * S(i, k) * 2.d0
end function
\end{minted}