\documentclass[../ising_model.tex]{subfiles} 

\begin{document} 
\section{Methods}\label{sec:methods}
\subsection{The Ising model}\label{subsec:ising_model}
The Ising model consist of a lattice of spins, which can be though of as atoms 
in a grid. In two dimensions, the length of the lattice is given by $L$, and the 
number of spins within a lattice is given by $N = L^{2}$. When we consider the 
entire lattice, the system spin configuration is represented as a matrix $L \times L$
\begin{align*}
    \mathbf{s} &= \begin{pmatrix}
        s_{1,1} & s_{1,2} & \dots & s_{1,L} \\
        s_{2,1} & s_{2,2} & \dots & s_{2,L} \\
        \vdots & \vdots & \ddots & \vdots \\
        s_{L,1} & s_{L,2} & \dots & s_{L,L}
    \end{pmatrix} \ .
\end{align*} % $\mathbf{s} = [s_{1}, s_{2}, \dots, s_{N}]$.
The total number of possible spin configurations, also called microstates, is $2^{N}$. 

A given spin $i$ can take one of two possible discrete values $s_{i} \in \{-1, +1\}$, 
where $-1$ represent down and $+1$ spin up. The spins interact with its nearest neighbors, 
and in a two-dimensional lattice each spin has up to four nearest neighbors. However, 
the model is not restricted to this dimentionality \cite[p. 3]{obermeyer:2020:ising}. 
In our experiment we will use periodic boundary conditions, meaning all spins 
have exactly four nearest neighbors. To find the analytical expressions necessary 
for validating our model implementation, we will assume a $2 \times 2$ lattice. 

The hamiltonian of the Ising model is given by
\begin{equation}
    E(\mathbf{s}) = -J \sum_{\langle k l \rangle}^{N} s_{k} s_{l} - B \sum_{k}^{N} s_{k}\ ,
    \label{eq:energy_hamiltonian}
\end{equation}
where $\langle k l \rangle$ denote a spin pair. $J$ is the coupling constant, and 
$B$ is the external magnetic field. For simplicity, we will consider the Ising model 
where $B = 0$, and find the total system energy given by 
\begin{equation}
    E(\mathbf{s}) = -J \sum_{\langle k l \rangle}^{N} s_{k} s_{l} \ .
    \label{eq:energy_total}
\end{equation}
To avoid counting duplicated, we count the neighboring spins using the pattern 
visualized in Figure \ref{fig:tikz_boundary}.
\begin{figure}
    \centering
    \begin{tikzpicture}
        \draw[help lines] (0, 0) grid (2, 2);
        % \node[inner] (s1) at (0.5, 0.5) {s};
        \node (s1) at (0.5, 1.5) {$s_1$};
        \node (s2) at (1.5, 1.5) {$s_2$};
        \node (s3) at (0.5, 0.5) {$s_3$};
        \node (s4) at (1.5, 0.5) {$s_4$};
        \node[gray] (s12) at (2.5, 1.5) {$s_1$};
        \node[gray] (s34) at (2.5, 0.5) {$s_3$};
        \node[gray] (s13) at (0.5, -0.5) {$s_1$};
        \node[gray] (s24) at (1.5, -0.5) {$s_2$};
    
        \draw[red, ->] (s1.east) -- (s2.west);
        \draw[red, ->] (s2.east) -- (s12.west);
        \draw[red, ->] (s1.south) -- (s3.north);
        \draw[red, ->] (s2.south) -- (s4.north);
        \draw[red, ->] (s3.east) -- (s4.west);
        \draw[red, ->] (s4.east) -- (s34.west);
        \draw[red, ->] (s3.south) -- (s13.north);
        \draw[red, ->] (s4.south) -- (s24.north);
      \end{tikzpicture}
      \caption{Visualization of counting rules for each pair of spins in a 
      $2 \times 2$ lattice, where periodic boundary conditions are applied.}
      \label{fig:tikz_boundary}
\end{figure} 
We also find the total magnetization of the system, which is given by
\begin{equation}
    M(\mathbf{s}) = \sum_{i}^{N} s_{i} \ .
    \label{eq:magnetization_total}
