Lab 7: Ising Model

For the class on Monday, February 5th

Tip

In this lab, questions are embedded in the implementation steps, but are always shown in bold. Don’t forget to answer the questions.

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.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

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

        indices = np.arange(self._size).reshape((self._side_size, self._side_size))
        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_size, self._side_size).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_size + 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\).

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?

---

// Write your answers to Part A here

---

# Do not edit this cell
m = IsingModel2D(side_size=3, initial_values=[1,-1,-1,-1,1,1,1,-1,1])
m.show()

# 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())

Demo: how the sweep method works

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

m.sweep(0.3, nsweeps=1)
m.show()

B. Critical Temperature

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

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

2. We will vary beta (\(\beta = 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\)?

---

// Write your answers to Part B here

---

# Include your implementation for Part B here

C. Experiment!

Choose at least one of the following to experiment, and 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_size when setting up the model)

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

Describe what you observe.

---

// Write your answers to Part C here

---

# Include your implementation for Part C here

Final question

Roughly how much time did you spend on this lab (not including the time you spent in class)?

---

// Write your answers here

---

Tip

How to submit this notebook on Canvas?

Once you complete this lab and have all your code, text answers, and all the desired results after running the notebook properly displayed in this notebook, please convert the notebook into HTML format by running the following:

jupyter nbconvert --to html /path/to/labs/07.ipynb 

Then, upload the resulting HTML file to Canvas for the corresponding assignment.