Cooking an Empanada with the Ising model. The Markov Chain Monte Carlo Recipe.

The Ising model is another example of a fascinating discrete model that can do a great job at modeling and simulating continuous systems. Among the many uses of the Ising model, it’s mainly known for simulating and studying ferromagnetic phase transitions which are notable changes in the properties of specific materials.

The most common shape of the model is a simple 2D grid that is divided into LxL cells. Each of these cells can occupy a binary number (+1 or -1). In the statistical mechanics' lingo, these correspond to magnetic dipole moments of atomic “spins” that can be either spin up or spin down (+1/-1). The energy of a unique configuration in the model is given by the following Hamiltonian function:

In that equation, i and j refer to coordinates in the lattice, and “s” refers to spins. So s_i is the spin at site i (can be +1 or -1). The sum over “<ij>” is shorthand for picking every pair of i and j that correspond to neighboring spins. The parameter J represents the strength of the intermolecular ferromagnetic interaction causing the spins to align, and the parameter h represents the effect of an external magnetic field that causes the spins to prefer aligning in one direction or the other. From that equation, we can conclude that any two neighboring spins that are aligned(having the same sign) decrease the energy, and any two neighboring spins that are anti-aligned (having opposite signs) increase the energy.

Since this is a model of ferromagnetism, the system is significantly affected by the externally applied temperature. At low temperature, the energy of the system prefers for neighboring spins to align, and below a critical temperature Tc​, the system will spontaneously acquire a net magnetization in one direction or the other. Watch this cute video to see the effect in practice. The reason that is happening is that at high temperature, the system is supplied with high energy and that makes all the states equally likely which means that the spins will not be aligned and the magnet will not work, while for low temperatures, the states with the low energy will be the most likely states of the system which means the system will favor the alignment of the spins in one direction which leads to the high performance of the magnet.

To model this interesting phase transition phenomenon, we can think of our system as a canonical ensemble in which the system is allowed to exchange heat with the external environment. This means that the equilibrium states of the systems should be statistically distributed according to the Boltzmann distribution:

Figure 2. PDF of the Boltzmann distribution. εi is the energy of a specific configuration, T is the temperature and k is a constant.

Then if we use Monte Carlo (MCMC) to get sample configurations for the Ising model from this distribution, then the phase transition should be simulated.

Now, I’m going to start this simulation in Python. I will use a random proposal distribution for MCMC and will define a class for the model:

# random proposal distribution for MCMC:
def proposal():
    return np.random.randint(sim.size, size=2)

J = 6.34369e-21  # Interaction constant for iron [Joule]
kB = 1.38065e-23  # Boltzmann constant [Joule / Kelvin]
class my_ising:
    def __init__(self,  T, size, magnetization=0  ,J= 6.34369e-21):
        self.J = J #Interaction Constant
        self.size = size # the size of the latttice
        self.T = T # the temperature
        self.config = np.zeros((self.size,self.size))
        self.steps=0
        self.magnetization=magnetization 
        for x in range(self.size): #starting with a radnom state
            for y in range(self.size):
                if self.magnetization == 0:
                    self.config[x, y] = 1 if random.random() < 0.5 else -1
               
    A function to visualize the model             
    def observe(self):
        pylab.cla()
        pylab.imshow(self.config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
        plt.title(f"after {self.steps} steps")
        plt.show()
        
    # A function to update based on the energy equation 
    def update(self, proposal):
        self.steps+=1
        x, y = proposal()
        energy = (
                    self.J * self.config[x, y]  *2 *(
                        self.config[(x + 1) % self.size, y] +
                        self.config[(x - 1) % self.size, y] +
                        self.config[x, (y + 1) % self.size] +
                        self.config[x, (y - 1) % self.size]))
if random.uniform(0, 1) < min(1,np.exp(-energy/ (T*kB) )):
            self.config[x, y] *= -1

In the code above, I started by defining a random proposal distribution. This would just randomly choose a cell in the lattice. I then wrote a class for the Ising model as a 2D lattice that is initialized randomly (cells have a random value of either +1 or -1). The update rule of that model works by taking a sample from the proposal distribution, which is just a coordinate of a cell, then it calculates the energy of the neighborhood that contains that cell using the hamiltonian equation that I described in the sections above. Once the energy is calculated, the Boltzmann distribution is used to decide whether the spin of that cell will change or not. According to our discussion above, this means that low energy states (aligned spins) are more likely to be sampled at low temperatures and random states will be sampled at high temperatures.

Let’s run this simulation:

Starting with a low temperature of 300 Kelvins:

states = []
T = 300 #kelvins
sim = my_ising(T = 300, size = 50)
pylab.cla()
pylab.imshow(sim.config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
plt.title(f"after {sim.steps} steps")
plt.show()
for i in range(50):
    for j in range(100000):
        sim.update(proposal)
    pylab.cla()
    pylab.imshow(sim.config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
    plt.title(f"after {sim.steps} steps")
    plt.show()
    states.append(sim.config)

We can see that at the low temperature, the configuration that led to the minimum energy was chosen. In this case, it was the one when all cells are aligned to +1.

Using a High temperature of 3000 Kelvins:

T = 3000
sim = my_ising(T = 300, size = 60)
pylab.cla()
pylab.imshow(sim.config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
plt.title(f"after {sim.steps} steps")
plt.show()
for i in range(100):
    for j in range(10000):
        sim.update(proposal)
    pylab.cla()
    pylab.imshow(sim.config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
    plt.title(f"after {sim.steps} steps")
    plt.show()

We see that at that high temperature, all the states are equally likely which means that the spins will be random with no alignment. This is the case when the magnet in the video didn’t work.

Finally, we can come to the fun part. I’m currently situated in Buenos Aires and everyone here loves to eat Empanadas. So I decided to make one with my Ising model. To make that Empanada, I will just change the Monte Carlo random proposal distribution to my new Empanada probability distribution:

def propose_empanda():
    x = random.randint(0, 60)
    y_lim = math.ceil(empanda(x))
    subt = math.ceil(0.2*y_lim)
    y = random.randint(0, y_lim-subt)
    return [60-x+5, 60-y-10 ]
def empanda(x, r= 30):
    return (np.sqrt(r**2-(x-30)**2))

Then I will run the Ising model simulation again with the high temperature 3000 Kelvins (certainly enough to cook an Empanada :D) but this time with the new fresh Empanada distribution:

T = 3000
sim = my_ising(T = 300, size = 70, magnetization = 1 )
pylab.cla()
pylab.imshow(sim.config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
plt.title(f"after {sim.steps} steps")
plt.show()
for i in range(10000):
    for j in range(100):
        sim.update(propose_empanda)
    pylab.cla()
    pylab.imshow(sim.config, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
    plt.title(f"after {sim.steps} steps")
    plt.show()

And now we can enjoy seeing our Empanada getting cooked :D

Hope you enjoyed this post. You can find the full code here.