{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical mechanics of fluids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General exprssion of probability distribution of fluids in phase space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The probability of a state in a classical fluid system is $f(x^N, p^N)$\n",
    "\n",
    "$${f(x^N, p^N) = \\frac{e^{-\\beta H(x^N p^N)}}{Z}}$$\n",
    "\n",
    "- Momentum and position coordinates are separated thanks to the form of kinetic and potential energy terms\n",
    "\n",
    "$$H(x^N, p^N) = K(p^N) + U(x^N)$$\n",
    "\n",
    "$$f(x^N, p^N) = \\Phi(p^N) P(r^N)$$\n",
    "\n",
    "\n",
    "- Kinetic energy part giveus us Maxwell-Botlzman distribution or the ideal gas part of the partition function $Z_p$ \n",
    "- Potential energy part gives us configruational parition function $Z_x$ evaluation of which is non-trivial and mostly can be done via simulaions:\n",
    "\n",
    "$$Z(\\beta, V, N) = Z_{p} \\cdot Z_x  = \\frac{1}{\\lambda^3 N!} Q(\\beta, V, N) $$\n",
    "\n",
    "$$Q = \\int dx^N e^{-\\beta U(x^N)}$$\n",
    "\n",
    "Pressure is related to Free energy as \n",
    "\n",
    "$$p= - \\frac{\\partial F}{\\partial V} = \\frac{\\partial }{\\partial V} log \\int dx^N e^{-\\beta U(x^N)}$$\n",
    "\n",
    "- **Volume dependence of partition function is in the integration limits!** As volume grows, so does partition function. Therefore p is always positive. We can thus conclude that **in equilibrium pressure is always a positive quantity**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reduced configruational distribution functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Key fact:** The full configurational probability $P(r^N)$ and the partition function $Q$ do *not* factorize due to interparticle interactions. Stronger interactions lead to stronger positional correlations.\n",
    "\n",
    "- **Joint probability of particle positions:**  \n",
    "To describe the probability of finding particles at positions $r_1, r_2, \\dots, r_N$:\n",
    "\n",
    "$${P(r^N) = P(r_1, r_2, \\dots, r_N)}$$\n",
    "\n",
    "- **Marginal (reduced) probability density:**  \n",
    "Probability of finding particles 1 and 2 at positions $r_1$ and $r_2$, regardless of others:\n",
    "\n",
    "$${\\rho^{(2)}(r_1, r_2) = \\int dr_3 \\dots dr_N \\, P(r^N)}$$\n",
    "\n",
    "- **Symmetrized reduced pair distribution:**  \n",
    "When particles are indistinguishable, the reduced two-particle distribution becomes:\n",
    "\n",
    "$${\\rho^{(2)}(r_1, r_2) = N(N-1) \\int dr_3 \\dots dr_N \\, P(r^N)}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Radial Distribution Function (RDF)\n",
    "\n",
    "For an **isotropic fluid**, the one-particle probability density is uniform:\n",
    "\n",
    "$$\n",
    "\\rho^{(1)} = N \\frac{\\int dr_2 \\dots dr_N}{\\int dr_1 \\dots dr_N} = N \\frac{V^{N-1}}{V^N} = \\frac{N}{V} = \\rho\n",
    "$$\n",
    "\n",
    "For an **ideal gas**, the two-particle (joint) probability density simplifies as:\n",
    "\n",
    "$$\n",
    "\\rho^{(2)} = N(N-1) \\frac{\\int dr_3 \\dots dr_N}{\\int dr_1 \\dots dr_N} = N(N-1) \\frac{V^{N-2}}{V^N} \\approx \\rho^2\n",
    "$$\n",
    "\n",
    "- To quantify spatial correlations between particles, we define the **Radial Distribution Function (RDF)** as:\n",
    "\n",
    "$$\n",
    "{g(r_1, r_2) = \\frac{\\rho^{(2)}(r_1, r_2)}{\\rho^2}}\n",
    "$$\n",
    "\n",
    "- For an **isotropic and homogeneous fluid**, RDF depends only on the **distance** between particles:\n",
    "\n",
    "$$\n",
    "{g(r) = g(|r_2 - r_1|)}\n",
    "$$\n",
    "\n",
    "\n",
    "### **Physical meaning of RDF**\n",
    "\n",
    "- Since $\\rho = \\rho^{(1)}$, the **conditional probability density** of finding a particle at distance $r$ from a tagged particle at the origin is:\n",
    "\n",
    "$$\n",
    "\\frac{\\rho^{(2)}(0, r)}{\\rho} = \\rho g(r)\n",
    "$$\n",
    "\n",
    "- Thus, $\\rho g(r)$ represents the **average local density** at distance $r$, given a particle is fixed at the origin.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# === System parameters ===\n",
    "N_side = 20                # particles per side for lattice\n",
    "N = N_side**2              # total particles\n",
    "L = 10.0                   # box length (L x L)\n",
    "r_max = L / 2              # max distance for RDF\n",
    "nbins = 100\n",
    "\n",
    "# === Make ideal gas configuration ===\n",
    "positions_ideal = np.random.uniform(0, L, size=(N, 2))\n",
    "\n",
    "# === Make square lattice configuration ===\n",
    "spacing = L / N_side\n",
    "x, y = np.meshgrid(np.linspace(0, L - spacing, N_side),\n",
    "                   np.linspace(0, L - spacing, N_side))\n",
    "positions_lattice = np.vstack([x.ravel(), y.ravel()]).T\n",
    "\n",
    "def compute_rdf(positions, L, r_max, nbins):\n",
    "    N = len(positions)\n",
    "    dists = []\n",
    "    for i in range(N):\n",
    "        for j in range(i + 1, N):\n",
    "            dx = positions[i, 0] - positions[j, 0]\n",
    "            dy = positions[i, 1] - positions[j, 1]\n",
    "            \n",
    "            # Periodic boundary conditions\n",
    "            dx -= L * np.round(dx / L)\n",
    "            dy -= L * np.round(dy / L)\n",
    "            r = np.sqrt(dx**2 + dy**2)\n",
    "            if r < r_max:\n",
    "                dists.append(r)\n",
    "    \n",
    "    dists = np.array(dists)\n",
    "    r_bins = np.linspace(0.0, r_max, nbins + 1)\n",
    "    r_centers = 0.5 * (r_bins[:-1] + r_bins[1:])\n",
    "    hist, _ = np.histogram(dists, bins=r_bins)\n",
    "    \n",
    "    dr = r_bins[1] - r_bins[0]\n",
    "    area = L**2\n",
    "    rho = N / area\n",
    "    shell_area = 2 * np.pi * r_centers * dr\n",
    "    norm = rho * shell_area * (N - 1) / 2  # expected number in each shell\n",
    "    g_r = hist / norm\n",
    "    \n",
    "    return r_centers, g_r\n",
    "\n",
    "# === Compute RDFs ===\n",
    "r_ideal, g_ideal = compute_rdf(positions_ideal, L, r_max, nbins)\n",
    "r_lattice, g_lattice = compute_rdf(positions_lattice, L, r_max, nbins)\n",
    "\n",
    "# === Plot ===\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot(r_ideal, g_ideal, label='Ideal Gas', lw=2)\n",
    "plt.plot(r_lattice, g_lattice, label='Square Lattice', lw=2)\n",
    "plt.axhline(1.0, color='gray', linestyle='--')\n",
    "plt.xlabel('Distance $r$')\n",
    "plt.ylabel('$g(r)$')\n",
    "plt.title('Radial Distribution Function Comparison')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coordination shells and structure in fluids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./figs/gr1_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./figs/gr2_2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./figs/gr3_3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./figs/gr4_4.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reversible work theorem and potential of mean force"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{admonition} **Reversible work theorem**\n",
    ":class: important\n",
    "\n",
    "$${g(r) = e^{-\\beta w(r)}}$$\n",
    "\n",
    "$${w(r) = - \\beta^{-1} log [g(r)]}$$\n",
    "\n",
    "- $w(r)$ Reversible work to bring two particles from infinity to distance r\n",
    "\n",
    "- $g(r)$ Radial distribution function\n",
    "\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\Big \\langle -\\frac{d  U(r^N)}{dr_1}  \\Big \\rangle_{r_1, r_2} = -\\frac{\\int dr_3...dr_N \\frac{dU(r^N)}{dr_1} e^{-\\beta U}}{\\int dr_3...r_N e^{-\\beta U}} =  \\frac{\\beta^{-1} \\frac{d}{d r_1} \\int dr_3...dr_N e^{-\\beta U}}{\\int dr_3...dr_N e^{-\\beta U}} = \\beta^{-1} \\frac{d}{d r_1} log \\int dr_3...dr_N e^{-\\beta U} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\Big \\langle -\\frac{d  U(r^N)}{dr_1}  \\Big \\rangle_{r_1, r_2} = \\beta^{-1} \\frac{d}{d r_1} log N (N-1)\\int dr_3...dr_N e^{-\\beta U} = \\beta^{-1} g(r_1, r_2)$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thermodynamic properites of $g(r)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\langle E \\rangle = N \\Big\\langle \\frac{p^2}{2m} \\Big\\rangle + \\Big\\langle \\sum_{j>i} u(|r_i - r_j|)\\Big \\rangle$$\n",
    "\n",
    "$${E/N  = \\frac{3}{2}k_B T +\\frac{1}{2}\\rho \\int dr g(r) u(r) }$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Low density approximation for $g(r)$\n",
    "\n",
    "$$w(r) =u(r) +\\Delta w(r) $$\n",
    "\n",
    "$$ g(r) = e^{-\\beta u(r)} \\Big (1 +O(\\rho) \\Big)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Density expanesion and virial coefficients\n",
    "\n",
    "$$\\beta p = \\rho + B_2(T) \\rho^2+ O(\\rho^3)$$\n",
    "\n",
    "$$B_2(T) = -\\frac{1}{2} \\int dr (e^{-\\beta u(r)}-1) $$"
   ]
  }
 ],
 "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.9.16"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}