phase_separation/phase_separation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "special-karen",
   "metadata": {},
   "source": [
    "| <img src=\"mixed.png\" alt=\"mixed.png\" width=\"300\"/> | <img src=\"chi=5_final.png\" alt=\"chi=5_final.png\" width=\"300\"/> |  <img src=\"chi=40_final.png\" alt=\"chi=40_final.png\" width=\"300\"/> | \n",
    "|:-:|:-:|:-:|\n",
    "| Randomly mixture, `0 ns` | &chi;=5 kJ/mol final structure, `0.1 ns` | &chi;=40 kJ/mol final structure, `0.1 ns` |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "parental-phase",
   "metadata": {},
   "source": [
    "### Example: Phase separation\n",
    "We setup a system of randomly dispersed `\"A\"` and `\"B\"` type particles, in equal proportions, run a short simulation, and examine the resulting trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "congressional-immune",
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def setup_system(box_size, number_density, out_path):\n",
    "    \"\"\"Setup a simple system of independent particles with no bonds\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    box_size : list of float\n",
    "        Simulation box size in nm.\n",
    "    number_density : float\n",
    "        Number density per nm squared of particles in the simulation box.\n",
    "    out_path : str\n",
    "        Path of the created structure file.\n",
    "    \"\"\"\n",
    "    box_size = np.asarray(box_size, dtype=np.float32)\n",
    "    volume = np.prod(box_size)\n",
    "    n_particles = int(number_density * volume)\n",
    "    n_particles -= np.mod(n_particles, 2)\n",
    "    \n",
    "    with h5py.File(out_path, \"w\") as out_file:\n",
    "        position_dataset = out_file.create_dataset(\n",
    "            \"coordinates\",\n",
    "            (1, n_particles, box_size.size,),\n",
    "            dtype=\"float32\",\n",
    "        )\n",
    "        types_dataset = out_file.create_dataset(\n",
    "            \"types\",\n",
    "            (n_particles,),\n",
    "            dtype=\"i\",\n",
    "        )\n",
    "        names_dataset = out_file.create_dataset(\n",
    "            \"names\",\n",
    "            (n_particles,),\n",
    "            dtype=\"S10\",\n",
    "        )\n",
    "        indices_dataset = out_file.create_dataset(\n",
    "            \"indices\",\n",
    "            (n_particles,),\n",
    "            dtype=\"i\",\n",
    "        )\n",
    "        \n",
    "        types_dataset[:n_particles//2] = 0\n",
    "        types_dataset[n_particles//2:] = 1\n",
    "        names_dataset[:n_particles//2] = np.string_(\"A\")\n",
    "        names_dataset[n_particles//2:] = np.string_(\"B\")\n",
    "        indices_dataset[...] = np.arange(n_particles)\n",
    "        position_dataset[...] = np.random.uniform(\n",
    "            low=0.0, high=box_size, size=(1, n_particles, box_size.size),\n",
    "        )\n",
    "\n",
    "# 1040 particles in a 5nm x 5nm x 5nm simulation box.\n",
    "setup_system([5, 5, 5], 8.333, \"phase_separation.HDF5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "scheduled-newfoundland",
   "metadata": {},
   "source": [
    "### Run the simulation\n",
    "We run one simulation with &chi; less than the critical point of phase separation, and one for &chi; well over: \n",
    "```bash\n",
    "# chi=5.0\n",
    "mpirun -n 6 python3 ${HYMD}                    \\\n",
    "                    phase_separation.toml      \\\n",
    "                    phase_separation.HDF5      \\\n",
    "                    --seed 1                   \\\n",
    "                    --verbose 1\n",
    "mv sim.H5 sim_chi=5.H5\n",
    "\n",
    "# chi=40.0\n",
    "gsed -i \"17s/^.*$/  ['A', 'B', 40.0],/\" phase_separation.toml\n",
    "mpirun -n 6 python3 ${HYMD}                    \\\n",
    "                    phase_separation.toml      \\\n",
    "                    phase_separation.HDF5      \\\n",
    "                    --seed 1                   \\\n",
    "                    --verbose 1\n",
    "mv sim.H5 sim_chi=40.H5\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "aboriginal-magnitude",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABBGUlEQVR4nO3dd3hUVf7H8feXFEIIBEhCL6EXQXq1YF1RXFi7+LN3d921rbu6xbrr7rrF1dW1rShYQMSyqKxYUVFEQKlBpPcSAiSQXs7vjxl0iIEUkrkzN5/X88yTzJ0zM98MyeFzzj33XnPOISIiIuInDbwuQERERKS2KeCIiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4EitMLPnzMxVcOtVhectqKTNbDObXsH2YWZWZGbJR1p/ZapSp4jUHTMbH+xTfvB3aGbJZvasme0xs2wze9HMUip5vfTg6515mDYnBNv0reCxmWb2YM1+mqqrSp1SsVivCxBf+Qa4oty29XX4fmOBOc657Dp8DxHxmJklAA8BOw7RZBrQA7gaKAP+ArwBHFdH9SQCJwJ/rovXl9qhgCO1Kdc590UY328sMCWM7yci3rgd2AKsAQ6aTTGzkcCPgNHOuU+C27YA88zsFOfc+3VQz0lAAfB5Hby21BLtovKJ4PTsbjN72MwsZPsrZrY5OAKKeGYWb2avmdlGM+t2mHatgUHA28H795jZLjMbbmYLzCzfzOaYWWcza2lmb5jZfjNbYWYnlXutmODzN5pZoZktN7OL6vYnFYl8kdCvmFlH4FfATYdocjqw40C4AXDOfQmsCz5Wnfc60cz2mdkDlTQdC8xyzpUEn+fM7BYz+7uZZQX7ol8GH7vMzNaa2V4zm1j+MzOzAWb2gZnlBXexvWhmrapTt1RMAcc/HgD+BvwC+AmAmQ0EzgHud84VHOqJZhZbhZsd6vkh+phZTjAkzDGz0dX5AYJ/+K8D/YHjnHOrD9P8DGCdc+6bkG2JwFMEprInAB2B5wnM8swBziYwCnwlOMV8wH3Ab4PPHQd8BrxoZhOqU7+ID0VCv/J3YJpz7qtDPN6LwO7x8lYEH6sSMzuNwIDpQefcbyppPjbYNtRtQBKBvucl4K/BNTqXE/j8fgP8H3BzyHumAbMJ9F0XAT8HRgPvmVl8VWuXQ3DO6eajGzAdeCP4/dsEpnTjKnmOq8Lt8kpe4ybgBgJ/nOcCc4EiYFglz3sOWEDgD/x9Ah1Vu3JtZgPTy217FXgk5P49wTpHh2z7aXDbXSHb+gS3nR683wLIBe4u9/ozgZXl6/T631c33by4edivnARkA62C93/wdwi8d6C2cttfAD4/zGunB2s4k8DApgD4Zbk2JwTb9A3Z1g8oBdLK/awfhdxvAGwD9gBNQ7ZPA+aF3P8zsLdcm+HB15tQvk6vfw+i7aY1OP4zCZgeXHF/BnCpc664kucMrcLrrjvcg865h0Pvm9lMYDmBUctPKnntxsA7QDMCAeVQCwkPvHYccCpwXrmHioBPQ+4fmAH6sIJt7YJf+xIIV6+Ue62XgefMLM05l1lJ/SJ+F/Z+xcxigUeAP1bWJxyhcwjMntzqnHusCu3HAl9W0C98cOAb51yZma0D8pxzOSFtVgOjQu4PA94NbeOcm2dm64Fj0RrDI6KA4z//IzDimQJkAC9W4TmLqtCmtDpFOOfygiHnx1Vo3hboTtU7suMJjJBml9u+zzlXFnK/KPh1b0hdRcFZ8QP7wdsEv5Z/3wP3WwAKOFLfedGvXAMkExhoNAtuiwdigvdzgyFrD5BWwfObBx+rzDhgN4Hd41VR0e4pCOlngooOsS10DU4bAgPB8nYQ6HvkCGgNjs+4wKK3lwnsC/59uf/wD6W4CrfLalJO8FaZVQQOL/+dmd1QhfZjgQ+cc4U1qKm8bcGvLcttP7DIb3ctvIdIVPOoX+kJtCfwn/2e4G0CMCD4/QXBdt9Q8VqbQ63NKe/nBNbmvWuVnzunOTCSigNOTWzjh30PBPof9T1HSDM4/jQk+HVJFdsf8S6q8sysEYEgsrAq7Z1zz5tZEvCome1zzr1wmOZjCSw8rA3LgDwCu7vuC9l+PvCtdk+JfCfc/cqjBM5lE+oOoDNwHYFFxBCYXfq9mR3rnJsDYGZDgC7BxyqTA5wGfAzMMrOTyu1WCjWGwBFbX1fhdatiHnCDmTVxzu0DMLOhBNbdzKml96i3FHB8JriPfASBP9pLgLsre45z7ojO0GuBMwm/RWBR32ogFbiFwK6n8utkDlfH48GQ86yZ7XfOvRH6cPC9uhE4odfMI6k55D13m9k/CcwelRBY8Hw2gXUGOopKBG/6FRc4ivKgIynN7HIg1Tk3O6TdXDN7F5gcPDT7wIn+5rgqngPHOZdlZqcSWMP3lpmNcc7lhTYJfh1LLfU9Qf8gcHDGLDP7C4EZsj8DSwkcSCFHQLuofCR4yOX9BE4+9Wfg4ioehnmkCgmsU/kdgT/+pwjsex5d3U7OOfdX4E/A1GCHA9CI79fTjAUWO+c210LdB9wVfM8bCAS144GLnXNTa/E9RKKSh/1KdVxAYAZmIjCZwMzxWdV5AefcNuBkArMnrwUP024UfLjIzBoQmMGprd1TBGeITyRwBNcU4DECIetU51zR4Z4rlbPgYWjiA2Z2HoHDEE8mMPJZD5zonPvYy7qORLBT2QZMdM7dGRypzXfO/dbj0kTqBT/2K1VlZjcCDwNNCKz9+QhIcc7t97IuqRoFHJ8IBoFlwHbn3EnBbR8QOBz6LmBmtP1RmtnxBM6p83NghHNunsclidQrfuxXqsLM2hM4H80/gAznXLXOiCyRQbuo/ONioDeB3UQH/BzYR2BtzCEvexDBPgBOAS5RuBHxhB/7laq4GngSmE/NjiCVCKAZHBEREfEdzeCIiIiI7/jiMPHU1FSXnp7udRki9crChQt3OecqOoOsb6mvEQm/mvY1vgg46enpLFhwRKdcEJFqMrMNXtcQbuprRMKvpn2NdlGJiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4IiIiIjvKOCIiIiI7yjgiIiIiO8o4IiIiIjvKODUsdIy53UJIiIi9Y4CTh0pLi3jXx+s4qi73+HmqV+zPbvA65JERETqjdhwvpmZTQTOBHY65/pW8LgBDwNnAHnA5c65r8JZY21YvjWb219ZQsa2HIZ1bsHMpdt5N2MHPzuxG1cd25mEuBivS6xz+wtLWLNzP6t37qdTSiJD0lt4XVK15ReVsm5XLut25bI2cz/rduWyPaeAlKSGtE1OoHVyAm2SE2id3Ig2yQmkJjUkpoF5XbZQf/oaETm0sAYc4DngUWDyIR4/HegevA0HHg9+jQqFJaU89uFq/j17Dc0S43ni4sGM6duajVl5/OHtDP46ayUvz9/Eb8f25kd9WhHoY6Nbdl4xqzP3sWrHflbtDNzW7NzPlr35B7U7plsKt57ag8GdIivolJY5tuzJZ+2u/azNDIaZXftZl5nL1nKzbgdCzdLNe3l3eQGFJWUHPR7bwGjVNNCmdXICbYLft2vWiE4pjUlPTSQxvnb/5JxzZO4vZENWHkUlZbRq2pBWTRNIahjri9+vI/AcPu5rRKRyYQ04zrlPzCz9ME3GA5Odcw74wsyamVkb59y28FRYc4s37eX26Yv5dsd+zh7Yjt+f2YfmjeMB6JiSyFOXDuHTVZnc+2YG1z2/kOO6p3LXmX3o3qqJx5VX3eY9ecxemcnqnftZtTMQanbuK/zu8YS4BnRNS2JoenMuatWRbi2T6JrWmNkrM3ni4zWc8/hcRvdI45ZTezCgQzNPfgbnHIs3Z/POsu3MXrmTtZm5FJV+H1SaJMTSJS2JEV1S6JzamC5pSXRObUzn1MY0io856HX25BWzLTuf7dkFbMsu+P5rTj4rtubw4Yqd5BeXHvT+rZo2JD2lMV3SGpOe0pj04Gt3bJF4yJk95xw79xWyblcuG7JyWZ+VF/i6K/A1t6j0B89JjI+hVdMEWjYJBJ4Dwadl0wRafbct4aCfyU/83NeISNWEewanMu2ATSH3Nwe3/aDTMbNrgWsBOnbsGJbiKlJQXMpD73/L05+sJa1JQ565bAgn925VYdvjuqfxv5uO4/m5G3jo/W8Z8/CnXDKiE7ec0oPkxLgwV151zjmmfLmJ+9/KIL+4lMbxMXRr1YTje6TRvWUS3Vsl0b1lE9o1a0SDCnbRdGvZhIuGd2Ty3A08+fEafvLYZ5zSuyU3n9KDvu2S67z+ktIy5q/fw6zl25m1fDvbsguIbWCM6JLC6J5pdE1NonNaIGikNI6v0syHmdGicTwtGsdzVNuKfwbnHDn5JWzem8f6XXmsz8r9bpfXu8t3kJVbFPJ60Da5EempiaSnNKZxw1g2ZOWyISuPDVl5BwWluBijQ/NE0lMbM7xLC9JTGtMpJZGGsTHs3FfAzpxCtucUsCMn8P3izXvZnv3DGSeACcM68MBZ/erjbE+V+xoRiU6RFnCqzDn3FPAUwJAhQzw5VGnhht3cPn0JazNzuWBIB34ztjfJjQ4fVOJiGnDlsZ0ZP6Atf3/vWybNXc+MxVu57Uc9uHBox4hbw5G5r5A7Xl3CB9/s5JhuKdw3vi9dUhtX+z/ExPhYrh/dlYtHdOK5z9bx1CdrOfNfcxhzVGtuPrU7vVo3rdW6C0tK+XxNFu8s3c57K3awO7eIhrENGN0jjdtP68nJvVrVeag0M5IT40hOTK4wBOUUFLN+1/ehZ/2uXNZl5fHWkm3kF5XSMSWR9JREjumWSnpqY9JTAuGnTXICsTHVOz7AOUdOQQk7cwrYkVPIjpwC5q7NYsqXm+jfvhkXDvNukBDpImUwJSLVY4EZ2jC+YWDa+K1DLPx7EpjtnJsSvL8SOKGyaeMhQ4a4BQsW1EW5FcovKuWvs1by7OfraJvciD+f04/juqfV6LWWb83m3hkZfLl+N33aNOWecUcxrHNkrFN5L2MHd7y6hH2FJdwxpheXj0qvcIamJnIKinnm03VMnLOOfYUljD26Dbec0p1uLWu+yy6vqISPV2byzvLtfLhiJ/sKS0hqGMtJvVoypm9rTuiZVutrYOqKc67OZ1XKyhyXTvyS+et3M+PGY+nZunqfvZktdM4NqaPyjpgf+hoRqXlfE2kBZyxwI4EjG4YDjzjnhlX2muHsdOauyeKO15awISuPS0Z04ten9yKp4ZH9p+mc460l23hg5gq2ZRdwet/WjOqa8t36j9ZNE2otWFRFbmEJ97+VwdT5m+jTpin/vHAAPepordDevCKe/nQtz362nvziUsb3b8tNp/Sgc2pjnHPkFpWSnV9MTn7xwV8LSr67n5NfTOb+Qr5ct5vCkjJaNI7n1N6tGNO3NaO6pdAw1p/rTGrDzn0FnPHwHJolxjHjxmOqFQCjPOBEfF8jIgFREXDMbApwApAK7ADuBuIAnHNPBA/dfBQYQ+DQzSucc5X2JuHodLbszeexj1bz0ryNdGyRyF/OOZqRXVNq9T3yikp44uO1TJyzjv2FJd9tT4hr8N3C1MCC10Dw6ZLa+LuFzLVl4YY93PLyIjbtyeP60V255ZQexMfW/emSsvYX8tQna5k0dz3FpY6mCbHkFJRUeqLEJgmxNE2II7lRHMM6t+C0o1ozNL15tXfh1GdzVu3ikonzOG9wex48t3+VnxfJASea+xoROVhUBJy6Upedzjfbc3jq47XMWLwVB1w2Mp1fntajTnd1OOfYkVMYOFx5Vy7rMr9fp7Fxdx4lIf/pN0+Mo3NqY7q1TGJ45xRGdk2hbbNG1X7P4tIyHvlgFY99tJq2zRrxj/MHeLKrLHNfIZPnric7v/i74NK0UWzI98GvCXEkJcRG3JqlaPX3d1fyrw9X89AF/TlrYPsqPSeSA05dUcARCb+a9jXRsSAhzJxzzF2bxZMfr+XjbzNJjI/h0pHpXHlsOu2bJ9b5+5vZd+dSGdU19aDHikvL2LQ77/sT0AUD0HsZO5i2YDMA6SmJjOyaysiuKYzskkJak4aHfb/VO/dzy8uLWLolm3MHt+fuH/ehSYI3R3WlNWnIbT/q6cl712c3ndydeWt389vXl9G/fTO6pCV5XZKIyBFRwAlRWuZ4Z9l2nvxkDUs2Z5OaFM/tp/Xk4uGdIuYw7riYBnRJS/rBf0BlZY6VO/bx+Zos5q7J4q3FW5ny5UYAurdMYlTXwOzOiC4pNEsM7NZyzjF57gYemLmCxPgYnrh4EGP6tgn7zyTei41pwMMTBnDGw5/ys5e+5vWfjqoXZ9wWEf/SLioCR0VNX7iJpz9dx8bdeXRObcw1x3Xh7EHtoraTLyktY/nWHOauzeLzNVnMX7eb/OJSzKB366aM6prCtzv388m3mYzukcZfzz2alk0TvC5bPPbhNzu48rkFXDKiE/f/5Adrcw+iXVQiEg7aRVUDe3KLmDx3A5Pmrmd3bhEDOzbjN2f05tQ+raJ+bUdsTAP6d2hG/w7NuH50V4pKyliyee93MzyTv9hAA4P7xx/FxSM61ccTvUkFTurVimuO68zTn65jVNcUTu+nGT0RiU71MuCUljn+9u5KngsemnxK75ZcN7orQzo19+1/9PGxDRiS3oIh6S34xcndKSgupaTMHfEh7uI/t5/Wiy/X7+FXry6hb7tkOrSo+3VnIiK1rd4dS1tUUsYvpnzN47PXMKZva9675Xj+c9lQhqa38G24qUhCXIzCjVQoPrYBj04YCMCNU76mqIJLPIiIRLp6FXDyi0q5ZvIC3l66jd+N7c1DFwyIqotdioRLhxaJPHjO0SzetJe/zvrG63JERKqt3gScnIJiLp04j09XZfKXc/px9XFdvC5JJKKd3q8Nl4zoxNOfruODFTu8LkdEpFrqRcDZtb+QCU99waJNe/nXhEFcMFQXzBOpit+O7U3vNk257ZXFbMvO97ocEZEq833A2bo3n/OfnMuazP08fekQxh6to0JEqiohLobHLhpIUUkZN01ZREmp1uOISHTwdcBZm7mf856YS2ZOIc9fNZwTerb0uiSRqNMlLYkHzurHl+t38/AHq7wuR0SkSnwbcDK25nD+k3MpKC5lyrUjGJoe/usqifjFTwa247zB7Xn0o9V8vmaX1+WIiI8558gvKj3i1/HlccILN+zmimfn07hhLC9cPZyuuq6OyBG7d/xRfLQykxfnbfzBNdJERGrKOcf6rDzmrsnii7WB2978YmbceAy9Wjet8ev6LuB8uiqTaycvpHVyAs9fNSwsF8cUqQ8S42M5vkcqH32zk7Ky6L/Ei4h4wznHxt15fLE2KxhqdrM9pwAIXHB5ZJcUPlmVyd3/Xc7Ua0fU+H18FXDeWbaNX0xZRNeWSUy+clilV9EWkeo5tlsqr321hYxtOV6XIiJRZNPuPOauzeKL4CzN1uxAoElNimdEl8CFoEd2TaFLamPMjBe+2MDv3ljG20u31fg9fRNwXlmwiV+/uoQBHZrx7OXDIubq3yJ+cky3wK6pz1ZrHY6IVG71zn1cM3kh63blAtCicTwjurTghmCg6ZqWVOFVBCYM68hL8zbyx7dX1Pi9fRFwdu0v5PbpSziueypPXjKYxHhf/FgiEadV0wR6tEpijgKOiFSiqKSMm6YuIju/mHvHHcWILin0aFVxoCkvpoFx7/ijOO+JuTV+f18cRbW/oITT+7bmP5cNUbgRqWPHdEvly3W7vS5DRCLcwx98y/KtOfzp7H5cNiqdnq2bVOuaj0PTW/CTAW1r/P6+CDidUhrzrwkDaRgb43UpIr53bLdUCnUBThE5jAXrd/P47DWcP6Q9px3Vusavc+cZvWv8XF8EHDOIjfHFjyIS8YZ3SSG2QdVHYSJSv+wvLOHWaYtp17wRd/34qCN6rVZNE2r8XKUCEamWpIaxDOzYzOsyRCRC3f9mBpv35PHQ+QNIaujdshEFHBGptgNHU4mIhHp3+XZeXrCJ60d3ZYjHVxBQwBGRajtWAUdEysncV8idry3lqLZNufmUHl6X44/DxEUkvPp3aOZ1CSISQZxz3PnaEvYVljD1ggHEx3o/f+J9BSISdeK0qF9EQkydv4n3V+zkjjG96N6qidflAAo4IiIicgTW78rl/rcyOKZbCpePSve6nO8o4IiIiEiNlJSWceu0RcQ2MP52Xn8aRNApJLQGR0RERGrkiY/X8NXGvTwyYSBtkht5Xc5BNIMjIiIi1bZk817++f4qxvVvy7j+Nb+kQl1RwBEREZFqyS8q5ZaXF5HWpCH3j+/rdTkV0i4qERERqZa/vPMNazJzefHq4SQnxnldToU0gyMiIiJV9sm3mTz3+XquPKZzRJ/VXAFHREREqmRvXhG3T19M95ZJ/GpMT6/LOSztohIREZFKOef47RvL2J1bxDOXDSUhLsbrkg5LMzgiIiJSqde/3sLbS7Zx8yk96Nsu2etyKqWAIyIiIoe1aXced/13OcPSW3D96K5el1MlCjgiIiJySCWlZdzy8iLM4B8X9Ccmgs5WfDhagyMiIiKH9PjsNSzYsIeHLxxA++aJXpdTZZrBERERkQp9vXEP//xgFeMHtGX8gHZel1MtCjgiIiLyA7mFJdzy8iJaN03gvgg9W/HhaBeViIiI/MB9b2awYXceU68ZQXKjyDxb8eFoBkdEREQO8s6y7by8YBM3jO7K8C4pXpdTIwo4IiIi8p0dOQXc8doSjm6fzM2n9PC6nBpTwBEREREAysocv3xlMYXFZfzzggHEx0ZvTIjeykVERKRWTfxsHZ+u2sXvz+xDl7Qkr8s5Igo4IiIiwoptOTz4zkpO7dOKCcM6eF3OEVPAERERqecKiku5eeoikhPj+PPZ/TCLjrMVH44OExcREann/vLON6zcsY/nrhhKSlJDr8upFZrBERERqcc++TaTZz9bz+Wj0jmhZ0uvy6k1CjgiIiL11O7cIm57ZTE9WiVxx+m9vC6nVmkXlYiISD3knOPXry4hO6+YSVcMIyEuxuuSapVmcEREROqhqfM38V7GDn41pid92jb1upxap4AjIiJSzyzbks29by7n2G6pXHlMZ6/LqRNhDzhmNsbMVprZajO7o4LHO5rZR2b2tZktMbMzwl2jiIiIX23PLuCqSfNJadyQhy4YQIMG0X9IeEXCGnDMLAZ4DDgd6ANMMLM+5Zr9DpjmnBsIXAj8O5w1iog/aDAl8kN5RSVcNWk++wtKeObyIaQ18cch4RUJ9wzOMGC1c26tc64ImAqML9fGAQd2BiYDW8NYn4j4gAZTIj9UVua4eeoiVmzL4dGLBtGrtf/W3YQKd8BpB2wKub85uC3UPcDFZrYZmAn8vKIXMrNrzWyBmS3IzMysi1pFJHppMCVSzl9mfcO7GTv4/Zl9OLGXf853cyiRuMh4AvCcc649cAbwvJn9oE7n3FPOuSHOuSFpaWlhL1JEIpoGUyIhps3fxJMfr+XiER25fFS61+WERbgDzhYg9Ape7YPbQl0FTANwzs0FEoDUsFQnIvWJBlNSL8xdk8VvXl/Kcd1TufvHR/niOlNVEe6AMx/obmadzSyewH7vGeXabAROBjCz3gQCjoZNIlIdGkyJAOt25XLDiwtJT23MoxcNIi4mEnfc1I2w/qTOuRLgRmAWsILAAr/lZnafmY0LNrsNuMbMFgNTgMudcy6cdYpI1NNgSuq9vXlFXPXcfAyYeNlQkhvFeV1SWIX9Ug3OuZkE9neHbrsr5PsM4Jhw1yUi/uGcKzGzA4OpGGDigcEUsMA5N4PAYOppM7uFwIJjDabEN4pLy7jhha/YvCefF68ZTseURK9LCjtdi0pEfEmDKamvnHP8/o1lzF2bxT/O78/Q9BZel+SJ+rMzTkREpB74z6frmDp/Ezee2I2zB7X3uhzPKOCIiIj4xLvLt/PA/1Ywtl8bbj21h9fleEoBR0RExAeWbcnmpqmLOLpdMn87r79vrzFVVQo4IiIiUW5HTgFXT1pA88Q4nr50CI3iY7wuyXNaZCwiIhLFikvLuO75heQUFDP9+lG0bJrgdUkRQQFHREQkiv3rg1Us2rSXxy4aRJ+2/r6AZnVoF5WIiEiUWrhhD49+tJqzB7Vj7NFtvC4noijgiIiIRKHcwhJunbaINsmNuGfcUV6XE3G0i0pERCQK/eHtDDbuzuPla0fSNKF+XYahKjSDIyIiEmXey9jBlC83cd3xXRnWuX6eqbgyCjgiIiJRJHNfIXe8uoTebZrW+5P5HY52UYmIiEQJ5xx3vraEfYUlTLlwAPGxmqc4FH0yIiIiUWLq/E28v2Invx7Tix6tmnhdTkRTwBEREYkC63flcv9bGRzTLYUrRqV7XU7EU8ARERGJcCWlZdz88iJiG5iuM1VFWoMjIiIS4f49ew2LNu3lkQkDaZPcyOtyooJmcERERCLY4k17efiDVYwf0JZx/dt6XU7UUMARERGJUHlFJdzy8iJaNWnIfeP7el1OVNEuKhERkQj1wMwVrN2Vy0vXDCe5kc5WXB2awREREYlAH63cyQtfbOTqYzszqmuq1+VEHQUcERGRCLM7t4hfTV9Cr9ZN+OVpPb0uJyppF5WIiEgEOXC24uy8YiZfOYyEuBivS4pKmsERERGJINMXbmbW8h388rQe9G7T1OtyopYCjoiISITYua+A+97KYFjnFlx9bBevy4lqCjgiIiIR4r43MygsKePPZ/fT2YqPkAKOiIhIBPho5U7eWrKNG0/sRpe0JK/LiXoKOCIiIh7LKyrhd68vo1vLJK4f3dXrcnxBR1GJiIh47J/vr2LL3nymXTeS+FjNPdQGfYoiIiIeWr41m2fmrOPCoR0Y1rmF1+X4hgKOiIiIR0rLHHe+tpTmiXHceXpvr8vxFQUcERERj0yeu54lm7P5/Zl9SE7UtaZqkwKOiIiIB7buzedvs1ZyfI80xvVv63U5vqOAIyIi4oF7Ziyn1Dn++JO+mOmcN7VNAUdERCTMZi3fzrsZO7j5lB50aJHodTm+pIAjIiISRvsKirn7v8vp1boJVx3b2etyfEvnwREREQmjv7/7LTv2FfD4xYOIi9E8Q13RJysiIhImizbtZdLc9VwyohMDOzb3uhxfU8AREREJg+LSMu58bSktmzTk9tN6el2O72kXlYiISBhMnLOOFdtyeOLiQTRJ0Dlv6ppmcEREROrYpt15PPT+t5zSuyWnHdXa63LqBQUcERGROuSc43dvLKOBGfeO1zlvwkUBR0REpA69tWQbH3+byW0/6km7Zo28LqfeUMARERGpI9l5xdz7Zgb92iVz+ah0r8upV7TIWEREpA6UlTnufH0Ju3MLee6KocQ00K6pcNIMjoiISC1zznHvm8uZuXQ7vx7Ti77tkr0uqd5RwBEREallj364mklzN3D1sZ259vguXpdTLyngiIiI1KIpX27k7+99y1kD2/GbM3rrqCmPKOCIiIjUkneWbee3ry/lhJ5pPHju0TTQuhvPKOCIiIjUgi/WZvGLqV9zdPtm/Pv/dCFNr+nTFxEROUIZW3O4ZtICOjRvxLOXDyUxXgcpe00BR0RE5Ahs2p3HZc9+SVJCLJOvGk7zxvFelyQo4IiIiNTYrv2FXPLMPIpKyph05TCdqTiChD3gmNkYM1tpZqvN7I5DtDnfzDLMbLmZvRTuGkVERCqzv7CEK56dz/acAiZePoQerZp4XZKECGvAMbMY4DHgdKAPMMHM+pRr0x24EzjGOXcUcHM4axSR6KeBlNS1wpJSrn9+IRnbcvj3/w1icKcWXpck5YR7FdQwYLVzbi2AmU0FxgMZIW2uAR5zzu0BcM7tDHONIhLFQgZSpwKbgflmNsM5lxHSJnQgtcfMWnpTrUSjsjLHbdMWM2f1Lv52Xn9O6tXK65KkAuHeRdUO2BRyf3NwW6geQA8z+8zMvjCzMRW9kJlda2YLzGxBZmZmHZUrIlHou4GUc64IODCQCqWBlNTIgUswvLVkG3ec3otzB7f3uiQ5hEhcZBwLdAdOACYAT5tZs/KNnHNPOeeGOOeGpKWlhbdCEYlktTaQAg2m5GCPffT9JRiu0yUYIlq4A84WoEPI/fbBbaE2AzOcc8XOuXXAtwQCj4hIbanSQAo0mJLv/XfRFv72ri7BEC3CHXDmA93NrLOZxQMXAjPKtXmDQKeDmaUSGGmtDWONIhLdNJCSWrdlbz6/e30ZQ9Ob6xIMUSKsAcc5VwLcCMwCVgDTnHPLzew+MxsXbDYLyDKzDOAj4HbnXFY46xSRqKaBlNSqsjLH7a8spsw5/n7eAF2CIUqE/VzSzrmZwMxy2+4K+d4BtwZvIiLV4pwrMbMDA6kYYOKBgRSwwDk3I/jYj4IDqVI0kJLDmDR3PZ+vyeLPZ/ejY0qi1+VIFeliGSLiOxpISW1ZvXM/f/7fN5zUqyUXDO1Q+RMkYmieTUREpAIlpWXcNm0RifEx/PmcflpUHGU0gyMiIlKBf89ew+LN2Tx20SBaNknwuhypJs3giIiIlLN0czaPfLCK8QPaMvboNl6XIzWggCMiIhKioLiUW6YtIiUpnvvG9fW6HKkh7aISEREJ8bdZK1m9cz+TrxxGcmKc1+VIDWkGR0REJGjumiye+Wwdl4zoxPE9dObqaFbtGZzgibPSgARg94GL1YmI1Cb1NRJu+wqK+eUri+nUIpE7z+jldTlyhKoUcMzsKOBS4BSgH4GTZx14LAv4HJgOvOqcy6+DOkWkHlBfI166/60MtmXn88r1o0iM1wqOaHfYXVRmdoyZfQQsAUYDHwPXAOOA04ALgAeBAuAhYKuZ/d7Mkuq0ahHxFfU14rX3M3YwbcFmbjihK4M7Nfe6HKkFlUXU14BHgEucc5sP19DMYgiMum4Obrr/iKsTkfpCfY14Jmt/IXe8toTebZpy08k9vC5HakllAaeTc66gKi/knCslcH2XWWamMyKJSHWorxFPOOf47evLyMkv4YWr+xMfq2Nv/OKw/5IHOhwzSzCzp81sRFVetKodlYgIqK8R77yxaAvvLN/OrT/qQa/WTb0uR2pRlaJqsBO5kMDRDCIidUJ9jYTT1r353PXf5Qzp1JxrjuvidTlSy6ozF/chcGJdFSIiEqS+RupcWZnj9umLKS1z/P38/sQ00IU0/aY6x8E9BvzHzBoDM4EdgAtt4JzLqMXaRKR+Ul8jde6JT9bw2eosHjirH51SGntdjtSB6gScd4Jfbw3eQjscC96PKf8kEZFqUl8jder9jB38ddZKfty/LROGdfC6HKkj1Qk4mjIWkXBQXyN15tsd+7hp6tf0bZvMg+ccjZl2TflVlQOOc+7juixERATU10jd2ZNbxNWTFpDYMJanLh1Mo3hNBPpZZWcyrtFZQs2sSc3KEZH6SH2N1LXi0jJ+9tJXbM8u4MlLBtMmuZHXJUkdq+woqo1m9gcz61rZC5lZQzM7x8w+4fszjIqIVIX6GqlTf3grg8/XZPHA2f0Y1FGXYqgPKttFdSqB06D/xswWE7jQ3TJgF1AINAM6A4MJXD8mH/gb8Ggd1Ssi/qS+RurMS/M2MmnuBq45rjPnDm7vdTkSJocNOM65hcAZZtYduAw4CbgSaBjSbCPwWXD7DOdccR3VKiI+pb5G6sq8tVnc9d9ljO6Rxh2n9/a6HAmjKi0yds6tAn534L6ZNSdwptHdzrnCOqpNROoZ9TVSmzbvyeOGF7+iY4tEHpkwUCfzq2eqfBSVmR1/iO0AZUAOsFKdkIgcCfU1UhtyC0u4etICikvLePqyISQ3ivO6JAmz6pwHZzYVn3Ar9H6+mf0HuDV4xV8RkeqajfoaOQJlZY5bpy3i2x37ePaKYXRNq9FBehLlqhNwTgGeIXDq9BlAJpAGjAfOAG4DegO/AfYDv63VSkWkvlBfI0fk4Q9WMWv5Dn43tjeje6R5XY54pDoB50ZgknPunnLbZ5nZPcDlzrkfm1kscDnqdESkZtTXSI3NXLqNhz9YxbmD23PVsZ29Lkc8VJ2rif8ImHOIxz7j+9OrfwK0OZKiRKReU18jNbJ8aza3TVvMoI7N+ONZfXUZhnquOgFnNzDuEI+NCz4OkAhkH0lRIlKvqa+Ratu1v5BrJy+kWWIcT1wymIaxugxDfVedXVQPAo+YWTrwJj/cL/7zYLsTgfm1WKOI1C/qa6RaCktKueGFhezaX8j060fRskmC1yVJBKjOxTYfNbMtwJ3Av4EYoBT4GjjHOfd6sOmfgKLaLlRE6gf1NVIdgXDzFfPX7+GRCQPp1z7Z65IkQlRnBodgx/K6mcUAqcCu8odoOud21WJ9IlIPqa+RqigqKePGl77mw2928sez+jKuf1uvS5IIUq2Ac0Cwo9lRy7WIiBxEfY0cSnFpGb+Y8jXvZezgvvFH8X/DO3ldkkSY6iwyFhER8VxJaRk3T13EO8u3c9eZfbh0ZLrXJUkEUsAREZGoUVJaxi3TFvP20m38bmxvrtS5buQQFHBERCQqlJY5fvnKYt5cvJU7Tu/F1cd18bokiWAKOCIiEvFKyxy3T1/MG4u2cvtpPbl+dFevS5IIp4AjIiIRrazMcedrS3jtqy3cemoPfnZiN69LkiiggCMiIhGrrMzx2zeWMm3BZn5xcnd+cXJ3r0uSKKGAIyIiEck5x10zljHly0387MSu3HKKwo1UnQKOiIhEHOcc98xYzgtfbOS60V345Y966uKZUi0KOCIiElGcc9z3VgaT5m7gmuM6c8eYXgo3Um0KOCIiEjGcc/zpf9/w7GfrueKYdH5zRm+FG6kRBRwREYkYf3/3W576ZC2XjuzEXWf2UbiRGlPAERGRiPCvD1bx6EermTCsA/f8+CiFGzkiCjgiIuK5pz5Zw9/f+5azB7bjjz/pR4MGCjdyZBRwRETEU5M+X88DM7/hzKPb8OC5RyvcSK1QwBEREc9M+XIjd89Yzo/6tOKhCwYQG6P/lqR26DdJREQ88erCzfzm9aWc0DONf100kDiFG6lF+m0SEZGwe3PxVm6fvphRXVN44uLBNIyN8bok8RkFHBERCatZy7dz88uLGNKpBU9fOoSEOIUbqX0KOCIiEjYffbOTG1/6in7tkpl4xVAS42O9Lkl8KuwBx8zGmNlKM1ttZnccpt05ZubMbEg46xMRkbrx2epdXPfCQnq2bsKkK4eR1FDhRupOWAOOmcUAjwGnA32ACWbWp4J2TYCbgHnhrE9E/EODqcgyb20WV02aT5fUxjx/5XCSG8V5XZL4XLhncIYBq51za51zRcBUYHwF7e4H/gIUhLM4EfEHDaYiy1cb93Dlc/Np16wRz181nOaN470uSeqBcAecdsCmkPubg9u+Y2aDgA7OubcP90Jmdq2ZLTCzBZmZmbVfqYhEMw2mIsSyLdlcNvFLUps05KVrRpDWpKHXJUk9EVGLjM2sAfAP4LbK2jrnnnLODXHODUlLS6v74kQkmmgwFQGWbcnm4mfm0TQhjpeuGUGrpglelyT1SLgDzhagQ8j99sFtBzQB+gKzzWw9MAKYoX3jIlKbNJiqe19v3MOEp7+gcXwsU64ZQbtmjbwuSeqZcAec+UB3M+tsZvHAhcCMAw8657Kdc6nOuXTnXDrwBTDOObcgzHWKSHTTYMpD89fv5pJnvqR5YjwvXzeCjimJXpck9VBYA45zrgS4EZgFrACmOeeWm9l9ZjYunLWIiK9pMOWRuWuyuGzil7Rs0pBp142kfXOFG/FG2E9C4JybCcwst+2uQ7Q9IRw1iYi/OOdKzOzAYCoGmHhgMAUscM7NOPwrSE3MWbWLqyfPp0PzRF68Zjgtm2jNjXhHZ1kSEV/SYCq8PvpmJ9e9sJAuqY154erhpCbpaCnxlgKOiIgckXeXb+dnL31Fz9ZNeP5KnedGIoMCjoiI1NjMpdv4xZSvOapdMpOvHKYzFEvEUMAREZEa+e+iLdw6bTEDOjTjuSuG0iRB4UYiR0Sd6E9ERKLD9IWbueXlRQzp1JzJVw5TuJGIoxkcERGplqlfbuTO15dyTNdUnr50CI3iY7wuSeQHNIMjIiJVNnnueu54bSmje6Txn8sUbiRyaQZHREQqlZ1XzH/mrOVfH67m1D6tePSigTSMVbiRyKWAIyIih7RsSzbPz93AfxdvoaC4jPED2vK38/oTF6MdABLZFHBEROQghSWlzFy6jclzN/D1xr00iovhrIHtuWREJ/q0bep1eSJVooAjIiIAbN6Tx4vzNvLy/E3szi2iS2pj7jqzD+cMbq/z20jUUcAREanHysocc1bvYvLcDXz4zQ4ATundiktHpjOqawoNGpjHFYrUjAKOiEg9lFNQzLT5m3hx3kbW7colNSmen57QjQnDO9KuWSOvyxM5Ygo4IiL1zNcb9/CzF79ia3YBgzs15+ZTujOmb2sdFSW+ooAjIlJPOOd4/osN3P9WBq2aJvDqDaMY3Km512WJ1AkFHBGReiC3sIQ7X1vKjMVbOblXS/5x/gCSE7VwWPxLAUdExOdW79zH9S98xdrM/dx+Wk9uGN1Vi4fF9xRwRER87M3FW/n1q0toFBfD81cN55huqV6XJBIWCjgiIj5UVFLGAzNX8Nzn6xncqTmPXTSI1skJXpclEjYKOCIiPrMtO5+fvfgVX23cy5XHdObOM3rp0gpS7yjgiIj4yJxVu/jF1K8pLC7lsYsGMfboNl6XJOIJBRwRER8oK3M89tFq/vH+t3RvmcTjFw+ma1qS12WJeEYBR0Qkyu3OLeK2aYv4aGUmPxnQlgfO7kdivLp3qd/0FyAiEqWy9hfynznrmPz5eopKy7j/J325eHhHzHQIuIgCjohIlMncV8jTn67l+bkbKCgp5cdHt+XnJ3Wje6smXpcmEjEUcEREosTOnAKe/GQtL87bQFFJGeMHtONnJ3ajW0uttREpTwFHRCTCbcvO58mP1/LSlxspLXOcNTAQbDqnNva6NJGIpYAjIhKhtuzN5/HZq5k2fzNlznHOoPb89MSudEpRsBGpjAKOiEgdKSwp5Z1l2ykudTSKiyEhrgGN4mJoGBfz/f34GBJiY2gUH0PD2AaYGZt25/Hv2auZvnAzAOcO7sBPT+hKhxaJHv9EItFDAUdEpI78/o1lTFuwuVrPSYhrQFFJGbENGnDh0I5cf0JX2jVrVEcViviXAo6ISB14ZcEmpi3YzPWju3LRsI7kF5dSUFz63dfvvy8jv6iUgpJSCopKKSgpIyEuhgnDOtAmWcFGpKYUcEREatnK7fv4/X+XMbJLCref1pOYBjovjUi46eprIiK1KLewhBteXEiThDgenjBA4UbEIwo4IiK1xDnHna8tZf2uXB65cCAtmyR4XZJIvaWAIyJSS16ct5EZi7dy2496MrJritfliNRrCjgiIrVg2ZZs7nszg9E90rhhdFevyxGp9xRwRESOUHZ+MT998StSkuJ56IIBNNC6GxHP6SgqEZEj4JzjV9MXs3VvPi9fN4IWjeO9LklE0AyOiMgRmfjZemYt38Edp/dicKcWXpcjIkEKOCIiNfTVxj38aeYKTu3TiquO7ex1OSISQgFHRKQG9uQWceOLX9GmWQJ/O68/Zlp3IxJJtAZHRKSaysoct0xbxK79Rbx6wyiSG8V5XZKIlKMZHBGRanr84zXMXpnJ73/ch37tk70uR0QqoIAjIlINX6zN4u/vruTH/dty8fCOXpcjIoeggCMiUkUlZY6fT/ma9JTG/Onsflp3IxLBFHBERKpo0+48cvKL+ffFg0hqqCWMIpFMf6EiIlWUW1jC4z/pS6/WTb0uRUQqoRkcEZEq6tGqCecP6eB1GSJSBQo4IiJVFB+rLlMkWuivVURERHxHAUdERER8RwFHREREfEcBR0RERHwn7AHHzMaY2UozW21md1Tw+K1mlmFmS8zsAzPrFO4aRUREJLqFNeCYWQzwGHA60AeYYGZ9yjX7GhjinDsamA48GM4aRST6aSAlIuGewRkGrHbOrXXOFQFTgfGhDZxzHznn8oJ3vwDah7lGEYliGkiJCIQ/4LQDNoXc3xzcdihXAf+r6AEzu9bMFpjZgszMzFosUUSinAZSIhK5i4zN7GJgCPDXih53zj3lnBvinBuSlpYW3uJEJJLV2kAKNJgSiVbhDjhbgNDznLcPbjuImZ0C/BYY55wrDFNtIlLPVDaQAg2mRKJVuAPOfKC7mXU2s3jgQmBGaAMzGwg8SSDc7AxzfSIS/TSQEpHwBhznXAlwIzALWAFMc84tN7P7zGxcsNlfgSTgFTNbZGYzDvFyIiIV0UBKRIgN9xs652YCM8ttuyvk+1PCXZOI+IdzrsTMDgykYoCJBwZSwALn3AwOHkgBbHTOjTvki4pI1Al7wBERqWsaSIlIxB5FJSIiIlJTCjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7CjgiIiLiOwo4IiIi4jsKOCIiIuI7YQ84ZjbGzFaa2Wozu6OCxxua2cvBx+eZWXq4axQREZHoFtaAY2YxwGPA6UAfYIKZ9SnX7Cpgj3OuG/AQ8Jdw1igi/qDBlEj9Fu4ZnGHAaufcWudcETAVGF+uzXhgUvD76cDJZmZhrFFEopwGUyIS7oDTDtgUcn9zcFuFbZxzJUA2kFL+hczsWjNbYGYLMjMz66hcEYlSGkyJ1HNRu8jYOfeUc26Ic25IWlqa1+WISGSptcGUiESncAecLUCHkPvtg9sqbGNmsUAykBWW6kREytFssUh0CnfAmQ90N7POZhYPXAjMKNdmBnBZ8PtzgQ+dcy6MNYpI9Ku1wZRmi0WiU1gDTnAa+EZgFrACmOacW25m95nZuGCzZ4AUM1sN3Ar84OgHEZFKaDAlUs/FhvsNnXMzgZnltt0V8n0BcF646xIR/3DOlZjZgcFUDDDxwGAKWOCcm0FgMPV8cDC1m0AIEhGfCHvAEREJBw2mROq3qD2KSkRERORQFHBERETEdxRwRERExHcUcERERMR3FHBERETEdxRwRERExHcUcERERMR3FHBERETEdxRwRERExHfMD5deMbNMYIPXdVRBKrDL6yKqSTWHRzTW3NM518TrIsIpSvqaaPxdUs3hEY01Qw37Gl9cqsE5FxWX+DWzBc65IV7XUR2qOTyitWavawi3aOhrovV3STXXvWisGWre12gXlYiIiPiOAo6IiIj4jgJOeD3ldQE1oJrDQzVLbYnGfxfVHB7RWDPUsG5fLDIWERERCaUZHBEREfEdBRwRERHxHQWcWmZmY8xspZmtNrM7Knj8cjPLNLNFwdvVXtRZrqaJZrbTzJYd4nEzs0eCP9MSMxsU7horqKmymk8ws+yQz/mucNdYQU0dzOwjM8sws+VmdlMFbSLqs65izRH3WdcH0dbXRGM/A9HX16ifCeGc062WbkAMsAboAsQDi4E+5dpcDjzqda3lajoeGAQsO8TjZwD/AwwYAcyLgppPAN7yus5yNbUBBgW/bwJ8W8HvR0R91lWsOeI+a7/forGvicZ+pop1R9Tvv/qZ72+awaldw4DVzrm1zrkiYCow3uOaKuWc+wTYfZgm44HJLuALoJmZtQlPdRWrQs0Rxzm3zTn3VfD7fcAKoF25ZhH1WVexZgm/qOtrorGfgejra9TPfE8Bp3a1AzaF3N9Mxf9I5wSnBaebWYfwlHZEqvpzRZqRZrbYzP5nZkd5XUwoM0sHBgLzyj0UsZ/1YWqGCP6sfcqPfU3E/u5XQUT+/tf3fkYBJ/zeBNKdc0cD7wGTPK7Hr74COjnn+gP/At7wtpzvmVkS8Cpws3Mux+t6qqKSmiP2s67n1NeER0T+/qufUcCpbVuA0FFS++C27zjnspxzhcG7/wEGh6m2I1HpzxVpnHM5zrn9we9nAnFmlupxWZhZHIE/4Bedc69V0CTiPuvKao7Uz9rn/NjXRNzvflVE4u+/+pkABZzaNR/obmadzSweuBCYEdqg3H7OcQT2NUa6GcClwZX3I4Bs59w2r4s6HDNrbWYW/H4Ygd/1LI9rMuAZYIVz7h+HaBZRn3VVao7Ez7oe8GNfE1G/+1UVab//6me+54uriUcK51yJmd0IzCJwlMNE59xyM7sPWOCcmwH8wszGASUEFq5d7lnBQWY2hcAK9VQz2wzcDcQBOOeeAGYSWHW/GsgDrvCm0u9VoeZzgRvMrATIBy50waX4HjoGuARYamaLgtt+A3SEiP2sq1JzJH7WvhaNfU009jMQlX2N+pkgXapBREREfEe7qERERMR3FHBERETEdxRwRERExHcUcERERMR3FHBERETEdxRwfMjM+pqZM7MTQra54GGlVX2NX4U+X0SkPPU1EskUcOqPkcAr1Wj/KwLnfhARqQ71NRIRdKK/eiJ4xVgRkTqlvkYihWZwfMDMfmpmm8ws18zeBH5w2fvy08ZmdqyZfWpmOcHbIjM7L/jYeiAFuDv4vO+moM3sNjObb2bZZrbDzN40s27l3mu2Ba5efJGZrQ6+/v/MrH25do3M7EEz22BmhWa2zsz+VK7N1Wa2PPj4BjP7VQ0+n/Tgz3C+mT0ZrH2zmd1rZg1C2t1jZrvMbLiZLTCzfDObEzwdfksze8PM9pvZCjM7qbp1iEQ79TWVfj7qayKIAk6UM7PxwGPAW8DZwFJgYiXPaRpsvxY4h8ApsJ8HmgWbnAVkE7g2yMjg7avgY+2BR4HxwDUEThP/uZkll3ub4cCNwG3AtcAg4KmQGgz4L3BDsP4zCJwCPTWkze3A4wSuGntm8Pv7rRr798t5ENgf/HlfAO4Kfh8qMVjnQ8AEAqcKfx6YAswh8BlvAV4xs8Qa1iESddTXVIv6mkjgnNMtim/Al8D/ym17GnDACSHbHHBj8PshwftNDvO6u4B7KnnvGKARsA+4NGT7bAKdVvOQbTcH37NR8P5pwfvjDvHaTQl0EHeX234fsB2IqcZnlB58r8nlti8CpobcvyfYbnTItp8Gt90Vsq1PcNvpXv/766ZbuG7qa6r0GamviaCbZnCimJnFEhit/LfcQz+41Hw5awj8Qb9kZuPNrFk13nOEmb1nZlkELuKXByQBPco1ne+c2xNyPyP4tV3w60nAbhe4KGBFRgKNCYxeYg/cgA+BVgRGd9X1brn7GRW8ThHwacj91cGvH1awrR0i9YD6mmpTXxMBFHCiWyqBkc3OctvL3z9IsDM4lcAVcacBmWb2tpl1OdzzzKwjgT9cA64jcAXYocH3SyjXfG+5+0XBrwfapQDbDvN2B6aPlwPFIbePgts7HK7WQ6iopvJ173POlZVrc9BznXPlfxYRv1NfUz0V1aS+Jsx0FFV02wWUAi3LbS9//wdc4EiHMWbWCDgF+AfwEjDiME8bQ2C/8XjnXC58N7JrUf3SyaKCBYohdge/ngnsqODxlTV4TxGpGfU1EnU0gxPFnHMlwNcEFuGFOrsar5HvnHuTwGLBPiEPVTTiaASUEZguPuB8ahaUPwBamNmZh3h8LpAPtHXOLajgtq8G7ykiNaC+RqKRZnCi3wPAa2b2OPA6MJrA6OeQzGwscCWBIwY2Eti/ex0H7/v9BhhrZu8Q2Ie+Mvh4DPCsmT0DHAX8kh9Ox1bFe8AsAvvm7yNw5EQb4Hjn3HXOub1mdg/wsJl1Aj4hEMh7ACc6584K/iwnEJhKPtE5N7sGdYhI1aivUV8TVTSDE+Wcc68DPwd+TKATGQhcVcnTVhNYmf8Agf3cDwLvEOiIDrgdyAXeBuYDg51zS4HLCRyW+RZwEXAegaMYqlu3I3CI6FMEjnr4H/AHAlPhB9o8SOCwz9MJLG6cAvwfBy/MO3D45GHXAojIkVFfo74m2ljg314kOpnZvQRGYid6XYuI+Jf6muijGRyJdqMILFoUEalL6muijGZwRERExHc0gyMiIiK+o4AjIiIivqOAIyIiIr6jgCMiIiK+o4AjIiIivqOAIyIiIr7z/6bJKA8rGNbVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os\n",
    "import h5py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from MDAnalysis.lib.distances import distance_array\n",
    "\n",
    "\n",
    "def compute_radial_distribution_function(simulation_file, \n",
    "                                         selection_1, \n",
    "                                         selection_2, \n",
    "                                         frames=None, \n",
    "                                         bins=30):\n",
    "    \"\"\"Compute the radial distribution function between selections of particles\n",
    "    \n",
    "    Borrows some fast distance calculation C code from MDAnalysis through \n",
    "    `MDAnalysis.lib.distances.distance_array`. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    simulation_file : str\n",
    "        File path of .H5MD simulation file.\n",
    "    selection_1 : str\n",
    "        Name of first selection type.\n",
    "    selection_2 : str\n",
    "        Name of second selection type.\n",
    "    frames : list, optional\n",
    "        List of frame numbers to calculate the RDF for.\n",
    "    bins : int, optional\n",
    "        Number of bins to use in the histogram binning of RDF data.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    edges : numpy.ndarray\n",
    "        Edges of the bins used to bin the RDF data.\n",
    "    histogram : numpy.ndarray\n",
    "        RDF values at the distances of the corresponding `edges`. \n",
    "    \"\"\"\n",
    "    with h5py.File(simulation_file, \"r\") as in_file:\n",
    "        box_size = in_file[\"/particles/all/box/edges\"][...]\n",
    "        positions = in_file[\"/particles/all/position/value\"][...]\n",
    "        types = in_file[\"/particles/all/species\"][...]\n",
    "        indexOfSpecies = in_file[\"/parameters/vmd_structure/indexOfSpecies\"][...]\n",
    "        names = in_file[\"/parameters/vmd_structure/name\"][...]\n",
    "    \n",
    "    volume = np.prod(box_size)\n",
    "    max_distance = np.sqrt(np.sum((0.5 * box_size)**2))\n",
    "    box = np.array([box_size[0], box_size[1], box_size[2], 90.0, 90.0, 90.0])\n",
    "    \n",
    "    index_1 = np.where(names == np.string_(selection_1))\n",
    "    index_2 = np.where(names == np.string_(selection_2))\n",
    "    \n",
    "    ind_1 = np.where(types == index_1[0])\n",
    "    ind_2 = np.where(types == index_2[0])\n",
    "    n_selection_1 = len(ind_1[0])\n",
    "    n_selection_2 = len(ind_2[0])\n",
    "        \n",
    "    if frames is None:\n",
    "        iter_range = list(range(positions.shape[0]))\n",
    "    else:\n",
    "        iter_range = frames\n",
    "    n_frames = len(iter_range)\n",
    "    \n",
    "    distances = np.zeros(shape=(n_selection_1, n_selection_2), dtype=np.float64)\n",
    "    histogram = np.zeros(shape=(bins,), dtype=np.float64)\n",
    "\n",
    "    for i in iter_range:\n",
    "        distance_array(\n",
    "            np.squeeze(positions[i, ind_1, :]), \n",
    "            np.squeeze(positions[i, ind_2, :]), \n",
    "            box=box, \n",
    "            result=distances,\n",
    "        )\n",
    "        count, edges = np.histogram(\n",
    "            distances, bins=bins, range=(0.05, max_distance)\n",
    "        )\n",
    "        histogram += count\n",
    "    \n",
    "    sphere_volumes = np.power(edges, 3)\n",
    "    fractional_volumes = 4/3 * np.pi * np.diff(sphere_volumes)\n",
    "    density = n_selection_1 * n_selection_2 / volume\n",
    "    return edges, histogram / (density * fractional_volumes * n_frames)\n",
    "    \n",
    "\n",
    "fig, axs = plt.subplots(1, 2)\n",
    "fig.set_figwidth(8)\n",
    "fig.set_figheight(6)\n",
    "font_opts = {\"fontsize\": 15}\n",
    "\n",
    "file_path_1 = \"sim_chi=5.H5\"\n",
    "file_path_2 = \"sim_chi=40.H5\"\n",
    "for ax, file_path in zip((axs[0], axs[1]), (file_path_1, file_path_2)):\n",
    "    edges, rdf = compute_radial_distribution_function(\n",
    "        file_path, \"A\", \"B\", range(2, 5), bins=30,\n",
    "    )\n",
    "    ax.plot(edges[:-1], rdf)\n",
    "    ax.set_xlim((0.04, 2.5))\n",
    "    ax.set_xlabel(\"distance, nm\", **font_opts)\n",
    "    ax.set_ylabel(\"g(r)\", **font_opts)\n",
    "    ax.set_ylim((None, 1.1))\n",
    "\n",
    "axs[0].set_title(\"χ = 5 kJ/mol\", **font_opts)\n",
    "axs[1].set_title(\"χ = 40 kJ/mol\", **font_opts)\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}