The Ising model is a model used to represent magnetic dipole moments in statistical physics. Physical details are on the Wikipedia page, but what is interesting is that it follows a complex probability distribution on a lattice, where each site can take the value +1 or -1.

# Mathematical definition

We have a lattice $\Lambda$ consisting of sites $k$. For each site,
there is a moment $\sigma_k \in \{ -1, +1 \}$. $\sigma =
(\sigma_k)_{k\in\Lambda}$ is called the *configuration* of the
lattice.

The total energy of the configuration is given by the *Hamiltonian*
\[
H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j,
\]
where $i\sim j$ denotes *neighbours*, and $J$ is the
*interaction matrix*.

The *configuration probability* is given by:
\[
\pi_\beta(\sigma) = \frac{e^{-\beta H(\sigma)}}{Z_\beta}
\]
where $\beta = (k_B T)^{-1}$ is the inverse temperature,
and $Z_\beta$ the normalisation constant.

For our simulation, we will use a constant interaction term $J > 0$. If $\sigma_i = \sigma_j$, the probability will be proportional to $\exp(\beta J)$, otherwise it would be $\exp(\beta J)$. Thus, adjacent spins will try to align themselves.

# Simulation

The Ising model is generally simulated using Markov Chain Monte Carlo (MCMC), with the Metropolis-Hastings algorithm.

The algorithm starts from a random configuration and runs as follows:

- Select a site $i$ at random and reverse its spin: $\sigma'_i = -\sigma_i$
- Compute the variation in energy (hamiltonian) $\Delta E = H(\sigma') - H(\sigma)$
- If the energy is lower, accept the new configuration
- Otherwise, draw a uniform random number $u \in ]0,1[$ and accept the new configuration if $u < \min(1, e^{-\beta \Delta E})$.

# Implementation

The simulation is in Clojure, using the Quil library (a Processing library for Clojure) to display the state of the system.

This post is "literate Clojure", and contains
`core.clj`

. The
complete project can be found on
GitHub.

```
(ns ising-model.core
(:require [quil.core :as q]
[quil.middleware :as m]))
```

The application works with Quil's functional mode, with each function taking a state and returning an updated state at each time step.

The `setup`

function generates the initial state, with random initial
spins. It also sets the frame rate. The matrix is a single vector in
row-major mode. The state also holds relevant parameters for the
simulation: $\beta$, $J$, and the iteration step.

```
(defn setup [size]
"Setup the display parameters and the initial state"
(q/frame-rate 300)
(q/color-mode :hsb)
(let [matrix (vec (repeatedly (* size size) #(- (* 2 (rand-int 2)) 1)))]
{:grid-size size
:matrix matrix
:beta 10
:intensity 10
:iteration 0}))
```

Given a site $i$, we reverse its spin to generate a new configuration state.

```
(defn toggle-state [state i]
"Compute the new state when we toggle a cell's value"
(let [matrix (:matrix state)]
(assoc state :matrix (assoc matrix i (* -1 (matrix i))))))
```

In order to decide whether to accept this new state, we compute the difference in energy introduced by reversing site $i$: \[ \Delta E = J\sigma_i \sum_{j\sim i} \sigma_j. \]

The `filter some?`

is required to eliminate sites outside of the
boundaries of the lattice.

```
(defn get-neighbours [state idx]
"Return the values of a cell's neighbours"
[(get (:matrix state) (- idx (:grid-size state)))
(get (:matrix state) (dec idx))
(get (:matrix state) (inc idx))
(get (:matrix state) (+ (:grid-size state) idx))])
(defn delta-e [state i]
"Compute the energy difference introduced by a particular cell"
(* (:intensity state) ((:matrix state) i)
(reduce + (filter some? (get-neighbours state i)))))
```

We also add a function to compute directly the hamiltonian for the entire configuration state. We can use it later to log its values across iterations.

```
(defn hamiltonian [state]
"Compute the Hamiltonian of a configuration state"
(- (reduce + (for [i (range (count (:matrix state)))
j (filter some? (get-neighbours state i))]
(* (:intensity state) ((:matrix state) i) j)))))
```

Finally, we put everything together in the `update-state`

function,
which will decide whether to accept or reject the new configuration.

```
(defn update-state [state]
"Accept or reject a new state based on energy
difference (Metropolis-Hastings)"
(let [i (rand-int (count (:matrix state)))
new-state (toggle-state state i)
alpha (q/exp (- (* (:beta state) (delta-e state i))))]
;;(println (hamiltonian new-state))
(update (if (< (rand) alpha) new-state state)
:iteration inc)))
```

The last thing to do is to draw the new configuration:

```
(defn draw-state [state]
"Draw a configuration state as a grid"
(q/background 255)
(let [cell-size (quot (q/width) (:grid-size state))]
(doseq [[i v] (map-indexed vector (:matrix state))]
(let [x (* cell-size (rem i (:grid-size state)))
y (* cell-size (quot i (:grid-size state)))]
(q/no-stroke)
(q/fill
(if (= 1 v) 0 255))
(q/rect x y cell-size cell-size))))
;;(when (zero? (mod (:iteration state) 50)) (q/save-frame "img/ising-######.jpg"))
)
```

And to reset the simulation when the user clicks anywhere on the screen:

```
(defn mouse-clicked [state event]
"When the mouse is clicked, reset the configuration to a random one"
(setup 100))
```

```
(q/defsketch ising-model
:title "Ising model"
:size [300 300]
:setup #(setup 100)
:update update-state
:draw draw-state
:mouse-clicked mouse-clicked
:features [:keep-on-top :no-bind-output]
:middleware [m/fun-mode])
```

# Conclusion

The Ising model is a really easy (and common) example use of MCMC and Metropolis-Hastings. It allows to easily and intuitively understand how the algorithm works, and to make nice visualizations!