{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating the Lennard-Jones fluid"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example will simulate a simple liquid that interacts *via* the Lennard-Jones pair potential:\n",
    "\n",
    "$$\n",
    "    u(r) = 4 \\varepsilon \\left[ \\left(\\frac{\\sigma}{r} \\right)^{12} - \\left(\\frac{\\sigma}{r}\\right)^6 \\right]\n",
    "$$\n",
    "\n",
    "The Lennard-Jones potential has two parameters: the energy $\\varepsilon$ and the size $\\sigma$.\n",
    "\n",
    "First, we import `relentless`. We don't need any simulation packages: `relentless` will take care of this for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import relentless"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model\n",
    "\n",
    "Next, we setup the model. We will assign the type \"1\" to the particles. We create an instance `lj` of the Lennard-Jones potential. Then, we specify the parameters for interactions between particles of type 1 by updating the pairwise coefficient matrix `lj.coeff`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "lj = relentless.model.potential.LennardJones(types=(\"1\",))\n",
    "lj.coeff[\"1\", \"1\"].update({\n",
    "    \"epsilon\": 1.0, \"sigma\": 1.0, \"rmax\": 3.0, \"shift\": True\n",
    "})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have chosen $\\varepsilon = 1$ and $\\sigma = 1$. We have also used `rmax` and `shift` options to truncate the potential to zero at $3\\sigma$, which is a common choice for the Lennard-Jones potential.\n",
    "\n",
    "To finish specifying the model, we should also define the thermodynamic state, e.g., the number of particles, the volume, and the temperature. We can conveniently collect this data in an `Ensemble`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "state = relentless.model.Ensemble(\n",
    "    T=1.5, V=relentless.model.extent.Cube(L=10.0), N={\"1\": 400}\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have chosen to have $N=400$ particles in a cubic box with edge length $L=10$ at temperature $T=1.5$. The units of $L$ are implicit to the model (i.e., they are *same* as $\\sigma$). The units of $T$ are also partially implicit to the model (i.e., they are *related* to the units of $\\varepsilon$) but can be manipulated through the Boltzmann constant $k_{\\rm B}$ that converts between energy and temperature units. By default, $k_{\\rm B} = 1$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation operations\n",
    "\n",
    "We now need to specify our simulation protocol. `relentless` runs a simulation as a sequence of operations.\n",
    "\n",
    "Every simulation protocol **must** have an initialization operation to set up the system. Here, we will initialize our system in a semi-random configuration. The `diameters` option will space out the particles to prevent overlap, so we set it equal to $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "init = relentless.simulate.InitializeRandomly(\n",
    "    seed=42,\n",
    "    T=state.T,\n",
    "    V=state.V,\n",
    "    N=state.N,\n",
    "    diameters={\"1\": lj.coeff[\"1\", \"1\"][\"sigma\"]}\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After initialization, we can run a sequence of simulation steps. Our first simulation step is a short equilibration using Langevin dynamics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq = relentless.simulate.RunLangevinDynamics(\n",
    "    steps=10000, timestep=0.005, T=state.T, friction=0.1, seed=2\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we will run a production simulation using Langevin dynamics. We attach an analysis operation to this step so that we can also write a trajectory `trajectory.gsd` of the particles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dump = relentless.simulate.WriteTrajectory(\n",
    "    filename=\"trajectory.gsd\", every=2000\n",
    ")\n",
    "\n",
    "prod = relentless.simulate.RunLangevinDynamics(\n",
    "    steps=50000, timestep=0.005, T=state.T, friction=0.1, seed=7, analyzers=dump\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running a simulation\n",
    "\n",
    "Now that we have defined our operations, we can combine them to be run by a specific simulation engine. Here we will use HOOMD-blue to perform the simulation. In order to run this cell, you need to make sure you have the necessary dependencies installed, or you will get an error. Always make sure to check relentless's documentation for the simulation engine you want to use for these!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = relentless.simulate.HOOMD(initializer=init, operations=[eq, prod])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then need to prepare the model's potentials for the run using `relentless.simulate.Potentials`. This object will turn one or more potentials, which are analytical functions in `relentless`, into a single *tabulated* potential that can be simulated. This means that you need to specify a few additional parameters for the tabulated potential (its starting point, stopping point, and number of subdivisions), along with the buffer that is used to construct pair neighbor lists in many simulation packages.\n",
    "\n",
    "Notes: The starting point should be greater than zero because the Lennard-Jones potential diverges at $r=0$. The stopping point should be sized based on where you would typically truncate interactions. Both the starting and stopping point must be fixed numbers so make sure that you choose them to be large enough or also remember to update them if you vary the potential parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pot = relentless.simulate.Potentials()\n",
    "pot.pair = relentless.simulate.PairPotentialTabulator(\n",
    "    potentials=lj, start=1e-6, stop=3.0, num=100, neighbor_buffer=0.5\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last, we can run the simulation. An output directory is required for all runs, so here we will use the working directory. This will be the slowest step of this example, but don't worry, it took less than a minute to run on an old laptop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = sim.run(potentials=pot, directory=\".\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `result` contains data from the run. Some of it can be useful for subsequent analysis, but some of it is just used internally by `relentless`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "To demonstrate that the simulation ran, we perform additional analysis of the trajectory file that was created. We will use the `gsd` package to read the file and the `freud` package to compute the radial distribution function $g(r)$. You need to make sure you have both installed, along with `matplotlib` to render the plot.\n",
    "\n",
    "Note that we use `result.directory` (which is a `relentless.data.Directory`) to get the path to `trajectory.gsd` created by the simulation!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHHCAYAAABeLEexAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOoUlEQVR4nO3deVzUdf4H8Nd3BmaGcwC5FcUTAxUUlTBNXVE0s6zdPDo0y67t/FHtSlvatUvXmh1utpVpa6W2ll2uWhSaiZko3pooiijDocDAcM98f3/gjEyCcszM9zszr+fj8X0kXz7znfdMOLz8fD+HIIqiCCIiIiI3ppC6ACIiIiKpMRARERGR22MgIiIiIrfHQERERERuj4GIiIiI3B4DEREREbk9BiIiIiJyewxERERE5PYYiIiIiMjtMRARERGR22MgIiKXsGLFCgiCYDk8PDzQvXt33HnnnThz5oxV23HjxlnaKRQK+Pv7IyYmBnfccQe+++67Vq8fHR1tdf2WR11dnSNeIhHZkYfUBRAR2dLzzz+P3r17o66uDjt27MCKFSuwbds2HDhwABqNxtKuR48eyMjIAAAYDAbk5eXh888/x6pVqzBjxgysWrUKnp6eVtdOSEjA448/fslzqlQq+74oIrI7BiIicilTpkzB8OHDAQDz589HcHAwXn75ZXz11VeYMWOGpZ1Wq8Xtt99u9diXXnoJjzzyCP71r38hOjoaL7/8stX3u3fvfsljiMg18JYZEbm0MWPGAACOHz9+xbZKpRJvvvkmYmNj8fbbb6OystLe5RGRTDAQEZFLO3nyJAAgMDCwXe2VSiVmz56NmpoabNu2zep7jY2NKCsrszpqampsXTIRSYCBiIhcSmVlJcrKylBYWIh169bhueeeg1qtxvXXX9/uawwaNAjApb1KmzdvRkhIiNXxyiuv2LR+IpIGxxARkUtJSUmx+jo6OhqrVq1Cjx492n0NX19fAEBVVZXV+aSkJLz44otW5/r06dPJSolIThiIiMilLF26FAMGDEBlZSWWL1+OrVu3Qq1Wd+ga1dXVAAA/Pz+r88HBwZcELiJyDQxERORSRo4caZllNn36dIwePRq33norjh49aun5uZIDBw4AAPr162e3OolIXjiGiIhcllKpREZGBs6ePYu33367XY8xGo345JNP4O3tjdGjR9u5QiKSCwYiInJp48aNw8iRI7FkyZIrrihtNBrxyCOP4PDhw3jkkUfg7+/voCqJSGq8ZUZELu/JJ5/ELbfcghUrVuD+++8H0DwbbdWqVQCAmpoay0rVx48fx6xZs/DCCy9IWTIRORgDERG5vJtvvhl9+/bFa6+9hnvuuQcAUFhYiDvuuANA86yyiIgIJCcn45133sHEiROlLJeIJCCIoihKXQQRERGRlDiGiIiIiNweAxERERG5PQYiIiIicnsMREREROT2GIiIiIjI7TEQERERkdvjOkTtYDKZcPbsWfj5+UEQBKnLISIionYQRRFVVVWIjIyEQnH5PiAGonY4e/YsoqKipC6DiIiIOuH06dPo0aPHZdswELWDn58fgOY3lHsbEREROQe9Xo+oqCjL7/HLYSBqB/NtMn9/fwYiIiIiJ9Oe4S4cVE1ERERuj4GIiIiI3B4DEREREbk9BiIiIiJyewxERERE5PYYiIiIiMjtMRARERGR22MgIiIiIrfHQERERERuj4GIiIiI3B4DEREREbk9BiIiIiJyewxERNSm2gYjRFGUugwiIrtjICKiVv1wpBgJz2/Gs18dlLoUIiK7YyAiokvklxnw6Opc1DeZ8PmeM2gymqQuiYjIrhiIiMiKob4J9360C1V1TQCAqrom5J6ukLYoIiI7YyAiIgtRFPHkf/fiWEk1Qv3UuKZfNwDAlt9KJa6MiMi+GIiIyOLdrSewYb8OnkoB79w+DDcN7QEA2MpAREQuzkPqAohIHn46VopXNh4BACyaFofEXkGICvQGAOw7U4lz1fXo5quWskQiIrthDxER4fT5Gjz86R6YRGDG8B64LaknACDUX4OrIvwhisC2vDKJqyQish8GIiI3V9tgxH3/yUFFTSPie2jx/I2DIAiC5ftjB4QAALYc5W0zInJdDEREbkwURTz1xX4cKtKjm48K79yeCI2n0qqNORBtPVYGk4mLNBKRa2IgInJjK7efxBd7zkCpEPD2rcMQGeB1SZvEXoHwUSlRVl2PQ0V6CaokIrI/BiIiNyWKIv65+TcAQPqUgUju263VdioPBZL7BgPg9Hsicl0MRERu6ryhAVX1zYsvzkmOvmzbsTEXbpsxEBGRi2IgInJTxfp6AECwrwoqj8t/FIzt3xyIck6Vo6qu0e61ERE5GgMRkZsq1tcBAEL9NFds27ObN3oH+6DJJGL78XP2Lo2IyOEYiIjclDkQhWuvHIiAFtPveduMiFwQAxGRm9JdCERh/u0LRNcOaB5YvfW3Uogip98TkWthICJyU+YxRGH+7duO4+o+3aBSKlBYXosTZQZ7lkZE5HAMRERuynLLrJ09RN4qD4zsHQSAq1YTkeuRVSDaunUrpk2bhsjISAiCgPXr11+2/Z133glBEC454uLiLG2effbZS74/cOBAO78SIvkr7uAtM6DFbbNjDERE5FpkFYgMBgPi4+OxdOnSdrV/4403UFRUZDlOnz6NoKAg3HLLLVbt4uLirNpt27bNHuUTOZXOBKKxA0IBADtOnENdo9EudRERScFD6gJamjJlCqZMmdLu9lqtFlqt1vL1+vXrUV5ejnnz5lm18/DwQHh4uM3qJHJ2DU0mlFU3AGj/GCIAGBDmi3B/DXT6OuzMP49rL8w8IyJydrLqIeqqDz74ACkpKejVq5fV+WPHjiEyMhJ9+vTBbbfdhoKCAokqJJKH0urmAdWeSgFBPqp2P04QBMttM06/JyJX4jKB6OzZs/jf//6H+fPnW51PSkrCihUrsHHjRrzzzjvIz8/HmDFjUFVV1ea16uvrodfrrQ4iV6KrvLgooyAIHXqs+bYZt/EgIlfiMoFo5cqVCAgIwPTp063OT5kyBbfccguGDBmC1NRUbNiwARUVFVi7dm2b18rIyLDcjtNqtYiKirJz9USOVdLBRRlbGt0vGAoBOFZSjTMVtbYujYhIEi4RiERRxPLly3HHHXdApbp8939AQAAGDBiAvLy8Ntukp6ejsrLScpw+fdrWJRNJ6uKijO0fP2Sm9fZEQlQAAPYSEZHrcIlAtGXLFuTl5eHuu+++Ytvq6mocP34cERERbbZRq9Xw9/e3OohcycVFGTveQwTwthkRuR5ZBaLq6mrk5uYiNzcXAJCfn4/c3FzLIOj09HTMmTPnksd98MEHSEpKwqBBgy753hNPPIEtW7bg5MmT2L59O2666SYolUrMnj3brq+FSM46uijj742NaZ5dtu1YGRqNJpvVRUQkFVlNu9+1axfGjx9v+TotLQ0AMHfuXKxYsQJFRUWXzBCrrKzEunXr8MYbb7R6zcLCQsyePRvnzp1DSEgIRo8ejR07diAkhNOFyX11Zg2ilgZ31yLA2xMVNY3YU1BhWcGaiMhZySoQjRs37rKbRq5YseKSc1qtFjU1NW0+ZvXq1bYojcildHRj199TKgRc0y8Y3+4rws78cwxEROT0ZHXLjIgco6SDG7u2JrFnIABgT0GFLUoiIpIUAxGRm6mub0J1fROAzvcQAcDQngEAgD2nKy7bs0tE5AwYiIjcjHn8kJ/aAz7qzt81j430h0qpwHlDAwrOt33bmojIGTAQEbmZ4gurVId1YlHGltQeSsR1b16SYndBeZfrIiKSEgMRkZvpyqKMvzc0iuOIiMg1MBARuZmuLsrYkmUcEQMRETk5BiIiN9PVNYhaGtaruYfocJEetQ3GLl+PiEgqDEREbqarq1S3FKnVINRPjSaTiANnK7t8PSIiqTAQEbmZri7K2JIgCC1um3FgNRE5LwYiIjdji0UZWxrKBRqJyAUwEBG5EZNJvHjLrIvT7s2GRgUAYCAiIufGQETkRs7XNKDJJEIQgGBf2/QQDe6hhVIhQKevw9mKWptck4jI0RiIiNyI7sKijMG+angqbfPX31vlgYHhfgDYS0REzouBiMiNlFTZblHGloZZxhFxYDUROScGIiI3oqtsHlBtiyn3LbXc6JWIyBkxEBG5EfOU+1CbB6LmHqL9ZyrR0GSy6bWJiByBgYjIjZTYcFHGlqK7eSPA2xMNTSYcLtLb9NpERI7AQETkRmy5sWtLgiC0mH7PcURE5HwYiIjciC03dv09ywKNHEdERE6IgYjIjdh6UcaWzAOrd7OHiIicEAMRkZuobzLivKEBABDmZ/tAFB8VAEEATp+vRWlVvc2vT0RkTwxERG7CvIeZykOBAG9Pm1/fX+OJ/qG+AIBc3jYjIifDQETkJlouyigIgl2eY2gUF2gkIufEQETkJuy1KGNLlgUauYUHETkZBiIiN1Fsp0UZWzLPNNtbWAGjSbTb8xAR2RoDEZGbKLbToowt9Qv1ha/aAzUNRvxWXGW35yEisjUGIiI3UWynRRlbUioExEdpAfC2GRE5FwYiIjdxcZVq+/UQARcHVnM9IiJyJgxERG7CnqtUtzSsVwAAzjQjIufCQETkBkRRdMgYIgBIuNBDdLzUgMqaRrs+FxGRrTAQEbmBqvom1DQYAdi/hyjIR4Xobt4AgNzCCrs+FxGRrTAQEbmBkgu9Q/4aD3iplHZ/PstGr7xtRkROgoGIyA1YFmW0w6aureECjUTkbBiIiNxAsYNmmJm13MLDxAUaicgJyCoQbd26FdOmTUNkZCQEQcD69esv2z4rKwuCIFxy6HQ6q3ZLly5FdHQ0NBoNkpKSsHPnTju+CiL5cdSUe7OBEX7wUSmhr2vCoSK9Q56TiKgrZBWIDAYD4uPjsXTp0g497ujRoygqKrIcoaGhlu+tWbMGaWlpWLRoEXbv3o34+HikpqaipKTE1uUTyVaJAxZlbMlTqUBy324AgJ+OlTnkOYmIukJWgWjKlCl48cUXcdNNN3XocaGhoQgPD7ccCsXFl7V48WLcc889mDdvHmJjY7Fs2TJ4e3tj+fLlti6fSLZ0Dppy39LofsEAgG15pQ57TiKizpJVIOqshIQEREREYOLEifj5558t5xsaGpCTk4OUlBTLOYVCgZSUFGRnZ0tRKpEkzIsy2nNj198bMyAEAPBrfjlqL0z5JyKSK6cORBEREVi2bBnWrVuHdevWISoqCuPGjcPu3bsBAGVlZTAajQgLC7N6XFhY2CXjjFqqr6+HXq+3OoicmaMWZWypT7APugd4ocFowi/55xz2vEREneHUgSgmJgb33XcfEhMTMWrUKCxfvhyjRo3C66+/3qXrZmRkQKvVWo6oqCgbVUzkeEaTiJIqx2zb0ZIgCBjT/8JtM44jIiKZc+pA1JqRI0ciLy8PABAcHAylUoni4mKrNsXFxQgPD2/zGunp6aisrLQcp0+ftmvNRPZ0zlAPo0mEQgCCfVUOfe7RFwIRB1YTkdy5XCDKzc1FREQEAEClUiExMRGZmZmW75tMJmRmZiI5ObnNa6jVavj7+1sdRM6q+MKijCF+angoHftX/pq+wRAE4GhxleW2HRGRHHlIXUBL1dXVlt4dAMjPz0dubi6CgoLQs2dPpKen48yZM/joo48AAEuWLEHv3r0RFxeHuro6vP/++/jhhx+wefNmyzXS0tIwd+5cDB8+HCNHjsSSJUtgMBgwb948h78+Iik4elHGlgJ9VBjSXYu9hZX46VgZ/pTYw+E1EBG1h6wC0a5duzB+/HjL12lpaQCAuXPnYsWKFSgqKkJBQYHl+w0NDXj88cdx5swZeHt7Y8iQIfj++++trjFz5kyUlpZi4cKF0Ol0SEhIwMaNGy8ZaE3kqhy9KOPvje4fjL2Fldh2rJSBiIhkSxBFkevqX4Fer4dWq0VlZSVvn5HTWbz5KN78IQ+3X90TL04f7PDn33HiHGb9eweCfVXY+VQKFArB4TUQkXvqyO9vlxtDRETWpFiUsaVhPQPhrVKirLoBh3VcwoKI5ImBiMjFSbEoY0sqDwWu7tO8jQen3xORXDEQEbk4KRZl/L0xnH5PRDLHQETk4qScZWY2pn/zNh47T57nNh5EJEsMREQurK7RiPKaRgDS9hD1DfFBpFaDhiYTdp48L1kdRERtYSAicmGlF7bsUHso4O8l3SobgiBYVq3edqxUsjqIiNrCQETkwiwzzLQaCIK0093Nt804joiI5IiBiMiF6SqlHz9kdk2/5m08juiqUMJtPIhIZhiIiFyYHAZUmwX5qDAoUgsA2JbHXiIikhcGIiIXdnHKvVriSppx+j0RyRUDEZELO1fdAAAI9pVLILo4jshk4q5BRCQfDERELkxf1zzl3t/LU+JKmg3rFXBhG496HNFVSV0OEZEFAxGRC6uqawIA+Gmkm3LfktpDiaTeQQCAbXmcfk9E8sFAROTCLgYiefQQAZx+T0TyxEBE5MKq6ptvmcmlhwgArh3QPLD6l/zzqGvkNh5EJA8MREQuzNxD5C+jQNQ3xBfh/he28cjnNh5EJA8MREQuShRFWd4yEwShxfR7jiMiInlgICJyUbWNRhgvTG2X0y0zALh2QPM4oswjJRJXQkTUjIGIyEWZe4eUCgFenkqJq7E2LiYEnkoBJ0oNyCuplrocIiIGIiJXVXVhDSJftYfkG7v+np/GE6P6Nt8223xIJ3E1REQMREQuSy+zNYh+b1JcGABg08FiiSshImIgInJZchxQ3dLE2DAIArD3dAV0lXVSl0NEbo6BiMhFmW+ZybWHKNRPg6FRAQCA7w6zl4iIpMVAROSiqmW4BtHvTYoLBwBsPshxREQkLQYiIhcl91tmADAptnkcUfbxc6isbZS4GiJyZwxERC5K7rfMAKBPiC/6h/qiySQi6yjXJCIi6TAQEbkouc8yMzPPNtvM2WZEJCEGIiIX5Qy3zABgUmzzOKKsoyXc7JWIJMNAROSinOGWGQAM7q5FuL8GhgYjth8vk7ocInJTDERELsrcQ+SrlncgUigE3jYjIskxEBG5qKr65h4if5nfMgMu3jb7/nCxZUNaIiJHYiAiclFVTjKoGgCS+gTBT+OBsuoG7Ckol7ocInJDDERELspZBlUDgKdSgQkDQwEAm7hIIxFJgIGIyAWJoug0g6rNUs2rVh8qhijythkROZasAtHWrVsxbdo0REZGQhAErF+//rLtP//8c0ycOBEhISHw9/dHcnIyNm3aZNXm2WefhSAIVsfAgQPt+CqIpFffZEKjsTlUOEsgunZACFQeCpw6V4PfiqulLoeI3IysApHBYEB8fDyWLl3arvZbt27FxIkTsWHDBuTk5GD8+PGYNm0a9uzZY9UuLi4ORUVFlmPbtm32KJ9INsy3ywQB8FE5RyDyUXtgTL9gANzbjIgcT1aflFOmTMGUKVPa3X7JkiVWX//jH//Al19+ia+//hpDhw61nPfw8EB4eLityiSSPfPtMl+1BxQKQeJq2m9SXBgyj5Rg0yEdHp7QX+pyiMiNyKqHqKtMJhOqqqoQFBRkdf7YsWOIjIxEnz59cNttt6GgoOCy16mvr4der7c6iJxJlWWne/kPqG4p5aowKATgwBk9zlTUSl0OEbkRlwpEr732GqqrqzFjxgzLuaSkJKxYsQIbN27EO++8g/z8fIwZMwZVVVVtXicjIwNardZyREVFOaJ8Iptxpin3LXXzVWN4r+Z/0HzH22ZE5EAuE4g++eQTPPfcc1i7di1CQ0Mt56dMmYJbbrkFQ4YMQWpqKjZs2ICKigqsXbu2zWulp6ejsrLScpw+fdoRL4HIZlreMnM2llWrD3HVaiJyHJcIRKtXr8b8+fOxdu1apKSkXLZtQEAABgwYgLy8vDbbqNVq+Pv7Wx1EzsRZe4iAi6tW/5J/HuWGBomrISJ34fSB6NNPP8W8efPw6aefYurUqVdsX11djePHjyMiIsIB1RFJQ29Zg8i5xhABQM9u3rgqwh9Gk4gXvj3ENYmIyCFkFYiqq6uRm5uL3NxcAEB+fj5yc3Mtg6DT09MxZ84cS/tPPvkEc+bMwT//+U8kJSVBp9NBp9OhsrLS0uaJJ57Ali1bcPLkSWzfvh033XQTlEolZs+e7dDXRuRIztxDBAB/u+4qKBUCPt99Bv/eekLqcojIDcgqEO3atQtDhw61TJlPS0vD0KFDsXDhQgBAUVGR1Qyxf//732hqasKDDz6IiIgIy/Hoo49a2hQWFmL27NmIiYnBjBkz0K1bN+zYsQMhISGOfXFEDuRM23a0ZnT/YDwz9SoAwEsbj+B7jiciIjsTRPZHX5Fer4dWq0VlZSXHE5FTePKzvfgspxBPpsbgwfH9pC6nU0RRxN/WH8AnvxTAR6XEuj+PwsBw/v0jovbryO9vWfUQEZFtVNeb1yFyzltmACAIAp67IQ7JfbrB0GDE/JW7cK66XuqyiMhFMRARuSBnv2Vm5qlU4F+3DUOvbt4oLK/F/atyUN9klLosInJBDERELsjZdrq/nEAfFT6YOxx+ag/8erIcT39xgDPPiMjmGIiIXJCr9BCZ9Qv1w1u3DoVCAD7LKcQH2/KlLomIXAwDEZEL0jv5tPvWjIsJxd+mxgIA/rHhMH48UiJxRUTkShiIiFyQM2/dcTl3XRONWSOiYBKBJz7biyajSeqSiMhFMBARuZiGJhPqm5qDgrPtdn8lgiDg+RsHIcDbE+cMDdhbWCF1SUTkIhiIiFyMuXcIAHxd6JaZmcpDgTH9mxdWzTpaKnE1ROQqGIiIXIx5QLWPSgmlQpC4GvsYO4CBiIhsi4GIyMW42gyz1lw7IBgAsP9MJcq4WCMR2QADEZGLcaU1iNoS6qdBXGTzMvxbf2MvERF1HQMRkYupqne9KfetGRfTfNtsCwMREdkAAxGRi3GHW2ZA87pEQHMPkdHElauJqGsYiIhcjDvcMgOAoVEB8NN4oLymEfs4/Z6IuoiBiMjFuEsPkYdSgTH9mwdX87YZEXUVAxGRi3GXHiKA0++JyHYYiIhcjKWHyMW27WjN2AHN44j2FlbgvKFB4mqIyJkxEBG5mCoX3Ni1LeFaDQaG+0EUgZ+OsZeIiDqPgYjIxegtt8xcewyR2Vjz9HveNiOiLmAgInIx7tRDBADjLtw223qsFCZOvyeiTmIgInIxVW7WQ5TYKxA+KiXKqhtw8Kxe6nKIyEkxEBG5GHfrIVJ5KHBNv+bp91lHSySuhoicFQMRkYupvrB1h7+b9BABF1et5npERNRZDERELqTJaEJNgxGA+/QQARcHVu8uKEdlTaPE1RCRM2IgInIh5t4hAPB1o0DUPcAL/UN9YRKBbXllUpdDDlLfZIQociA92QYDEZELMY8f8vJUwlPpXn+9L65azXFE7mD9njMY/sL3GPtqFj7fXcgNfqnL3OsTk8jFmdcgcqfeIbOW44jYa+C6ahqa8MRne/HYmlxU1Teh4HwN0tbuReqSrdiwv4hLL1CnMRARuRB3m2HW0ojegfDyVKKkqh6Hi6qkLofs4HCRHte/tQ3/zSmEQgAemdAff508EFovT+SVVOPPH+/GtLe34ccjJQzF1GHu96lJ5MLcZaf71qg9lBjVtxsyj5Qg67cSxEb6S10SXYYoith0UIf/5hSiVzcfpFwVhhHRgfBo5VavKIpY9UsBXvjmEBqaTAjzV2PJzKFI7tsNAHDb1T3x/k/5+OCnEzh4Vo95K35FYq9APD5xAJL6dINSITj65ZETYiAiciHmRRn93bCHCADGxYQg80gJthwtxZ/H9ZO6HGrD2YpaLPzyIL4/XGw598G2fGi9PDE+JgQpsWEYOyAEfhpPVNY0YsHn+/C/AzoAwB8GhuK1W+IR5KOyPNZf44m0iQNw56hoLNtyHCu3n0TOqXLc+v4vUCkV6NnNG72DfdAn2Ae9zUeID0J81RAEhiVq5p6fmkQuyp1vmQHA2AGhAA4i51Q5quoa3bKnTM6MJhErt5/EPzcfhaHBCA+FgLmjolFZ24gfjpTgvKEB63PPYn3uWXgqBVzdpxtOlBpwpqIWnkoBf508EHeP7t1miAnyUeGp667C3aN7Y+mPefhsVyFqG43IK6lGXkn1Je191R7oG+KDPiG+6Bvig74hvugT4ote3byh8VTa++0gmXHPT00iF2XZtkPtnkGgZzdv9An2wYkyA37OK8PkQRFSl0QXHDhTiae+2I99hZUAmrdcybh5MAaE+QFoDkt7Csrx3eFifHeoGCdKDfjpWPMSCj2DvPHW7KGIjwpo13OF+Wvw/I2D8Oy0OJytrEV+mQH5ZQacKDVY/lxYXoPq+ibsLazE3gs1mSkEoEegN7r5quCtUsLL0wPeKmXzny/8N9RPgz8m9oCvmr9GXQX/TxK5EHfvIQKaF2k8UWbAu1tPYFxMKP+lb2fnDQ3IOloCtYcSPmol/DQe8FV7wlfjAV+1BxQC8GbmMSz/+SSMJhF+Gg8smDIQs0f0hKLF2B6lQsDw6CAMjw5C+pSrcLy0GpmHi2E0NY8R6szK6wqFgB6B3ugR6I0x/UOsvlffZETBuRocLzXgeGk1jpdW48SFP1fVNc9eKzhfc9nrf7qzAMvvHIHIAK8O10by476fmkQuqKrefQdVm905Khrrcgqxp6ACC9btw+szEzhOxE6K9XW4aenPOFtZ1672UwdHYNG0WIT6a67Ytm+IL/qG+Ha1xDapPZToH+aH/hd6qMxEUURpdT3ySw2oqG1EbYMRNQ1G1DQ0Nf+50YjaBiO+3V+EI7oq3PSvn/HB3BEY1F1rt1rJMWQ17X7r1q2YNm0aIiMjIQgC1q9ff8XHZGVlYdiwYVCr1ejXrx9WrFhxSZulS5ciOjoaGo0GSUlJ2Llzp+2LJ5IB9hABvbr54J3bE+GhELA+9yze/iFP6pJcUnV9E+Z9+CvOVtYh3F+Dkb2DEBvhj55B3gjyUcFTeTGEdg/wwvI7h2PpbcPaFYakJAgCQv00SOrTDalx4Zg+tDtuTeqJ+WP64OEL0/yfvSEOX/x5FAaE+aJYX48Z72bjhyPFV744yZqsPjUNBgPi4+Nx11134eabb75i+/z8fEydOhX3338/Pv74Y2RmZmL+/PmIiIhAamoqAGDNmjVIS0vDsmXLkJSUhCVLliA1NRVHjx5FaGiovV8SkUNZxhC5cSACgGv6BeOF6YOQ/vl+/PO739A7xAfXD4mUuiyX0Wg04cGPd+NQkR7dfFRYe18yenbzvqRdfZMR1XVNCPRWWd0ecwU9Ar3x3wdG4c+rdmNbXhnmr9yF526Iwx3J0VKXRp0kiDJdvUoQBHzxxReYPn16m23++te/4ttvv8WBAwcs52bNmoWKigps3LgRAJCUlIQRI0bg7bffBgCYTCZERUXh4YcfxoIFC9pVi16vh1arRWVlJfz9ubYJydcf39mOnFPlWHb7MA4oBvDiN4fw/rZ8qD0UWHNfMhLaOSiX2iaKIp76Yj8+3XkaGk8FVt/r3u9ro9GEv32xH2t3FQIA5o/ujaeuu6rNAFhWXY/8MgP6h/oiwFvVahs50lXW4btDOvQJ8cXI3kFOszVQR35/O/U/I7Ozs5GSkmJ1LjU1FY899hgAoKGhATk5OUhPT7d8X6FQICUlBdnZ2W1et76+HvX19Zav9Xq9bQsnspOLPUTuO4aopfTrrkJ+mQGZR0owf+UufPnQNejOAbBd8q+s4/h052kIAvDmrKFuHYYAwFOpwMt/HIKeQd54bfNveH9bPgrLa/HPGfEoqqzDoSI9Dhfpcehs839Lqpp/t2g8FfhTYg/cdU1v9LHBWKlGowk7888j+/g5eKmU6B7ghQitBpEBXgjXajodYOoajXj/pxNY+uNx1DYaAQBaL09MuCoUqXHhuLZ/CLxUrjFxwakDkU6nQ1hYmNW5sLAw6PV61NbWory8HEajsdU2R44cafO6GRkZeO655+xSM5E9cQyRNaVCwBuzh+JP72zHEV0V7l7xK/77wChOle6kL3PP4NVNRwEAi66PxaS4cIkrkgdBEPDQH/ojKsgbT362DxsP6rDpWR1au/8iCEA3HxXKqhuwakcBPv6lABMGhmH+mN5I6h3UoQkAdY1G/HSsDBsP6JB5pBgVNY1t1AeE+qkRofVC72AfTIwNw/iY0MsGGVEU8b8DOvxjw2EUltcCAAaG+6Gkqh7nDQ34fPcZfL77DDSeCowdEILUuHCMiwm1WjCzPQz1TTheWo1jxdW4KsJf0hXmu/Sp0NjYCJ1Oh5qaGoSEhCAoKMhWdUkqPT0daWlplq/1ej2ioqIkrIiofdx56462+Ko98P7c4Zi+tDkUPbZ6D969Yzi3c+igHSfO4cnP9gEA7h7dG3de01viiuTnxoTuiNB64d7/7EJFTSO8PJWICfdDbKR/8y/7CH8MDPeDt0qJ7BPn8MFP+cg8UoLvDxfj+8PFGNTdH/NH98HUIRGWHh1RFNFkElHfZELdhRluuwvKsemgDllHS1HTYLQ8f5CPCuNiQiBAwNmKWhRV1uJsZR0amkwo1tejWF+P3NMV+GLPGXirlEi5KgzXD4nAtQNCrJanOHRWj+e+Pohf8s8DACK0GiyYMhA3xEfCJAK7Tp7HpoPF2HRQhzMVtRf+3DyoXOvliaggL/QM8kZUkDeiAr3RM8gbPQK9oK9rwm/FVcgrqcax4iocK6m2hC2geW86pwpEVVVVWLVqFVavXo2dO3eioaEBoihCEAT06NEDkyZNwr333osRI0bYo14r4eHhKC62HtlfXFwMf39/eHl5QalUQqlUttomPLztf9mo1Wqo1Wq71ExkL0aTiOp69hC1pkegN96bk4iZ/96B7w+X4PmvD+KhP/RHiB//nrdHXkkV7v1oFxqMJkyOC8ffrrtK6pJka2TvIOxInwBdZR2igrzbDN6j+gZjVN9gHC+txvJt+Vi3uxAHzujx2JpcPLP+ADyUgiUEmS4z0jdSq8GkuHBMHhSO4b0u3QtOFEWcMzTgbEUtzlbUYc/pcnyztwhnKmrx1d6z+GrvWfipPTAxNgypg8Kx5bdSrN5ZAJMIqD0UuH9sX9w3tg+8Vc2fKUoBSOrTDUl9uuGZ66/CoSI9Nh0sxuaDOhzRVaGythGVZxpx4Ez7h5oE+6rQL9QXkVppZyB2aFD14sWL8fe//x19+/bFtGnTMHLkSERGRsLLywvnz5/HgQMH8NNPP2H9+vVISkrCW2+9hf79+3eusHYOqt6wYQP2799vOXfrrbfi/PnzVoOqR44cibfeegtA86Dqnj174qGHHuKganIplbWNiH9uMwDg6IuTofZwjfv6tvT13rN4+NM9lq8jtBoM7q7FkB5aDO4RgMHdtR3u8nd1ZypqMfPdbBSW12JozwB8es/VXOzSDs4bGvDJL6ewMvsUSqvq22yn8lCgZ5A3UuPCkBoXjsHdtR1eZ0sUReSersA3+4rw7b4i6PSXriN1/ZAILJgyED0CL5092JaahiYUltei4FyNZWHLwvLm/54+XwtfjQcGhPmif6gf+oX6YkBY83/t+XeuI7+/OxSIZs+ejaeffhoDBw7E119/jQkTJsDPz++SdvX19fjwww+hUqlw1113tbvw6upq5OU1rxkydOhQLF68GOPHj0dQUBB69uyJ9PR0nDlzBh999BGA5mn3gwYNwoMPPoi77roLP/zwAx555BF8++23VtPu586di3fffRcjR47EkiVLsHbtWhw5cuSSsUVtYSAiZ1BYXoPRL/8IlYcCv704RepyZGv1zgJ8sC0feaXVrY7x6BHohZuGdscjE/o7zUyaKxFFEUaT2OpO8peTdbQEj63JRUVNI3p188bnD4xCN1/2qtlTfZMRx0sM8FQKUHsoofZUQHPhvyqlwubLF5hMInIKyvHtviJ8f7gYoX5qLJhyFUb2do0hMHYLRC15eXnh4MGD6NOnT6eKbE1WVhbGjx9/yfm5c+dixYoVuPPOO3Hy5ElkZWVZPeb//u//cOjQIfTo0QPPPPMM7rzzTqvHv/3223j11Veh0+mQkJCAN998E0lJSe2ui4GInMHhIj2mvPETgn1V2PX0RKnLkT1DfRMOntVjX2EF9p+pxP7CSpwoM1i+PzI6CG/fOlT2Cwm2R8aGw3h/Wz5mDI9C2sQBV7xVaDSJeCPzGN764RhEERjcXYt3bh/Wod4CIjlwSCC69tprsWjRIkyYMKFTRToTBiJyBr+ePI9blmWjd7APfnxinNTlOKXK2kZkHi7Gwi8Porq+CSF+aiy9dZhT/2t5f2Elbli6zdIb5qv2wJ/H98Vd1/Ru9dbXuep6PLo6F9vymjdWvS2pJ565Ppa3ycgpdeT3d6f7gx9++GE89dRTOH36dGcvQUQ2xFWqu07r5Ymbh/XAVw9dg5gwP5RW1WP2ezvw/k8nINM1bC9LFEU89/VBiCIwul8whvTQorq+Ca9sPIqUxVvwzb6zVq8r59R5TH1zG7bllcHLU4nXZ8bj7zcNZhgit9DpT86ZM2cCAOLi4nDDDTdg3LhxGDp0KAYPHgyVioMSiRzNPOWea+x0XZ8QX3zx4Cikf74fX+aexYvfHsbugnK88qd4p3p/v9p7FrtOlcPLU4lXbxmCMD8N1ueewSsbj6KwvBYPfbIHH/Y6iaenXoWcU+V46X9H0GQS0SfEB8tuT8SAsEvHiBK5qk7/zc7Pz8fevXuRm5uLvXv3IiMjAydPnoSHhwdiYmKwb98+W9ZJRFeg56KMNuWt8sCSmQkY3isQz39zCBv2N08rfvf2xEt2SJcjQ30TMjY0L0D753F9EaFtXqH75mE9MGVQBP699QSWbTmOnFPluOlf2y2Pu35IBF764xCnCn5EttDpn/hevXqhV69euOGGGyznqqqqkJubyzBEJAFu22F7giDgjuRoxHXX4sGPd+NEqQE3Lv0Zb84aipTY9s1Slcq/svKg09chKsgL91xrPfnFS6XEoyn9MXNEFF7ddBTrdhfCUyng6amxmJPcq8PTuIlcgWw3d5UTDqomZ/DyxiN4J+s45l0TjUXT4qQux+W0HGzsoRDw5uyhuG6wPDfQLThXg5TXt6ChyYRltydi8qDLb7FxvLQaSkFAdLCPgyokcgy7DaouKCjoUCFnzpzpUHsi6jz2ENlXN181VswbgRsTItFkEvHQJ7uxfo88P+Ne/PYQGppMuKZfN6TGXbknq2+IL8MQub0OBaIRI0bgvvvuw6+//tpmm8rKSrz33nsYNGgQ1q1b1+UCiah9zIOq/TmGyG48lAosnpGAPyX2gEkE/m9tLtbuktdM25+OlWLzoWIoFQIWTYvj7S+idurQJ+ehQ4fw97//HRMnToRGo0FiYiIiIyOh0WhQXl6OQ4cO4eDBgxg2bBheeeUVXHfddfaqm4h+hzvdO4ZSIeCVPw6BykOBT34pwF/+uw+NRhNuS+oldWloNJrw3NeHAAB3XN2Ls8SIOqBDPUTdunXD4sWLUVRUhLfffhv9+/dHWVkZjh07BgC47bbbkJOTg+zsbIYhIgfjLTPHUSgE/H36IMy7JhoA8LcvDmD5tny7PV99kxFbfivFd4eKYbiwgW9r/pN9Cnkl1Qj09sT/pQywWz1ErqhT/5T08vJCbGwsbr75ZigUrrHXD5GzYw+RYwmCgIXXx0LlocC7W07g+W8OocFowv1j+9rk+vq6RmQdLcWmgzpkHSmBocEIoHkH8jH9Q5AaF4aUq8IQeGFjzHPV9Xj9+98AAE+kxkDrzWBM1BGd/uQcNGgQNBoNYmNjER8fb3UEBATYsEQiao+LgYi/CB1FEAQsmDwQag8l3sw8hpf+dwTVdU14cHw/eKk6vrpzSVUdvjtUjM0Hi7H9eBkajRcnAYf5q6HyUOD0+Vp8f7gY3x9uHieU1DsIqXHhyD1dgaq6JsRG+GPWiJ62fJlEbqHTgWjLli245ZZb0L17d1RVVeG9997DwYMHIQgC+vXrhxkzZuDxxx9nOCJyEG7dIQ1BEJA2cQDUHgq8uuko3v4xDx9ln8TNw3rg1qSeVxzHc666HhsO6PD13rP49eR5tFwIpW+IDybFhSM1LhxDumshCMARXRU2HtBh08HmhSK3Hz+H7cfPWR7z7A1xUNp4R3Qid9DpdYiGDRuGZ555BjfddJPlXGZmJu677z7MnTsX33//PQoKCrBz506EhITYrGApcB0ikjtRFNH3qQ0wicDOpya4xA7tzmjtr6fx5g/HUFheazk3vFcgbk3qiesGR1j2BNPXNWLzwWJ8vfcstuWVwWi6+DEcHxWA1LgwTIoNR79Q38s+36lzBmw+WIyNB3XYXVCOmcOj8NIfh9jnxRE5IYfsdu/t7Y3c3FwMGGA9cO/rr7/GypUr8dlnn2HGjBkICAjAe++915mnkA0GIpK76vomDFq0CQBw+PnJnbpdQ7ZhMon4Ka8Mn/xyCt8fLrGEHa2XJ6YnREKnr8OPR0vR0GSyPGZwdy2mxUfg+iGRiAzw6tTz1jUaofZQcJo9UQsd+f3d6b71xMREfPzxx3juueeszg8aNAibN2+GIAh48sknLZvAEpH9mG+XeSgEaDw50UFKCoWAsQNCMHZACIr1dVj762ms/vU0zlTUYmX2KUu7fqG+uCE+EtcPiUCfkMv3BLUHd6Qn6ppOB6LXXnsNKSkpOHHiBP72t79h4MCBaGhowOuvv46goCAAQEhICIqLi21WLBG1ruUMM/YQyEeYvwYPT+iPP4/vh62/leLb/UUI8VPjhvhIDAz34/8rIhnpdCBKSkpCdnY2Hn30UcTGxkKtVqOpqQkeHh748MMPAQB79uxBZGSkzYolotZxDSJ5UyoEjB8YivEDQ6UuhYja0KXpKIMGDUJmZiYKCgqQm5sLpVKJxMREhIc3byQYEhKCl156ySaFElHb9FyDiIioS2zy6dmzZ0/07HnpuhdjxoyxxeWJ6Aq4KCMRUddw9CWRC+AtMyKirmEgInIB7CEiIuoaBiIiF1B9IRD5s4eIiKhTGIiIXID5lpmvmj1ERESdwUBE5AJ4y4yIqGsYiIhcgJ473RMRdQkDEZEL4E73RERdw0BE5AJ4y4yIqGsYiIhcQFU91yEiIuoKBiIiF1BlmXbPHiIios5gICJycqIotrhlxh4iIqLOYCAicnK1jUYYTSIAjiEiIuosBiIiJ2fuHVIqBHirlBJXQ0TknBiIiJycORD5qj0gCILE1RAROScGIiInx207iIi6TpaBaOnSpYiOjoZGo0FSUhJ27tzZZttx48ZBEIRLjqlTp1ra3HnnnZd8f/LkyY54KUR2xzWIiIi6TnafoGvWrEFaWhqWLVuGpKQkLFmyBKmpqTh69ChCQ0Mvaf/555+joaHB8vW5c+cQHx+PW265xard5MmT8eGHH1q+VqvV9nsRRA5UxZ3uiYi6THY9RIsXL8Y999yDefPmITY2FsuWLYO3tzeWL1/eavugoCCEh4dbju+++w7e3t6XBCK1Wm3VLjAw0BEvh8juuG0HEVHXySoQNTQ0ICcnBykpKZZzCoUCKSkpyM7Obtc1PvjgA8yaNQs+Pj5W57OyshAaGoqYmBg88MADOHfuXJvXqK+vh16vtzqI5Iq3zIiIuk5WgaisrAxGoxFhYWFW58PCwqDT6a74+J07d+LAgQOYP3++1fnJkyfjo48+QmZmJl5++WVs2bIFU6ZMgdFobPU6GRkZ0Gq1liMqKqrzL4rIzi72EPGWGRFRZ7nUPyk/+OADDB48GCNHjrQ6P2vWLMufBw8ejCFDhqBv377IysrChAkTLrlOeno60tLSLF/r9XqGIpItPXuIiIi6TFY9RMHBwVAqlSguLrY6X1xcjPDw8Ms+1mAwYPXq1bj77ruv+Dx9+vRBcHAw8vLyWv2+Wq2Gv7+/1UEkV9y2g4io62QViFQqFRITE5GZmWk5ZzKZkJmZieTk5Ms+9rPPPkN9fT1uv/32Kz5PYWEhzp07h4iIiC7XTCQ1DqomIuo6WQUiAEhLS8N7772HlStX4vDhw3jggQdgMBgwb948AMCcOXOQnp5+yeM++OADTJ8+Hd26dbM6X11djSeffBI7duzAyZMnkZmZiRtvvBH9+vVDamqqQ14TkT1xUDURUdfJ7hN05syZKC0txcKFC6HT6ZCQkICNGzdaBloXFBRAobDOcUePHsW2bduwefPmS66nVCqxb98+rFy5EhUVFYiMjMSkSZPwwgsvcC0icgnV9QxERERdJYiiKEpdhNzp9XpotVpUVlZyPBHJzrhXf8TJczX47P5kjIgOkrocIiLZ6Mjvb9ndMiOijuEtMyKirmMgInJynGVGRNR1DERETqyu0YgGowkAe4iIiLqCgYjIiZl7hwQB8FUxEBERdRYDEZETM69B5KvygEIhSFwNEZHzYiAicmIcUE1EZBsMREROjAOqiYhsg4GIyImV1zQAAAK8GYiIiLqCgYjIiZkDUaC3SuJKiIicGwMRkRMrNzQPqg70YSAiIuoKBiIiJ3axh4i3zIiIuoKBiMiJmQNREHuIiIi6hIGIyImdN5gHVTMQERF1BQMRkROrqGkeQxTkw1tmRERdwUBE5MQuTrtnDxERUVcwEBE5sfILt8yCGIiIiLqEgYjISdU3GWFoMALgOkRERF3FQETkpMzjh5QKgXuZERF1EQMRkZOyjB/y8uRO90REXcRAROSkzFPuuUo1EVHXMRAROSnzLTOuUk1E1HUMREROiosyEhHZDgMRkZOqqOGUeyIiW2EgInJS5y/sdB/AVaqJiLqMgYjISbGHiIjIdhiIiJzU+QuBiIsyEhF1HQMRkZMqN88y47R7IqIuYyAiclLmfcw47Z6IqOsYiIiclHmlavYQERF1HQMRkRNqNJpQVdcEgGOIiIhsgYGIyAmZV6kWBEDrxVtmRERdxUBE5ITMU+61Xp5QcmNXIqIuYyAickLmbTu4BhERkW0wEBE5IfOU+wDOMCMisglZBqKlS5ciOjoaGo0GSUlJ2LlzZ5ttV6xYAUEQrA6NRmPVRhRFLFy4EBEREfDy8kJKSgqOHTtm75dBZDfmGWZBnGFGRGQTsgtEa9asQVpaGhYtWoTdu3cjPj4eqampKCkpafMx/v7+KCoqshynTp2y+v4rr7yCN998E8uWLcMvv/wCHx8fpKamoq6uzt4vh8guzIGIO90TEdmG7ALR4sWLcc8992DevHmIjY3FsmXL4O3tjeXLl7f5GEEQEB4ebjnCwsIs3xNFEUuWLMHTTz+NG2+8EUOGDMFHH32Es2fPYv369Q54RUS2x0UZiYhsS1aBqKGhATk5OUhJSbGcUygUSElJQXZ2dpuPq66uRq9evRAVFYUbb7wRBw8etHwvPz8fOp3O6pparRZJSUltXrO+vh56vd7qIJITbttBRGRbsgpEZWVlMBqNVj08ABAWFgadTtfqY2JiYrB8+XJ8+eWXWLVqFUwmE0aNGoXCwkIAsDyuI9fMyMiAVqu1HFFRUV19aUQ2dbGHiIGIiMgWZBWIOiM5ORlz5sxBQkICxo4di88//xwhISF49913O33N9PR0VFZWWo7Tp0/bsGKirivnTvdERDYlq0AUHBwMpVKJ4uJiq/PFxcUIDw9v1zU8PT0xdOhQ5OXlAYDlcR25plqthr+/v9VBJCeWW2YcQ0REZBOyCkQqlQqJiYnIzMy0nDOZTMjMzERycnK7rmE0GrF//35EREQAAHr37o3w8HCra+r1evzyyy/tviaR3HDaPRGRbXlIXcDvpaWlYe7cuRg+fDhGjhyJJUuWwGAwYN68eQCAOXPmoHv37sjIyAAAPP/887j66qvRr18/VFRU4NVXX8WpU6cwf/58AM0z0B577DG8+OKL6N+/P3r37o1nnnkGkZGRmD59ulQvk6jTjCYRlbXmhRkZiIiIbEF2gWjmzJkoLS3FwoULodPpkJCQgI0bN1oGRRcUFEChuNixVV5ejnvuuQc6nQ6BgYFITEzE9u3bERsba2nzl7/8BQaDAffeey8qKiowevRobNy48ZIFHImcQWVtI0Sx+c9cqZqIyDYEUTR/tFJb9Ho9tFotKisrOZ6IJJdXUo2UxVvgp/HA/mdTpS6HiEi2OvL7W1ZjiIjoyio4foiIyOYYiIiczMWNXRmIiIhshYGIyMmYF2UM4vghIiKbYSAicjJclJGIyPYYiIiczHnudE9EZHMMREROpsLQPIYoyIe3zIiIbIWBiMjJsIeIiMj2GIiInAyn3RMR2R4DEZGTOW8w9xDxlhkRka0wEBE5mYoa8xgi9hAREdkKAxGREzGZRE67JyKyAwYiIidSVdcEEzd2JSKyOQYiIidinmHmo1JC7aGUuBoiItfBQETkRCy3yzh+iIjIphiIiJyIeR8zjh8iIrItBiIiJ2Le6Z49REREtsVAROREKiwzzDigmojIlhiIiJzIed4yIyKyCwYiIidiuWXGQEREZFMMREROxDKomjvdExHZFAMRkRPhKtVERPbBQETkRBiIiIjsg4GIyIlcnHbPW2ZERLbEQETkJERR5MKMRER2wkBE5CSq65vQdGFnVwYiIiLbYiAichLlhubbZRpPBbxU3NiViMiWGIiInIR5QHUQe4eIiGyOgYjISZy/EIgCGIiIiGyOgYjISZj3MQvixq5ERDbHQETkJM5fGEMUwI1diYhsjoGIyEmwh4iIyH4YiIicRDnHEBER2Q0DEZGTME+7D+ItMyIim2MgInISln3MeMuMiMjmZBmIli5diujoaGg0GiQlJWHnzp1ttn3vvfcwZswYBAYGIjAwECkpKZe0v/POOyEIgtUxefJke78MIps6b+AtMyIie5FdIFqzZg3S0tKwaNEi7N69G/Hx8UhNTUVJSUmr7bOysjB79mz8+OOPyM7ORlRUFCZNmoQzZ85YtZs8eTKKioosx6effuqIl0NkMxU15ltmDERERLYmu0C0ePFi3HPPPZg3bx5iY2OxbNkyeHt7Y/ny5a22//jjj/HnP/8ZCQkJGDhwIN5//32YTCZkZmZatVOr1QgPD7ccgYGBjng5RDYhimKLhRk5hoiIyNZkFYgaGhqQk5ODlJQUyzmFQoGUlBRkZ2e36xo1NTVobGxEUFCQ1fmsrCyEhoYiJiYGDzzwAM6dO9fmNerr66HX660OIinVNhrR0GQCwGn3RET2IKtAVFZWBqPRiLCwMKvzYWFh0Ol07brGX//6V0RGRlqFqsmTJ+Ojjz5CZmYmXn75ZWzZsgVTpkyB0Whs9RoZGRnQarWWIyoqqvMvisgGzOOHVEoFvLmxKxGRzXlIXYAtvfTSS1i9ejWysrKg0Wgs52fNmmX58+DBgzFkyBD07dsXWVlZmDBhwiXXSU9PR1pamuVrvV7PUESSMo8fCvTxhCAIEldDROR6ZNVDFBwcDKVSieLiYqvzxcXFCA8Pv+xjX3vtNbz00kvYvHkzhgwZctm2ffr0QXBwMPLy8lr9vlqthr+/v9VBJCVzD1EgB1QTEdmFrAKRSqVCYmKi1YBo8wDp5OTkNh/3yiuv4IUXXsDGjRsxfPjwKz5PYWEhzp07h4iICJvUTWRvljWIGIiIiOxCVoEIANLS0vDee+9h5cqVOHz4MB544AEYDAbMmzcPADBnzhykp6db2r/88st45plnsHz5ckRHR0On00Gn06G6uhoAUF1djSeffBI7duzAyZMnkZmZiRtvvBH9+vVDamqqJK+RqKPKzT1EPpxhRkRkD7IbQzRz5kyUlpZi4cKF0Ol0SEhIwMaNGy0DrQsKCqBQXMxx77zzDhoaGvCnP/3J6jqLFi3Cs88+C6VSiX379mHlypWoqKhAZGQkJk2ahBdeeAFqtdqhr42os8rNY4jYQ0REZBeCKIqi1EXInV6vh1arRWVlJccTkSQWfnkAH2WfwkPj++GJ1BipyyEicgod+f0tu1tmRHQpSw8R1yAiIrILBiIiJ1BhGVTNMURERPbAQETkBDjtnojIvhiIiJxABW+ZERHZFQMRkRO42EPEW2ZERPbAQEQkc3WNRtQ2Nu+7xx4iIiL7YCAikjnzKtUeCgF+atktHUZE5BIYiIhkrtzQPH4owFvFjV2JiOyEgYhI5so55Z6IyO4YiIhkzhKIOH6IiMhuGIiIZK6cM8yIiOyOgYhI5szbdgSxh4iIyG4YiIhkzrwGUQBXqSYishsGIiKZM+9jFsRARERkNwxERDJ3vsY87Z5jiIiI7IWBiEjmLD1EHENERGQ3DEREMscxRERE9sdARCRjJ0qrUVheC4UA9OrmLXU5REQui4GISMa+2HMGADCmfwiCfdUSV0NE5LoYiIhkymQS8fnu5kD0x8QeEldDROTaGIiIZOqX/PM4U1ELP7UHJsWGSV0OEZFLYyAikql1uwsBAFOHREDjqZS4GiIi18ZARCRDNQ1N+N/+IgC8XUZE5AgMREQytOmgDoYGI3oGeWN4r0CpyyEicnkMREQyZB5MffOw7hAEQeJqiIhcHwMRkcwUVdZiW14ZAODmobxdRkTkCAxERDKzfs9ZiCIwMjoIPbkYIxGRQzAQEcmIKIqW2WU3D+sucTVERO6DgYhIRvafqUReSTXUHgpcNyRC6nKIiNwGAxGRjKzLae4dmhQXDn+Np8TVEBG5DwYiIploaDLhq71nAQB/5O0yIiKHYiAikokfj5agvKYRoX5qjO4XLHU5RERuhYGISCY+vzCYevrQ7vBQ8q8mEZEj8VOXSAbKDQ344UgJAOCPw7j2EBGRo8kyEC1duhTR0dHQaDRISkrCzp07L9v+s88+w8CBA6HRaDB48GBs2LDB6vuiKGLhwoWIiIiAl5cXUlJScOzYMXu+BKIO+XrfWTQaRcRF+iMm3E/qcoiI3I7sAtGaNWuQlpaGRYsWYffu3YiPj0dqaipKSkpabb99+3bMnj0bd999N/bs2YPp06dj+vTpOHDggKXNK6+8gjfffBPLli3DL7/8Ah8fH6SmpqKurs5RL4vossyzy9g7REQkDUEURVHqIlpKSkrCiBEj8PbbbwMATCYToqKi8PDDD2PBggWXtJ85cyYMBgO++eYby7mrr74aCQkJWLZsGURRRGRkJB5//HE88cQTAIDKykqEhYVhxYoVmDVr1hVr0uv10Gq1qKyshL+/v41eKaCva4S+ttFm1yP5EkXAJIpoNJrQ0NT830ajCQ1GE0qr6vHo6lx4KATseGoCgn3VUpdLROQSOvL728NBNbVLQ0MDcnJykJ6ebjmnUCiQkpKC7OzsVh+TnZ2NtLQ0q3OpqalYv349ACA/Px86nQ4pKSmW72u1WiQlJSE7O7vVQFRfX4/6+nrL13q9visvq02rdpzCKxuP2uXa5HzGxYQwDBERSURWgaisrAxGoxFhYWFW58PCwnDkyJFWH6PT6Vptr9PpLN83n2urze9lZGTgueee69Rr6AgPhQC1h+zuWpKdeCgEeHoo4KlUQKVUwFMpwFPZ/LWfxgOPpQyQukQiIrclq0AkF+np6Va9Tnq9HlFRUTZ/nnuv7Yt7r+1r8+sSERFRx8iqeyI4OBhKpRLFxcVW54uLixEeHt7qY8LDwy/b3vzfjlxTrVbD39/f6iAiIiLXJatApFKpkJiYiMzMTMs5k8mEzMxMJCcnt/qY5ORkq/YA8N1331na9+7dG+Hh4VZt9Ho9fvnllzavSURERO5FdrfM0tLSMHfuXAwfPhwjR47EkiVLYDAYMG/ePADAnDlz0L17d2RkZAAAHn30UYwdOxb//Oc/MXXqVKxevRq7du3Cv//9bwCAIAh47LHH8OKLL6J///7o3bs3nnnmGURGRmL69OlSvUwiIiKSEdkFopkzZ6K0tBQLFy6ETqdDQkICNm7caBkUXVBQAIXiYsfWqFGj8Mknn+Dpp5/GU089hf79+2P9+vUYNGiQpc1f/vIXGAwG3HvvvaioqMDo0aOxceNGaDQah78+IiIikh/ZrUMkR/Zah4iIiIjspyO/v2U1hoiIiIhICgxERERE5PYYiIiIiMjtMRARERGR22MgIiIiIrfHQERERERuj4GIiIiI3B4DEREREbk9BiIiIiJye7LbukOOzIt56/V6iSshIiKi9jL/3m7PphwMRO1QVVUFAIiKipK4EiIiIuqoqqoqaLXay7bhXmbtYDKZcPbsWfj5+UEQhC5dS6/XIyoqCqdPn+a+aO3A96vj+J51DN+vjuH71XF8zzrGlu+XKIqoqqpCZGSk1cbwrWEPUTsoFAr06NHDptf09/fnX4wO4PvVcXzPOobvV8fw/eo4vmcdY6v360o9Q2YcVE1ERERuj4GIiIiI3B4DkYOp1WosWrQIarVa6lKcAt+vjuN71jF8vzqG71fH8T3rGKneLw6qJiIiIrfHHiIiIiJyewxERERE5PYYiIiIiMjtMRARERGR22MgsoOlS5ciOjoaGo0GSUlJ2Llz52Xbf/bZZxg4cCA0Gg0GDx6MDRs2OKhSeejI+7VixQoIgmB1aDQaB1Yrra1bt2LatGmIjIyEIAhYv379FR+TlZWFYcOGQa1Wo1+/flixYoXd65STjr5nWVlZl/yMCYIAnU7nmIIllpGRgREjRsDPzw+hoaGYPn06jh49esXHuevnWGfeL3f+HHvnnXcwZMgQy6KLycnJ+N///nfZxzjqZ4uByMbWrFmDtLQ0LFq0CLt370Z8fDxSU1NRUlLSavvt27dj9uzZuPvuu7Fnzx5Mnz4d06dPx4EDBxxcuTQ6+n4BzauXFhUVWY5Tp045sGJpGQwGxMfHY+nSpe1qn5+fj6lTp2L8+PHIzc3FY489hvnz52PTpk12rlQ+OvqemR09etTq5yw0NNROFcrLli1b8OCDD2LHjh347rvv0NjYiEmTJsFgMLT5GHf+HOvM+wW47+dYjx498NJLLyEnJwe7du3CH/7wB9x44404ePBgq+0d+rMlkk2NHDlSfPDBBy1fG41GMTIyUszIyGi1/YwZM8SpU6danUtKShLvu+8+u9YpFx19vz788ENRq9U6qDp5AyB+8cUXl23zl7/8RYyLi7M6N3PmTDE1NdWOlclXe96zH3/8UQQglpeXO6QmuSspKREBiFu2bGmzjbt/jrXUnveLn2PWAgMDxffff7/V7znyZ4s9RDbU0NCAnJwcpKSkWM4pFAqkpKQgOzu71cdkZ2dbtQeA1NTUNtu7ks68XwBQXV2NXr16ISoq6rL/siD3/vnqqoSEBERERGDixIn4+eefpS5HMpWVlQCAoKCgNtvw5+yi9rxfAD/HAMBoNGL16tUwGAxITk5utY0jf7YYiGyorKwMRqMRYWFhVufDwsLaHH+g0+k61N6VdOb9iomJwfLly/Hll19i1apVMJlMGDVqFAoLCx1RstNp6+dLr9ejtrZWoqrkLSIiAsuWLcO6deuwbt06REVFYdy4cdi9e7fUpTmcyWTCY489hmuuuQaDBg1qs507f4611N73y90/x/bv3w9fX1+o1Wrcf//9+OKLLxAbG9tqW0f+bHG3e3IqycnJVv+SGDVqFK666iq8++67eOGFFySsjFxFTEwMYmJiLF+PGjUKx48fx+uvv47//Oc/ElbmeA8++CAOHDiAbdu2SV2KU2jv++Xun2MxMTHIzc1FZWUl/vvf/2Lu3LnYsmVLm6HIUdhDZEPBwcFQKpUoLi62Ol9cXIzw8PBWHxMeHt6h9q6kM+/X73l6emLo0KHIy8uzR4lOr62fL39/f3h5eUlUlfMZOXKk2/2MPfTQQ/jmm2/w448/okePHpdt686fY2Ydeb9+z90+x1QqFfr164fExERkZGQgPj4eb7zxRqttHfmzxUBkQyqVComJicjMzLScM5lMyMzMbPP+aHJyslV7APjuu+/abO9KOvN+/Z7RaMT+/fsRERFhrzKdmjv/fNlSbm6u2/yMiaKIhx56CF988QV++OEH9O7d+4qPceefs868X7/n7p9jJpMJ9fX1rX7PoT9bNh+m7eZWr14tqtVqccWKFeKhQ4fEe++9VwwICBB1Op0oiqJ4xx13iAsWLLC0//nnn0UPDw/xtddeEw8fPiwuWrRI9PT0FPfv3y/VS3Cojr5fzz33nLhp0ybx+PHjYk5Ojjhr1ixRo9GIBw8elOolOFRVVZW4Z88ecc+ePSIAcfHixeKePXvEU6dOiaIoigsWLBDvuOMOS/sTJ06I3t7e4pNPPikePnxYXLp0qahUKsWNGzdK9RIcrqPv2euvvy6uX79ePHbsmLh//37x0UcfFRUKhfj9999L9RIc6oEHHhC1Wq2YlZUlFhUVWY6amhpLG36OXdSZ98udP8cWLFggbtmyRczPzxf37dsnLliwQBQEQdy8ebMoitL+bDEQ2cFbb70l9uzZU1SpVOLIkSPFHTt2WL43duxYce7cuVbt165dKw4YMEBUqVRiXFyc+O233zq4Yml15P167LHHLG3DwsLE6667Tty9e7cEVUvDPCX894f5PZo7d644duzYSx6TkJAgqlQqsU+fPuKHH37o8Lql1NH37OWXXxb79u0rajQaMSgoSBw3bpz4ww8/SFO8BFp7rwBY/dzwc+yizrxf7vw5dtddd4m9evUSVSqVGBISIk6YMMEShkRR2p8tQRRF0fb9TkRERETOg2OIiIiIyO0xEBEREZHbYyAiIiIit8dARERERG6PgYiIiIjcHgMRERERuT0GIiIiInJ7DERERETk9hiIiIiIyO0xEBERtdDU1CR1CUQkAQYiInJbJ0+ehCAIWLt2LcaMGQO1Wo2vvvpK6rKISAIeUhdARCSVvXv3AgBeffVV/OMf/0Dv3r0REhIicVVEJAUGIiJyW7m5ufDx8cFnn32G6OhoqcshIgnxlhkRua29e/fihhtuYBgiIgYiInJfubm5GDdunNRlEJEMMBARkVvS6/U4efIkhg4dKnUpRCQDDERE5Jb27t0LpVKJwYMHS10KEckAAxERuaW9e/ciJiYGGo1G6lKISAYEURRFqYsgIiIikhJ7iIiIiMjtMRARERGR22MgIiIiIrfHQERERERuj4GIiIiI3B4DEREREbk9BiIiIiJyewxERERE5PYYiIiIiMjtMRARERGR22MgIiIiIrfHQERERERu7/8BjUGjQ++6lpQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import freud\n",
    "import gsd.hoomd\n",
    "\n",
    "rdf = freud.density.RDF(bins=60, r_max=3.0)\n",
    "with gsd.hoomd.open(result.directory.file(\"trajectory.gsd\")) as traj:\n",
    "    for snap in traj:\n",
    "        rdf.compute(\n",
    "            (snap.configuration.box, snap.particles.position), reset=False\n",
    "        )\n",
    "rdf.plot();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks like a typical correlation function for a liquid!"
   ]
  }
 ],
 "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.11.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "6834e2eb4702651dab439b72249722ae87da91a96342472ed2409ca0cf3240b1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}