From 1905d853f7b2e967a74ee95bcb94c7d53db33308 Mon Sep 17 00:00:00 2001 From: Anders Hafreager Date: Wed, 8 Jul 2026 21:45:39 +0200 Subject: [PATCH] Add md-basics tutorial series; test the DriveFS mount; fix truncation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Six new notebooks under md-basics/ — statistical mechanics you can watch, each validated headlessly in a real browser: - 01-maxwell-boltzmann: every atom starts at the same speed; the M-B distribution assembles itself (WCA gas, arcsine → Gaussian vx). - 02-condensation: quench a 2D vapor below coexistence; droplets nucleate and coarsen, colored by pe/atom. Optional KOKKOS "go big" cell with the basics/04 isolation check + fallback. - 03-real-vs-ideal-gas: measure Z = P/ρT across (ρ, T), overlay the second virial coefficient computed in numpy from the pair potential. - 04-reversibility-and-chaos: Loschmidt's paradox live — flip every velocity and un-melt a crystal (RMS return error ~1e-12), then find the chaos horizon where the return fails (positions matched by atom id, since LAMMPS re-sorts atoms in memory). - 05-thermal-conductivity: Müller-Plathe flux + ave/chunk temperature profile read from a file → κ ≈ 6 for the 2D LJ liquid. - 06-polymer-chains: numpy-generated FENE bead-spring data files, Rg ~ N^ν with the Flory exponent from a log-log fit. tests/drivefs-mount.spec.ts unit-tests the DriveFS mount snippet from the Python bindings against the real wasm filesystem with a mocked JupyterLite ContentsAPI — every reachable branch of toAPIPath (all four cwd cases, re-evaluated per call), node_ops, stream_ops, and mounting, plus LAMMPS end-to-end cases (write_dump, ave/time + close, ave/chunk + unfix, read_data of Python-written files, re-run overwrite). The tests exposed a real bug, now fixed: truncation arrives as a node-level setattr size change (Emscripten strips O_TRUNC before stream open), which the proxy ignored — so re-opening an existing output file with "w" kept the stale tail of the previous, longer content. setattr now honors attr.size (shrink, zero, and extend). Co-Authored-By: Claude Fable 5 --- NOTEBOOK_TUTORIALS.md | 15 +- examples/notebook/content/index.ipynb | 14 + .../md-basics/01-maxwell-boltzmann.ipynb | 193 ++++++++ .../content/md-basics/02-condensation.ipynb | 220 ++++++++ .../md-basics/03-real-vs-ideal-gas.ipynb | 153 ++++++ .../04-reversibility-and-chaos.ipynb | 200 ++++++++ .../md-basics/05-thermal-conductivity.ipynb | 158 ++++++ .../content/md-basics/06-polymer-chains.ipynb | 183 +++++++ python/lammps/__init__.py | 13 + tests/drivefs-mount.spec.ts | 468 ++++++++++++++++++ 10 files changed, 1616 insertions(+), 1 deletion(-) create mode 100644 examples/notebook/content/md-basics/01-maxwell-boltzmann.ipynb create mode 100644 examples/notebook/content/md-basics/02-condensation.ipynb create mode 100644 examples/notebook/content/md-basics/03-real-vs-ideal-gas.ipynb create mode 100644 examples/notebook/content/md-basics/04-reversibility-and-chaos.ipynb create mode 100644 examples/notebook/content/md-basics/05-thermal-conductivity.ipynb create mode 100644 examples/notebook/content/md-basics/06-polymer-chains.ipynb create mode 100644 tests/drivefs-mount.spec.ts diff --git a/NOTEBOOK_TUTORIALS.md b/NOTEBOOK_TUTORIALS.md index 2ba67a2..9f98a7f 100644 --- a/NOTEBOOK_TUTORIALS.md +++ b/NOTEBOOK_TUTORIALS.md @@ -90,7 +90,13 @@ examples/notebook/content/ │ ├── 04-grain-boundary.ipynb │ ├── 05-fracture.ipynb │ └── 06-nanoindentation.ipynb -├── md-basics/ # molecular dynamics concepts (LJ systems) +├── md-basics/ # SHIPPED: statistical mechanics you can watch +│ ├── 01-maxwell-boltzmann.ipynb +│ ├── 02-condensation.ipynb +│ ├── 03-real-vs-ideal-gas.ipynb +│ ├── 04-reversibility-and-chaos.ipynb +│ ├── 05-thermal-conductivity.ipynb +│ └── 06-polymer-chains.ipynb ├── materials/ # matsci-tutorials.pdf ports (needs MANYBODY) ├── soft-matter/ # LiveCoMS soft-matter ports (needs KSPACE/MOLECULE) ├── showcase/ # atomify example ports @@ -112,6 +118,13 @@ examples/notebook/content/ ### Series 1 — MD fundamentals with LJ systems (default wasm ✓) +Partially SHIPPED as `md-basics/`: Maxwell–Boltzmann relaxation from a +delta-function start, vapor→droplet condensation (with an optional KOKKOS +"go big" cell), real-vs-ideal gas with a numpy $B_2$ virial overlay, +Loschmidt velocity-reversal + chaos horizon, Müller-Plathe thermal +conductivity, and FENE bead-spring polymers with the Flory exponent +(data files generated by numpy). Remaining ideas below. + 1. **Melting an fcc crystal** — energy/temperature vs time, spot the phase transition. 2. **Radial distribution function** — `compute rdf`, plot g(r) for solid vs liquid. 3. **Diffusion & mean-squared displacement** — Einstein relation, extract D from MSD slope (atomify has a `diffusion` example to port). diff --git a/examples/notebook/content/index.ipynb b/examples/notebook/content/index.ipynb index 58d2d7d..dc968c6 100644 --- a/examples/notebook/content/index.ipynb +++ b/examples/notebook/content/index.ipynb @@ -95,6 +95,20 @@ "| [basics/04-multithreading-kokkos](basics/04-multithreading-kokkos.ipynb) | The multithreaded KOKKOS build |\n", "| [basics/05-thermo-output-to-file](basics/05-thermo-output-to-file.ipynb) | The classic workflow: `fix ave/time` writes thermo to a file, a separate cell plots it |\n", "\n", + "## Molecular dynamics fundamentals\n", + "\n", + "Statistical mechanics you can watch: each notebook takes one famous piece\n", + "of physics and makes it happen live, in a system small enough to see.\n", + "\n", + "| Notebook | The physics |\n", + "|---|---|\n", + "| [md-basics/01-maxwell-boltzmann](md-basics/01-maxwell-boltzmann.ipynb) | Start every atom at the same speed, watch the Maxwell–Boltzmann distribution assemble itself |\n", + "| [md-basics/02-condensation](md-basics/02-condensation.ipynb) | Quench a vapor below its condensation point: droplets nucleate and coarsen |\n", + "| [md-basics/03-real-vs-ideal-gas](md-basics/03-real-vs-ideal-gas.ipynb) | Measure Z = P/ρT, then predict it from the pair potential with a numpy virial integral |\n", + "| [md-basics/04-reversibility-and-chaos](md-basics/04-reversibility-and-chaos.ipynb) | Flip every velocity and un-melt a crystal — then find the chaos horizon where it fails |\n", + "| [md-basics/05-thermal-conductivity](md-basics/05-thermal-conductivity.ipynb) | A real transport measurement: Müller-Plathe flux, temperature gradient, κ |\n", + "| [md-basics/06-polymer-chains](md-basics/06-polymer-chains.ipynb) | Build bead-spring polymers with numpy, measure the Flory exponent |\n", + "\n", "## Materials science: from a perfect crystal to failure\n", "\n", "A six-part course following the LiveCoMS\n", diff --git a/examples/notebook/content/md-basics/01-maxwell-boltzmann.ipynb b/examples/notebook/content/md-basics/01-maxwell-boltzmann.ipynb new file mode 100644 index 0000000..92c336e --- /dev/null +++ b/examples/notebook/content/md-basics/01-maxwell-boltzmann.ipynb @@ -0,0 +1,193 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 01 — The Maxwell–Boltzmann distribution, live\n", + "\n", + "Give every atom in a gas *exactly the same speed*, pointing in random\n", + "directions — about as far from thermal equilibrium as a velocity\n", + "distribution gets. Then let them collide. Within a few hundred collisions\n", + "the most famous distribution in statistical mechanics assembles itself in\n", + "front of you: velocity components become Gaussian, speeds become\n", + "Maxwell–Boltzmann. No thermostat, no randomness added — just Newton's laws\n", + "(this is Boltzmann's H-theorem happening in real time).\n", + "\n", + "We use a 2D gas with purely repulsive interactions (the WCA potential — a\n", + "Lennard-Jones potential cut at its minimum), so collisions redistribute\n", + "energy without any of it hiding as potential energy." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install lammps-js matplotlib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A gas where every atom has the same speed\n", + "\n", + "`velocity create` draws random directions; two atom-style variables then\n", + "rescale every velocity to the same magnitude $v_0 = 1.5$. In LJ units\n", + "($m = k_B = 1$) the temperature in 2D is $T = \\langle v^2 \rangle / 2$,\n", + "so it starts at $v_0^2/2 = 1.125$ and stays within a few percent of it\n", + "(during a collision a little energy briefly sits in the repulsive\n", + "potential) while the distribution's shape changes completely." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from lammps import lammps, LMP_VAR_ATOM\n", + "\n", + "lmp = await lammps(output=None)\n", + "lmp.commands_string(\"\"\"\n", + "units lj\n", + "dimension 2\n", + "lattice sq 0.3\n", + "region box block 0 32 0 32 -0.1 0.1\n", + "create_box 1 box\n", + "create_atoms 1 box\n", + "mass 1 1.0\n", + "pair_style lj/cut 1.122462\n", + "pair_coeff 1 1 1.0 1.0\n", + "pair_modify shift yes\n", + "\n", + "velocity all create 1.0 777 dist gaussian\n", + "variable v0 equal 1.5\n", + "variable s atom v_v0/sqrt(vx*vx+vy*vy)\n", + "variable ux atom vx*v_s\n", + "variable uy atom vy*v_s\n", + "velocity all set v_ux v_uy NULL\n", + "\n", + "variable vxa atom vx\n", + "variable spd atom sqrt(vx*vx+vy*vy)\n", + "\n", + "fix 1 all nve\n", + "fix 2d all enforce2d\n", + "run 0\n", + "\"\"\")\n", + "print(lmp.get_natoms(), \"atoms, all at speed 1.5, T =\", round(lmp.get_thermo(\"temp\"), 4))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Watch equilibrium assemble\n", + "\n", + "Sample the velocity components at a few moments. At $t=0$ every speed is\n", + "$v_0$, so the distribution of $v_x = v_0\\cos\\theta$ has spikes at\n", + "$\\pm v_0$ (an arcsine distribution — most of a circle is near its\n", + "edges). Then collisions take over:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "snapshots = {}\n", + "for label, steps in [(\"t = 0\", 0), (\"after 100 steps\", 100),\n", + " (\"after 400 steps\", 300), (\"after 4000 steps\", 3600)]:\n", + " if steps:\n", + " lmp.command(f\"run {steps}\")\n", + " snapshots[label] = (lmp.extract_variable(\"vxa\", vartype=LMP_VAR_ATOM).copy(),\n", + " lmp.extract_variable(\"spd\", vartype=LMP_VAR_ATOM).copy())\n", + "T = lmp.get_thermo(\"temp\")\n", + "print(\"temperature after equilibration:\", round(T, 4), \"(started at v0²/2 = 1.125)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(2, 2, figsize=(9, 6), sharex=True)\n", + "v = np.linspace(-3.5, 3.5, 300)\n", + "gauss = np.exp(-v**2 / (2 * T)) / np.sqrt(2 * np.pi * T)\n", + "for ax, (label, (vx, _)) in zip(axes.flat, snapshots.items()):\n", + " ax.hist(vx, bins=45, range=(-3.5, 3.5), density=True, alpha=0.7)\n", + " ax.plot(v, gauss, \"k--\", lw=1.2, label=\"Gaussian, measured T\")\n", + " ax.set_title(label)\n", + " ax.legend(fontsize=8)\n", + "for ax in axes[1]:\n", + " ax.set_xlabel(\"$v_x$\")\n", + "for ax in axes[:, 0]:\n", + " ax.set_ylabel(\"probability density\")\n", + "fig.suptitle(\"A Gaussian velocity distribution assembles itself\")\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Speeds: the 2D Maxwell–Boltzmann distribution\n", + "\n", + "Same data, but now the speed $|v|$. It starts as a delta function at\n", + "$v_0$ and relaxes to the 2D Maxwell–Boltzmann form\n", + "$f(v) = (v/T)\\,e^{-v^2/2T}$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "v = np.linspace(0, 4, 300)\n", + "mb2d = (v / T) * np.exp(-v**2 / (2 * T))\n", + "\n", + "plt.figure(figsize=(6, 3.6))\n", + "plt.hist(snapshots[\"t = 0\"][1], bins=45, range=(0, 4), density=True,\n", + " alpha=0.55, label=\"t = 0 (all speeds = 1.5)\")\n", + "plt.hist(snapshots[\"after 4000 steps\"][1], bins=45, range=(0, 4), density=True,\n", + " alpha=0.55, label=\"after 4000 steps\")\n", + "plt.plot(v, mb2d, \"k--\", lw=1.2, label=\"2D Maxwell–Boltzmann\")\n", + "plt.xlabel(\"speed $|v|$\"); plt.ylabel(\"probability density\")\n", + "plt.legend(); plt.tight_layout(); plt.show()\n", + "\n", + "lmp.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Every collision conserved energy and momentum exactly, yet the initial\n", + "condition was forgotten completely — only its total energy survives, as the\n", + "temperature setting the width of the Gaussian.\n", + "\n", + "Next: [02 — Condensation](02-condensation.ipynb), where the same gas is\n", + "quenched below its condensation point and droplets rain out." + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python", + "display_name": "Python (Pyodide)" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/notebook/content/md-basics/02-condensation.ipynb b/examples/notebook/content/md-basics/02-condensation.ipynb new file mode 100644 index 0000000..48f45c0 --- /dev/null +++ b/examples/notebook/content/md-basics/02-condensation.ipynb @@ -0,0 +1,220 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 02 — Condensation: quench a vapor, watch droplets form\n", + "\n", + "Take a hot 2D Lennard-Jones vapor and suddenly cool it below its\n", + "condensation point. The uniform gas becomes unstable: density fluctuations\n", + "grow, droplets nucleate everywhere, then bigger droplets eat smaller ones\n", + "(coarsening / Ostwald ripening). This is a genuine first-order phase\n", + "transition happening in your browser tab." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install lammps-js matplotlib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A hot, uniform vapor\n", + "\n", + "Density $\\rho = 0.3$, temperature $T = 1.0$ — comfortably above the 2D LJ\n", + "critical temperature ($T_c \\approx 0.5$), so the gas stays uniform. A weak\n", + "Langevin thermostat plays the role of the environment that carries heat\n", + "away. The `pe/atom` compute lets us color each atom by how tightly it is\n", + "bound — gas atoms ≈ 0, atoms inside a droplet strongly negative:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from lammps import lammps, LMP_STYLE_ATOM, LMP_TYPE_VECTOR\n", + "\n", + "lmp = await lammps(output=None)\n", + "lmp.commands_string(\"\"\"\n", + "units lj\n", + "dimension 2\n", + "lattice sq 0.3\n", + "region box block 0 40 0 40 -0.1 0.1\n", + "create_box 1 box\n", + "create_atoms 1 box\n", + "mass 1 1.0\n", + "velocity all create 1.0 8712 dist gaussian\n", + "pair_style lj/cut 2.5\n", + "pair_coeff 1 1 1.0 1.0 2.5\n", + "compute peatom all pe/atom\n", + "fix 1 all nve\n", + "fix lang all langevin 1.0 1.0 1.0 3216\n", + "fix 2d all enforce2d\n", + "thermo 1000\n", + "run 3000\n", + "\"\"\")\n", + "\n", + "def snapshot():\n", + " xy = lmp.extract_atom(\"x\")[:, :2].copy()\n", + " pe = lmp.extract_compute(\"peatom\", LMP_STYLE_ATOM, LMP_TYPE_VECTOR).copy()\n", + " return xy, pe\n", + "\n", + "frames = [(\"hot vapor, T = 1.0\", *snapshot())]\n", + "print(lmp.get_natoms(), \"atoms equilibrated\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Quench\n", + "\n", + "Drop the thermostat target to $T = 0.4$ and keep integrating, saving\n", + "snapshots as the system condenses:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lmp.command(\"fix lang all langevin 0.4 0.4 1.0 3216\")\n", + "pe_trace, steps_trace = [], []\n", + "for label, steps in [(\"just after the quench\", 2000),\n", + " (\"droplets nucleate\", 8000),\n", + " (\"droplets coarsen\", 20000)]:\n", + " for _ in range(steps // 1000):\n", + " lmp.command(\"run 1000\")\n", + " steps_trace.append(lmp.extract_global(\"ntimestep\"))\n", + " pe_trace.append(lmp.get_thermo(\"pe\"))\n", + " frames.append((label, *snapshot()))\n", + "print(\"final potential energy per atom:\", round(pe_trace[-1], 3))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "L = 40 / np.sqrt(0.3) # box edge in LJ units (lattice sq 0.3, 40 cells)\n", + "fig, axes = plt.subplots(2, 2, figsize=(9, 9))\n", + "for ax, (label, xy, pe) in zip(axes.flat, frames):\n", + " sc = ax.scatter(xy[:, 0], xy[:, 1], c=pe, s=4, cmap=\"viridis_r\",\n", + " vmin=-3.5, vmax=0.5)\n", + " ax.set_title(label, fontsize=10)\n", + " ax.set_aspect(\"equal\"); ax.set_xticks([]); ax.set_yticks([])\n", + "fig.colorbar(sc, ax=axes, shrink=0.8, label=\"potential energy / atom\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Dark atoms are bound inside droplets; pale ones are still vapor. The\n", + "system's potential energy tracks the transition — flat while the vapor is\n", + "stable, then falling as bound pairs form:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(6, 3.2))\n", + "plt.plot(steps_trace, pe_trace)\n", + "plt.xlabel(\"timestep since quench\"); plt.ylabel(\"potential energy / atom\")\n", + "plt.title(\"Condensation releases latent heat\")\n", + "plt.tight_layout(); plt.show()\n", + "\n", + "lmp.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Go bigger (optional)\n", + "\n", + "More atoms means more droplets and better coarsening statistics. On a\n", + "cross-origin-isolated host this cell uses the multithreaded KOKKOS build\n", + "(see [basics/04](../basics/04-multithreading-kokkos.ipynb)); on the public\n", + "site it falls back to the single-threaded engine with a smaller box — still\n", + "entirely CPU-bound WebAssembly, so be patient with big boxes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import js\n", + "from lammps import lammps as lammps2\n", + "\n", + "isolated = bool(getattr(js, \"crossOriginIsolated\", False))\n", + "n = 100 if isolated else 60\n", + "args = [\"-k\", \"on\", \"t\", \"4\", \"-sf\", \"kk\"] if isolated else None\n", + "big = await lammps2(cmdargs=args, output=None)\n", + "big.commands_string(f\"\"\"\n", + "units lj\n", + "dimension 2\n", + "lattice sq 0.3\n", + "region box block 0 {n} 0 {n} -0.1 0.1\n", + "create_box 1 box\n", + "create_atoms 1 box\n", + "mass 1 1.0\n", + "velocity all create 1.0 8712 dist gaussian\n", + "pair_style lj/cut 2.5\n", + "pair_coeff 1 1 1.0 1.0 2.5\n", + "fix 1 all nve\n", + "fix lang all langevin 0.4 0.4 1.0 3216\n", + "fix 2d all enforce2d\n", + "thermo 2000\n", + "run 12000\n", + "\"\"\")\n", + "xy = big.extract_atom(\"x\")[:, :2]\n", + "plt.figure(figsize=(6, 6))\n", + "plt.scatter(xy[:, 0], xy[:, 1], s=2)\n", + "plt.gca().set_aspect(\"equal\"); plt.xticks([]); plt.yticks([])\n", + "plt.title(f\"{big.get_natoms()} atoms, \"\n", + " f\"{'KOKKOS ×4 threads' if isolated else 'single-threaded'}\")\n", + "plt.tight_layout(); plt.show()\n", + "big.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next: [03 — Real vs ideal gas](03-real-vs-ideal-gas.ipynb) — putting a\n", + "number on how non-ideal this gas is." + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python", + "display_name": "Python (Pyodide)" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/notebook/content/md-basics/03-real-vs-ideal-gas.ipynb b/examples/notebook/content/md-basics/03-real-vs-ideal-gas.ipynb new file mode 100644 index 0000000..987bbff --- /dev/null +++ b/examples/notebook/content/md-basics/03-real-vs-ideal-gas.ipynb @@ -0,0 +1,153 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 03 — How ideal is a real gas?\n", + "\n", + "The ideal gas law says $PV = Nk_BT$, i.e. the *compressibility factor*\n", + "\n", + "$$Z = \\frac{P}{\\rho\\,k_B T} = 1$$\n", + "\n", + "for every density and temperature. A Lennard-Jones gas disagrees: at low\n", + "temperature attraction pulls the pressure **down** ($Z < 1$), at high\n", + "density the atoms' own volume pushes it **up** ($Z > 1$). We measure\n", + "$Z(\\rho, T)$ directly, then predict the low-density behavior from the pair\n", + "potential alone with the second virial coefficient — computed in numpy, no\n", + "simulation required." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install lammps-js matplotlib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Measure Z across densities and temperatures\n", + "\n", + "One short NVT run per state point (500 atoms, 3D). Each takes a couple of\n", + "seconds — the loop prints as it goes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from lammps import lammps\n", + "\n", + "async def measure_Z(rho, T):\n", + " L = (500 / rho) ** (1 / 3)\n", + " lmp = await lammps(output=None)\n", + " lmp.commands_string(f\"\"\"\n", + "units lj\n", + "region box block 0 {L:.4f} 0 {L:.4f} 0 {L:.4f}\n", + "create_box 1 box\n", + "create_atoms 1 random 500 {1000 + int(rho * 100) + int(T * 7)} box overlap 0.9 maxtry 200\n", + "mass 1 1.0\n", + "pair_style lj/cut 2.5\n", + "pair_coeff 1 1 1.0 1.0 2.5\n", + "minimize 1e-4 1e-6 200 2000\n", + "reset_timestep 0\n", + "velocity all create {T} 4928459 dist gaussian\n", + "fix 1 all nvt temp {T} {T} 0.5\n", + "thermo 1000\n", + "run 2500\n", + "\"\"\")\n", + " # average P over a further run\n", + " samples = []\n", + " for _ in range(6):\n", + " lmp.command(\"run 250\")\n", + " samples.append(lmp.get_thermo(\"press\"))\n", + " Tm = lmp.get_thermo(\"temp\")\n", + " lmp.close()\n", + " return np.mean(samples) / (rho * Tm)\n", + "\n", + "temps = [1.5, 3.0, 5.0]\n", + "rhos = [0.05, 0.1, 0.2, 0.3, 0.45, 0.6]\n", + "Z = {}\n", + "for T in temps:\n", + " Z[T] = [await measure_Z(rho, T) for rho in rhos]\n", + " print(f\"T = {T}:\", \" \".join(f\"{z:.3f}\" for z in Z[T]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The second virial coefficient, from numpy\n", + "\n", + "To lowest order in density, $Z = 1 + B_2(T)\\,\\rho$ with\n", + "\n", + "$$B_2(T) = -2\\pi \\int_0^\\infty \\left(e^{-u(r)/k_BT} - 1\\right) r^2\\,dr,$$\n", + "\n", + "a one-line numerical integral over the pair potential. $B_2 < 0$ means\n", + "attraction wins (pressure deficit); it crosses zero at the **Boyle\n", + "temperature** $T_B \\approx 3.4$, where the gas *looks* ideal far beyond\n", + "its right to:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r = np.linspace(0.01, 8, 4000)\n", + "u = 4 * (r**-12 - r**-6)\n", + "\n", + "def B2(T):\n", + " return -2 * np.pi * np.trapezoid((np.exp(-u / T) - 1) * r**2, r)\n", + "\n", + "print(\"Boyle temperature check: B2(3.0) =\", round(B2(3.0), 3),\n", + " \" B2(3.42) =\", round(B2(3.42), 3), \" B2(4.0) =\", round(B2(4.0), 3))\n", + "\n", + "rho_line = np.linspace(0, 0.65, 50)\n", + "plt.figure(figsize=(6.5, 4))\n", + "for T, color in zip(temps, [\"tab:blue\", \"tab:orange\", \"tab:green\"]):\n", + " plt.plot(rhos, Z[T], \"o\", color=color, label=f\"measured, T = {T}\")\n", + " plt.plot(rho_line, 1 + B2(T) * rho_line, \"--\", color=color, lw=1,\n", + " label=f\"virial: 1 + B₂ρ (B₂ = {B2(T):.2f})\")\n", + "plt.axhline(1, color=\"k\", lw=0.8)\n", + "plt.xlabel(\"density ρ\"); plt.ylabel(\"Z = P / ρT\")\n", + "plt.legend(fontsize=8); plt.tight_layout(); plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The virial line nails the low-density limit at every temperature — that's\n", + "the whole content of $B_2$ — and the measured points peel away from it as\n", + "higher-order terms kick in. At $T = 1.5$ attraction dominates ($Z < 1$\n", + "everywhere shown); by $T = 5$ the atoms act like hard spheres and $Z$\n", + "climbs above 1 at all densities.\n", + "\n", + "Next: [04 — Reversibility and chaos](04-reversibility-and-chaos.ipynb) —\n", + "running Newton's laws backwards." + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python", + "display_name": "Python (Pyodide)" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/notebook/content/md-basics/04-reversibility-and-chaos.ipynb b/examples/notebook/content/md-basics/04-reversibility-and-chaos.ipynb new file mode 100644 index 0000000..f61cd69 --- /dev/null +++ b/examples/notebook/content/md-basics/04-reversibility-and-chaos.ipynb @@ -0,0 +1,200 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 04 — Running time backwards: reversibility and chaos\n", + "\n", + "Newton's laws don't care which way time runs — so if a crystal melts, then\n", + "at some instant every velocity is flipped ($v \\to -v$), the atoms must\n", + "retrace their paths exactly and **re-assemble the crystal**. This is\n", + "Loschmidt's paradox: microscopic physics is reversible, yet you have never\n", + "seen a puddle spontaneously freeze into a snowflake lattice.\n", + "\n", + "The velocity Verlet integrator is *also* time-reversible, so we can\n", + "actually run this experiment — and discover the resolution of the paradox\n", + "in the process." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install lammps-js matplotlib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Melt, flip, un-melt\n", + "\n", + "A small 2D crystal given enough kinetic energy to melt, integrated for 500\n", + "steps with pure NVE (no thermostat — it would inject randomness and break\n", + "reversibility). Then two atom-style variables flip every velocity and we\n", + "keep integrating:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from lammps import lammps\n", + "\n", + "SETUP = \"\"\"\n", + "units lj\n", + "dimension 2\n", + "lattice hex 0.9\n", + "region box block 0 12 0 12 -0.1 0.1\n", + "create_box 1 box\n", + "create_atoms 1 box\n", + "mass 1 1.0\n", + "velocity all create 2.6 4321 dist gaussian\n", + "pair_style lj/cut 2.5\n", + "pair_coeff 1 1 1.0 1.0 2.5\n", + "fix 1 all nve\n", + "fix 2d all enforce2d\n", + "run 0\n", + "\"\"\"\n", + "FLIP = \"\"\"\n", + "variable nvx atom -vx\n", + "variable nvy atom -vy\n", + "velocity all set v_nvx v_nvy NULL\n", + "\"\"\"\n", + "\n", + "lmp = await lammps(output=None)\n", + "lmp.commands_string(SETUP)\n", + "\n", + "def positions_by_id(lmp):\n", + " # LAMMPS re-sorts atoms in memory as they move; order by atom id\n", + " # so positions stay comparable across time.\n", + " order = np.argsort(lmp.extract_atom(\"id\"))\n", + " return lmp.extract_atom(\"x\")[order]\n", + "\n", + "x0 = positions_by_id(lmp)\n", + "box = lmp.extract_box()\n", + "\n", + "steps, pes = [0], [lmp.get_thermo(\"pe\")]\n", + "def advance(n, chunk=50):\n", + " for _ in range(n // chunk):\n", + " lmp.command(f\"run {chunk}\")\n", + " steps.append(lmp.extract_global(\"ntimestep\"))\n", + " pes.append(lmp.get_thermo(\"pe\"))\n", + "\n", + "advance(500, chunk=25) # melt\n", + "lmp.commands_string(FLIP) # time reversal\n", + "advance(500, chunk=25) # ...and back\n", + "\n", + "d = positions_by_id(lmp) - x0\n", + "L = np.array([box[1][i] - box[0][i] for i in range(3)])\n", + "d -= L * np.round(d / L) # minimum image convention\n", + "rmsd = float(np.sqrt((d**2).sum(axis=1).mean()))\n", + "print(f\"RMS distance from the initial lattice: {rmsd:.2e}\")\n", + "lmp.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(7, 3.4))\n", + "plt.plot(steps, pes)\n", + "plt.axvline(500, color=\"k\", ls=\"--\", lw=1)\n", + "plt.text(510, min(pes) + 0.05, \"velocities flipped\", fontsize=9)\n", + "plt.xlabel(\"timestep\"); plt.ylabel(\"potential energy / atom\")\n", + "plt.title(\"Melting, played forward and then exactly backwards\")\n", + "plt.tight_layout(); plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The energy trace is a perfect mirror image and the atoms return to their\n", + "lattice sites to ~12 decimal places. Entropy decreased for 500 straight\n", + "steps — Loschmidt wins this round.\n", + "\n", + "## Chaos: why you still can't unscramble an egg\n", + "\n", + "Now repeat the experiment, waiting longer and longer before flipping.\n", + "Molecular dynamics is *chaotic*: the tiny floating-point roundoff of the\n", + "flip operation grows exponentially (a positive Lyapunov exponent), and\n", + "past some horizon the return trip no longer finds the crystal:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results = []\n", + "for nf in [250, 1000, 2000, 4000, 8000, 16000]:\n", + " lmp = await lammps(output=None)\n", + " lmp.commands_string(SETUP)\n", + " x0 = positions_by_id(lmp)\n", + " box = lmp.extract_box()\n", + " lmp.command(f\"run {nf}\")\n", + " lmp.commands_string(FLIP)\n", + " lmp.command(f\"run {nf}\")\n", + " d = positions_by_id(lmp) - x0\n", + " L = np.array([box[1][i] - box[0][i] for i in range(3)])\n", + " d -= L * np.round(d / L)\n", + " err = float(np.sqrt((d**2).sum(axis=1).mean()))\n", + " results.append((nf, err))\n", + " print(f\"flip after {nf:>6} steps → RMS return error {err:.2e}\")\n", + " lmp.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nf, err = zip(*results)\n", + "plt.figure(figsize=(6, 3.6))\n", + "plt.semilogy(nf, err, \"o-\")\n", + "plt.xlabel(\"steps before the velocity flip\")\n", + "plt.ylabel(\"RMS distance from initial lattice\")\n", + "plt.title(\"The reversibility horizon\")\n", + "plt.tight_layout(); plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There is the resolution of the paradox: reversibility is exact, but it\n", + "demands *exact* information. Any error — one part in $10^{16}$ from\n", + "roundoff, or one photon's kick in a real fluid — doubles every few\n", + "collision times, and beyond that horizon the backward movie degenerates\n", + "into just another forward one. Macroscopic irreversibility is microscopic\n", + "chaos plus finite information.\n", + "\n", + "Next: [05 — Thermal conductivity](05-thermal-conductivity.ipynb), a real\n", + "transport-coefficient measurement." + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python", + "display_name": "Python (Pyodide)" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/notebook/content/md-basics/05-thermal-conductivity.ipynb b/examples/notebook/content/md-basics/05-thermal-conductivity.ipynb new file mode 100644 index 0000000..271505e --- /dev/null +++ b/examples/notebook/content/md-basics/05-thermal-conductivity.ipynb @@ -0,0 +1,158 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 05 — Measuring thermal conductivity\n", + "\n", + "Fourier's law says a heat flux $J$ through a material sustains a\n", + "temperature gradient: $J = -\\kappa\\, \\nabla T$. Here we measure the\n", + "thermal conductivity $\\kappa$ of a 2D Lennard-Jones liquid with the\n", + "**Müller-Plathe** trick (`fix thermal/conductivity`): periodically swap the\n", + "kinetic energy of the hottest atom in one slab with the coldest atom in\n", + "another. That imposes a *known* heat flux; the fluid answers with a\n", + "temperature gradient, which we read from a `fix ave/chunk` profile written\n", + "to a file — the same file-based workflow as\n", + "[basics/05](../basics/05-thermo-output-to-file.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install lammps-js matplotlib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Impose the flux\n", + "\n", + "A wide, short box (aspect ratio 4:1). The swap region is slab 0 (cold,\n", + "left edge) and slab 10 of 20 (hot, middle), so heat flows through the two\n", + "halves in opposite directions. `f_mp` accumulates the total kinetic energy\n", + "moved by the swaps. At the end, `unfix tp` closes the profile file (an\n", + "averaging fix keeps its output file open while it is active), and we grab\n", + "the numbers the analysis needs before closing the session:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from lammps import lammps\n", + "\n", + "lmp = await lammps(output=None)\n", + "lmp.commands_string(\"\"\"\n", + "units lj\n", + "dimension 2\n", + "lattice hex 0.7\n", + "region box block 0 40 0 10 -0.1 0.1\n", + "create_box 1 box\n", + "create_atoms 1 box\n", + "mass 1 1.0\n", + "velocity all create 0.9 87287 dist gaussian\n", + "pair_style lj/cut 2.5\n", + "pair_coeff 1 1 1.0 1.0 2.5\n", + "fix 1 all nve\n", + "fix 2d all enforce2d\n", + "thermo 2000\n", + "run 4000\n", + "reset_timestep 0\n", + "\n", + "compute ck all chunk/atom bin/1d x lower 0.05 units reduced\n", + "fix tp all ave/chunk 10 1000 20000 ck temp file tprofile.txt\n", + "fix mp all thermal/conductivity 20 x 20\n", + "variable swapped equal f_mp\n", + "run 20000\n", + "\"\"\")\n", + "e_swapped = lmp.extract_variable(\"swapped\")\n", + "lmp.command(\"unfix tp\") # closes tprofile.txt so it lands in the notebook filesystem\n", + "t_run = 20000 * lmp.extract_global(\"dt\")\n", + "box = lmp.extract_box()\n", + "Lx = box[1][0] - box[0][0]\n", + "Ly = box[1][1] - box[0][1]\n", + "lmp.close()\n", + "print(f\"kinetic energy transported by swaps: {e_swapped:.1f} (LJ units) \"\n", + " f\"over t = {t_run:.0f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Read the temperature profile and fit the gradient\n", + "\n", + "`ave/chunk` wrote one block of 20 bins to `tprofile.txt` (check the file\n", + "browser — same file-based workflow as\n", + "[basics/05](../basics/05-thermo-output-to-file.ipynb)), time-averaged over\n", + "the second half of the run; the first half is the transient while the\n", + "gradient builds up. The profile is tent-shaped — cold at the swap slab on\n", + "the left edge, hot in the middle — so we fit a straight line to each half\n", + "and average the slopes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "raw = np.loadtxt(\"tprofile.txt\", skiprows=4) # bin, coord, natoms, temp\n", + "xfrac, temp = raw[:, 1], raw[:, 3]\n", + "x = xfrac * Lx\n", + "\n", + "up, down = slice(1, 10), slice(11, 20) # skip the swap bins\n", + "s1 = np.polyfit(x[up], temp[up], 1)\n", + "s2 = np.polyfit(x[down], temp[down], 1)\n", + "dTdx = (abs(s1[0]) + abs(s2[0])) / 2\n", + "\n", + "# Fourier: J = kappa dT/dx. The swapped energy splits over two directions\n", + "# and crosses area Ly (a length, in 2D).\n", + "J = e_swapped / (2 * t_run * Ly)\n", + "kappa = J / dTdx\n", + "print(f\"dT/dx = {dTdx:.4f} J = {J:.4f} → κ ≈ {kappa:.2f} (LJ units)\")\n", + "\n", + "plt.figure(figsize=(6.5, 3.6))\n", + "plt.plot(x, temp, \"o\", ms=4)\n", + "plt.plot(x[up], np.polyval(s1, x[up]), \"--\", color=\"tab:red\")\n", + "plt.plot(x[down], np.polyval(s2, x[down]), \"--\", color=\"tab:red\")\n", + "plt.xlabel(\"x\"); plt.ylabel(\"temperature\")\n", + "plt.title(f\"Steady-state gradient → κ ≈ {kappa:.2f}\")\n", + "plt.tight_layout(); plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A real transport coefficient from a 30-second browser simulation. Things\n", + "worth trying: double the run length (the profile gets smoother), change the\n", + "density, or measure at several temperatures — for a dense LJ liquid κ\n", + "falls as it warms.\n", + "\n", + "Next: [06 — Polymer chains](06-polymer-chains.ipynb), where numpy builds\n", + "the molecules." + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python", + "display_name": "Python (Pyodide)" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/notebook/content/md-basics/06-polymer-chains.ipynb b/examples/notebook/content/md-basics/06-polymer-chains.ipynb new file mode 100644 index 0000000..f36af5e --- /dev/null +++ b/examples/notebook/content/md-basics/06-polymer-chains.ipynb @@ -0,0 +1,183 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 06 — Polymer chains: build molecules with numpy\n", + "\n", + "The bead-spring (Kremer–Grest) polymer is *the* model of polymer physics:\n", + "beads connected by FENE bonds, purely repulsive LJ between everything, a\n", + "Langevin thermostat as implicit solvent. We build the chains ourselves —\n", + "a LAMMPS data file is just text, and thanks to the shared filesystem\n", + "LAMMPS reads what Python writes with plain `open()`.\n", + "\n", + "The physics payoff: a self-avoiding chain of $N$ beads swells as\n", + "$R_g \\sim N^{\\nu}$ with the Flory exponent $\\nu \\approx 0.588$ —\n", + "measurably bigger than the random-walk value $1/2$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install lammps-js matplotlib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A data file from numpy\n", + "\n", + "`make_chain(N)` writes an `N`-bead zig-zag chain (bond lengths near the\n", + "FENE minimum) with `N−1` bonds. Note the file is created with ordinary\n", + "Python `open()` — check the file browser after running:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from lammps import lammps\n", + "\n", + "def make_chain(N, path):\n", + " lines = [f\"bead-spring chain, N = {N}\", \"\",\n", + " f\"{N} atoms\", f\"{N - 1} bonds\",\n", + " \"1 atom types\", \"1 bond types\", \"\",\n", + " \"-60 60 xlo xhi\", \"-60 60 ylo yhi\", \"-60 60 zlo zhi\", \"\",\n", + " \"Masses\", \"\", \"1 1.0\", \"\", \"Atoms\", \"\"]\n", + " for i in range(N):\n", + " lines.append(f\"{i + 1} 1 1 {i * 0.93:.3f} {0.15 * (i % 2):.3f} 0.0\")\n", + " lines += [\"\", \"Bonds\", \"\"]\n", + " for i in range(N - 1):\n", + " lines.append(f\"{i + 1} 1 {i + 1} {i + 2}\")\n", + " with open(path, \"w\") as f:\n", + " f.write(\"\\n\".join(lines) + \"\\n\")\n", + "\n", + "make_chain(16, \"chain.data\")\n", + "print(open(\"chain.data\").read()[:180], \"…\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulate one chain length\n", + "\n", + "Standard Kremer–Grest parameters: FENE with $k = 30$, $R_0 = 1.5$, and a\n", + "WCA (purely repulsive, cut at $2^{1/6}$) pair interaction. `compute\n", + "gyration` reports $R_g$; we sample it every 500 steps after an\n", + "equilibration proportional to the chain's relaxation time. The longer\n", + "chains take a minute — progress prints as it goes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "async def gyration(N):\n", + " make_chain(N, \"chain.data\")\n", + " lmp = await lammps(output=None)\n", + " lmp.commands_string(f\"\"\"\n", + "units lj\n", + "atom_style bond\n", + "special_bonds fene\n", + "read_data chain.data\n", + "bond_style fene\n", + "bond_coeff 1 30.0 1.5 1.0 1.0\n", + "pair_style lj/cut 1.122462\n", + "pair_coeff 1 1 1.0 1.0\n", + "pair_modify shift yes\n", + "neigh_modify every 1 delay 0\n", + "timestep 0.005\n", + "fix 1 all nve\n", + "fix lang all langevin 1.0 1.0 2.0 {90429 + N}\n", + "compute g all gyration\n", + "variable rg equal c_g\n", + "run {1000 * N}\n", + "\"\"\")\n", + " samples = []\n", + " for _ in range(4 * N):\n", + " lmp.command(\"run 500\")\n", + " samples.append(lmp.extract_variable(\"rg\"))\n", + " lmp.close()\n", + " return float(np.mean(samples)), float(np.std(samples))\n", + "\n", + "Ns = [8, 16, 32, 48]\n", + "rg = {}\n", + "for N in Ns:\n", + " rg[N] = await gyration(N)\n", + " print(f\"N = {N:>3}: Rg = {rg[N][0]:.3f} ± {rg[N][1]:.3f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The Flory exponent\n", + "\n", + "Fit $\\log R_g$ against $\\log(N-1)$ (the number of bonds). Self-avoiding\n", + "walks in 3D give $\\nu \\approx 0.588$; an ideal random walk would give\n", + "$1/2$. Short chains overshoot the asymptotic exponent a little — real\n", + "(and simulated) polymer physics has finite-size corrections:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_bonds = np.array([N - 1 for N in Ns], dtype=float)\n", + "rg_mean = np.array([rg[N][0] for N in Ns])\n", + "nu, log_pref = np.polyfit(np.log(n_bonds), np.log(rg_mean), 1)\n", + "\n", + "n_line = np.linspace(n_bonds[0], n_bonds[-1], 50)\n", + "plt.figure(figsize=(6, 4))\n", + "plt.errorbar(n_bonds, rg_mean, yerr=[rg[N][1] for N in Ns], fmt=\"o\",\n", + " capsize=3, label=\"measured\")\n", + "plt.plot(n_line, np.exp(log_pref) * n_line**nu, \"-\",\n", + " label=f\"fit: ν = {nu:.3f}\")\n", + "plt.plot(n_line, rg_mean[0] * (n_line / n_bonds[0])**0.588, \"--\",\n", + " label=\"SAW theory: ν = 0.588\")\n", + "plt.plot(n_line, rg_mean[0] * (n_line / n_bonds[0])**0.5, \":\",\n", + " label=\"ideal chain: ν = 0.5\")\n", + "plt.xscale(\"log\"); plt.yscale(\"log\")\n", + "plt.xlabel(\"number of bonds N − 1\"); plt.ylabel(\"radius of gyration $R_g$\")\n", + "plt.legend(fontsize=9); plt.tight_layout(); plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The excluded-volume repulsion visibly beats the random walk — the chain\n", + "*swells*. From here it's a short path to real soft-matter simulations:\n", + "many chains in a melt (where, beautifully, ν drops back to ½), grafted\n", + "brushes, or pulling a single chain by its ends with `fix spring`.\n", + "\n", + "That's the end of this series — back to [the index](../index.ipynb)." + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python", + "display_name": "Python (Pyodide)" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/python/lammps/__init__.py b/python/lammps/__init__.py index af7e228..5485e20 100644 --- a/python/lammps/__init__.py +++ b/python/lammps/__init__.py @@ -175,6 +175,19 @@ const node = nodeOrStream.node ?? nodeOrStream; if (attr.mode !== undefined) node.mode = attr.mode; if (attr.timestamp !== undefined) node.timestamp = attr.timestamp; + if (attr.size !== undefined && FS.isFile(node.mode)) { + // Truncation (e.g. fopen "w" on an existing file) arrives here as + // a node-level size change, before any stream is opened. + const path = toAPIPath(realPath(node)); + let data = new Uint8Array(attr.size); + if (attr.size > 0) { + try { + const old = pyAPI.get(path).data || new Uint8Array(0); + data.set(old.subarray(0, Math.min(old.length, attr.size))); + } catch {} + } + pyAPI.put(path, { data }); + } }, lookup(parentOrStream, name) { const parent = parentOrStream.node ?? parentOrStream; diff --git a/tests/drivefs-mount.spec.ts b/tests/drivefs-mount.spec.ts new file mode 100644 index 0000000..7f8ce98 --- /dev/null +++ b/tests/drivefs-mount.spec.ts @@ -0,0 +1,468 @@ +import { beforeAll, beforeEach, describe, expect, it } from "vitest"; +import { readFileSync } from "node:fs"; +import { join } from "node:path"; + +import type { LammpsModule } from "../types"; +import { loadModule } from "./helpers/lammps"; + +// The DriveFS mount lives as a JS snippet inside the Python bindings +// (python/lammps/__init__.py, _MOUNT_DRIVEFS_JS). These tests extract that +// exact snippet and exercise every branch of it against the real wasm +// filesystem, with the JupyterLite ContentsAPI replaced by an in-memory +// mock. Branch checklist: +// +// toAPIPath : cwd inside /drive (prefix), cwd == /drive (no prefix), +// cwd outside /drive, cwd() throwing +// node_ops : lookup hit/miss, mknod (file + dir), getattr, setattr +// (mode, timestamp), rename, unlink, rmdir, readdir, +// symlink (EPERM), readlink (EINVAL) +// stream_ops : open existing / missing+write (mknod) / missing+read +// (ENOENT) / O_TRUNC, close with and without write flags, +// read (normal, past-EOF, zero-length), write (grow, +// in-place), llseek (CUR, END, negative) +// mount : first mount, re-mount (unmount + mkdirTree branches) +// end-to-end : LAMMPS write_dump / fix ave-time file / fix ave-chunk + +// unfix / read_data + include of Python-written files / +// re-run overwriting a previous output + +function extractMountJs(): string { + const py = readFileSync( + join(process.cwd(), "python", "lammps", "__init__.py"), + "utf8", + ); + const match = py.match(/_MOUNT_DRIVEFS_JS = """\n([\s\S]*?)\n"""/); + if (!match) throw new Error("_MOUNT_DRIVEFS_JS not found in bindings"); + // Interpret the python string escapes (only \\ occurs in the snippet). + return match[1].replace(/\\\\/g, "\\"); +} + +interface StoredFile { + data: Uint8Array; +} + +/** In-memory stand-in for JupyterLite's DriveFS ContentsAPI. */ +class MockContentsAPI { + files = new Map(); + dirs = new Set(["/"]); + putCalls: string[] = []; + + seed(path: string, text: string) { + const parts = path.split("/").slice(1, -1); + let prefix = ""; + for (const part of parts) { + prefix += "/" + part; + this.dirs.add(prefix); + } + this.files.set(path, { data: new TextEncoder().encode(text) }); + } + + text(path: string): string | undefined { + const f = this.files.get(path); + return f ? new TextDecoder().decode(f.data) : undefined; + } + + lookup(path: string) { + if (path === "/" || this.dirs.has(path)) return { ok: true, mode: 16895 }; + if (this.files.has(path)) return { ok: true, mode: 33206 }; + return { ok: false }; + } + + getattr(path: string) { + const size = this.files.get(path)?.data.length ?? 0; + const now = new Date(0); + return { + dev: 1, ino: 0, mode: this.files.has(path) ? 33206 : 16895, nlink: 1, + uid: 0, gid: 0, rdev: 0, size, atime: now, mtime: now, ctime: now, + blksize: 4096, blocks: Math.ceil(size / 4096), + }; + } + + get(path: string): StoredFile { + const f = this.files.get(path); + if (!f) throw new Error(`get: ${path} not found`); + // Return a copy-holder so proxy writes don't mutate the store until put. + return { data: f.data.slice() }; + } + + put(path: string, file: StoredFile) { + this.putCalls.push(path); + this.files.set(path, { data: file.data.slice() }); + } + + mknod(path: string, mode: number) { + if ((mode & 61440) === 16384) this.dirs.add(path); + else this.files.set(path, { data: new Uint8Array() }); + } + + rename(oldPath: string, newPath: string) { + const f = this.files.get(oldPath); + if (f) { + this.files.delete(oldPath); + this.files.set(newPath, f); + } else if (this.dirs.has(oldPath)) { + this.dirs.delete(oldPath); + this.dirs.add(newPath); + } + } + + rmdir(path: string) { + this.files.delete(path); + this.dirs.delete(path); + } + + readdir(path: string): string[] { + // The proxy hands us the mount root as "/", subdirs as "/name/" or + // "/name" — normalize away a trailing slash like the real API does. + if (path.length > 1 && path.endsWith("/")) path = path.slice(0, -1); + const prefix = path === "/" ? "/" : path + "/"; + const names = new Set(); + for (const p of [...this.files.keys(), ...this.dirs]) { + if (p !== path && p.startsWith(prefix)) { + names.add(p.slice(prefix.length).split("/")[0]); + } + } + return [".", "..", ...names]; + } +} + +let wasm: LammpsModule; +let api: MockContentsAPI; +let cwd: () => string; +let mount: (mod: LammpsModule, pyFS: unknown) => boolean; + +function pyodideFsMock() { + return { + lookupPath: () => ({ node: { mount: { type: { API: api } } } }), + cwd: () => cwd(), + }; +} + +function remount() { + api = new MockContentsAPI(); + cwd = () => "/drive/nb"; + expect(mount(wasm, pyodideFsMock())).toBe(true); +} + +beforeAll(async () => { + wasm = await loadModule(); + // eslint-disable-next-line no-eval + mount = (0, eval)(extractMountJs()); +}); + +beforeEach(() => { + remount(); +}); + +describe("mounting", () => { + it("mounts /work on the contents API and chdirs into it", () => { + expect(wasm.FS.cwd()).toBe("/work"); + }); + + it("can re-mount over an existing mount (unmount + mkdirTree branches)", () => { + remount(); // beforeEach already mounted once; this exercises unmount + wasm.FS.writeFile("/work/x.txt", "hello"); + expect(api.text("/nb/x.txt")).toBe("hello"); + }); +}); + +describe("toAPIPath cwd prefix", () => { + it("prefixes with the kernel cwd inside /drive", () => { + cwd = () => "/drive/sub/dir"; + wasm.FS.writeFile("/work/a.txt", "A"); + expect(api.text("/sub/dir/a.txt")).toBe("A"); + }); + + it("uses the drive root when cwd is exactly /drive", () => { + cwd = () => "/drive"; + wasm.FS.writeFile("/work/b.txt", "B"); + expect(api.text("/b.txt")).toBe("B"); + }); + + it("uses the drive root when cwd is outside /drive", () => { + cwd = () => "/tmp"; + wasm.FS.writeFile("/work/c.txt", "C"); + expect(api.text("/c.txt")).toBe("C"); + }); + + it("falls back to the drive root when cwd() throws", () => { + cwd = () => { + throw new Error("no cwd"); + }; + wasm.FS.writeFile("/work/d.txt", "D"); + expect(api.text("/d.txt")).toBe("D"); + }); + + it("re-evaluates the cwd on every operation", () => { + cwd = () => "/drive/one"; + wasm.FS.writeFile("/work/e.txt", "E1"); + cwd = () => "/drive/two"; + wasm.FS.writeFile("/work/e.txt", "E2"); + expect(api.text("/one/e.txt")).toBe("E1"); + expect(api.text("/two/e.txt")).toBe("E2"); + }); +}); + +describe("node_ops", () => { + it("lookup finds files seeded in the contents store", () => { + api.seed("/nb/seeded.txt", "from python"); + expect(wasm.FS.readFile("/work/seeded.txt", { encoding: "utf8" })).toBe( + "from python", + ); + }); + + it("lookup of a missing file throws ENOENT", () => { + expect(() => wasm.FS.readFile("/work/missing.txt")).toThrow(); + }); + + it("mknod creates directories in the store", () => { + wasm.FS.mkdir("/work/newdir"); + expect(api.dirs.has("/nb/newdir")).toBe(true); + wasm.FS.writeFile("/work/newdir/f.txt", "nested"); + expect(api.text("/nb/newdir/f.txt")).toBe("nested"); + }); + + it("getattr reports the stored size", () => { + api.seed("/nb/sized.txt", "12345"); + const st = wasm.FS.stat("/work/sized.txt"); + expect(st.size).toBe(5); + }); + + it("setattr accepts mode and timestamp changes", () => { + api.seed("/nb/attr.txt", "x"); + wasm.FS.chmod("/work/attr.txt", 0o600); // mode branch + wasm.FS.utime("/work/attr.txt", 1000, 1000); // timestamp branch + expect(api.text("/nb/attr.txt")).toBe("x"); + }); + + it("rename moves the file in the store", () => { + wasm.FS.writeFile("/work/old.txt", "payload"); + wasm.FS.rename("/work/old.txt", "/work/new.txt"); + expect(api.text("/nb/new.txt")).toBe("payload"); + expect(api.text("/nb/old.txt")).toBeUndefined(); + }); + + it("unlink removes the file from the store", () => { + wasm.FS.writeFile("/work/gone.txt", "bye"); + wasm.FS.unlink("/work/gone.txt"); + expect(api.text("/nb/gone.txt")).toBeUndefined(); + }); + + it("rmdir removes a directory from the store", () => { + wasm.FS.mkdir("/work/tmpdir"); + wasm.FS.rmdir("/work/tmpdir"); + expect(api.dirs.has("/nb/tmpdir")).toBe(false); + }); + + it("readdir lists store entries", () => { + wasm.FS.writeFile("/work/l1.txt", "1"); + api.seed("/nb/l2.txt", "2"); + const names = wasm.FS.readdir("/work"); + expect(names).toContain("l1.txt"); + expect(names).toContain("l2.txt"); + }); + + it("symlink and readlink are rejected", () => { + expect(() => wasm.FS.symlink("/work/l1.txt", "/work/link")).toThrow(); + expect(() => wasm.FS.readlink("/work/l1.txt")).toThrow(); + }); +}); + +describe("stream_ops", () => { + it("open of a missing file in read mode throws ENOENT", () => { + expect(() => wasm.FS.open("/work/nope.txt", "r")).toThrow(); + }); + + it("write-mode open of a missing file creates it (mknod branch)", () => { + const s = wasm.FS.open("/work/created.txt", "w"); + wasm.FS.write(s, new TextEncoder().encode("made"), 0, 4); + wasm.FS.close(s); + expect(api.text("/nb/created.txt")).toBe("made"); + }); + + it("read-only close does not put", () => { + api.seed("/nb/ro.txt", "read me"); + api.putCalls = []; + const s = wasm.FS.open("/work/ro.txt", "r"); + const buf = new Uint8Array(7); + wasm.FS.read(s, buf, 0, 7); + wasm.FS.close(s); + expect(new TextDecoder().decode(buf)).toBe("read me"); + expect(api.putCalls).toEqual([]); + }); + + it("read past EOF returns 0 bytes", () => { + api.seed("/nb/short.txt", "ab"); + const s = wasm.FS.open("/work/short.txt", "r"); + const buf = new Uint8Array(10); + const n1 = wasm.FS.read(s, buf, 0, 10); + const n2 = wasm.FS.read(s, buf, 0, 10); + wasm.FS.close(s); + expect(n1).toBe(2); + expect(n2).toBe(0); + }); + + it("writes grow the buffer and in-place writes preserve the rest", () => { + wasm.FS.writeFile("/work/grow.txt", "0123456789"); + const s = wasm.FS.open("/work/grow.txt", "r+"); + wasm.FS.write(s, new TextEncoder().encode("AB"), 0, 2, 3); // in place + wasm.FS.close(s); + expect(api.text("/nb/grow.txt")).toBe("012AB56789"); + }); + + it("append mode seeks to the end (SEEK_END branch)", () => { + wasm.FS.writeFile("/work/log.txt", "line1\n"); + const s = wasm.FS.open("/work/log.txt", "a"); + wasm.FS.write(s, new TextEncoder().encode("line2\n"), 0, 6); + wasm.FS.close(s); + expect(api.text("/nb/log.txt")).toBe("line1\nline2\n"); + }); + + it("llseek with SEEK_CUR advances relative to the position", () => { + api.seed("/nb/seek.txt", "abcdef"); + const s = wasm.FS.open("/work/seek.txt", "r"); + const buf = new Uint8Array(2); + wasm.FS.read(s, buf, 0, 2); + wasm.FS.llseek(s, 2, 1); // SEEK_CUR: 2 + 2 = 4 + wasm.FS.read(s, buf, 0, 2); + wasm.FS.close(s); + expect(new TextDecoder().decode(buf)).toBe("ef"); + }); + + it("llseek to a negative position throws EINVAL", () => { + api.seed("/nb/neg.txt", "abc"); + const s = wasm.FS.open("/work/neg.txt", "r"); + expect(() => wasm.FS.llseek(s, -10, 0)).toThrow(); + wasm.FS.close(s); + }); + + it("re-opening an existing file with 'w' truncates the old content", () => { + wasm.FS.writeFile("/work/trunc.txt", "a much longer previous content"); + wasm.FS.writeFile("/work/trunc.txt", "short"); + expect(api.text("/nb/trunc.txt")).toBe("short"); + }); + + it("truncate to a nonzero size keeps the leading bytes", () => { + wasm.FS.writeFile("/work/part.txt", "abcdefgh"); + wasm.FS.truncate("/work/part.txt", 3); + expect(api.text("/nb/part.txt")).toBe("abc"); + }); + + it("truncate can also extend a file with zero padding", () => { + wasm.FS.writeFile("/work/pad.txt", "ab"); + wasm.FS.truncate("/work/pad.txt", 4); + expect(api.files.get("/nb/pad.txt")!.data).toEqual( + new Uint8Array([97, 98, 0, 0]), + ); + }); +}); + +describe("LAMMPS end-to-end through the mount", () => { + const SETUP = ` +units lj +atom_style atomic +lattice fcc 0.8442 +region box block 0 2 0 2 0 2 +create_box 1 box +create_atoms 1 box +mass 1 1.0 +velocity all create 1.44 87287 +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 +fix 1 all nve +`; + + function session() { + const lmp = new wasm.LAMMPSWeb(); + lmp.start(); + return lmp; + } + + it("write_dump lands next to the notebook and parses", () => { + const lmp = session(); + lmp.runScript(SETUP + "run 10\nwrite_dump all atom out.dump\n"); + lmp.stop(); + const dump = api.text("/nb/out.dump"); + expect(dump).toBeDefined(); + const lines = dump!.trim().split("\n"); + expect(lines[0]).toBe("ITEM: TIMESTEP"); + expect(lines.length).toBe(9 + 32); // header + one row per atom + }); + + it("fix ave/time output is complete after the session closes", () => { + const lmp = session(); + lmp.runScript( + SETUP + + "variable t equal temp\n" + + "fix out all ave/time 5 1 5 v_t file thermo_test.txt\n" + + "run 50\n", + ); + lmp.stop(); // closes the fix's file -> flushed through put + const rows = api + .text("/nb/thermo_test.txt")! + .split("\n") + .filter((l) => l && !l.startsWith("#")); + expect(rows.length).toBe(11); // steps 0, 5, …, 50 + for (const row of rows) { + const [step, value] = row.trim().split(/\s+/).map(Number); + expect(step % 5).toBe(0); + expect(Number.isFinite(value)).toBe(true); + } + }); + + it("fix ave/chunk output is complete after unfix", () => { + const lmp = session(); + lmp.runScript( + SETUP + + "compute ck all chunk/atom bin/1d x lower 0.25 units reduced\n" + + "fix tp all ave/chunk 5 2 10 ck temp file profile_test.txt\n" + + "run 10\n" + + "unfix tp\n", + ); + const rows = api + .text("/nb/profile_test.txt")! + .split("\n") + .filter((l) => l && !l.startsWith("#")); + lmp.stop(); + expect(rows.length).toBe(1 + 4); // block header + 4 bins + }); + + it("LAMMPS reads scripts and data files written to the store", () => { + api.seed("/nb/tiny.in", SETUP + "run 5\nwrite_data tiny.data\n"); + const lmp = session(); + lmp.runFile("tiny.in"); + lmp.stop(); + expect(api.text("/nb/tiny.data")).toContain("32 atoms"); + + const lmp2 = session(); + lmp2.runScript( + "units lj\natom_style atomic\npair_style lj/cut 2.5\nread_data tiny.data\n", + ); + expect(lmp2.getThermo("atoms")).toBe(32); + lmp2.stop(); + }); + + it("re-running a script overwrites the previous, longer output cleanly", () => { + const lmp = session(); + lmp.runScript( + SETUP + + "variable t equal temp\n" + + "fix out all ave/time 5 1 5 v_t file rerun_test.txt\n" + + "run 100\n", + ); + lmp.stop(); + const first = api.text("/nb/rerun_test.txt")!; + + const lmp2 = session(); + lmp2.runScript( + SETUP + + "variable t equal temp\n" + + "fix out all ave/time 5 1 5 v_t file rerun_test.txt\n" + + "run 10\n", + ); + lmp2.stop(); + const second = api.text("/nb/rerun_test.txt")!; + expect(second.length).toBeLessThan(first.length); + const rows = second.split("\n").filter((l) => l && !l.startsWith("#")); + expect(rows.length).toBe(3); // steps 0, 5, 10 — no stale tail from the first run + }); +});