# Ferromagnetic Ising model on a square lattice

## No external field
Ising model on square lattice are $L \times L$ spin that can take 2 values, either $S_i=+1$ or $S_i=-1$. Each spin interact with neighbors and with external field. The Hamiltonian of the system is
\begin{equation}
H = -\sum_{<ij>} J_{ij}S_iS_j - \sum_i H_iS_i
\end{equation}
A particular solution for no external field are obtained by Lars Onsager where the second-order phase transition is happened at
$\dfrac{J}{k_BT_c} = \dfrac{ln(1+\sqrt{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 \times 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.

First of all we can measure the mean energy at different temperature after reaching equilibrium by using following formula.
\begin{equation}
<E>=-\dfrac{1}{N}<-\sum_{<ij>} J_{ij}S_iS_j - \sum_i H_iS_i>
\end{equation}

After equilibrium we can measure some physical quantity by averaging over time. For calculating the magnetization we use following formula.
\begin{equation}
\langle m \rangle = \dfrac{1}{N}\sum_{i}S_i
\end{equation}
No surprise to see that the magnetization become zero when temperature are greater than $T_c$. The simulation result are agreed with Onsager solution of magnetization.
\begin{equation}
m(T) = (1-[\sinh 2\beta J]^{-4})^{1/8} \quad, T \lt T_c
\end{equation}


Heat capacity and susceptibility are best represent a phase transition at $T_c$, 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.
\begin{equation}
C_v = k_b\beta^2[\langle (E-\langle E \rangle)^2 \rangle]
\end{equation}
and
\begin{equation}
\chi = \beta B[\langle (m-\langle m\rangle)^2 \rangle]
\end{equation}


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.

\begin{equation}
{\large T=1 \ [J/k_B]}
\end{equation}



\begin{equation}
{\large T=2 \ [J/k_B]}
\end{equation}


\begin{equation}
{\large T=3 \ [J/k_B]}
\end{equation}


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.
\begin{equation}
P(E)=\dfrac{e^{-\beta E}}{\sum_E e^{-\beta E}}
\end{equation}

## With external field

Then for the external field aren't zero, we can observe the first-order transition for $\beta J > \dfrac{1}{2}ln(1+\sqrt{2})$ when H changes from negative to positive. Mean energy show a transition at $H=0$ and has the minimum when $H=\pm 1$

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.
\begin{equation}
{\large H=-1}
\end{equation}

\begin{figure}[H]
\centering
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bT1Pm32} }}%
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bT1Pm64} }}%
\end{figure}

\begin{center}
{\large $H=1$}
\end{center}
\begin{figure}[H]
\centering
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bT2Pm8} }}%
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bT2Pm16} }}%
\end{figure}
\begin{figure}[H]
\centering
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bT2Pm32} }}%
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bT2Pm64} }}%
\end{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.
\begin{figure}[H]
\centering
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bPE8} }}%
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bPE16} }}%
\end{figure}
\begin{figure}[H]
\centering
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bPE32} }}%
\subfloat{{\includegraphics[width=0.5\textwidth]{../figures/bPE64} }}%
\end{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}