\end{equation}
In addition, we have to consider the state degeneracy, the number of different 
microstates sharing the same value of total magnetization. In the case where we 
have two spins oriented up, the total energy have two possible values, as shown 
in Appendix \ref{sec:energy_special}.
\begin{table}[H]
    \centering
    \begin{tabular}{cccc} 
        \hline
        Spins up & $E(\mathbf{s})$ & $M(\mathbf{s})$ & Degeneracy \\
        \hline 
        4 & -8J & 4 & 1 \\
        3 & 0 & 2 & 4 \\
        \multirow{2}*{2} & 0 & 0 & 4 \\
            & 8J & 0 & 2 \\ 
        1 & 0 & -2 & 4 \\
        0 & -8J & -4 & 1 \\
        \hline
    \end{tabular}
    \caption{Values of total energy, total magnetization, and degeneracy for the 
    possible states of a system for a $2 \times 2$ lattice, with periodic boundary 
    conditions.}
    \label{tab:lattice_config}
\end{table}
We calculate the analytical values for a $2 \times 2$ lattice, found in Table \ref{tab:lattice_config}.
However, we scale the total energy and total magnetization by number of spins, 
to compare these quantities for lattices where $L > 2$. Energy per spin is given by
\begin{equation}
    \epsilon (\mathbf{s}) = \frac{E(\mathbf{s})}{N} \ ,
    \label{eq:energy_spin}
\end{equation}
and magnetization per spin given by 
\begin{equation}
    m (\mathbf{s}) = \frac{M(\mathbf{s})}{N} \ .
    \label{eq:magnetization_spin}
\end{equation}


\subsection{Statistical mechanics}\label{subsec:statistical_mechanics}
When we study the behavior of the Ising model, we have to consider the probability 
of a microstate $\mathbf{s}$ at a fixed temperature $T$. The probability distribution 
function (pdf) is given by 
\begin{equation}
    p(\mathbf{s} \ | \ T) = \frac{1}{Z} \exp^{-\beta E(\mathbf{s})},
    \label{eq:boltzmann_distribution}
\end{equation}
known as the Boltzmann distribution. This is an exponential distribution, where 
$\beta$ is given by 
\begin{equation}
    \beta = \frac{1}{k_{B}} \ ,
    \label{eq:beta}
\end{equation}
where and $k_{B}$ is the Boltzmann constant. $Z$ is a normalizing factor of the 
pdf, given by
\begin{equation}
    Z = \sum_{\text{all possible } \mathbf{s}} e^{-\beta E(\mathbf{s})} \ , 
    \label{eq:partition}
\end{equation}
and is known as the partition function. We derive $Z$ in Appendix \ref{sec:partition_function},
which gives us
\begin{equation*}
    Z = 4 \cosh (8 \beta J) + 12 \ .
\end{equation*}
Using the partition function and Equation \eqref{eq:boltzmann_distribution}, the 
probability of a microstate at a fixed temperature is given by 
\begin{equation}
    p(\mathbf{s} \ | \ T) = \frac{1}{4 \cosh (8 \beta J) + 12} e^{-\beta E(\mathbf{s})} \ .
    \label{eq:pdf}
\end{equation}
% Add something about why we use the expectation values?
We find the expected total energy
\begin{equation}\label{eq:energy_total_result}
	\langle E \rangle = -\frac{8J \sinh(8 \beta J)}{\cosh(8 \beta J) + 3} \ ,
\end{equation}
and the expected energy per spin
\begin{equation}\label{eq:energy_spin_result}
	\langle \epsilon \rangle = \frac{-2J \sinh(8 \beta J)}{ \cosh(8 \beta J) + 3} \ .
\end{equation}
We find the expected absolute total magnetization
\begin{equation}\label{eq:magnetization_total_result}
	\langle |M| \rangle = \frac{2(e^{8 \beta J} + 2)}{\cosh(8 \beta J) + 3} \ ,
\end{equation}
and the expected magnetization per spin
\begin{equation}\label{eq:magnetization_spin_result}
	\langle |m| \rangle = \frac{e^{8 \beta J} + 1}{2( \cosh(8 \beta J) + 3)} \ .
\end{equation}
We derive the analytical expressions for expectation values in Appendix 
\ref{sec:expectation_values}. 

We also need to determine the heat capacity 
\begin{equation}
    C_{V} = \frac{1}{k_{B} T^{2}} (\mathbb{E}(E^{2}) - [\mathbb{E}(E)]^{2}) \ ,
    \label{eq:heat_capacity}
