{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bbec08f5-2ccb-4e82-8184-afb1e9e507ab",
   "metadata": {
    "id": "bbec08f5-2ccb-4e82-8184-afb1e9e507ab"
   },
   "source": [
    "# 2D Ising model Monte Carlo\n",
    "\n",
    "MC code to simulate Ising model in 2D\n",
    "$$H = -J \\sum_{\\langle ij \\rangle} s_i s_j - h \\sum_i s_i$$\n",
    "This model, for $J>0$ has a Curie temperature of\n",
    "$$T_\\text{c} = \\frac{2J}{\\ln\\left(\\sqrt{2}+1\\right)}\\approx 2.269J$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af02609f-a296-492d-9211-e577d213da41",
   "metadata": {
    "id": "af02609f-a296-492d-9211-e577d213da41"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "rng = np.random.default_rng()\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [6, 6]\n",
    "from numba import jit\n",
    "from tqdm.auto import trange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "682ad713-6811-441f-ba38-16c2e8b4bbfb",
   "metadata": {
    "id": "682ad713-6811-441f-ba38-16c2e8b4bbfb"
   },
   "outputs": [],
   "source": [
    "# simple ferromagnetic coupling\n",
    "param = {'J': 1.0, 'h': 0.0}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "884541d6-2477-4a81-9b92-4295b9e8181c",
   "metadata": {
    "id": "682ad713-6811-441f-ba38-16c2e8b4bbfb"
   },
   "outputs": [],
   "source": [
    "# set up cell and PBC\n",
    "Nx, Ny = 16, 16\n",
    "N = Nx*Ny\n",
    "\n",
    "d, z = 2, 4 # dimensionality and coordination number\n",
    "site = np.zeros((N, d), dtype=int)\n",
    "neigh = np.zeros((N, z), dtype=int)\n",
    "\n",
    "n = 0\n",
    "for i in range(Nx):\n",
    "    for j in range(Ny):\n",
    "        site[n] = i, j\n",
    "        n += 1\n",
    "\n",
    "for n in range(N):\n",
    "    i, j = site[n]\n",
    "    neigh[n, 0] = ((i+1)%Nx)*Ny + j\n",
    "    neigh[n, 1] = ((i-1)%Nx)*Ny + j\n",
    "    neigh[n, 2] = i*Ny + (j+1)%Ny\n",
    "    neigh[n, 3] = i*Ny + (j-1)%Ny"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a84fb7d-e46f-411d-bfb8-88e082b59e4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# neighbor list check\n",
    "if (False):\n",
    "    for n in range(N):\n",
    "        print(f'{site[n]}:')\n",
    "        for m in neigh[n]:\n",
    "            print(f'  {site[m]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7777ff0f-9bf7-4cd8-8b7c-9215e5e45fcb",
   "metadata": {
    "id": "7777ff0f-9bf7-4cd8-8b7c-9215e5e45fcb"
   },
   "outputs": [],
   "source": [
    "# initialize our spins on a square lattice\n",
    "spin = rng.choice((-1, 1), size=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f1acbfa-7e22-4c73-8522-57637664dbb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize\n",
    "plt.imshow(spin.reshape(Nx, Ny), cmap=\"grey\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0139e74e-4f1c-4ef9-99db-e3c5ae710ebf",
   "metadata": {
    "id": "0139e74e-4f1c-4ef9-99db-e3c5ae710ebf"
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def ene(spin, neigh, J=param['J'], h=param['h']):\n",
    "    \"\"\"Computes the energy for a given configuration `spin` using `neigh` as a neighbor\n",
    "    list, with `J` and `h` as our parameters\"\"\"\n",
    "    N = neigh.shape[0]\n",
    "    z = neigh.shape[1]\n",
    "    E_J = 0\n",
    "    for n in range(N):\n",
    "        neigh_list = neigh[n]\n",
    "        s = spin[n]\n",
    "        for i in range(z):\n",
    "            if s==spin[neigh_list[i]]:\n",
    "                E_J += 1\n",
    "            else:\n",
    "                E_J += -1\n",
    "    return -0.5*J*E_J - h*sum(spin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f81430a-47fc-4f9f-9fef-861fae16d6f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "ene(spin, neigh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e707f176-f500-41ba-b6bc-e3a59e9725d1",
   "metadata": {
    "id": "0139e74e-4f1c-4ef9-99db-e3c5ae710ebf"
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def MC(spin, kT, urand, neigh, J=param['J'], h=param['h']):\n",
    "    \"\"\"Run MC for Ising spins where we attempt to flip each spin in order.\n",
    "    Returns acceptance ratio and (relative) energy, and overwrites spin in place.\"\"\"\n",
    "    N = neigh.shape[0]\n",
    "    z = neigh.shape[1]\n",
    "    \n",
    "    accept = 0\n",
    "    E = 0\n",
    "    for i in range(len(urand)):\n",
    "        u = urand[i]\n",
    "        n = i%N # spin to attempt to flip\n",
    "        s = spin[n]\n",
    "        dE = 0\n",
    "        for j in range(z):\n",
    "            m = neigh[n, j]\n",
    "            if s == spin[m]:\n",
    "                dE += -2\n",
    "            else:\n",
    "                dE += 2\n",
    "        dE = -dE*J + 2*h*s\n",
    "        if u < np.exp(-dE/kT):\n",
    "            spin[n] = -s\n",
    "            E += dE\n",
    "            accept += 1\n",
    "    return accept/len(urand), E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af05a73c-3e06-4ba3-a828-301a1fc2fd3f",
   "metadata": {
    "id": "af05a73c-3e06-4ba3-a828-301a1fc2fd3f"
   },
   "outputs": [],
   "source": [
    "N_MC = N\n",
    "\n",
    "kT = 2.5\n",
    "E = ene(spin, neigh)\n",
    "urand = rng.uniform(size=(N_MC))\n",
    "accept, dE = MC(spin, kT, urand, neigh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a316999-e3e7-4177-9ffb-0656851fe331",
   "metadata": {
    "id": "6a316999-e3e7-4177-9ffb-0656851fe331"
   },
   "outputs": [],
   "source": [
    "accept, dE, np.average(spin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "636aecd5-e3b9-4743-b4bd-409708ffab8b",
   "metadata": {
    "id": "636aecd5-e3b9-4743-b4bd-409708ffab8b"
   },
   "outputs": [],
   "source": [
    "# visualize\n",
    "plt.imshow(spin.reshape(Nx, Ny), cmap=\"grey\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "xMcfuTzxdBHY",
   "metadata": {
    "id": "xMcfuTzxdBHY"
   },
   "outputs": [],
   "source": [
    "ene(spin, neigh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "TiLmWr31kbjb",
   "metadata": {
    "id": "TiLmWr31kbjb"
   },
   "outputs": [],
   "source": [
    "N_sweeps = 1024\n",
    "E_list, accept_list, m_list = [], [], []\n",
    "E = ene(spin, neigh)\n",
    "for _ in trange(N_sweeps):\n",
    "    accept, dE = MC(spin, kT, rng.uniform(size=N), neigh)\n",
    "    E += dE\n",
    "    E_list.append(E/N)\n",
    "    accept_list.append(accept)\n",
    "    m_list.append(np.average(spin))\n",
    "E_list = np.array(E_list)\n",
    "accept_list = np.array(accept_list)\n",
    "m_list = np.array(m_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8fade96-0a98-4de0-8714-1d26a1165c52",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.average(accept_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49fbd63b-9f1a-4160-b66e-b7fc7394a3cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(E_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03013af7-c18b-4ffa-a7f2-e6bd9a0b7af0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(spin.reshape(Nx, Ny), cmap='grey');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba1305ec-71fb-4aad-b5b0-05cfece6242a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01f40d8e-8fc7-42db-99ac-a50b571c79da",
   "metadata": {},
   "outputs": [],
   "source": [
    "kT = 2.0\n",
    "\n",
    "N_sweeps = 1024\n",
    "E_list, accept_list, m_list = [], [], []\n",
    "E = ene(spin, neigh)\n",
    "for _ in trange(N_sweeps):\n",
    "    accept, dE = MC(spin, kT, rng.uniform(size=N), neigh)\n",
    "    E += dE\n",
    "    E_list.append(E/N)\n",
    "    accept_list.append(accept)\n",
    "    m_list.append(np.average(spin))\n",
    "E_list = np.array(E_list)\n",
    "accept_list = np.array(accept_list)\n",
    "m_list = np.array(m_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9740055d-45e7-47a6-a9ea-65c7d5fffb2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.average(accept_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd21fa6f-5353-4c3c-bf5f-588f9187b336",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(E_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "915a5bb0-0245-4cb1-8f4f-7dded065ed22",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(spin.reshape(Nx, Ny), cmap='grey');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8435a90a-ce76-4c19-a438-7f37d941e1a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a40bcbfe-91bb-44f2-9460-65d9e29c868e",
   "metadata": {},
   "outputs": [],
   "source": [
    "kT = 4.0\n",
    "\n",
    "N_sweeps = 1024\n",
    "E_list, accept_list, m_list = [], [], []\n",
    "E = ene(spin, neigh)\n",
    "for _ in trange(N_sweeps):\n",
    "    accept, dE = MC(spin, kT, rng.uniform(size=N), neigh)\n",
    "    E += dE\n",
    "    E_list.append(E/N)\n",
    "    accept_list.append(accept)\n",
    "    m_list.append(np.average(spin))\n",
    "E_list = np.array(E_list)\n",
    "accept_list = np.array(accept_list)\n",
    "m_list = np.array(m_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8816747a-7fa5-48c5-b201-ae4a3d61e7ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.average(accept_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d02ce677-34b6-47c2-995e-f460a94bef7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(E_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "835e0c2b-7ad6-46f8-97fe-0708cfa0e7d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(spin.reshape(Nx, Ny), cmap='grey');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23c037f2-f4e9-4d7e-9ffc-3fbcbd245036",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_list);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e98193b5-d5fd-4d77-b1f9-b2718cab51d2",
   "metadata": {},
   "source": [
    "### Magnetization with temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e7afaa7-5ef3-44e0-b3dc-cd8a54b734fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "kT_list = np.linspace(3.0, 2.0, num=11)\n",
    "N_sweeps = 4096\n",
    "m_T = []\n",
    "for kT in kT_list:\n",
    "    m_list = []\n",
    "    for _ in trange(N_sweeps):\n",
    "        # just capture magnetization\n",
    "        MC(spin, kT, rng.uniform(size=N), neigh)\n",
    "        m_list.append(np.average(spin))\n",
    "    m_T.append(np.array(m_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c38e24a4-69e9-4e05-84e3-d551962dfb8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ave_m_T = np.array([np.average(m) for m in m_T])\n",
    "err_m_T = np.array([np.std(m)/np.sqrt(N_sweeps) for m in m_T])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef4ee08-35e7-4b19-8795-33cf5a1b67e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "err_m_T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "927c8a40-7759-4f56-83b4-1c39d66461fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(kT_list, ave_m_T, 30*err_m_T, fmt='o-');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3347b2df-7b6b-4940-b315-c92040bcd7d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argwhere(kT_list==2.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "272aec99-7f7a-4299-ac8f-9c91bc151748",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_T[8]);\n",
    "kT_list[8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bff9d468-e73b-455d-90d4-11595fdd29cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_T[7]);\n",
    "kT_list[7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ce6cbeb-206c-4111-af18-827a6863cc8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_T[6]);\n",
    "kT_list[6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ab110e9-1877-4161-b42c-d2bf3bbb8504",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_T[0]);\n",
    "kT_list[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "277973cf-87ff-498a-8b22-34bc94e0b552",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_T[-1]);\n",
    "kT_list[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4381c35d-5514-4686-8e0e-827d42ec1d6f",
   "metadata": {},
   "source": [
    "### Long run near Tc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4d9dc92-2850-468b-bd7e-c4e7e794fc8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "kT = 2.3\n",
    "\n",
    "N_sweeps = 16384\n",
    "E_list, accept_list, m_list = [], [], []\n",
    "E = ene(spin, neigh)\n",
    "for _ in trange(N_sweeps):\n",
    "    accept, dE = MC(spin, kT, rng.uniform(size=N), neigh)\n",
    "    E += dE\n",
    "    E_list.append(E/N)\n",
    "    accept_list.append(accept)\n",
    "    m_list.append(np.average(spin))\n",
    "E_list = np.array(E_list)\n",
    "accept_list = np.array(accept_list)\n",
    "m_list = np.array(m_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa4a2c77-eb25-45e1-b23b-b1599db683e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.average(accept_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "084e9f95-d72c-4114-9af0-43f9a9af1356",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "512299e2-4c17-40d3-b2ee-94b8d2bd4bde",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(spin.reshape(Nx, Ny), cmap='grey');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc11ebb2-7182-49a1-a28c-92abb2c544e1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "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.13.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}