Ising model simulation

Posted on February 5, 2018 by Dimitri Lozeve

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.


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:

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


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.

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.

Given a site \(i\), we reverse its spin to generate a new configuration state.

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.

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.

Finally, we put everything together in the update-state function, which will decide whether to accept or reject the new configuration.

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

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


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!