Lab 7: Ising Model#

For the class on Wednesday, February 5th

Code for simulating a 2D Ising Model#

The following cell defines the IsingModel2D class, which implements a numerical representation of a 2D Ising Model and the associated MCMC simulation methods.

You need to run this cell, but you don’t need to edit anything in it. You are welcome to take a look at the code of the class, but it is not necessary to understand the code in detail to complete the lab.

Hide code cell source
import numpy as np
import matplotlib.pyplot as plt

class IsingModel2D:
    def __init__(self, side_length: int, boundaries=None, initial_values="cold", antiferro=False):
        self._side_length = int(side_length)
        self._size = self._side_length * self._side_length

        indices = np.arange(self._size).reshape((self._side_length, self._side_length))
        neighbor_indices = np.stack([
            np.roll(indices, 1, 0),   # up
            np.roll(indices, -1, 1),  # right
            np.roll(indices, -1, 0),  # down
            np.roll(indices, 1, 1),   # left
        ])

        if boundaries not in [None, "periodic"]:
            try:
                bu, br, bd, bl = boundaries
            except TypeError:
                bu = br = bd = bl = boundaries
            idx = np.arange(self._size, self._size + 3)
            for bv, s in zip([bu, br, bd, bl], [(0,0), (1,Ellipsis,-1), (2,-1), (3,Ellipsis,0)]):
                neighbor_indices[s] = idx[int(np.sign(bv))]

        self._neighbor_indices = neighbor_indices.reshape(4, self._size).T

        self._data = np.empty(self._size + 3, dtype=np.int32)
        self._data[-3:] = [0, 1, -1]
        if isinstance(initial_values, str):
            if initial_values == "hot":
                self.data = np.random.default_rng().random(self._size) >= 0.5
            elif initial_values == "cold":
                self.data = 1
        else:
            self.data = initial_values

        self._J = float(-1 if antiferro else 1)

    @property
    def size(self):
        return self._size

    @property
    def data(self):
        return self._data[:-3].reshape(self._side_length, self._side_length).copy()

    @data.setter
    def data(self, value):
        value = np.ravel(value)
        if len(value) not in [1, self._size]:
            raise ValueError(f"`value` must be a single value or an array of size {self._size}")
        new_value = np.ones(len(value), dtype=np.int32)
        new_value[~(value > 0)] = -1
        if len(new_value) == 1:
            new_value = new_value[0]
        self._data[:-3] = new_value

    def show(self, **kwargs):
        plt.matshow(self.data, vmax=1, vmin=-1, **kwargs)

    def _check_site_index(self, site_index: int):
        if not isinstance(site_index, int):
            i, j = site_index
            site_index = int(i) * self._side_length + int(j)
        if site_index >= 0 and site_index < self._size:
            return site_index
        raise ValueError(f"`site_index` must be between 0 and {self._size - 1} (inclusive)")

    def calc_magnetization(self):
        return self._data[:self._size].sum(dtype=np.float64) / self._size

    def calc_total_energy(self, ext_field=0.0):
        energy = (
            self._data[self._neighbor_indices].sum(axis=1, dtype=np.float64)
            * self._data[:self._size]
        ).sum(dtype=np.float64) * self._J * (-0.5)
        if ext_field:
            energy -= self._data[:self._size].sum(dtype=np.float64) * float(ext_field)
        return energy

    def calc_local_energy(self, site_index: int, ext_field=0.0):
        site_index = self._check_site_index(site_index)
        energy = (
            self._data[self._neighbor_indices[site_index]].sum(dtype=np.float64)
            * self._data[site_index]
            * self._J
            * (-1.0)
        )
        if ext_field:
            energy -= self._data[site_index] * float(ext_field)
        return energy

    def flip(self, site_index):
        self._data[self._check_site_index(site_index)] *= -1

    def transition_prob_metropolis(self, site_index: int, beta: float, ext_field=0.0):
        return np.exp(2.0 * beta * self.calc_local_energy(site_index, ext_field))

    def transition_prob_gibbs(self, site_index: int, beta: float, ext_field=0.0):
        return 1.0 / (1.0 + np.exp(-2.0 * beta * self.calc_local_energy(site_index, ext_field)))

    def sweep(self, beta, ext_field=0.0, nsweeps=1, seed=None, transition_prob="metropolis"):
        rng = np.random.default_rng(seed)
        transition_prob_func = getattr(self, f"transition_prob_{transition_prob}")
        accept = 0
        for _  in range(nsweeps):
            for i in range(self._size):
                if rng.random() < transition_prob_func(i, beta, ext_field):
                    self.flip(i)
                    accept += 1
        return accept / self._size / nsweeps