\end{equation}
and the magnetic susceptibility 
\begin{align*}
    \chi = \frac{1}{k_{\text{B} T}} (\mathbb{E}(M^{2}) - [\mathbb{E}(M)]^{2}) \ .
    \label{eq:susceptibility}
\end{align*}
In Appendix \ref{sec:heat_susceptibility} we derive expressions for the heat 
capacity and the susceptibility, and find the heat capacity 
\begin{align*}
    \frac{C_{V}}{N} &= \frac{64J^{2} }{N k_{B} T^{2}} \bigg( \frac{3\cosh(8 \beta J) + 1}{(\cosh(8 \beta J) + 3)^{2}} \bigg) \ ,
\end{align*}
and magnetic susceptibility
\begin{align*}
    \frac{\chi}{N} &= \frac{4}{N k_{B} T} \bigg( \frac{3e^{8 \beta J} + e^{-8 \beta J} + 3}{(\cosh(8 \beta J) + 3)^{2}} \bigg) \ .
\end{align*}

For scalability of our model, we want to work with unitless spins. To obtain this, 
we introduce the coupling constant $J$ as the base unit for energy. With the 
Boltzmann constant we derive the remaining units, which can be found in Table 
\ref{tab:units}.
\begin{table}[H]
    \centering
    \begin{tabular}{ccc} % @{\extracolsep{\fill}}
        \hline
        Value & Unit &  \\
        \hline 
        $[E]$ & $J$ &  \\
        $[M]$ & $1$ & (unitless) \\
        $[T]$ & $J / k_{B}$ &  \\
        $[C_{V}]$ & $k_{B}$ &  \\
        $[\chi]$ & $1 / J$ &  \\
        \hline
    \end{tabular}
    \caption{Values of total energy, total magnetization, and degeneracy for the possible states of a system for a $2 \times 2$ lattice, with periodic boundary conditions.}
    \label{tab:units}
\end{table}


\subsection{Phase transition and critical temperature}\label{subsec:phase_critical}
We consider the Ising model in two dimensions, with no external external magnetic 
field. At temperatures below the critical temperature $T_{c}$, the Ising model will 
magnetize spontaneous. When increasing the temperature of the external field, 
the Ising model transition from an ordered to an unordered phase. The spins become 
more correlated, and we can measure the discontinous behavior as an increase in 
correlation length $\xi (T)$ \cite[p. 432]{hj:2015:comp_phys}. At $T_{c}$, the 
correlation length is proportional with the lattice size, resulting in the critical 
temperature scaling relation
\begin{equation}
    T_{c}(L) - T_{c}(L = \infty) = aL^{-1} \ ,
    \label{eq:critical_infinite}
\end{equation}
where $a$ is a constant. For the Ising model in two dimensions, with an lattice of 
infinite size, the analytical solution of the critical temperature is 
\begin{equation}
    T_{c}(L = \infty) = \frac{2}{\ln (1 + \sqrt{2})} J/k_{B} \approx 2.269 J/k_{B}
    \label{eq:critical_solution}
\end{equation}
This result was found analytically by Lars Onsager in 1944. We can estimate the 
critical temperature of an infinite lattice, using finite lattices critical temperature 
and linear regression. 


\subsection{The Markov chain Monte Carlo method}\label{subsec:mcmc_method}
Markov chains consist of a sequence of samples, where the probability of the next 
sample depend on the probability of the current sample. Whereas the Monte Carlo 
method introduces randomness to the sampling, which allows us to approximate 
statistical quantities \cite{tds:2021:mcmc} without using the pdf directly. We 
generate samples of spin configuration, to approximate $\langle \epsilon \rangle$, 
$\langle |m| \rangle$, $C_{V}$, and $\chi$, using the Markov chain Monte Carlo 
(MCMC) method. 

The Ising model is initialized in a random state, and the Markov chain is evolved 
until the model reaches an equilibrium state. However, generating new random states 
require ergodicity and detailed balance. A Markov chain is ergoditic when all system 
states can be reached at every current state, whereas detaild balance implies no 
net flux of probability. To satisfy these criterias we use the Metropolis-Hastings 
algorithm, found in Figure \ref{algo:metropolis}, to generate samples of microstates. 

