notebook/statistical-mechanics/ising_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Ising Model in 2D**\n",
    "\n",
    "<i class=\"fa fa-home fa-2x\"></i><a href=\"../index.ipynb\" style=\"font-size: 20px\"> Go back to index</a>\n",
    "\n",
    "**Source code:** https://github.com/osscar-org/quantum-mechanics/blob/master/notebook/statistical-mechanics/ising_model.ipynb\n",
    "\n",
    "The Cython computing kernel is adapted from the code presented in the \n",
    "[Pythonic Perambulations](https://jakevdp.github.io/blog/2017/12/11/live-coding-cython-ising-model) blog.\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goals\n",
    "\n",
    "* Understand how the Ising model can predict the ferromagnetic or antiferromagnetic behavior of a square (2D) spin lattice.\n",
    "* Compute the average magnetization of the Ising model using a Monte Carlo algorithm.\n",
    "* Examine the effects of the interaction parameter $J$ (ferromagnetic vs. antiferromagnetic behavior).\n",
    "* Identify the critical temperature of the system from the simulation results.\n",
    "* Investigate finite-size effects and their effect on fluctuations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Background Theory**\n",
    "\n",
    "[More in the background theory](./theory/theory_ising_model.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## **Tasks and exercises**\n",
    "\n",
    "1. Run the simulation with the default parameters (and $J>0$) and check the plot of the magnetization per spin as a function of the simulation step. What happens at very low temperatures? And at intermediate temperatures? And at very high temperatures? And what happens for $J<0$?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "  For a positive $J$ value, parallel neighboring spins are energetically favored.\n",
    "  Therefore, at $T=0$ we expect all spins to be aligned: the system is ferromagnetic.\n",
    "  At low temperature, spin will still tend to be aligned with their neighbors, with some\n",
    "  local fluctuations as temperature increases.\n",
    "  Fluctuations become very large close to $T_c$, and above $T_c$\n",
    "  they become predominant and, on average, we obtain zero net magnetization.\n",
    "  On the other hand, $J<0$ will lead to a checkerboard pattern for the final spin \n",
    "  configuration at low temperatures, i.e., an antiferromagnetic configuration.\n",
    "    When $J=0$, there is no interaction between spins and we always obtain a random spin configuration.<br><br>\n",
    "    </details>\n",
    "\n",
    "2. Simulate a large system (200x200) with, e.g., $J=1$ and $T=2.5$. Is the simulation converging with the selected number of steps? Which values can you inspect to determine the convergence of the simulation?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "  A large number of simulation steps can be needed before properties converge, especially when we are close to $T_C$.<br>\n",
    "  \n",
    "   You can check the plot of the integrated quantities (total energy, magnetization per spin) to see if the value is converged (their values are constant save for small fluctuations with simulation step), or if you need to increase the number of steps.<br><br>\n",
    "    </details>\n",
    "\n",
    "3. Set $J=1$ and start from all spin up (i.e., disable a random starting configuration). Run multiple simulations at various temperatures. Can you approximately identify the critical (Curie) temperature using the plot below the interactive controls? \n",
    "\n",
    "  <details>\n",
    "   <summary style=\"color: red\">Solution</summary>\n",
    "  We suggest that you start from a small size system (50x50) to avoid excessively long simulations, and choose a long number of simulation steps (1000) to be approximately converged (see previous question).<br>\n",
    "\n",
    "  Theoretically, $T_C$ for the 2D Ising model (and $J=1$) is about 2.27 (see the <a href=\"https://en.wikipedia.org/wiki/Square_lattice_Ising_model\">Wikipedia page of the Square lattice Ising model</a>), as visualized by the analytical exact solution plotted below the interactive controls.<br>\n",
    "    \n",
    "  You should see that the simulation results follow relatively closely the analytical curve for $T \\ll T_c$ or for $T \\gg T_C$, but large fluctuations occur close to the Curie temperature. We suggest that you run multiple runs for each temperature to see the fluctuations. You should be able to reproduce the curve and thus approximately identify the transition temperature.<br>\n",
    "    \n",
    "  We suggest that you try the same with different values of $J$.<br><br>\n",
    "  </details>\n",
    "\n",
    "4. Set $J=1$, and simulate the system close to $T_C$. Investigate the magnitude of fluctuations as a function of system size.\n",
    "\n",
    "  <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "  Consider for instance both a small (50x50) and a large system (200x200), $J=1$, and $T=2.5$. Run a few simulations (at least 5 or 10) for each system size (note that each simulation for the large system will take many seconds). Verify that fluctuations are less pronounced for the larger system.<br>\n",
    "    \n",
    "  Note that the difference in fluctuations might not be very large, since the two systems are still relatively similar in size. The difference in the magnitude of the fluctuation can also be visualized by inspecting the fluctuations of the averaged quantities in the two plots of the magnetization per spin and total energy as a function of the simulation step.\n",
    "  </details>\n",
    "\n",
    "5. Consider a 100x100 system, set $J=1$ and $T=1.7$. Use a large number of simulation steps, and enable the randomization of the initial spin configuration. Investigate the formation of domains. Are they more stable for large or small systems?\n",
    "\n",
    "  <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "  If you run the simulation multiple times, you will notice that often the simulation does not reach the expected (positive or negative) analytical exact result, but will have intermediate magnetization values due to the formation of domains. These might disappear during the simulation, but sometimes they will remain even after 1000 simulation steps.<br>\n",
    "    \n",
    "  If you consider a smaller system (e.g. 50x50), you should notice that there is a higher probability that the domains disappear during the simulation. In fact, the cost of a domain scales as the length of its boundary, which (for very large systems) becomes a negligible cost with respect to the total energy of the system (that scales with the number of spins, i.e. quadratically with respect to the system size) and therefore can exist for a long time; in addition, domains will have a higher probability to merge for a small system.\n",
    "  </details>\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext Cython\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ipywidgets import interact, FloatSlider, Button, Output, IntSlider, VBox\n",
    "from ipywidgets import HBox, Checkbox, IntProgress, HTML\n",
    "import matplotlib.pyplot as plt\n",
    "from time import sleep\n",
    "import matplotlib.gridspec as gridspec\n",
    "from scipy.ndimage import convolve\n",
    "from matplotlib.animation import FuncAnimation\n",
    "import base64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_spin_field(N, M):\n",
    "    \"\"\"Randomize the initial spin configuration.\"\"\"\n",
    "    return np.random.choice([-1, 1], size=(N, M))\n",
    "\n",
    "def all_up_spin_field(N, M):\n",
    "    \"\"\"Set all spin up.\"\"\"\n",
    "    return np.ones((N, M), dtype=int)\n",
    "\n",
    "run_button = Button(description='Run simulation')\n",
    "run_button.style.button_color = 'lightgreen'\n",
    "play_button = Button(description='Play', disabled=True)\n",
    "\n",
    "random_checkbox = Checkbox(value=False, description=\"Randomize initial spin configuration\", style={'description_width': 'initial'})\n",
    "\n",
    "jvalue_slider = FloatSlider(value = 1.0, min = -2.0, max = 2.0, description = 'Exchange interaction J',\n",
    "                    style={'description_width': 'initial'}, continuous_update=False)\n",
    "num_slider = IntSlider(value=100, min=50, max=200, step=10, description=\"Size\", continuous_update=False)\n",
    "temp_slider = FloatSlider(value=2, min=0.5, max=4, step=0.1, description=\"Temperature\", continuous_update=False) # Units of J/k\n",
    "step_slider = IntSlider(value=100, min=100, max=1000, step=50, description=\"Num steps\", continuous_update=False)\n",
    "frame_slider = IntSlider(value=0, min=0, max=step_slider.value, description=\"Frame\", layout={'width':'800px'}, disabled=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "cimport cython\n",
    "\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "\n",
    "from libc.math cimport exp\n",
    "from libc.stdlib cimport rand\n",
    "cdef extern from \"limits.h\":\n",
    "    int RAND_MAX\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "def cy_ising_step(np.int64_t[:, :] field, float beta, float J):\n",
    "    \"\"\"Update the Ising step, each step actually contains N*M steps.\n",
    "       \n",
    "    Args:\n",
    "        field: matrix representation for the spin configuration.\n",
    "        beta: 1/kbT\n",
    "        J: the strength of exchange interaction.\n",
    "        \n",
    "    Returns:\n",
    "        New spin configuration and the total energy change.\n",
    "    \"\"\"\n",
    "    cdef int N = field.shape[0]\n",
    "    cdef int M = field.shape[1]\n",
    "    cdef int i\n",
    "    cdef np.ndarray x, y\n",
    "    cdef float dE = 0.0\n",
    "    \n",
    "    x = np.random.randint(N, size=(N*M))\n",
    "    y = np.random.randint(M, size=(N*M))\n",
    "    \n",
    "    for i in range(N*M):\n",
    "        dE += _cy_ising_update(field, x[i], y[i], beta, J)\n",
    "        \n",
    "    return np.array(field), dE\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "cdef _cy_ising_update(np.int64_t[:, :] field, int n, int m, float beta, float J):\n",
    "    \"\"\"Monte Carlo simulation using the Metropolis algorithm.\n",
    "    \n",
    "    Args:\n",
    "        field: matrix representation for the spin configuration.\n",
    "        n: chosen row index.\n",
    "        m: chosen column index.\n",
    "        beta: 1/kbT\n",
    "        J: the strength of exchange interaction.\n",
    "        \n",
    "    Returns:\n",
    "        The total energy change.\n",
    "    \"\"\"\n",
    "    cdef int total\n",
    "    cdef int N = field.shape[0]\n",
    "    cdef int M = field.shape[1]\n",
    "\n",
    "    total = field[(n+1)%N, m] + field[n, (m+1)%M] + field[(n-1)%N, m] + field[n, (m-1)%M]\n",
    "    cdef float dE = 2.0 * J * field[n, m] * total\n",
    "    if dE <= 0:\n",
    "        field[n, m] *= -1\n",
    "        return dE\n",
    "    elif exp(-dE * beta) * RAND_MAX > rand():\n",
    "        field[n, m] *= -1\n",
    "        return dE\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pause = True;\n",
    "\n",
    "def on_frame_change(b):\n",
    "    \"\"\"Update the plot for playing the animation.\"\"\"\n",
    "    global fig, v1\n",
    "    fig.set_data(images[frame_slider.value])\n",
    "    v1.set_data([frame_slider.value, frame_slider.value],[-1.1, 1.1])\n",
    "    v2.set_data([frame_slider.value, frame_slider.value],[-5000, 5000])\n",
    "\n",
    "    \n",
    "def compute_total_energy(M, J):\n",
    "    \"\"\"Compute the total energy of the given spin configuration.\"\"\"\n",
    "    a = np.ones(np.shape(M));\n",
    "    c = convolve(M, a, mode='constant')\n",
    "    c = (c-M)*M*J\n",
    "    return c.sum()\n",
    "\n",
    "\n",
    "def update(frame):\n",
    "    \"\"\"Update function for the animation.\"\"\"\n",
    "    global pause\n",
    "    \n",
    "    if pause:\n",
    "        ani.event_source.stop()\n",
    "    else:\n",
    "        frame_slider.value = frame\n",
    "        return (fig)\n",
    "    \n",
    "def play_animation(event):\n",
    "    \"\"\"OnClick function the 'Play' button.\"\"\"\n",
    "    global pause\n",
    "    \n",
    "    pause ^= True\n",
    "    if play_button.description == \"Pause\":\n",
    "        play_button.description = \"Play\"\n",
    "        ani.event_source.stop()\n",
    "    else:\n",
    "        play_button.description = \"Pause\"\n",
    "        ani.event_source.start()\n",
    "    \n",
    "\n",
    "frame_slider.observe(on_frame_change, names='value')\n",
    "\n",
    "images = [all_up_spin_field(num_slider.value, num_slider.value)]\n",
    "\n",
    "img = plt.figure(tight_layout=True, figsize=(8,5))\n",
    "img.canvas.header_visible = False\n",
    "gs = gridspec.GridSpec(4, 2)\n",
    "\n",
    "ax1 = img.add_subplot(gs[:, 0])\n",
    "ax2 = img.add_subplot(gs[0:2, 1])\n",
    "ax3 = img.add_subplot(gs[2:4, 1])\n",
    "\n",
    "fig = ax1.imshow(images[0], vmin=-1, vmax=1)\n",
    "ax1.axes.xaxis.set_ticklabels([])\n",
    "ax1.axes.yaxis.set_ticklabels([])\n",
    "ax1.set_title('Spin up (yellow), spin down (purple)', fontsize=12)\n",
    "\n",
    "line1, = ax2.plot([0], [0], 'r-')\n",
    "line2, = ax3.plot([0], [0], 'r-')\n",
    "v1 = ax2.axvline(x=0, c='black')\n",
    "v2 = ax3.axvline(x=0, c='black')\n",
    "\n",
    "ax2.set_xlim([0, step_slider.value])\n",
    "ax2.set_ylim([-1.1, 1.1])\n",
    "ax2.set_title('Magnetization per spin', fontsize=12)\n",
    "\n",
    "ax3.set_xlim([0, step_slider.value])\n",
    "ax3.set_ylim([-1, 1])\n",
    "ax3.set_title(r'Total energy (E$_{init}$=0)', fontsize=12)\n",
    "ax3.set_xlabel('Step', fontsize=12)\n",
    "\n",
    "ani = FuncAnimation(img, update, interval= 20, frames=np.arange(0, step_slider.value+1), blit=True)\n",
    "\n",
    "display(frame_slider,\n",
    "    HBox([num_slider, step_slider]),\n",
    "    HBox([jvalue_slider, temp_slider]),\n",
    "    HBox([random_checkbox, run_button, play_button]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_analytical_plot_data(J, k_B=1):\n",
    "    \"\"\"Exact solution of the 2D square Ising model.\n",
    "    \n",
    "    See e.g.: https://en.wikipedia.org/wiki/Square_lattice_Ising_model\n",
    "    \"\"\"\n",
    "    temperature = np.linspace(0.01, 4, 500, dtype=np.float128)\n",
    "    magnetization = np.zeros(len(temperature))\n",
    "    Tc = 2 * J / k_B / np.log(1 + np.sqrt(2))\n",
    "    magnetization[temperature < Tc] = (1-1./(np.sinh(2 * J / k_B / temperature[temperature < Tc])**4))**(1/8)\n",
    "\n",
    "    return temperature, magnetization\n",
    "\n",
    "def reset_mag_plot(ax):\n",
    "    \"\"\"Reset the magnetization plot.\n",
    "    \n",
    "    Clear the axes and re-draw the analytical exact solution.\n",
    "    \"\"\"\n",
    "    global ax_mag, random_checkbox\n",
    "    temps, magnetization = get_analytical_plot_data(J=jvalue_slider.value, k_B=1)\n",
    "    \n",
    "    ax.clear()\n",
    "    ax.plot(temps, magnetization, 'b-')\n",
    "    if random_checkbox.value:\n",
    "        ax.plot(temps, -magnetization, 'r-')        \n",
    "    ax.set_xlabel(\"$T$\")\n",
    "    ax.set_ylabel(r\"$\\langle \\sigma \\rangle$\")\n",
    "    ax.set_title(\"Magnetization per site vs. temperature\")\n",
    "    ax.set_xlim((0, 4))\n",
    "\n",
    "def create_mag_plot():\n",
    "    \"\"\"Create the magnetization plot.\n",
    "    \n",
    "    To be called only one to create it and return the axes object.\n",
    "    \"\"\"\n",
    "    img_mag = plt.figure(tight_layout=True, figsize=(5,3.5))\n",
    "    img_mag.canvas.header_visible = False\n",
    "    gs_mag = gridspec.GridSpec(1, 1)\n",
    "    ax_mag = img_mag.add_subplot(gs_mag[:, :])\n",
    "\n",
    "    reset_mag_plot(ax_mag)\n",
    "    \n",
    "    return ax_mag    \n",
    "\n",
    "# Create and display the magnetization plot.\n",
    "ax_mag = create_mag_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(b):\n",
    "    \"\"\"Callback to be called when the 'Run simulation' button is clicked.\n",
    "    \n",
    "    Takes care of temporarily disabling the button, running the simulation and updating\n",
    "    the various plots, and re-enabling the button again.\n",
    "    \"\"\"\n",
    "    global ax_mag\n",
    "\n",
    "    play_button.disabled = True\n",
    "    run_button.disabled = True\n",
    "    frame_slider.disabled = True\n",
    "    run_button.style.button_color = 'red'\n",
    "    global images, fig\n",
    "    if random_checkbox.value:\n",
    "        images = [random_spin_field(num_slider.value, num_slider.value)]\n",
    "    else:\n",
    "        images = [all_up_spin_field(num_slider.value, num_slider.value)]\n",
    "    \n",
    "    x = np.arange(step_slider.value + 1)\n",
    "    y1 = []\n",
    "    y2 = [0]\n",
    "\n",
    "    for i in range(step_slider.value):\n",
    "        imag, dE = cy_ising_step(images[-1].copy(), beta=1.0/temp_slider.value, J=jvalue_slider.value)\n",
    "        images.append(imag)\n",
    "        y2.append(dE+y2[-1])\n",
    "\n",
    "    frame_slider.max = step_slider.value\n",
    "    ax2.set_xlim([0, step_slider.value])\n",
    "\n",
    "    fig.set_data(images[frame_slider.max - 1])\n",
    "\n",
    "    for i in images:\n",
    "        y1.append(i.sum()*1.0/(num_slider.value * num_slider.value))\n",
    "\n",
    "    y1 = np.array(y1)\n",
    "    y2 = np.array(y2)\n",
    "\n",
    "    line1.set_data(x, y1)\n",
    "    line2.set_data(x, y2)\n",
    "\n",
    "    ax3.set_ylim([y2.min(), y2.max()])\n",
    "    ax3.set_xlim([0, step_slider.value])\n",
    "\n",
    "    frame_slider.value = frame_slider.max\n",
    "\n",
    "    ani.frames = np.arange(0, step_slider.value+1)\n",
    "    ani._iter_gen = lambda: iter(ani.frames)\n",
    "    ani.save_count = len(ani.frames)\n",
    "    ani.frame_seq = ani.new_frame_seq()\n",
    "\n",
    "    ax_mag.errorbar([temp_slider.value], [y1[-1]], fmt='.g')\n",
    "    \n",
    "    frame_slider.disabled = False\n",
    "    play_button.disabled = False\n",
    "    run_button.disabled = False\n",
    "    run_button.style.button_color = 'lightgreen'\n",
    "    return y1\n",
    "\n",
    "def on_needs_reset(b):\n",
    "    \"\"\"Callback to be called when the magnetization plot needs to be cleared and reset.\"\"\"\n",
    "    global ax_mag\n",
    "    reset_mag_plot(ax_mag)\n",
    "\n",
    "# Attach actions to the buttons\n",
    "run_button.on_click(run_simulation)\n",
    "play_button.on_click(play_animation)\n",
    "# Attach reset actions (for the magnetization plot) when some of the\n",
    "# simulation parameters are changed (all except the temperature)\n",
    "jvalue_slider.observe(on_needs_reset, names='value')\n",
    "num_slider.observe(on_needs_reset, names='value')\n",
    "step_slider.observe(on_needs_reset, names='value')\n",
    "random_checkbox.observe(on_needs_reset, names='value')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Legend\n",
    "\n",
    "(How to use the interactive visualization)\n",
    "\n",
    "### Controls\n",
    "\n",
    "The \"size\" slider defines the number of spins along each dimension.\n",
    "You can also adapt the number of simulation steps, the value of the\n",
    "exchange interaction parameter $J$, and the temperature of the simulation\n",
    "(note: units have been chosen so that the Boltzmann constant $k_B=1$).\n",
    "\n",
    "By default, we set all spin up for the initial configuration.\n",
    "This induces a bias for short simulations (but can help avoid the formation of\n",
    "domains); you can instead randomize the initial configuration by ticking\n",
    "the checkbox \"Randomize initial spin configuration\".\n",
    "\n",
    "### Running the simulation and output plots\n",
    "Click the \"Run simulation\" button to start the simulation with the selected parameters.\n",
    "\n",
    "In particular, the figures above the controls display the time evolution of the last simulation\n",
    "that was executed.\n",
    "The figure on the left shows the spin configuration. A yellow pixel represents \n",
    "a spin up and blue represents a spin down.\n",
    "After the simulation finished, you can click \"Play\" to view the evolution of the\n",
    "simulation step by step (or drag manually the \"Frame\" slider).\n",
    "The figure on the top right shows the evolution of the magnetization per spin over\n",
    "the simulation as a function of the Monte Carlo simulation step number.\n",
    "The figure on the bottom right shows the total energy as a function of the step number.\n",
    "\n",
    "In addition, at every new simulation, a new green point is added to the bottom plot,\n",
    "showing the magnetization per spin as a function of temperature $T$.\n",
    "The analytical exact solution for the chosen value of $J$ is also shown\n",
    "(see [Wikipedia page on the square lattice Ising model](https://en.wikipedia.org/wiki/Square_lattice_Ising_model)\n",
    "for more details on the exact solution). Note that when a random initial configuration is selected instead of all spin up, there shall be both blue and red curves for the analytical solution due to the possibility of having negative magnetization in that case.\n",
    "\n",
    "Changing any slider (except for the temperature slider)\n",
    "will reset the bottom plot.\n",
    "You can change the temperature and run the simulation multiple times\n",
    "to compare the analytical solution with the results of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "voila": {
   "authors": "Dou Du and Giovanni Pizzi"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}