{
 "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$$\n",
    "Now coded with more efficient lookup tables for Arrhenius factors."
   ]
  },
  {
   "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 = 256, 256\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": "407b47a8-1451-4598-a9f9-0a5bb97ca0f2",
   "metadata": {
    "id": "682ad713-6811-441f-ba38-16c2e8b4bbfb"
   },
   "outputs": [],
   "source": [
    "def arrhenius(kT, J=param['J'], h=param['h'], z=z):\n",
    "    \"\"\"Returns a matrix arr[si.sj,si] = exp(-(2J*(si.sj) + 2h*si)/kT) that is the Arrhenius factor where\n",
    "    si.sj = sum_j si*sj and si. Since we flip si -> -si, this is the Arrhenius factor for the change in \n",
    "    energy with such a flip.\"\"\"\n",
    "    arr = np.zeros((2*z+1, 3))\n",
    "    for sisj in range(-z, z+1):\n",
    "        for si in (-1,1):\n",
    "            arr[sisj, si] = np.exp(-2*(J*sisj + h*si)/kT)\n",
    "    return arr"
   ]
  },
  {
   "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",
    "def MC(spin, kT, urand, neigh, arr):\n",
    "    \"\"\"Run MC for Ising spins where we attempt to flip each spin in order.\n",
    "    Returns acceptance ratio only, and overwrites spin in place.\"\"\"\n",
    "    N = neigh.shape[0]\n",
    "    z = neigh.shape[1]\n",
    "    \n",
    "    accept = 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",
    "        sisj = 0\n",
    "        for j in range(z):\n",
    "            m = neigh[n, j]\n",
    "            if s == spin[m]:\n",
    "                sisj += 1\n",
    "            else:\n",
    "                sisj += -1\n",
    "        if u < arr[sisj, s]:\n",
    "            spin[n] = -s\n",
    "            accept += 1\n",
    "    return accept/len(urand)"
   ]
  },
  {
   "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",
    "arr = arrhenius(kT)\n",
    "E = ene(spin, neigh)\n",
    "urand = rng.uniform(size=(N_MC))\n",
    "accept = MC(spin, kT, urand, neigh, arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a316999-e3e7-4177-9ffb-0656851fe331",
   "metadata": {
    "id": "6a316999-e3e7-4177-9ffb-0656851fe331"
   },
   "outputs": [],
   "source": [
    "accept, 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",
    "for _ in trange(N_sweeps):\n",
    "    accept = MC(spin, kT, rng.uniform(size=N), neigh, arr)\n",
    "    E_list.append(ene(spin, neigh)/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",
    "arr = arrhenius(kT)\n",
    "\n",
    "N_sweeps = 1024\n",
    "E_list, accept_list, m_list = [], [], []\n",
    "for _ in trange(N_sweeps):\n",
    "    accept = MC(spin, kT, rng.uniform(size=N), neigh, arr)\n",
    "    E_list.append(ene(spin, neigh)/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",
    "arr = arrhenius(kT)\n",
    "\n",
    "N_sweeps = 1024\n",
    "E_list, accept_list, m_list = [], [], []\n",
    "for _ in trange(N_sweeps):\n",
    "    accept = MC(spin, kT, rng.uniform(size=N), neigh, arr)\n",
    "    E_list.append(ene(spin, neigh)/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": "code",
   "execution_count": null,
   "id": "7e53f4e6-5a7b-4595-87fe-a50ab55efc9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(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 = 2**14\n",
    "m_T = []\n",
    "for kT in kT_list:\n",
    "    arr = arrhenius(kT)\n",
    "    m_list = []\n",
    "    for _ in trange(N_sweeps):\n",
    "        # just capture magnetization\n",
    "        MC(spin, kT, rng.uniform(size=N), neigh, arr)\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.0/(np.log(np.sqrt(2)+1))*param['J']\n",
    "kT = 2.3\n",
    "arr = arrhenius(kT)\n",
    "\n",
    "N_sweeps = 2**16\n",
    "E_list, accept_list, m_list = [], [], []\n",
    "for _ in trange(N_sweeps):\n",
    "    accept = MC(spin, kT, rng.uniform(size=N), neigh, arr)\n",
    "    E_list.append(ene(spin, neigh)/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": [
    "plt.plot(m_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa9d3292-039c-4c64-8f41-9f3f0b3fe44c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def acf_calc(v):\n",
    "    \"\"\"Returns the auto-correlation function for the sequence v\"\"\"\n",
    "    v_ave = np.mean(v)\n",
    "    v_len = len(v)\n",
    "    return np.correlate(v-v_ave, v-v_ave, \"full\")[(v_len-1):] / \\\n",
    "        (np.var(v)*np.array([v_len - x for x in range(v_len)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04ffe616-2ece-46f0-9b19-9b2e55dbcb53",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(acf_calc(m_list)[:10000]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "258aee60-8719-48dc-aa04-c69561f88d2c",
   "metadata": {},
   "source": [
    "### Cluster algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0ece6a5-b2d3-42ed-a696-a66008fb000f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import sparse # graph analysis for clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a2ccad5-a929-464e-9813-dd8a82b6b701",
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def spin__make_graph(row, col, spin, neigh_sorted, connect):\n",
    "    \"\"\"Makes a sparse matrix of connected bonds. `spin` is the state, `neigh_sorted` is our\n",
    "    *sorted* neighbor list, `row` and `col` are the outputs, and `connect` is a random matrix\n",
    "    of booleans deciding if we connect a given bond. Returns the number of connections.\"\"\"\n",
    "    Nconnect = 0\n",
    "    counter = 0\n",
    "    for i, s, nn in zip(range(len(spin)), spin, neigh_sorted):\n",
    "        for j in nn:\n",
    "            if i<j and s==spin[j] and connect[counter]:\n",
    "                row[Nconnect] = i\n",
    "                col[Nconnect] = j\n",
    "                Nconnect += 1\n",
    "            counter += 1\n",
    "    return Nconnect\n",
    "\n",
    "neigh_sorted = neigh.copy()\n",
    "for i in range(N):\n",
    "    neigh_sorted[i] = sorted(neigh[i])\n",
    "\n",
    "row = np.zeros(N*z, dtype=int)\n",
    "col = np.zeros(N*z, dtype=int)\n",
    "data = np.ones(N*z, dtype=int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4ca4b62-6626-4657-8298-7b9e0d370fb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cluster algorithm...\n",
    "NclusterMC = 1024\n",
    "E_list, m_list = [], []\n",
    "\n",
    "for _ in trange(NclusterMC):\n",
    "    # make our connections:\n",
    "    Nconnect = spin__make_graph(row, col, spin, neigh_sorted, rng.uniform(size=(N*z)) < (1-np.exp(-2*param['J']/kT)))\n",
    "    # find the clusters:\n",
    "    connectivity = sparse.csr_matrix((data[:Nconnect], (row[:Nconnect], col[:Nconnect])), shape=(N,N))\n",
    "    nsub, labels = sparse.csgraph.connected_components(connectivity, directed=False)\n",
    "    # label[i] = which cluster site i belongs to; we need to find the clusters of size > 1 to flip\n",
    "    label_count = np.zeros(nsub, dtype=int)\n",
    "    for i in labels:\n",
    "        label_count[i] += 1\n",
    "    # now we make a list of which sites we're going to flip; we can either *always* flip\n",
    "    # or *randomly* choose to flip\n",
    "    # label_flip = set(label for label, count in enumerate(label_count) if count > 1)  # always flip\n",
    "    label_flip = set(label for label, count, flip in \n",
    "                     zip(range(len(label_count)), label_count, rng.choice((True, False), size=nsub)) \n",
    "                     if count > 1 and flip)\n",
    "    for i, label in enumerate(labels):\n",
    "        if label in label_flip:\n",
    "            spin[i] = -spin[i]\n",
    "    m_list.append(np.average(spin))\n",
    "    E_list.append(ene(spin, neigh_sorted))\n",
    "m_list = np.array(m_list)\n",
    "E_list = np.array(E_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5055c6f6-d74e-486b-96a1-968ea49d89f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(m_list[:256]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fea9bf7f-b6ca-43a6-8ca3-c760edbfc608",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(acf_calc(m_list)[:256]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53537c68-1cbd-48dc-a959-7ccd809ee234",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_list);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24151303-884e-4967-a469-4e9c90d9bfb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(E_list/N);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b25ed26-c4b9-40af-8435-73f2df0b9d82",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(spin.reshape(Nx, Ny), cmap='grey');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d535fa-d53b-4329-a902-3924a34a5a57",
   "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
}