One Monte Carlo cycle consist of changing the initial configuration of the lattice,
by randomly flipping a spin. When a spin is flipped, the change in energy is evaluated as
\begin{align*}
    \Delta E &= E_{\text{after}} - E_{\text{before}} \ ,
\end{align*}
and accepted if $\Delta E \leq 0$. However, if $\Delta E > 0$ we have to compute
the Boltzmann factor
\begin{equation}
    p(\mathbf{s} | T) = e^{-\beta \Delta E}
    \label{eq:boltzmann_factor}
\end{equation} 
Since the total system energy only takes three different values, the change in 
energy can take $3^{2}$ values. However, there are only five distinct values 
$\Delta E = \{-16J, -8J, 0, 8J, 16J\}$, we derive these values in Appendix \ref{sec:delta_energy}.
We can avoid computing the Boltzmann factor, by using a look up table (LUT). 
We use $\Delta E$ as an index to access the resulting value of the exponential function, 
in an array. 
\begin{figure}[H]
    \begin{algorithm}[H]
    \caption{Metropolis-Hastings Algorithm}
    \label{algo:metropolis}
        \begin{algorithmic}
            \Procedure{Metropolis}{$lattice, energy, magnetization$}
            \For{ Each spin in the lattice }
            \State $ri \leftarrow \text{random index of the lattice}$
            \State $rj \leftarrow \text{random index of the lattice}$
            \State $dE \leftarrow 2 \cdot lattice_{ri,rj} \cdot (lattice_{ri,rj-1} + lattice_{ri,rj+1} 
            + lattice_{ri-1,rj} + lattice_{ri+1,rj} )$

            \State $rn \leftarrow \text{random number} \in [0,1)$
            \If{$rn \leq e^{- \frac{dE}{T}}$}
                \State $lattice_{ri,rj} = -1 \cdot lattice_{ri,rj}$
                \State $magnetization = magnetization + 2 \cdot lattice_{ri,rj}$
                \State $energy = energy + dE$
            \EndIf
            \EndFor
            \EndProcedure
        \end{algorithmic}
    \end{algorithm}
\end{figure}
The Markov process reach an equilibrium after a certain number of Monte Carlo cycles, 
where the system state is reflecting the state of a real system. After this burn-in time, 
given by number of Monte Carlo cycles, we can start sampling microstates.
The probability distribution of the samples will tend toward the actual probability 
distribution of the system. 


\subsection{Implementation and testing}\label{subsec:implementation_test}
We implemented a test suite, and compared the numerical estimates to the analytical 
results from Section \ref{subsec:statistical_mechanics}. In addition, we set a tolerance to 
verify convergence, and a maximum number of Monte Carlo cycles to avoid longer 
runtimes during implementation.

We used a pattern to access all neighboring spins, where all indices of the neighboring 
spins are put in an $L \times 2$ matrix. The indices are accessed using pre-defined 
constants, where the first column contain the indices for neighbors to the left 
and up, and the second column right and down. This method avoids the use of if-tests,
and takes advantage of the parallel optimization. 

We parallelize our code using both a message passing interface (OpenMPI) and 
multi-threading (OpenMP). First, we divide the temperatures into smaller ranges, 
and each MPI process receives a set of temperatures. Every MPI process spawn a 
set of threads, which initialize an Ising model and performs the Metropolis-Hastings 
algorithm. We limit the number of times threads are spawned and joined, by using 
single parallel regions, reducing parallel overhead. We used Fox \footnote{Technical 
specifications for Fox can be found at \url{https://www.uio.no/english/services/it/research/platforms/edu-research/help/fox/system-overview.md}}, 
a high-performance computing cluster, to run our program. 



\subsection{Tools}\label{subsec:tools}
The Ising model and MCMC methods are implemented in C++, and parallelized using 
both \verb|OpenMPI| \cite{gabriel:2004:open_mpi} and \verb|OpenMP| \cite{openmp:2018}. We used the Python library 
\verb|matplotlib| \cite{hunter:2007:matplotlib} to produce all the plots, and 
\verb|seaborn| \cite{waskom:2021:seaborn} to set the theme in the figures. In 
addition, we used \verb|Scalasca| \cite{scalasca} and \verb|Score-P| \cite{scorep} 
to profile and optimize our implementation.
% Add profiler

\end{document}