{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bbec08f5-2ccb-4e82-8184-afb1e9e507ab",
   "metadata": {
    "id": "bbec08f5-2ccb-4e82-8184-afb1e9e507ab"
   },
   "source": [
    "# 2D hard disk Monte Carlo\n",
    "\n",
    "MC code, similar to Metropolis *et al.* to simulate hard disks in 2D."
   ]
  },
  {
   "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": [
    "rc = 0.5\n",
    "L, Nx, Ny = 18, 14, 16\n",
    "N = Nx*Ny\n",
    "V = L*L/(N*np.pi*rc**2)\n",
    "rho = 1.0/V # fill fraction\n",
    "rho, V, 0.5*np.sqrt(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7777ff0f-9bf7-4cd8-8b7c-9215e5e45fcb",
   "metadata": {
    "id": "7777ff0f-9bf7-4cd8-8b7c-9215e5e45fcb"
   },
   "outputs": [],
   "source": [
    "# initialize our particles on a square lattice, staggered\n",
    "r_i = np.zeros((N, 2))\n",
    "n = 0\n",
    "for i in range(Nx):\n",
    "    for j in range(Ny):\n",
    "        r_i[n] = (i+0.5*(j%2))*(L/Nx), j*(L/Ny)\n",
    "        n += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb985a2e-64c4-48a6-ba41-12636bb6de28",
   "metadata": {
    "id": "fb985a2e-64c4-48a6-ba41-12636bb6de28"
   },
   "outputs": [],
   "source": [
    "plt.plot(r_i[:,0], r_i[:,1], 'o');"
   ]
  },
  {
   "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 MC(pos, rc, i_rand, delta, L):\n",
    "    \"\"\"Run MC for hard disks where we attempt to move atom i_rand by delta.\n",
    "    Returns acceptance ratio, and overwrites pos in place. PBC with\n",
    "    sidelength L; uses N^2 method for checking distance\"\"\"\n",
    "    N = len(pos)\n",
    "    accept = 0\n",
    "    rc2 = 4*rc*rc\n",
    "    L2 = L*0.5\n",
    "    for n in range(len(i_rand)):\n",
    "        i = i_rand[n]\n",
    "        rp = pos[i] + delta[n]\n",
    "        collide = False\n",
    "        for j in range(N):\n",
    "            if j == i: continue\n",
    "            r_ij = pos[j] - rp\n",
    "            if (r_ij[0] < -L2): r_ij[0] += L\n",
    "            if (r_ij[0] >  L2): r_ij[0] -= L\n",
    "            if (r_ij[1] < -L2): r_ij[1] += L\n",
    "            if (r_ij[1] >  L2): r_ij[1] -= L\n",
    "            if (r_ij[0]*r_ij[0] + r_ij[1]*r_ij[1]) < rc2:\n",
    "                collide = True\n",
    "                break\n",
    "        if not collide:\n",
    "            accept += 1\n",
    "            if (rp[0] < 0): rp[0] += L\n",
    "            if (rp[0] > L): rp[0] -= L\n",
    "            if (rp[1] < 0): rp[1] += L\n",
    "            if (rp[1] > L): rp[1] -= L\n",
    "            pos[i,0] = rp[0]\n",
    "            pos[i,1] = rp[1]\n",
    "    return accept/len(i_rand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af05a73c-3e06-4ba3-a828-301a1fc2fd3f",
   "metadata": {
    "id": "af05a73c-3e06-4ba3-a828-301a1fc2fd3f"
   },
   "outputs": [],
   "source": [
    "N_MC = 1024\n",
    "Delta = 0.5*np.sqrt(V)-rc\n",
    "i_rand = rng.choice(N, size=N_MC)\n",
    "delta = rng.uniform(-Delta, Delta, size=(N_MC, 2))\n",
    "accept = MC(r_i, rc, i_rand, delta, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a316999-e3e7-4177-9ffb-0656851fe331",
   "metadata": {
    "id": "6a316999-e3e7-4177-9ffb-0656851fe331"
   },
   "outputs": [],
   "source": [
    "plt.plot(r_i[:,0], r_i[:,1], 'o');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "636aecd5-e3b9-4743-b4bd-409708ffab8b",
   "metadata": {
    "id": "636aecd5-e3b9-4743-b4bd-409708ffab8b"
   },
   "outputs": [],
   "source": [
    "accept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8665000-65dc-49c4-891b-56abcac4ff8d",
   "metadata": {
    "id": "c8665000-65dc-49c4-891b-56abcac4ff8d"
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def g_r(pos, L, g_bins, r2_cut):\n",
    "    \"\"\"Take in positions, add to histograms\"\"\"\n",
    "    N = len(pos)\n",
    "    Nr = len(r2_cut)\n",
    "    L2 = L/2\n",
    "    for i in range(N):\n",
    "        for j in range(i):\n",
    "            r_ij = pos[i] - pos[j]\n",
    "            if (r_ij[0] < -L2): r_ij[0] += L\n",
    "            if (r_ij[0] >  L2): r_ij[0] -= L\n",
    "            if (r_ij[1] < -L2): r_ij[1] += L\n",
    "            if (r_ij[1] >  L2): r_ij[1] -= L\n",
    "            r2 = (r_ij[0]*r_ij[0] + r_ij[1]*r_ij[1])\n",
    "            for n in range(1,Nr):\n",
    "                if r2_cut[n-1] < r2 < r2_cut[n]:\n",
    "                    g_bins[n] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49963e55-258a-4202-8c7a-9e97fb9782f8",
   "metadata": {
    "id": "49963e55-258a-4202-8c7a-9e97fb9782f8"
   },
   "outputs": [],
   "source": [
    "rhist = np.linspace(0,L/2,num=100)\n",
    "r2_cut = rhist*rhist\n",
    "g_bins = np.zeros_like(rhist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10b4a857-4050-4e6b-98de-76f5a6277051",
   "metadata": {
    "id": "10b4a857-4050-4e6b-98de-76f5a6277051"
   },
   "outputs": [],
   "source": [
    "g_r(r_i, L, g_bins, r2_cut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c73f725-d025-4bcb-8493-b10b1381171b",
   "metadata": {
    "id": "7c73f725-d025-4bcb-8493-b10b1381171b"
   },
   "outputs": [],
   "source": [
    "g_bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cef74136-7072-4f44-a5ce-437d8b3c40ed",
   "metadata": {
    "id": "cef74136-7072-4f44-a5ce-437d8b3c40ed"
   },
   "outputs": [],
   "source": [
    "plt.plot(g_bins);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b3ac513-e53c-440d-a2c9-63466d39a073",
   "metadata": {
    "id": "8b3ac513-e53c-440d-a2c9-63466d39a073"
   },
   "outputs": [],
   "source": [
    "N_sweeps = 1024\n",
    "N_MC = 1024\n",
    "Delta = 0.5*np.sqrt(V)-rc\n",
    "\n",
    "rhist = np.linspace(0,L/2,num=100)\n",
    "r2_cut = rhist*rhist\n",
    "g_bins = np.zeros_like(rhist)\n",
    "\n",
    "accept_ave = 0\n",
    "for _ in trange(N_sweeps):\n",
    "    i_rand = rng.choice(N, size=N_MC)\n",
    "    delta = rng.uniform(-Delta, Delta, size=(N_MC, 2))\n",
    "    accept_ave += MC(r_i, rc, i_rand, delta, L)\n",
    "    g_r(r_i, L, g_bins, r2_cut)\n",
    "g_bins /= N_sweeps\n",
    "accept_ave /= N_sweeps\n",
    "accept_ave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff8a8312-d90e-410a-9c91-d75802d1afd7",
   "metadata": {
    "id": "ff8a8312-d90e-410a-9c91-d75802d1afd7"
   },
   "outputs": [],
   "source": [
    "# normalize g_bins:\n",
    "for n in range(len(g_bins)-1):\n",
    "    g_bins[n+1] /= 0.5*N*(N-1)*np.pi*(r2_cut[n+1]-r2_cut[n])/(L*L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42897cf5-99ec-498f-a5e6-fbe1fb080852",
   "metadata": {
    "id": "42897cf5-99ec-498f-a5e6-fbe1fb080852"
   },
   "outputs": [],
   "source": [
    "plt.plot(rhist, g_bins, '.-');\n",
    "plt.axhline(1.0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c47230d-64df-4afe-83ba-6cd599f894ed",
   "metadata": {
    "id": "3c47230d-64df-4afe-83ba-6cd599f894ed"
   },
   "outputs": [],
   "source": [
    "plt.plot(r_i[:,0], r_i[:,1],'o');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7L80TZkagydK",
   "metadata": {
    "id": "7L80TZkagydK"
   },
   "source": [
    "Now, let's scan through some different different supercell sizes to see how $g(r)$ changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "EwJ-ZtiNgx5-",
   "metadata": {
    "id": "EwJ-ZtiNgx5-"
   },
   "outputs": [],
   "source": [
    "g_r_dict = {18.0: g_bins.copy()}\n",
    "L_range = (17.5, 17.0, 16.5, 16.0, 15.5, 15.0, 14.5)\n",
    "L_old = L\n",
    "for L in L_range:\n",
    "  # scale positions down to fit in box\n",
    "  r_i *= L/L_old\n",
    "  Delta *= L/L_old # do smaller steps too\n",
    "  g_bins = np.zeros_like(rhist)\n",
    "\n",
    "  accept_ave = 0\n",
    "  for _ in trange(N_sweeps):\n",
    "      i_rand = rng.choice(N, size=N_MC)\n",
    "      delta = rng.uniform(-Delta, Delta, size=(N_MC, 2))\n",
    "      accept_ave += MC(r_i, rc, i_rand, delta, L)\n",
    "      g_r(r_i, L, g_bins, r2_cut)\n",
    "  g_bins /= N_sweeps\n",
    "  accept_ave /= N_sweeps\n",
    "  # normalize g_bins:\n",
    "  for n in range(len(g_bins)-1):\n",
    "      g_bins[n+1] /= 0.5*N*(N-1)*np.pi*(r2_cut[n+1]-r2_cut[n])/(L*L)\n",
    "  g_r_dict[L] = g_bins.copy()\n",
    "  L_old = L\n",
    "  print(L, accept_ave)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "NAy2wOlybwmm",
   "metadata": {
    "id": "NAy2wOlybwmm"
   },
   "outputs": [],
   "source": [
    "for L in L_range[:]:\n",
    "  plt.plot(rhist, g_r_dict[L], label=f\"L={L}\");\n",
    "plt.legend();\n",
    "plt.axhline(1.0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "xMcfuTzxdBHY",
   "metadata": {
    "id": "xMcfuTzxdBHY"
   },
   "outputs": [],
   "source": [
    "plt.plot(r_i[:,0], r_i[:,1], 'o');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "TiLmWr31kbjb",
   "metadata": {
    "id": "TiLmWr31kbjb"
   },
   "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
}