Add sphere and isotropic rotational diffusion models

docs/docs/tutorials/diffusion_in_sphere.ipynb

@@ -0,0 +1,279 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# Diffusion inside a sphere\n",
+    "The model<sup>1,2</sup> describes a particle that diffuses with diffusion coefficient $D$ inside a sphere of radius $R$ that it cannot leave, for example an atom of a side chain moving within a limited volume. The sphere is isotropic, so the scattering depends only on the magnitude of $Q$:\n",
+    "\n",
+    "$$ S(Q, E) = A_0(Q)\\,\\delta(E) + \\sum_{(l, n)} A_{nl}(Q)\\, L(E, \\Gamma_{nl}) $$\n",
+    "\n",
+    "with\n",
+    "\n",
+    "$$ A_0(Q) = \\left(\\frac{3 j_1(QR)}{QR}\\right)^2,\\quad\n",
+    "A_{nl}(Q) = (2l+1) \\frac{6 x_{nl}^2}{x_{nl}^2 - l(l+1)}\n",
+    "\\left(\\frac{QR\\, j_{l+1}(QR) - l\\, j_l(QR)}{(QR)^2 - x_{nl}^2}\\right)^2,\\quad\n",
+    "\\Gamma_{nl} = \\frac{\\hbar x_{nl}^2 D}{R^2}, $$\n",
+    "\n",
+    "where $j_l$ are the spherical Bessel functions and $x_{nl}$ the roots of $j_l'(x) = 0$, which express that no particles flow through the surface. $A_0$ is the Elastic Incoherent Structure Factor (EISF) and $A_{nl}$ the Quasi-Elastic Incoherent Structure Factor (QISF) of the root $x_{nl}$. The widths do not depend on $Q$; only the QISF do.\n",
+    "\n",
+    "The sum is infinite. The model uses the 98 roots below 26.84, the same roots as Mantid<sup>3</sup>, and by default only the lowest 10 of them, because every term adds Parameters that slow down a fit."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "1",
+   "metadata": {},
+   "source": [
+    "import warnings\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "import easydynamics as edyn\n",
+    "\n",
+    "%matplotlib widget"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2",
+   "metadata": {},
+   "source": [
+    "## The model at different Q values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "3",
+   "metadata": {},
+   "source": [
+    "# Q is in Angstrom^-1, the radius in Angstrom and D in Angstrom^2/ps.\n",
+    "# With R = 1 Angstrom, QR stays below 1.5, where 10 terms are enough.\n",
+    "Q = np.linspace(0.3, 1.5, 6)\n",
+    "\n",
+    "sphere_model = edyn.DiffusionInSphere(\n",
+    "    scale=1.0,\n",
+    "    radius=1.0,\n",
+    "    diffusion_coefficient=0.1,\n",
+    "    Q=Q,\n",
+    ")\n",
+    "print(sphere_model)"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "code",
+   "id": "4",
+   "metadata": {},
+   "source": [
+    "# Plot the quasi-elastic part (all Lorentzians) at every Q. The delta function\n",
+    "# would be a spike at E = 0; its weight, the EISF, is plotted further down.\n",
+    "energy = np.linspace(-0.5, 0.5, 501)\n",
+    "component_collections = sphere_model.get_component_collections()\n",
+    "\n",
+    "cmap = plt.cm.jet\n",
+    "nQ = len(component_collections)\n",
+    "plt.figure()\n",
+    "for Q_index, collection in enumerate(component_collections):\n",
+    "    quasielastic = sum(\n",
+    "        component.evaluate(energy)\n",
+    "        for component in collection\n",
+    "        if isinstance(component, edyn.Lorentzian)\n",
+    "    )\n",
+    "    color = cmap(Q_index / (nQ - 1))\n",
+    "    plt.plot(energy, quasielastic, label=f'Q={Q[Q_index]:.2f} Å^-1', color=color)\n",
+    "\n",
+    "plt.legend()\n",
+    "plt.xlabel('Energy (meV)')\n",
+    "plt.ylabel('Intensity (arb. units)')\n",
+    "plt.title('Diffusion inside a sphere, R = 1 Å')\n",
+    "plt.show()"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5",
+   "metadata": {},
+   "source": [
+    "## Roots and widths"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "metadata": {},
+   "source": [
+    "The Lorentzian of order $j$ belongs to the $j$-th smallest root. `roots` gives each root as $(x_{nl}, l, n)$, and the width grows with $x_{nl}^2$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "for order, (x, angular_order, n) in enumerate(sphere_model.roots, start=1):\n",
+    "    width = sphere_model.calculate_width(order=order)[0]\n",
+    "    print(f'order {order:2d}: x = {x:6.3f}, l = {angular_order}, n = {n}, HWHM = {width:.4f} meV')"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8",
+   "metadata": {},
+   "source": [
+    "## EISF and QISF as a function of Q"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9",
+   "metadata": {},
+   "source": [
+    "At small $QR$ almost all intensity is elastic. As $QR$ grows, the weight moves into the Lorentzians, first into the one with the smallest root ($l = 1$). The EISF and all QISF together sum to one, minus the weight of the terms that are left out."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "10",
+   "metadata": {},
+   "source": [
+    "Q_fine = np.linspace(0.01, 1.5, 200)\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.plot(Q_fine, sphere_model.calculate_EISF(Q_fine), color='k', label='EISF')\n",
+    "for order in range(1, 5):\n",
+    "    _, angular_order, n = sphere_model.roots[order - 1]\n",
+    "    plt.plot(\n",
+    "        Q_fine,\n",
+    "        sphere_model.calculate_QISF_order(Q_fine, order=order),\n",
+    "        label=f'QISF of order {order} (l = {angular_order}, n = {n})',\n",
+    "    )\n",
+    "plt.xlabel('Q (Å$^{-1}$)')\n",
+    "plt.ylabel('Weight')\n",
+    "plt.title('EISF and QISF for R = 1 Å')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11",
+   "metadata": {},
+   "source": [
+    "## Truncation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "12",
+   "metadata": {},
+   "source": [
+    "Leaving out terms loses intensity at large $QR$. A warning is raised when more than $10^{-3}$ of the weight is lost at any of the model's $Q$ values, and `calculate_missing_weight` returns the weight that is left out. With 10 terms this starts at $QR \\approx 1.8$, and with all 98 at $QR \\approx 9$. A model without $Q$ values builds no components, so it is quick to create for this check."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "13",
+   "metadata": {},
+   "source": [
+    "QR = np.linspace(0.01, 12, 300)\n",
+    "\n",
+    "plt.figure()\n",
+    "for number_of_lorentzians in [10, 30, 98]:\n",
+    "    model = edyn.DiffusionInSphere(radius=1.0, number_of_lorentzians=number_of_lorentzians)\n",
+    "    plt.semilogy(\n",
+    "        QR,\n",
+    "        np.maximum(model.calculate_missing_weight(QR), 1e-16),\n",
+    "        label=f'{number_of_lorentzians} terms',\n",
+    "    )\n",
+    "plt.axhline(1e-3, color='k', linestyle='dashed', label='warning threshold')\n",
+    "plt.xlabel('QR')\n",
+    "plt.ylabel('Missing weight')\n",
+    "plt.title('Weight left out by the truncated sum')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "14",
+   "metadata": {},
+   "source": [
+    "A larger sphere reaches the threshold at smaller $Q$. With $R = 3$ Å and the default 10 terms, $Q$ up to 2 Å$^{-1}$ is already too far:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "15",
+   "metadata": {},
+   "source": [
+    "with warnings.catch_warnings(record=True) as caught:\n",
+    "    warnings.simplefilter('always')\n",
+    "    large_sphere = edyn.DiffusionInSphere(radius=3.0, Q=np.linspace(0.3, 2, 4))\n",
+    "for warning in caught:\n",
+    "    print(warning.message)\n",
+    "\n",
+    "print('Missing weight per Q:', np.round(large_sphere.calculate_missing_weight(), 4))"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "16",
+   "metadata": {},
+   "source": [
+    "More terms remove the warning but make every update of the model slower. The spherical Bessel functions are faster with the optional numba dependency, `pip install easydynamics[fast]`; without it SciPy gives the same values."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "17",
+   "metadata": {},
+   "source": [
+    "## References\n",
+    "\n",
+    "1. F. Volino and A. J. Dianoux, *Mol. Phys.* **41**, 271 (1980), [doi:10.1080/00268978000102761](https://doi.org/10.1080/00268978000102761).\n",
+    "\n",
+    "2. D. J. Bicout, *Phys. Rev. E* **62**, 261 (2000), [doi:10.1103/PhysRevE.62.261](https://doi.org/10.1103/PhysRevE.62.261).\n",
+    "\n",
+    "3. Mantid, [`InelasticDiffSphere`](https://github.com/mantidproject/mantid/blob/6ec66a6f524a7ec9d4258351fe3001c66c63dbe5/Framework/CurveFitting/src/Functions/InelasticDiffSphere.cpp#L53-L69): the table of the 98 roots."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "default",
+   "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.14.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}