A. Energy and Magnetization Calculations#

The following cell shows a 3-by-3 two-dimensional Ising model with a specific initial condition.

The yellow/light color in the image represents a spin of \(+1\), and the purple/dark color represents a spin of \(-1\).

📝 Questions:

  1. Calculate the average magnetization \(\frac{1}{N} \sum_i \sigma_i\).

  2. Calculate the energy from pairs associated with the central cell [with indices (1,1)]:

    \[ - \sum_{j \in D} \sigma_i \sigma_j \text{, where } i=(1,1) \text{ and } D = \text{neighbors of } i \]
  3. If the spin of the central cell flips from \(+1\) to \(-1\), how much will the total energy of the system change?

  4. How do your answers to (2) and (3) relate to each other? Briefly explain.

  5. Verify your answers with the code below. Does the result agree with your calculations?

// Add your answers to Part A here

# Do not edit this cell
m = IsingModel2D(side_length=3, initial_values=[1,-1,-1,-1,1,1,1,-1,1])
m.show()
../_images/48db1751d6a23a38a4ece4beac0967b3e8525a2de40a6e6c7b33ca6770d4c226.png
# Uncomment and run the following codes to verify your answers

#print("1. Average magnetization:", m.calc_magnetization())

#print("2. Local energy of site (1,1):", m.calc_local_energy(site_index=(1,1)))

#print("3. Total energy before flip:", m.calc_total_energy())
#m.flip(site_index=(1,1))
#print("3. Total energy after flip:", m.calc_total_energy())

B. Critical Temperature#

Let’s first take a look at how the sweep method works. When calling sweep, the argument beta sets the temperature of the system with \(\beta = 1/kT\). The argument nsweeps sets the number of sweeps to perform. Each sweep updates all sites in the lattice once.

m = IsingModel2D(side_length=3, initial_values="hot")
m.show()

m.sweep(beta=0.3, nsweeps=1)
m.show()
../_images/34d2d7ee3d42f3e77c3d8015eb3c9bba17d83f490db706756e4122fbb98915bc.png ../_images/6764820542593d493b14707700fd9a7f48acb421d46e383c5b2beb5a84323b56.png

We will run a simulation and gradually lower the temperature, in order to observe how the magnetization changes with temperature.

🚀 Tasks and Questions:

  1. Set up a 25-by-25 lattice with the “hot” initial condition.

  2. We will vary beta (\(\beta = 1/kT\)) from 0.1 to 0.9, with a step of 0.02.

  3. For each value of beta, we will run 20 “sweeps”.

  4. We will measure the average magnetization (averaged over sites) for each value of beta, after the 20 sweeps ran.

  5. Make a plot of the average magnetization as a function of beta (\(m\)-\(\beta\) relation).

  6. Describe what you observe. Does the magnetization show a sudden change at a certain value of \(\beta\)?

# Complete "..." in the the follow code
from tqdm.notebook import tqdm

m = IsingModel2D(...)

beta_values = np.arange(0.1, 0.91, 0.02)
magnetization_values = []

for beta in tqdm(beta_values):
    m.sweep(...)
    magnetization_values.append(m.calc_magnetization())
fig, ax = plt.subplots(dpi=150)
ax.plot(beta_values, magnetization_values)
ax.grid(True)
ax.set_xlabel("$\\beta = 1/kT$")
ax.set_ylabel("$\\langle m \\rangle$ Average magnetization")

// Write your answers to Part B here

C. Your Experiment!#

🚀 Tasks and Questions:

Choose at least one of the following conditions to experiment. Observe how the \(m\)-\(\beta\) relation changes:

  • as the external field strength changes (choose a few values from 0 to 2; set ext_field when sweep is called);

  • as the lattice size changes (choose a few numbers between 5 and 100; set side_length when setting up the model); or

  • as the boundary condition changes (try 1, -1, [1,1,-1,-1], [1,-1,1,-1]; set boundaries when setting up the model)

and describe what you observe.

Use the earlier example as a reference to for how to implement you experiment. If you have time, you are welcome to do more than one experiment, or come up with your own ideas.

# Include your implementation for Part C here

// Write your answers to Part C here

Tip

Submit your notebook

Follow these steps when you complete this lab and are ready to submit your work to Canvas:

  1. Check that all your text answers, plots, and code are all properly displayed in this notebook.

  2. Run the cell below.

  3. Download the resulting HTML file 07.html and then upload it to the corresponding assignment on Canvas.

!jupyter nbconvert --to html --embed-images 07.ipynb