alignment_distribution_function.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mishael Derla, FAU Erlangen-Nürnberg, winter term 2024/2025"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os import listdir\n",
    "\n",
    "import numpy\n",
    "from scipy.special import gamma as euler_gamma_function\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from xyz_file_io import XYZ, active_particle_metadata\n",
    "from physics import pbc_iterator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our $n$-body alignment distribution function reads\n",
    "$$\\boldsymbol{a}^{(n)}(\\vec{r}_1,\\dots,\\vec{r}_n)=\\left\\langle\\;\\sum_{i_1}\\sum_{i_1\\ne i_2}\\cdots\\sum_{i_n\\ne i_1,i_2,\\dots,i_{n-1}}\\frac{\\vec{F}^\\text{(self)}_{i_1}}{F^\\text{(self)}}\\frac{\\delta(\\vec{r}_1-\\vec{x}_{i_1})}{\\rho(\\vec{r}_1)}\\otimes\\cdots\\otimes\\frac{\\vec{F}^\\text{(self)}_{i_n}}{F^\\text{(self)}}\\frac{\\delta(\\vec{r}_1-\\vec{x}_{i_n})}{\\rho(\\vec{r}_n)}\\;\\right\\rangle$$\n",
    "When the system is homogeneous and isotropic (which our system would be for $N,V\\to\\infty$ as $N/V=:\\rho$ remains constant), this reduces to a scalar $a{(2)}(r)$, which is approximated by the following binning\n",
    "$$a^{(2)}(r)\\approx\\frac{1}{Nd}\\cdot\\frac{1}{S_dr^{d-1}\\Delta r\\,\\rho}\\cdot\\underbrace{\\left\\langle\\sum_i\\sum_{j\\ne i}\\frac{\\vec{F}_i^\\text{(self)}}{F^\\text{(self)}}\\cdot\\frac{\\vec{F}_j^\\text{(self)}}{F^\\text{(self)}}\\,\\boldsymbol{1}_{|\\vec{x}_i-\\vec{x}_j|\\in[r,r+\\Delta r]}\\right\\rangle}_{\\substack{\\text{average over the sum of all pair alignments}\\\\ \\vec{F}_i^\\text{(self)}\\vec{F}_j^\\text{(self)}/(F^\\text{(self)})^2\\text{ such that $|\\vec{x}_i-\\vec{x}_j|\\in[r,r+\\Delta r]$}}}$$\n",
    "which we have implemented below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: pair correlation tensor? but how justify mere r dependence?\n",
    "def pair_alignment_radial_binning(\n",
    "\t\tpositions: numpy.ndarray,\n",
    "\t\tself_propulsions: numpy.ndarray,\n",
    "\t\tbox_sidelength:float,\n",
    "\t\tbins:int=10\n",
    "\t\t) -> tuple[numpy.ndarray, list[float], list[float], list[float]]:\n",
    "\t'''\n",
    "\tComputes Alignment(r) and its error; also computes a\n",
    "\tdistance histogtam Distances(r) for reference\n",
    "\n",
    "\tArgs:\n",
    "\t\tpositions:\n",
    "\t\t\tnumpy.ndarray of shape (N,d) containing the\n",
    "\t\t\tpositions of N particles in d spatial dimensions\n",
    "\t\tself_propulsions:\n",
    "\t\t\tnumpy.ndarray of shape (N,d) containint the\n",
    "\t\t\tcomponents of the self propulsion force of\n",
    "\t\t\tN particles in d spatial dimensions\n",
    "\t\tbox_sidelength:\n",
    "\t\t\tdomain sidelength, the distance traveled before\n",
    "\t\t\twrapping around in periodic boundary conditions\n",
    "\t\tbins:\n",
    "\t\t\tthe number of radial bins (indexing the 'bags',\n",
    "\t\t\ti.e. arbitrarily sized lists of alignment values)\n",
    "\t\n",
    "\tReturns:\n",
    "\t\ttuple containing, in the following order:\n",
    "\t\t\t0. the radial bin lower edges\n",
    "\t\t\t1. a distance histogram (for reference)\n",
    "\t\t\t2. a list of alignment values\n",
    "\t\t\t3. a list of alignment value errors\n",
    "\t'''\n",
    "\tparticle_number, spatial_dimensions = positions.shape\n",
    "\n",
    "\t# normalized self-propulsions\n",
    "\tnormalized_self_propulsions = numpy.array([\n",
    "\t\tpropulsion / numpy.linalg.norm(propulsion)\n",
    "\t\tfor propulsion in self_propulsions\n",
    "\t])\n",
    "\n",
    "\t# binning from 0 to box_sidelength\n",
    "\thighest_meaningful_distance = box_sidelength / 2\n",
    "\tbin_lower_edges = numpy.linspace(1e-5, highest_meaningful_distance, bins)\n",
    "\talignment_bags = [[] for _ in range(bins)]\n",
    "\tdistance_histogram = [0] * bins\n",
    "\t\t\n",
    "\t# looping over all unique ...\n",
    "\tfor i in range(particle_number):\n",
    "\t\t# ... pairs of particles ...\n",
    "\t\tfor j in range(particle_number):\n",
    "\n",
    "\t\t\t# no measuring of self alignment\n",
    "\t\t\tif i == j:\n",
    "\t\t\t\tcontinue\n",
    "\n",
    "\t\t\t# for all periodic images\n",
    "\t\t\tfor pbc_image in pbc_iterator(spatial_dimensions):\n",
    "\n",
    "\t\t\t\t# distance of i and j squared, including periodic images\n",
    "\t\t\t\tseparation_ij = positions[j] - positions[i] + box_sidelength * numpy.array(pbc_image)\n",
    "\t\t\t\tr_ij = numpy.linalg.norm(separation_ij)\n",
    "\n",
    "\t\t\t\t# only consider alignment up until periodic boundary ...\n",
    "\t\t\t\t# ... box ends\n",
    "\t\t\t\tif r_ij > highest_meaningful_distance:\n",
    "\t\t\t\t\tcontinue\n",
    "\n",
    "\t\t\t\t# bin index\n",
    "\t\t\t\tr_ij_bin_index = int(bins * r_ij / highest_meaningful_distance)\n",
    "\n",
    "\t\t\t\t# appending the alignment \n",
    "\t\t\t\ta_ij = numpy.dot(normalized_self_propulsions[i], normalized_self_propulsions[j])\n",
    "\t\t\t\talignment_bags[r_ij_bin_index].append(a_ij)\n",
    "\n",
    "\t\t\t\t# increasing the counter of the distance histogram\n",
    "\t\t\t\tdistance_histogram[r_ij_bin_index] += 1\n",
    "\n",
    "\treturn (\n",
    "\t\tbin_lower_edges,\n",
    "\t\tdistance_histogram,\n",
    "\t\t[\n",
    "\t\t\tnumpy.sum(bag) if bag else 0 # only take average if the bag has content\n",
    "\t\t\tfor bag in alignment_bags \n",
    "\t\t]\n",
    "\t)\n",
    "\n",
    "\n",
    "def hypershell_volume(\n",
    "\t\tradii: numpy.ndarray,\n",
    "\t\tspatial_dimensions: float,\n",
    "\t\tshell_thickness: float\n",
    "\t):\n",
    "\t'''\n",
    "\tComputes the volume of a hypershell of given thickness\n",
    "\n",
    "\tArgs:\n",
    "\t\tradii:\n",
    "\t\t\tThe radii that hypershell volume should be\n",
    "\t\t\tcomputed\n",
    "\t\tspatial_dimensions:\n",
    "\t\t\tnumber of spatial dimensions\n",
    "\t\tshell_thickness:\n",
    "\t\t\tthe shell thickness dr\n",
    "\n",
    "\tReturns:\n",
    "\t\tnumpy.ndarray of hypershell volumes at the passed radii\n",
    "\t'''\n",
    "\t# d-dimensional surface\n",
    "\td = spatial_dimensions\n",
    "\tdr = shell_thickness\n",
    "\tunit_shell = 2*numpy.pi**(d/2) / euler_gamma_function(d/2)\n",
    "\treturn unit_shell * radii**(d-1) * dr\n",
    "\n",
    "\n",
    "def alignment_correlation_function(\n",
    "\t\tpositions: numpy.ndarray,\n",
    "\t\tself_propulsions: numpy.ndarray,\n",
    "\t\tbox_sidelength:float,\n",
    "\t\tbins:int=10\n",
    "\t):\n",
    "\t'''\n",
    "\tComputes the alignment correlation function a^(2)(r) and its error.\n",
    "\n",
    "\tArgs:\n",
    "\t\tpositions:\n",
    "\t\t\tnumpy.ndarray of shape (N,d) containing the\n",
    "\t\t\tpositions of N particles in d spatial dimensions\n",
    "\t\tself_propulsions:\n",
    "\t\t\tnumpy.ndarray of shape (N,d) containint the\n",
    "\t\t\tcomponents of the self propulsion force of\n",
    "\t\t\tN particles in d spatial dimensions\n",
    "\t\tbox_sidelength:\n",
    "\t\t\tdomain sidelength, the distance traveled before\n",
    "\t\t\twrapping around in periodic boundary conditions\n",
    "\t\tbins:\n",
    "\t\t\tthe number of radial bins (indexing the 'bags',\n",
    "\t\t\ti.e. arbitrarily sized lists of alignment values)\n",
    "\n",
    "\tReturns:\n",
    "\t\t\n",
    "\t'''\n",
    "\n",
    "\tbin_lower_edges, distance_histogram, alignment_sums = pair_alignment_radial_binning(\n",
    "\t\tpositions,\n",
    "\t\tself_propulsions,\n",
    "\t\tbox_sidelength,\n",
    "\t\tbins=bins\n",
    "\t)\n",
    "\n",
    "\t# bin size\n",
    "\tshell_thickness = bin_lower_edges[1] - bin_lower_edges[0]\n",
    "\n",
    "\t# unpacking particle number and number of ...\n",
    "\t# ... spatial dimensions of the domain\n",
    "\tN, d = positions.shape\n",
    "\n",
    "\thyper_shell_volumes = hypershell_volume(\n",
    "\t\tbin_lower_edges,\n",
    "\t\td,\n",
    "\t\tshell_thickness\n",
    "\t)\n",
    "\n",
    "\t# computing the correct prefactor\n",
    "\tL = box_sidelength\n",
    "\trho = N / L**d\n",
    "\n",
    "\t# number of tuples (i,j) considered\n",
    "\tpair_count = N * (N - 1)\n",
    "\n",
    "\treturn (\n",
    "\t\tbin_lower_edges,\n",
    "\t\tdistance_histogram / hyper_shell_volumes / pair_count,\n",
    "\t\talignment_sums / (N*d*rho * hyper_shell_volumes)\n",
    "\t)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dependence of $a^{(2)}(r)$ on packing fraction $\\phi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# which packing fraction to render\n",
    "PHI_STUDY = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_directory = f'alignment_study_datasets/phi={PHI_STUDY}'\n",
    "\n",
    "xyz_files = {\n",
    "\t# simulation index : file path\n",
    "\tint(file_name.split('_')[2]) : f'{data_directory}/{file_name}'\n",
    "\tfor file_name in listdir(data_directory)\n",
    "\tif file_name.endswith('.xyz')\n",
    "}\n",
    "\n",
    "alignment_correlations = []\n",
    "radial_distributions = []\n",
    "\n",
    "for i, path in xyz_files.items():\n",
    "\n",
    "\tif i == 0:\n",
    "\t\tcontinue\n",
    "    \n",
    "\tprint(f'({i}) Analyzing {path}')\n",
    "\n",
    "\tmetadata = active_particle_metadata(path)\n",
    "\n",
    "\tfor index, block in enumerate(XYZ.file_iterator(path)):\n",
    "\n",
    "\t\t# try not to get time correlation too involved. Also wait\n",
    "\t\t# ... until the simulation has become stationary.\n",
    "\t\tsampling_rate = 100\n",
    "\t\tburn_in_time = 1000\n",
    "\t\tif not (burn_in_time < index and index % sampling_rate == 0):\n",
    "\t\t\tcontinue\n",
    "\n",
    "\t\t# data unpacking. NOTE: assumes d=2\n",
    "\t\tpositions, self_propulsions = block.xyz_data[:,(0,1)], block.xyz_data[:,(2,3)]\n",
    "\n",
    "\t\tN, d = positions.shape\n",
    "\n",
    "\t\tr, g, a = alignment_correlation_function(\n",
    "\t\t\tpositions,\n",
    "\t\t\tself_propulsions,\n",
    "\t\t\tmetadata['box_sidelength'],\n",
    "\t\t\tbins=100\n",
    "\t\t)\n",
    "\t\t\n",
    "\t\talignment_correlations.append(a)\n",
    "\t\tradial_distributions.append(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a, g = numpy.mean(alignment_correlations, axis=0), numpy.mean(radial_distributions, axis=0)\n",
    "sigma = metadata['particle_diameter']\n",
    "\n",
    "def triangular_lattice_vector_length(j, k):\n",
    "\t'''\n",
    "\tReturns the length of triangular lattice vectors\n",
    "\tj * a_1 + k * a_2, where a_1 = (1,0) and a_2 is\n",
    "\trotated counter-clockwise by 120deg\n",
    "\t'''\n",
    "\ttri_a_1 = numpy.array([1, 0])\n",
    "\ttri_a_2 = numpy.array([numpy.cos(2*numpy.pi * 120/360), numpy.sin(2*numpy.pi * 120/360)])\n",
    "\n",
    "\treturn numpy.linalg.norm(j * tri_a_1 + k * tri_a_2)\n",
    "\n",
    "if PHI_STUDY >= 0.8:\n",
    "\tfor spike_index, spike in enumerate([\n",
    "\t\ttriangular_lattice_vector_length(1,0),\n",
    "\t\ttriangular_lattice_vector_length(1,1),\n",
    "\t\ttriangular_lattice_vector_length(2,1),\n",
    "\t\ttriangular_lattice_vector_length(2,2),\n",
    "\t\ttriangular_lattice_vector_length(2,3),\n",
    "\t\ttriangular_lattice_vector_length(3,3)\n",
    "\t]):\n",
    "\t\tplt.axvline(\n",
    "\t\t\tx=spike,\n",
    "\t\t\tcolor='gray',\n",
    "\t\t\tlinestyle='-.',\n",
    "\t\t\tlabel=r'triangular lattice $g^{(2)}(r)$ peak locations' if spike_index == 0 else None\n",
    "\t\t)\n",
    "\n",
    "plt.plot(r / sigma, (N*d)*a, label=r'$Nd\\cdot a^{(2)}(r)$', color='black')\n",
    "plt.plot(r / sigma, g, label=r'$g^{(2)}(r)$', linestyle='dashed', color='black')\n",
    "\n",
    "plt.axhline(y=1, color='black', linestyle='dotted')\n",
    "plt.axhline(y=0, color='black', linestyle='dotted')\n",
    "plt.axhline(y=-1, color='black', linestyle='dotted', label='y-value $-1$, $0$ and $1$ respectively')\n",
    "#plt.axvline(x=1, color='black', linestyle='dotted')\n",
    "\n",
    "plt.xlabel(r'$r/\\sigma$')\n",
    "\n",
    "plt.title(fr'Alignment distribution and radial distribution in $d=2$ at $\\phi={PHI_STUDY}$')\n",
    "#plt.grid()\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig(f'alignment_distribution_phi={PHI_STUDY}.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dependence of $a^{(2)}(r)$ on mass $m$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "masses = [0.1, 1.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alignment_correlations = {}\n",
    "radial_distributions = {}\n",
    "\n",
    "for m in masses:\n",
    "\n",
    "\tdata_directory = f'alignment_study_datasets/m={m}_phi=0.1'\n",
    "\n",
    "\txyz_files = {\n",
    "\t\t# simulation index : file path\n",
    "\t\tint(file_name.split('_')[3]) : f'{data_directory}/{file_name}'\n",
    "\t\tfor file_name in listdir(data_directory)\n",
    "\t\tif file_name.endswith('.xyz')\n",
    "\t}\n",
    "\n",
    "\talignment_correlations[m] = []\n",
    "\tradial_distributions[m] = []\n",
    "\n",
    "\tfor i, path in xyz_files.items():\n",
    "\n",
    "\t\tif i == 0:\n",
    "\t\t\tcontinue\n",
    "\t\t\n",
    "\t\tprint(f'({i}) Analyzing {path}')\n",
    "\n",
    "\t\tmetadata = active_particle_metadata(path)\n",
    "\n",
    "\t\tfor index, block in enumerate(XYZ.file_iterator(path)):\n",
    "\n",
    "\t\t\t# try not to get time correlation too involved. Also wait\n",
    "\t\t\t# ... until the simulation has become stationary.\n",
    "\t\t\tsampling_rate = 25\n",
    "\t\t\tburn_in_time = 1000\n",
    "\t\t\tif not (burn_in_time < index and index % sampling_rate == 0):\n",
    "\t\t\t\tcontinue\n",
    "\n",
    "\t\t\t# data unpacking. NOTE: assumes d=2\n",
    "\t\t\tpositions, self_propulsions = block.xyz_data[:,(0,1)], block.xyz_data[:,(2,3)]\n",
    "\n",
    "\t\t\tN, d = positions.shape\n",
    "\n",
    "\t\t\tr, g, a = alignment_correlation_function(\n",
    "\t\t\t\tpositions,\n",
    "\t\t\t\tself_propulsions,\n",
    "\t\t\t\tmetadata['box_sidelength'],\n",
    "\t\t\t\tbins=100\n",
    "\t\t\t)\n",
    "\t\t\t\n",
    "\t\t\talignment_correlations[m].append(a)\n",
    "\t\t\tradial_distributions[m].append(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for m, alignment_correlation_functions in alignment_correlations.items():\n",
    "\n",
    "\ta = numpy.mean(alignment_correlation_functions, axis=0)\n",
    "\t\n",
    "\tsigma = metadata['particle_diameter']\n",
    "\n",
    "\tplt.plot(\n",
    "\t\tr / sigma,\n",
    "\t\t(N*d)*a,\n",
    "\t\tlabel=fr'$m={m}$',\n",
    "\t\tcolor='Red' if m == 1 else 'Blue'\n",
    "\t)\n",
    "\t#plt.plot(r / sigma, g, label=r'$g^{(2)}(r)$', linestyle='dashed', color='black')\n",
    "\n",
    "plt.axhline(y=0, color='black', linestyle='dotted')\n",
    "plt.axhline(y=-1, color='black', linestyle='dotted', label='y-value $-1$ and $0$ respectively')\n",
    "\n",
    "#plt.ylim(-6, 5)\n",
    "\n",
    "plt.xlabel(r'$r/\\sigma$')\n",
    "\n",
    "plt.title(r'Alignment distribution function $Nd\\cdot a^{(2)}(r)$ in $d=2$ at $\\phi=0.1$')\n",
    "#plt.grid()\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig(f'alignment_distribution_for_various_masses_phi=0.1.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}