{
 "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": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHFCAYAAAAT5Oa6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMs0lEQVR4nO3de1xUdf4/8NcZLjOgzAioXAIRTUXUCEEFTctMTEvtpq5tZq3m+tvaMr+7W3Tbtd2ydiu1vJS7Fl+3Qkq8tekqflPxQpkGVpp3FEQQkctwEYaZOb8/cEZGLnKZmTNnzuv5eMyD5viZw3umI7z8nM9FEEVRBBEREZGCqaQugIiIiEhqDERERESkeAxEREREpHgMRERERKR4DERERESkeAxEREREpHgMRERERKR4DERERESkeAxEREREpHgMRETkFlJSUiAIgvXh6emJkJAQ/OpXv8KpU6ds2t51113WdiqVCn5+frj11lsxbdo0rF+/Hmazucn5e/fubXP+xo+qqipnvU0ichBPqQsgIrKnTz75BFFRUaitrcX+/fvxxhtvYNeuXTh+/Dj8/f2t7fr06YPPPvsMAFBdXY3c3Fxs2rQJ06ZNw+jRo/HVV19Bp9PZnHvUqFF45513mnxPX19fx74pInI4BiIiciuDBw9GfHw8gIaeIJPJhD//+c/YtGkTnnzySWs7Hx8fJCQk2Lx27ty5+OSTT/Cb3/wG8+bNQ1pams2fd+vWrclriMg98JYZEbk1Szi6dOlSm9o/+eSTmDRpEr788kucP3/ekaURkQthICIit5abmwsA6N+/f5tfM2XKFIiiiL1799ocF0URRqPR5tHceCMikh8GIiJyKyaTCUajEVVVVdi+fTv+9re/YcyYMZgyZUqbzxEREQEAuHjxos3xrVu3wsvLy+bx2muv2bV+IpIGxxARkVu5cYzPwIEDsXnzZnh6tv3HnSiKzR6/4447sGTJEptjoaGh7S+SiFwOAxERuZW1a9di4MCBqKysRFpaGj766CPMnDkT27Zta/M5LGOHbgw7Op3OOiaJiNwLAxERuZWBAwdaQ8vYsWNhMpnwr3/9C+vXr8cjjzzSpnNs2bIFgiBgzJgxjiyViFwIxxARkVv7+9//Dn9/f7z22mttGgD9ySefYNu2bZg5cyZ69erlhAqJyBWwh4iI3Jq/vz+Sk5Pxpz/9CZ9//jkee+wxAMDVq1fx7bffWv/77Nmz2LRpE/7zn//gzjvvxIcffihl2UTkZAxEROT2fv/732P58uV4/fXXMXPmTADA2bNnkZiYCADo0qULgoKCMHToUHz55Zd46KGHoFKxA51ISQSxpekURERERArBfwIRERGR4jEQERERkeIxEBEREZHiMRARERGR4jEQERERkeIxEBEREZHicR2iNjCbzbh48SL8/PwgCILU5RAREVEbiKKIyspKhIaG3nRtMQaiNrh48SLCw8OlLoOIiIg6ID8/H2FhYa22YSBqAz8/PwANH6hWq5W4GiIiImoLvV6P8PBw6+/x1jAQtYHlNplWq2UgIiIikpm2DHfhoGoiIiJSPAYiIiIiUjwGIiIiIlI8BiIiIiJSPAYiIiIiUjwGIiIiIlI8BiIiIiJSPAYiIiIiUjwGIiIiIlI8BiIiIiJSPAYiIiIiUjwGIiIiIlI8BiIiapHZLKLeZJa6DCIih2MgIqJm/XihHCPf+gYPrNiP2nqT1OUQETkUAxERNZF15gpmrv4WRfpaHL2ox+ff5UldEhGRQzEQEZGNjGOXMPuTg6g2mBCq0wAAVuw6jeo6o8SVERE5DgMREVmlH76A+Z8ehsFoRlJ0EDIW3onegb64Um3AJ/tzpS6PiMhhGIiICADwyf5c/M+XR2Ayi3h4aBhW/noouqg98fz4/gCAjzLPorzGIHGVRESOwUBEpHCiKGJJxkks+uoYAODJUb3xj0dug6dHw4+HybeFIirYD5W1RnyUeVbKUomIHIaBiEjBzGYRi746hmX/dwoAsHB8f7x2fzRUKsHaRqUS8IekAQAaepGKK2slqZWIyJEYiIgUbPmu00g5cA4A8JfJ0Xh2XD8IgtCk3biBPRHbqxtq681YueuMk6skInI8BiIiBdt1ohgA8NKkKDwxKrLFdoIg4I8TGnqJPvvuPC6U1TilPiIiZ2EgIlKwy5V1AIC4iICbth3ZtztG3RqIepOIZTtPObo0IiKnYiAiUihRFFF8LRD19FO36TWWsUTpP1zA6eIqh9VGRORsDERECqW/aoTB2LBPWY82BqLYXv4YHx0EswgsyTjpyPKIiJyKgYhIoS5XNcwW02o8ofHyaPPr/iepPwQB+PqnQvxcUOGo8oiInMqlAlFmZiYmT56M0NBQCIKATZs2tdr+iSeegCAITR6DBg2ytklJSWm2TW0tpw6TshXrr90u02ra9bqoYC2mxoQCAN7dccLudRERScGlAlF1dTViYmKwfPnyNrVftmwZCgsLrY/8/HwEBARg2rRpNu20Wq1Nu8LCQmg07fslQORuLOOHenRt2+2yxhbc0x+eKgG7TlzG9+dK7V0aEZHTeUpdQGMTJ07ExIkT29xep9NBp9NZn2/atAllZWV48sknbdoJgoDg4GC71UnkDiwzzHpq2x+Ienfvgmnx4Ug9mIfPvj2PYb1vPkuNiMiVuVQPUWetWbMG99xzDyIiImyOV1VVISIiAmFhYbj//vuRnZ3d6nnq6uqg1+ttHkTuxrLidFtnmN3onoE9AQDHiyrtVhMRkVTcJhAVFhZi27ZtmDt3rs3xqKgopKSkYMuWLUhNTYVGo8GoUaNw6lTL66gsXrzY2vuk0+kQHh7u6PKJnM7SQ9TWGWY36h/kBwA4e7kaRpPZbnUREUnBbQJRSkoKunXrhgceeMDmeEJCAh577DHExMRg9OjR+OKLL9C/f3988MEHLZ4rOTkZFRUV1kd+fr6DqydyvutrEHVsPN0t3Xzg4+UBg8mMc1e4cjURyZtbBCJRFPHxxx9j1qxZ8Pb2brWtSqXCsGHDWu0hUqvV0Gq1Ng8id9PeRRlvpFIJ6BfUFQBw6hJvmxGRvLlFINqzZw9Onz6NOXPm3LStKIrIyclBSEiIEyojcl2dvWUGAP16Ntw2O3mJq1YTkby51CyzqqoqnD592vo8NzcXOTk5CAgIQK9evZCcnIyCggKsXbvW5nVr1qzBiBEjMHjw4CbnXLRoERISEtCvXz/o9Xq8//77yMnJwYoVKxz+fohcVW29CRVX6wF0/JYZAPS/1kN0spg9REQkby4ViA4dOoSxY8dany9cuBAAMHv2bKSkpKCwsBB5eXk2r6moqEB6ejqWLVvW7DnLy8sxb948FBUVQafTITY2FpmZmRg+fLjj3giRi7P0Dnl7qqD16fiPAcvAat4yIyK5c6lAdNddd0EUxRb/PCUlpckxnU6HmpqWB3QuWbIES5YssUd5RG7jctX1RRkFQejweSxjiHJLqlFvMsPLwy3uwhORAvGnF5ECXd+2o+Pjh4CGmWZdvD1QbxJxrqTaHqUREUmCgYhIgRr3EHWGIAjoF8SB1UQkfwxERAp0WX9tlepO9hABjQZWcxwREckYAxGRAnV2UcbGrAOrOdOMiGSMgYhIgeyxBpEFb5kRkTtgICJSoM6uUt2Y5ZbZuZJq1BlNnT4fEZEUGIiIFMiePUTBWg381J4wmkXkcqYZEckUAxGRwpjNIkqq7DeGqGGmmWVgNW+bEZE8MRARKUxpjQFGswhBAAK7tr4ZcltxxWoikjsGIiKFsdwuC/D1ttvK0tcHVjMQEZE8MRARKUyxHccPWVgGVp/iLTMikikGIiKFKbYuytj58UMWlltm565Uo7aeM82ISH4YiIgUxl7bdjTW008NrcYTZhE4e5kzzYhIfhiIiBTGXhu7NiYIAlesJiJZYyAiUhhH9BABHFhNRPLGQESkMJcd0EMENN7klQOriUh+GIiIFKa48tqgajssytgY1yIiIjljICJSGHtu29GYZbXq86U1nGlGRLLDQESkINV1RlQbGsKKPTZ2baxHVzW6+XpBFIHTxbxtRkTywkBEpCCWRRm7eHugi9rTrucWBAH9e3KmGRHJEwMRkYI46naZBTd5JSK5YiAiUhBHDai24MBqIpIrBiIiBWEPERFR8xiIiBTEERu7NjbgWg9RflkNrho404yI5IOBiEhBHLFtR2OBXdUI7OLNmWZEJDsMREQK4qhtOxq7ftuM44iISD4YiIgUpFh/bVC11jGDqoHrA6tPcuo9EckIAxGRglgGVdt7UcbG+llnmvGWGRHJBwMRkULUm8worTEAcNygagDo37PhltmJIvYQEZF8MBARKcSVKgNEEfBQCQjw9XbY97HcMisov4rqOqPDvg8RkT0xEBEphOV2Wfeu3lCpBId9H/8u3uh+bdD2Kc40IyKZYCAiUghHr1LdWH/ONCMimWEgIlKIYicMqLbgFh5EJDcMREQK4ehtOxrjFh5EJDcuFYgyMzMxefJkhIaGQhAEbNq0qdX2u3fvhiAITR7Hjx+3aZeeno7o6Gio1WpER0dj48aNDnwXRK7p+i0z9hAREd3IpQJRdXU1YmJisHz58na97sSJEygsLLQ++vXrZ/2zrKwszJgxA7NmzcKRI0cwa9YsTJ8+Hd999529yydyadYeIgcuymjRv2dDILpYUYvK2nqHfz8ios7ylLqAxiZOnIiJEye2+3U9e/ZEt27dmv2zpUuXYvz48UhOTgYAJCcnY8+ePVi6dClSU1M7Uy6RrFg3dnXgth0WOl8v9PRTo7iyDicvVSEuwt/h35OIqDNcqoeoo2JjYxESEoJx48Zh165dNn+WlZWFpKQkm2MTJkzAgQMHnFkikeQcvbHrjQaFagEAOfnlTvl+RESdIetAFBISgtWrVyM9PR0bNmzAgAEDMG7cOGRmZlrbFBUVISgoyOZ1QUFBKCoqavG8dXV10Ov1Ng8iORNF0SkbuzY2PDIQAHAw94pTvh8RUWe41C2z9howYAAGDBhgfZ6YmIj8/Hy88847GDNmjPW4INguQieKYpNjjS1evBiLFi2yf8FEEtFfNcJgNANwziwzABge2XCb7PtzZTf9O0dEJDVZ9xA1JyEhAadOnbI+Dw4ObtIbVFxc3KTXqLHk5GRUVFRYH/n5+Q6rl8gZLDPMdD5e0Hh5OOV7DrmlG9SeKpRWG3DmMqffE5Frc7tAlJ2djZCQEOvzxMREZGRk2LTZsWMHRo4c2eI51Go1tFqtzYNIzpy5BpGFt6cKsb26AQAO5pY57fsSEXWES90yq6qqwunTp63Pc3NzkZOTg4CAAPTq1QvJyckoKCjA2rVrATTMIOvduzcGDRoEg8GATz/9FOnp6UhPT7ee47nnnsOYMWPw9ttvY+rUqdi8eTN27tyJffv2Of39EUnFmatUNzY8MhDfni3FwdwreHREL6d+byKi9nCpQHTo0CGMHTvW+nzhwoUAgNmzZyMlJQWFhYXIy8uz/rnBYMAf/vAHFBQUwMfHB4MGDcLXX3+NSZMmWduMHDkS69atwyuvvIJXX30Vffv2RVpaGkaMGOG8N0YksctSBaLeAQAaxhEREbkyQRRFUeoiXJ1er4dOp0NFRQVvn5EsvfH1Mfxzby6eGh2Jl++Ldtr3rTEYcdtfdsBoFrHvhbEI8/d12vcmImrP72+3G0NERE1dv2Xm+FWqG/P19sSgW3QAgO/PlTr1exMRtQcDEZECWG+ZOWlRxsZGRDbcNjuYy0BERK6LgYhIAZy5bceNhvVmICIi18dARKQAxfprO91L0EM0rHfDAo1nLlej5Npq2UREroaBiMjN1daboK81AgB6dHXuGCIA6ObrjahgPwDAIY4jIiIXxUBE5OYs44e8PVXQ+kiz0oblttl3vG1GRC6KgYjIzVk2de3pp5ZsP7HhHFhNRC6OgYjIzRXrnb9tx40sgeiXQj30tfWS1UFE1BIGIiI3d/naxq7OXqW6sSCtBhGBvjCLwOHzXLWaiFwPAxGRm7ss0aKMN7KMI/qet82IyAUxEBG5uWIJdrpvDscREZErYyAicnNSbex6I8tGrz9eqEBtvUnSWoiIbsRAROTmXKWHKCLQFz391DCYzMjJL5e0FiKiGzEQEbm5YuugamnHEAmCgGGRHEdERK6JgYjIjZnNIkqqDACk2bbjRtaNXrliNRG5GAYiIjdWWmOAySxCEIDALt5Sl2OdaXb4fBmMJrPE1RARXcdAROTGLIsyBnbxhqeH9H/dBwT5QavxRI3BhKMX9VKXQ0RkJf1PSCJyGMu2HT0kHj9koVIJ1un33/O2GRG5EAYiIjdWcm2GWfeu0t8us+BGr0TkihiIiNyYZd8wnY+XxJVc17iHyGwWJa6GiKgBAxGRG9NfNQIAtC4UiAbfooOPlwfKa+px+nKV1OUQEQFgICJya5XXeoj8NJ4SV3Kdl4cKQyO6AeBtMyJyHQxERG6ssvZaD5HGdXqIAG70SkSuh4GIyI1ZxhBpXaiHCLg+jijr7BWYOI6IiFwAAxGRG7P0EPm5WA9RXIQ/uvl64XJlHTJPXpa6HCIiBiIid6Z3wTFEAKD29MDDQ8MAAJ8fzJO4GiIiBiIit2YdQ+RCs8wsZg4PBwB8c7wYl/S1EldDRErHQETkxlxxlpnFrT39MKy3P0xmEV8eype6HCJSOAYiIjdmWYfI1cYQWfxqWC8AwLrv87lIIxFJioGIyE3V1ptguLajvKvNMrO477YQaDWeuFB2FftOl0hdDhEpGAMRkZuyjB8SBKCLt2sGIo2XBx6MvQUAsO57Dq4mIukwEBG5KcsMs65qT6hUgsTVtGzmiIbbZjuOXsLla5vREhE5GwMRkZty1VWqbxQVrMXt4d1gNItI/+GC1OUQkUIxEBG5KVeeYXajR4dfG1x9MI+Dq4lIEgxERG7KutO9i/cQAcD9MSHoqvbEuSs1+PbsFanLISIFcqlAlJmZicmTJyM0NBSCIGDTpk2ttt+wYQPGjx+PHj16QKvVIjExEdu3b7dpk5KSAkEQmjxqa7kQHLk3Sw+R1sf1e4h8vT0x9fZQAEDq91yTiIicz6UCUXV1NWJiYrB8+fI2tc/MzMT48eOxdetWHD58GGPHjsXkyZORnZ1t006r1aKwsNDmodFoHPEWiFyGq+5j1pKZ126bbf+5CKXVBomrISKlcal/Ok6cOBETJ05sc/ulS5faPH/zzTexefNmfPXVV4iNjbUeFwQBwcHB9iqTSBZcdR+zlgy+RYcht+jwU0EFNvxwAXNH95G6JCJSEJfqIeoss9mMyspKBAQE2ByvqqpCREQEwsLCcP/99zfpQbpRXV0d9Hq9zYNIbuQyy6wxSy/R5wfzIIocXE1EzuNWgejdd99FdXU1pk+fbj0WFRWFlJQUbNmyBampqdBoNBg1ahROnTrV4nkWL14MnU5nfYSHhzujfCK7klsPEQBMuT0Uvt4eOHu5Gt+fK5O6HCJSELcJRKmpqfjLX/6CtLQ09OzZ03o8ISEBjz32GGJiYjB69Gh88cUX6N+/Pz744IMWz5WcnIyKigrrIz+fgzxJfqyzzFxwp/uWdFV7YkrMtcHVB7lyNRE5j1sEorS0NMyZMwdffPEF7rnnnlbbqlQqDBs2rNUeIrVaDa1Wa/Mgkhs5rUPUmOW22dc/FaK8hoOricg5ZB+IUlNT8cQTT+Dzzz/Hfffdd9P2oigiJycHISEhTqiOSDpym2VmcVuYDgNDtDAYzUg9yN5ZInIOlwpEVVVVyMnJQU5ODgAgNzcXOTk5yMtr6DpPTk7G448/bm2fmpqKxx9/HO+++y4SEhJQVFSEoqIiVFRUWNssWrQI27dvx9mzZ5GTk4M5c+YgJycH8+fPd+p7I3I2yxgiV93pviWCIODJkb0BAEt3nsSJokppCyIiRXCpQHTo0CHExsZap8wvXLgQsbGxeO211wAAhYWF1nAEAB999BGMRiOefvpphISEWB/PPfectU15eTnmzZuHgQMHIikpCQUFBcjMzMTw4cOd++aInEyuPUQAMC0+DHcN6IE6oxnPpmajtt4kdUlE5OYEkXNbb0qv10On06GiooLjiUgWRFFE35e2wiwCB18ah55a+S1EWlJVh3uX7kVJVR0eT4zA61MHS10SEclMe35/u1QPERHZR7XBBMseqXKaZdZY965qvDs9BgCwNus8Mo5dkrgiInJnDEREbsgyw8zLQ4DaU75/ze/s3wNz74gEAPxp/REUVXAPQiJyDPn+pCSiFlnWIPLTeEEQBImr6Zw/3jsAg0K1KKupx8IvcmAy8y4/EdkfAxGRG6qU6Qyz5qg9PfD+zFj4eHngwJkr+CjzjNQlEZEbYiAickNynmHWnL49umLRlEEAgPd2nEROfrm0BRGR22EgInJDctzH7GamxYfhvttCYDSLeDY129oLRkRkDwxERG5IL8Od7m9GEAS8+eAQ3NLNB3mlNfjz5qNSl0REboSBiMgNyXUfs5vR+Xhh2a9uh0oANmQX4GBuqdQlEZGbYCAickONZ5m5m/jeAZgxLBwAsDbrnLTFEJHbYCAickPWWWY+7tVDZPFYQgQA4L8/F6G4kmsTEVHnMRARuSF3m2V2o0GhOsRF+MNoFpF2MF/qcojIDTAQEbkhd5xldqNZ13qJPj+YB6PJLHE1RCR3DEREbqjSDWeZ3WjikGAEdvFGYUUtdv5SLHU5RCRzDEREbsidVqpuidrTwzq4+tNvz0tcDRHJHQMRkRty51lmjT06ohcEAdh3ugRnLldJXQ4RyRgDEZEbcvdZZhZh/r4YF9UTAPDZt3kSV0NEcsZARORmjCYzqg0mAO7fQwRcn4L/5eF81BiMEldDRHLFQETkZqrqrocCd55lZjGmXw9EBPqistaILTkXpS6HiGSKgYjIzVhmmPl4ecDLw/3/iqtUAh4b0dBLtDbrPERRlLgiIpIj9/9pSaQwFVfdfw2iGz0SFwa1pwrHCvXIzi+XuhwikiEGIiI3c32VauUEIv8u3pgcEwoA+DSLU/CJqP0YiIjczPUZZu4/oLoxy8rV//mxEFeq6iSuhojkhoGIyM3o3Xwfs5bEhHfDbWE6GExmfHHogtTlEJHMMBARuRklrFLdEssU/M++Ow+TmYOriajtGIiI3Iy773TfmikxodD5eOFC2VXsOcn9zYio7RiIiNyM/qpye4g0Xh6YHh8GAPhk/zlpiyEiWWEgInIz1p3uFTao2mJWQm94qgTsPVWCXSfYS0REbcNARORmKuuUtw5RY70CffHkqN4AgL9+dQwGo1nagohIFpT5E5PIjV3f6V65f72fHdcPG7Mv4mxJNT7en4v5d/aVuiSys7wrNfj79uM4eakSQVoNQnQaBOt8rn1teB6i9YFKBRiMZtQZzTAYzTCYzNeem2AWAX9fb/ToqobWxxOCIEj9tkhCyv2JSeSmrs8yU+YtM6BhQPmLE6Pwhy+P4IP/O4UHY29BkFYjdVnUDMtWK20NIzUGI1buOoPVe89ae/9OXqrqdB2eKgGBXb0R2EWN7n5qdO/ijfjeAZg5PJxBSSEYiIjcjJJnmTX2UOwt+Oy788jOK8db245jyYzbpS6JrtHX1uP/frmErT8VYc/Jy+jRVY3JMaGYensoooL9mg0goijiqx8LsXjrLyisqAUA3HFrdzwxsjfKagwoqqhFob624WtFLYoqrqKspt76em8PFbw9Gx7qa18FAKXVBuhrjTCaRVzS1+GSvg4obHjNhuwCfHv2Cv7+yG3QeHk446MhCTEQEbkZfa2yxxBZqFQCFk0ZhKkr9mNjdgF+PaIX4nsHSF2WYlXU1GPHsSJs+7kI+06VwGC6PraroPwqPtxzBh/uOYN+Pbti6u2hmBJzC3oF+gIAjl3U4y9bjuLguVIAQJi/D165LxoTBgW12ntTW28C0BCGVKqW29UZTSitNuBKlQGXq+pwpcqAs5ersDrzLLYcuYiC8qtYPSsOgV3V9vgoyEUJIreGvim9Xg+dToeKigpotVqpyyFqVf9XtsFgNGP/i3fjlm4+UpcjuRfTf8S67/MxKFSLLc/cAY9WfjGS/e05eRkf78vF/tMlMDZaLLNvjy64b0gIkgYFI6+0BltyLuKb48U2Qen28G6I7N4Fm3MKYBYBjZcKv7vrVswb08cpPTYHTpdg/qeHoa81IjzABx/PHoZ+QX4O/75kP+35/c1A1AYMRCQXtfUmRL36XwDAj39JUvQ4IouSqjqMfWc3KmuNeOPBwfj1iAipS1IEURTxz71nsXjbcVh+ywwI8sPEIcG4b0hIs8Gi4mo9th8twldHLmL/6RI0Xmz8viEhSJ4UhTB/Xye9gwani6sw53+/x/krNfDTeGLVr+NwR7/u7T6P0WTGzxf1OHCmBN+dLUVk9y549f5oBnQHYyCyMwYikovLlXUY9sZOCAJw5o1Jrd4mUJJP9udi0VfH4O/rhV1/uAvdfL2lLsmt1ZvMeG3zUaQezAMAzIgPx1Nj+uDWnl3bfI7iylps/bEQJy5VYvJtoRh5a/tDiL2UVhswb+0hHDpfBg+VgL9OHYxHR/Rq9TVms4iTxZU4cPoKDpy5gu/OXkFlndGmza+GhWPxQ0M4aNuB2vP726XWIcrMzMTkyZMRGhoKQRCwadOmm75mz549iIuLg0ajQZ8+ffDhhx82aZOeno7o6Gio1WpER0dj48aNDqieSHqWGWZd1Z4MQ43MSohA/6CuKKupx3sZJ6Uux63pa+vxm5TvkXowD4IAvHp/NN56eEi7whAA9PTT4IlRkVj80G2ShiEACOjijc+eGoEHbg+FySzipY0/4Y2vj+HkpUrsPXUZ6YcvYNXuM1j01VE8/fkPmP5hFoa9sRP3Lt2L1/9zDDt/uYTKOiO0Gk8kRQdh/p19oRKAdd/n4+/bT0j63ug6lxp1WV1djZiYGDz55JN4+OGHb9o+NzcXkyZNwlNPPYVPP/0U+/fvx+9+9zv06NHD+vqsrCzMmDEDf/3rX/Hggw9i48aNmD59Ovbt24cRI0Y4+i0ROZV1lWreKrPh6aHCX6YMwqP//A6ffnsevxrWC9Gh7O21t/zSGsz53+9x8lIVfLw88P7MWIyPDpK6LLtQe3pgyYzbEdm9K5bsPIl/7s3FP/fmtvoaHy8PDI8MwMi+gRjZtzuiQ7XWW2SR3X3xQvpPWLX7DPx9vTBvDNfKqjeZ4eUhXT+Ny94yEwQBGzduxAMPPNBimxdeeAFbtmzBL7/8Yj02f/58HDlyBFlZWQCAGTNmQK/XY9u2bdY29957L/z9/ZGamtqmWnjLjORi76nLmLXmIKKC/fDfBWOkLsflPP3ZD/j6p0IM7x2AtN8m8FaFHeXkl2Pu/36PkioDgrRqrJk9DINv0UldlkNszinAX//zC4xmM3r6qdHTT9PwVWv5qkZoNx8MDtXB27PlX/Af7jmDt7YdBwD8/ZHbMD0+3FlvQVJVdUaculSJU5eqcOJSJU5ee8SEdcPqx+Pt+r3a8/vbpXqI2isrKwtJSUk2xyZMmIA1a9agvr4eXl5eyMrKwvPPP9+kzdKlS1s8b11dHerq6qzP9Xq9XesmchT2ELXupfsG4v+OX8LBc6V44+tf8IcJA7i+TBsYjGZU1RmhEhqWM/AQBKgEASoV4CEIyDh2CQvSclBnNGNgiBYfPxGPEJ37znCcevstmHr7LZ0+z/w7+6Ks2oCPMs/ixfQfofPxwoRBwXao0P7yS2uwJOMkrlQb0EXtAV9vT3Tx9oCv+tpXb09ovDxQ32glcMvq4HXXnl/S1+FEUSUKyq82+z1OFXd+gc3OkHUgKioqQlCQbXdsUFAQjEYjSkpKEBIS0mKboqKiFs+7ePFiLFq0yCE1EzmSZad7pa9B1JJbuvlg4fj+eHPrcfxrXy6+OV6Mtx6+DcMjuT5RS/JLa/DAiv24Um24adu7o3ri/Zmx6Krm9ddWL06MQlmNAV8cuoDfp2bjf58cjsS+gU3aiaKIc1dqsO/UZeSV1qCbrzf8fb0R0MUbgV2vfe3iDa3Gy67jB0VRxPrDF7Doq2OoumFQeGf08FNjQJAf+gV1vfa14b+lJPur9sYu7+aWgW+uTWtd5cnJyVi4cKH1uV6vR3i4MroySd6ur1It+7/aDjNvTF/0CuiCVzf/jLMl1Zj+URZmJUTgT/cOUPzq3jcymsx4bl32TcOQSgBmj+yNlycNhKeEY0DkSBAEvPngEJTX1GPHsUt4au0hrJuXgMG36FBWbcD+MyXYd6oEe0+VtNiz0piHSsDIvoFYNGUQ+vToXMC4UlWHlzb+hO1HLwEA4iL8MWNYOGrrTaiuM6HGYLz+1WDCVYMJ6kYrgas9VVB7eViP6Xy90b9nV/QP8oN/F9eb6Snrn5rBwcFNenqKi4vh6emJwMDAVtvc2GvUmFqthlrNFUlJfqz7mPnwF3tr7h0cjMQ+gXhz6y9IO5SPf397Hjt/uYQ3HhyMu6PcYxCwPbz/zWn8kFcOP40ntj47GsE6DcyiCLMZMIsiTKII0Qx4egjowl6hDvP0UOH9mbF44pOD+PZsKR7/+CBu6eaDny9WoPEoXy8PAXER/hgUqoP+an3D6trVBpRee1TVGWEyi9h7qgT3LtuL34+9Fb+9s2+r45ha8s3xS/jT+p9QUlUHT5WA58f3x/w7+7r1ukmyvoITExPx1Vdf2RzbsWMH4uPj4eXlZW2TkZFhM45ox44dGDlypFNrJXIGPXuI2kzn64W3H7kNU24PRfKGn5BXWoPfpBzC1NtD8dr90YrfpuFgbimWf3MKAPDGg0MQHuDcBRGVRuPlgX8+Ho+Z//wWPxfoUXqtV25AkB/u6Ncdd/TrjhGRAfD1bvnvdp3RhLwrNfjb179gz8nLeDfjJLYcuYjFDw1p87Y11XVGvLH1F3z+XcMaUv16dsWSGbe77QD5xjr1U7O+vh5FRUWoqalBjx49EBDQufvwVVVVOH36tPV5bm4ucnJyEBAQgF69eiE5ORkFBQVYu3YtgIYZZcuXL8fChQvx1FNPISsrC2vWrLGZPfbcc89hzJgxePvttzF16lRs3rwZO3fuxL59+zpVK5Erur6PGXuI2mrUrd2xfcEYvJdxAmv25WJzzkXsP30F6+cnonf3LlKXJ4mKq/V4Pi0HZhF4eGgYpsSESl2SIvhpvLD2NyOQcuAcegf64o5bu6OnVtPm16s9PdAvyA8pTw7DliMX8df/HMOp4io88mEWHh3RCy/cGwVdM73H9SYzckuq8UuhHksyTuLclRoAwJw7IvFHBU08aPe0+6qqKnz22WdITU3FwYMHbWZjhYWFISkpCfPmzcOwYcPaXczu3bsxduzYJsdnz56NlJQUPPHEEzh37hx2795t/bM9e/bg+eefx9GjRxEaGooXXngB8+fPt3n9+vXr8corr+Ds2bPo27cv3njjDTz00ENtrovT7kkunlp7CBnHLuHNB4fcdCVdaupIfjn+58sjOF1chd6Bvtjwu1EIcMGxDo4kiiKeSc3G1z8WIiLQF18/O5qDpGWqvMaAxVuPI+1QPoCGgcwvTYqCr7cnThZV4mRxFU4WVeJsSRXqTdejQIhOg3emxWCUxAti2oPDtu5YsmQJ3njjDURGRmLy5MkYPnw4brnlFvj4+KC0tBQ///wz9u7di40bNyIhIQEffPAB+vXr1+k3JDUGIpKLGR9l4bvcUnwwMxaT+a/6DimurMWDKw6goPwqhvbqhs+fSpD9v5BrDEZsP1qEYb0DbroX2JeH8vHH9T/CUyVg/f8bidvDuzmnSHKY785eQfLGn3D2cnWLbbqqPdE/qCviIvzxzN39mu1JkiOHBaJp06bhtddeQ3R0NDZt2oSkpCT4+TXdoK+urg5r1qyBt7c35s6d2/534GIYiEguJi3bi2OFeqQ8OQx3DegpdTmydbq4Eg+tPAB9rRETBwdjxaNDZbsVyrGLevw+9QecuVwNLw8Bvx4RgafH3ooefk3HSJ0rqcak9/eixmDCHycMwNNjb5WgYnKEOqMJK3edwReH8hHY1Rv9g/zQP8gPA4L80D/YD6E6jVsuVOqUzV19fHxw9OhR9OnTp0NFygkDEcnF6L9/g/zSq9jwu5EY2stf6nJk7duzV/D4moMwmMyYe0ckXrk/WuqS2kUURazNOo83tv4Cg9EMHy8PXK03AWjYUuI3d/TGvNF9ofNt6AkwGM145MMD+PFCBRL6BOCzuQluPaOIlMEpm7sOHz4cubmt7+NCRM6lv2pZqZpjPjoroU8g/jHtNgDAv/blImW/fH7elVUbMO/fh/HnLUdhMJoxLqon9r94Nz6bOwIx4d1wtd6EFbvOYPTfv8GKXadRYzBiyc6T+PFCBXQ+Xlgy43aGIVKcDv/UfPbZZ/HSSy9h/fr1XLSQyAWIomhdSZazzOxj6u234ELZVfxj+wks+s8xhHbzQZJEWyvU1ptwvKgS4f4+rS4J8N3ZK1iQloPCilp4e6iQPCkKT4zsDUEQMOrW7tjUNxAZxy7hnR0ncPJSFf6x/QQ+3peL0pqGad5vPTTErbfdIGpJhwPRtGnTAACDBg3ClClTcNdddyE2NhZDhgyBt7eyZmUQuYIagwkmc8MdcO5lZj+/u6svLpRdRerBPDy7LhupTyUg1sm3I/NLazD7k4PWQbHBWg2iQ7WIDtEiOlSLQaFahHbzwYpdp/H+/52CWQQiu3fBBzNjm6wfIwgCkgYFY9zAIHx15CLeyziJvNKGadYzh4dj4pAQp743IlfR4TFE58+fR05ODo4cOWL9eu7cOXh4eCAqKgo//vijvWuVDMcQkRwUVlxF4uJv4KkScOqNiW45QFIqRpMZT609hF0nLiOwizc2/G4kIgKds0bRTxcq8GTK9yipqoPGS4XaenOz7bw8BOvU6YeHhuH1qYPatHq0wWhG+g8XkF9ag9/f3Q8+3vKeUUfUmFN2u4+IiEBERASmTp1qPVZZWYmcnBy3CkNEcmHd6d7Hi2HIzjw9VFj+6FBM/ygLRy/q8fCqA3h3+u24s3+Pdp3neJEedfVmxLRxKvvuE8X43Wc/oMZgQlSwH/73N8PRRe2J44V6HCvU42hBw9cTlyphMJrRxdsDf3twMB6MDWtzTd6eKswczjWriDrcQ6Qk7CEiOTh8vhQPr8pCRKAv9vyx6QKn1HnF+lo8/vFBHC+qBAD8dkwf/E/SgJvuFVVWbcDb/z2Odd83LJA3sm8gFtzTH8MjW17d/4vv85G88SeYzCLuuLU7Vj02tMWxYfUmM86VVKOnVuM268cQ2YPDZpnl5eW1q5CCgoJ2tSeijrPMMOM+Zo7TU6vBpqdH4fHECADAR5lnMe3DA8i7ttXBjcxmEesO5mHsu7utYchTJeDAmSuY/lEWfv2vb/H9uVKb14iiiKU7T+JP6T/CZBbxUOwt+PiJYa0OlPfyUKFfkB/DEFEntCsQDRs2DE899RQOHjzYYpuKigr885//xODBg7Fhw4ZOF0hEbWPZx4wDqh1L4+WB16cOxoePxUHn44UjFyow6f292Jxj+w/Anwsq8PCHB/Dihp9QXlOPqGA/rJ+fiN1/vAuPjugFLw8B+09fwbQPrwcjo8mM5A0/YenOhk1Vf3dXX7w7PaZDu5UTUfu065+Sv/zyC958803ce++98PLyQnx8PEJDQ6HRaFBWVoZjx47h6NGjiI+Pxz/+8Q9MnDjRUXUT0Q0qudO9U907OBhDwnRYsC4b358rw3PrcrD/dAn+J2kAVu46jX9/ex5msWFLhOfH98fsxAh4ejQEmzcfHILf3dUXK3efwZeH8rH/9BXsP52FEJ0GhRW1UAnAoqmDMSshQuJ3SaQcHRpDVFtbi61bt2Lv3r04d+4crl69iu7duyM2NhYTJkzA4MGDHVGrZDiGiORg5e7T+Pt/T+CRuDC8My1G6nIUw2gy4/1vTuODb05BFAFBACw/VSfHhOKV+wYiqJUdyy+U1ViDUb1JhMZLhQ9mDsX46CAnvQMi9+XwWWYajQb9+vXDlClT4OnJf40SuQLrLDPeMnMqTw8VFo7vj8Q+gViQlo1L+jr06dEFf506uE27hYf5+1p7jDb+UICxUT2brB1ERI7X4TQTExMDb29vREdHIyYmBrfffrv1a7du3exYIhG1hf5qwxgi3jKTRmLfQOxYcCd+yCvDqFu7t3vcT5i/L34/rp+DqiOim+nwSL19+/YhICAAkZGRqKurQ0pKCu6++24EBgZiwIABePXVV1FeXm7HUomoNRxDJD2drxfGRvXkIGgiGerw39pnnnkGK1euRHp6Oj7//HNkZ2dj165d6NOnD2bPno29e/ciNjYWly9ftme9RNSCSsssM069JiJqtw4HouPHjyM6Otrm2J133oklS5bghx9+wK5duxAfH4+XXnqp00US0c3pa7nTPRFRR3U4EA0bNgyffvppk+ODBg3Cjh07IAgC/vjHP2Lnzp2dKpCI2sbSQ8Sd7omI2q/DgWjlypVYunQpHn30URw/fhwAYDAYsGTJEgQENCxH36NHD1y6dMk+lRJRqzjLjIio4zrctz5o0CBkZWXhmWeeQXR0NNRqNYxGIzw9PfHJJ58AALKzsxEaGmq3YomoZZxlRkTUcZ36yTlo0CDs2rUL58+fx5EjR+Dh4YG4uDgEBwcDaOgheuutt+xSKBG1zGQWUW0wAWAgIiLqCLv85IyIiEBERNMl5kePHm2P0xPRTVRdu10GcAwREVFHcLEMIjdg2dhV46XiGjhERB3An5xEbkDPGWZERJ3CQETkBiq5BhERUacwEBG5geszzNhDRETUEQxERG6A+5gREXUOAxGRG+A+ZkREncNAROQGuI8ZEVHnMBARuQHuY0ZE1DkMRERugLPMiIg6h4GIyA1wHSIios5hICJyA5xlRkTUOQxERG7Asg6Rlj1EREQdwkBE5AbYQ0RE1DkuF4hWrlyJyMhIaDQaxMXFYe/evS22feKJJyAIQpPHoEGDrG1SUlKabVNbW+uMt0PkFNZp91yHiIioQ1wqEKWlpWHBggV4+eWXkZ2djdGjR2PixInIy8trtv2yZctQWFhofeTn5yMgIADTpk2zaafVam3aFRYWQqPROOMtETnF9UHV7CEiIuoIlwpE7733HubMmYO5c+di4MCBWLp0KcLDw7Fq1apm2+t0OgQHB1sfhw4dQllZGZ588kmbdoIg2LQLDg52xtshcoo6owkGoxkAZ5kREXWUywQig8GAw4cPIykpyeZ4UlISDhw40KZzrFmzBvfccw8iIiJsjldVVSEiIgJhYWG4//77kZ2d3ep56urqoNfrbR5ErsoyfkgQAD81e4iIiDrCZQJRSUkJTCYTgoKCbI4HBQWhqKjopq8vLCzEtm3bMHfuXJvjUVFRSElJwZYtW5CamgqNRoNRo0bh1KlTLZ5r8eLF0Ol01kd4eHjH3hSRE1hmmHX19oRKJUhcDRGRPLlMILIQBNsf6KIoNjnWnJSUFHTr1g0PPPCAzfGEhAQ89thjiImJwejRo/HFF1+gf//++OCDD1o8V3JyMioqKqyP/Pz8Dr0XImfgDDMios5zmZ+g3bt3h4eHR5PeoOLi4ia9RjcSRREff/wxZs2aBW9v71bbqlQqDBs2rNUeIrVaDbVa3fbiiSRUyRlmRESd5jI9RN7e3oiLi0NGRobN8YyMDIwcObLV1+7ZswenT5/GnDlzbvp9RFFETk4OQkJCOlUvkavgDDMios5zqZ+gCxcuxKxZsxAfH4/ExESsXr0aeXl5mD9/PoCGW1kFBQVYu3atzevWrFmDESNGYPDgwU3OuWjRIiQkJKBfv37Q6/V4//33kZOTgxUrVjjlPRE5Gne6JyLqPJcKRDNmzMCVK1fw+uuvo7CwEIMHD8bWrVuts8YKCwubrElUUVGB9PR0LFu2rNlzlpeXY968eSgqKoJOp0NsbCwyMzMxfPhwh78fImfgTvdERJ0niKIoSl2Eq9Pr9dDpdKioqIBWq5W6HCIb7+04gfe/OY1ZCRH46wNNe0mJiJSqPb+/XWYMERF1jJ6zzIiIOo2BiEjmOMuMiKjzGIiIZI6zzIiIOo+BiEjmOMuMiKjzGIiIZI6zzIiIOo+BiEjm9OwhIiLqNAYiIpljDxERUecxEBHJmCiKnGVGRGQHDEREMlZVZ4TJ3LC2qo6BiIiowxiIiGSsrLph/JDGSwWNl4fE1RARyRcDEZGMldUYAAABvt4SV0JEJG8MREQyVnotEHVjICIi6hQGIiIZK7f0EHVhICIi6gwGIiIZs4wh6ubLAdVERJ3BQEQkY2XsISIisgsGIiIZK+MYIiIiu2AgIpIxyy0zf94yIyLqFAYiIhnjLTMiIvtgICKSsdJq3jIjIrIHBiIiGSuvabhlxoUZiYg6h4GISKZEUWw0qJpjiIiIOoOBiEimrtabUGc0A+AYIiKizmIgIpKpsmu3y7w9VPD15sauRESdwUBEJFNl1ddvlwmCIHE1RETyxkBEJFOcck9EZD8MREQyVVrNAdVERPbCQEQkU9Yp9+whIiLqNAYiIpniooxERPbDQEQkU+WWMUQMREREncZARCRTlmn3HENERNR5DEREMmWZZebPHiIiok5jICKSKU67JyKyHwYiIpkqq+YtMyIie2EgIpIp9hAREdmPywWilStXIjIyEhqNBnFxcdi7d2+LbXfv3g1BEJo8jh8/btMuPT0d0dHRUKvViI6OxsaNGx39NogcqrbehBqDCQCn3RMR2YNLBaK0tDQsWLAAL7/8MrKzszF69GhMnDgReXl5rb7uxIkTKCwstD769etn/bOsrCzMmDEDs2bNwpEjRzBr1ixMnz4d3333naPfDpHDWBZl9FAJ0Go8Ja6GiEj+BFEURamLsBgxYgSGDh2KVatWWY8NHDgQDzzwABYvXtyk/e7duzF27FiUlZWhW7duzZ5zxowZ0Ov12LZtm/XYvffeC39/f6SmprapLr1eD51Oh4qKCmi12va9KSIH+KVQj4nL9qJ7V28cemW81OUQEbmk9vz+dpkeIoPBgMOHDyMpKcnmeFJSEg4cONDqa2NjYxESEoJx48Zh165dNn+WlZXV5JwTJky46TmJXFkZV6kmIrIrl+lrLykpgclkQlBQkM3xoKAgFBUVNfuakJAQrF69GnFxcairq8O///1vjBs3Drt378aYMWMAAEVFRe06JwDU1dWhrq7O+lyv13f0bRE5hGVRRq5STURkHy4TiCwEQbB5Lopik2MWAwYMwIABA6zPExMTkZ+fj3feeccaiNp7TgBYvHgxFi1a1JHyiZyitIY73RMR2ZPL3DLr3r07PDw8mvTcFBcXN+nhaU1CQgJOnTplfR4cHNzucyYnJ6OiosL6yM/Pb/P3J3KG8mpOuScisieXCUTe3t6Ii4tDRkaGzfGMjAyMHDmyzefJzs5GSEiI9XliYmKTc+7YsaPVc6rVami1WpsHkSu53kPEQEREZA8udcts4cKFmDVrFuLj45GYmIjVq1cjLy8P8+fPB9DQc1NQUIC1a9cCAJYuXYrevXtj0KBBMBgM+PTTT5Geno709HTrOZ977jmMGTMGb7/9NqZOnYrNmzdj586d2LdvnyTvkcgeLNPuA7rwlhkRkT24VCCaMWMGrly5gtdffx2FhYUYPHgwtm7dioiICABAYWGhzZpEBoMBf/jDH1BQUAAfHx8MGjQIX3/9NSZNmmRtM3LkSKxbtw6vvPIKXn31VfTt2xdpaWkYMWKE098fkb2UsYeIiMiuXGodIlfFdYjI1Uxdvg9HLlTgX4/H457oto+xIyJSElmuQ0REbWeZdu/PW2ZERHbBQEQkQ1yYkYjIvhiIiGSm3mRGZZ0RABdmJCKyFwYiIpmxDKgWBEDrw1tmRET2wEBEJDOWKffdfLzgoWp5xXUiImo7BiIimbGMH/Ln7TIiIrthICKSGcstM39u20FEZDcMREQyY51yz41diYjshoGISGZKOeWeiMjuGIiIZKa8hjvdExHZGwMRkcyUVl+bZcZbZkREdsNARCQz1h4i3jIjIrIbBiIimSnlTvdERHbHQEQkM5aFGTmGiIjIfhiIiGTGug4RxxAREdkNAxGRjJjMIiquWgZVs4eIiMheGIiIZKTiaj1EseG/OcuMiMh+GIiIZMSyKKOfxhNeHvzrS0RkL/yJSiQjXJSRiMgxGIiIZITbdhAROQYDEZGMWKfcc/wQEZFdMRARycj1KffsISIisicGIiIZ4SrVRESOwUBEJCPl1ZZVqnnLjIjInhiIiGSEPURERI7BQEQkI5x2T0TkGAxERDJyfdo9b5kREdkTAxGRjHCneyIix2AgIpIJs1lE+bWNXTntnojIvhiIiGSistYIk7lhZ1feMiMisi8GIiKZsCzK2MXbA2pPD4mrISJyLwxERDLBKfdERI7DQEQkE5xyT0TkOAxERDJRem2Vao4fIiKyPwYiIplgDxERkeO4XCBauXIlIiMjodFoEBcXh71797bYdsOGDRg/fjx69OgBrVaLxMREbN++3aZNSkoKBEFo8qitrXX0WyGyK8uijJxyT0Rkfy4ViNLS0rBgwQK8/PLLyM7OxujRozFx4kTk5eU12z4zMxPjx4/H1q1bcfjwYYwdOxaTJ09Gdna2TTutVovCwkKbh0ajccZbIrKbshquQURE5CieUhfQ2HvvvYc5c+Zg7ty5AIClS5di+/btWLVqFRYvXtyk/dKlS22ev/nmm9i8eTO++uorxMbGWo8LgoDg4GCH1k7kaJZbZv7c6Z6IyO5cpofIYDDg8OHDSEpKsjmelJSEAwcOtOkcZrMZlZWVCAgIsDleVVWFiIgIhIWF4f7772/Sg3Sjuro66PV6mweR1K7vY8YeIiIie3OZQFRSUgKTyYSgoCCb40FBQSgqKmrTOd59911UV1dj+vTp1mNRUVFISUnBli1bkJqaCo1Gg1GjRuHUqVMtnmfx4sXQ6XTWR3h4eMfeFJEdWfcxYyAiIrI7lwlEFoIg2DwXRbHJseakpqbiL3/5C9LS0tCzZ0/r8YSEBDz22GOIiYnB6NGj8cUXX6B///744IMPWjxXcnIyKioqrI/8/PyOvyEiO7m+MCNvmRER2ZvLjCHq3r07PDw8mvQGFRcXN+k1ulFaWhrmzJmDL7/8Evfcc0+rbVUqFYYNG9ZqD5FarYZarW578UQOJooip90TETmQy/QQeXt7Iy4uDhkZGTbHMzIyMHLkyBZfl5qaiieeeAKff/457rvvvpt+H1EUkZOTg5CQkE7XTOQsVXVG1JsaNnblLDMiIvtzmR4iAFi4cCFmzZqF+Ph4JCYmYvXq1cjLy8P8+fMBNNzKKigowNq1awE0hKHHH38cy5YtQ0JCgrV3ycfHBzqdDgCwaNEiJCQkoF+/ftDr9Xj//feRk5ODFStWSPMmiTrAMn5I46WCjzc3diUisjeXCkQzZszAlStX8Prrr6OwsBCDBw/G1q1bERERAQAoLCy0WZPoo48+gtFoxNNPP42nn37aenz27NlISUkBAJSXl2PevHkoKiqCTqdDbGwsMjMzMXz4cKe+N6LOsOx0z94hIiLHEERRFKUuwtXp9XrodDpUVFRAq9VKXQ4p0O4TxXjik+8xMESLbc+NlrocIiJZaM/vb5cZQ0RELbNOueeijEREDsFARCQDXJSRiMixGIiIZMA65Z6BiIjIIRiIiGSg1DqomrfMiIgcgYGISAasO91zUUYiIodgICKSgXJOuycicigGIiIZKK1u6CHiPmZERI7BQEQkA9zHjIjIsRiIiGTAMu2et8yIiByDgYjIxV01mFBnNAPgoGoiIkdhICJycZYp914eArpwY1ciIodgICJycWWNbpcJgiBxNURE7omBiMjFcad7IiLHYyAicnGWRRk55Z6IyHEYiIhcHKfcExE5HgMRkYvjTvdERI7HQETk4sqv3TIL6MJbZkREjsJAROTiLpRdBcBB1UREjsRAROTCKmrqkXnyMgAgoU+gxNUQEbkvBiIiF/bVjxdhMJkRFeyHQaFaqcshInJbDERELiz9hwsAgIeHhnFRRiIiB2IgInJRZy5XITuvHB4qAVNjQ6Uuh4jIrTEQEbmo9MMNvUN39u+Bnn4aiashInJvDERELshkFrExuwAA8EhcmMTVEBG5PwYiIhd04EwJCitqofPxwriBPaUuh4jI7TEQEbkgy+2yyTEhUHt6SFwNEZH7YyAicjGVtfX479EiAMAjceESV0NEpAwMREQuZutPhaitN6Nvjy6ICdNJXQ4RkSIwEBG5mPTDDYOpH47j2kNERM7CQETkQs5fqcbBc6VQCcBDsZxdRkTkLAxERC4k/YeG3qFRt3ZHsI5rDxEROQsDEZGLMJtFbLi2VQfXHiIici4GIiIX8V1uKS6UXYWf2hNJ0cFSl0NEpCgMREQuwrKR6323hcDHm2sPERE5EwMRkQuorjNi60+FAHi7jIhICi4XiFauXInIyEhoNBrExcVh7969rbbfs2cP4uLioNFo0KdPH3z44YdN2qSnpyM6OhpqtRrR0dHYuHGjo8on6pD//lyEGoMJvQN9ERfhL3U5RESK41KBKC0tDQsWLMDLL7+M7OxsjB49GhMnTkReXl6z7XNzczFp0iSMHj0a2dnZeOmll/Dss88iPT3d2iYrKwszZszArFmzcOTIEcyaNQvTp0/Hd99956y3RXRTlttlDw3l2kNERFIQRFEUpS7CYsSIERg6dChWrVplPTZw4EA88MADWLx4cZP2L7zwArZs2YJffvnFemz+/Pk4cuQIsrKyAAAzZsyAXq/Htm3brG3uvfde+Pv7IzU1tU116fV66HQ6VFRUQKvVdvTtNVFnNOFyZZ3dzkeurcbQ8P/7cmUdiitrG/13HQ6cuQIA2PfCWIT5+0pcKRGRe2jP729PJ9V0UwaDAYcPH8aLL75oczwpKQkHDhxo9jVZWVlISkqyOTZhwgSsWbMG9fX18PLyQlZWFp5//vkmbZYuXdpiLXV1dairux5U9Hp9O99N2xy9qMdDK5t/b6Q8Ywf0YBgiIpKIywSikpISmEwmBAUF2RwPCgpCUVFRs68pKipqtr3RaERJSQlCQkJabNPSOQFg8eLFWLRoUQffSdsJANSeLnXXkhxI7alCT60GPbqq0VOrvv7VT42efhoM7cWxQ0REUnGZQGRx4/gJURRbHVPRXPsbj7f3nMnJyVi4cKH1uV6vR3i4/Xcdj+3ljxN/m2j38xIREVH7uEwg6t69Ozw8PJr03BQXFzfp4bEIDg5utr2npycCAwNbbdPSOQFArVZDrVZ35G0QERGRDLnM/Rpvb2/ExcUhIyPD5nhGRgZGjhzZ7GsSExObtN+xYwfi4+Ph5eXVapuWzklERETK4zI9RACwcOFCzJo1C/Hx8UhMTMTq1auRl5eH+fPnA2i4lVVQUIC1a9cCaJhRtnz5cixcuBBPPfUUsrKysGbNGpvZY8899xzGjBmDt99+G1OnTsXmzZuxc+dO7Nu3T5L3SERERK7HpQLRjBkzcOXKFbz++usoLCzE4MGDsXXrVkRERAAACgsLbdYkioyMxNatW/H8889jxYoVCA0Nxfvvv4+HH37Y2mbkyJFYt24dXnnlFbz66qvo27cv0tLSMGLECKe/PyIiInJNLrUOkaty1DpERERE5Djt+f3tMmOIiIiIiKTCQERERESKx0BEREREisdARERERIrHQERERESKx0BEREREisdARERERIrHQERERESKx0BEREREiudSW3e4Ksti3nq9XuJKiIiIqK0sv7fbsikHA1EbVFZWAgDCw8MlroSIiIjaq7KyEjqdrtU23MusDcxmMy5evAg/Pz8IgtCpc+n1eoSHhyM/P5/7orUBP6/242fWPvy82o+fWfvw82o/e31moiiisrISoaGhUKlaHyXEHqI2UKlUCAsLs+s5tVot/2K0Az+v9uNn1j78vNqPn1n78PNqP3t8ZjfrGbLgoGoiIiJSPAYiIiIiUjwGIidTq9X485//DLVaLXUpssDPq/34mbUPP6/242fWPvy82k+Kz4yDqomIiEjx2ENEREREisdARERERIrHQERERESKx0BEREREisdA5AArV65EZGQkNBoN4uLisHfv3lbb79mzB3FxcdBoNOjTpw8+/PBDJ1XqGtrzee3evRuCIDR5HD9+3IkVSyczMxOTJ09GaGgoBEHApk2bbvoapV9f7f3MlH6NLV68GMOGDYOfnx969uyJBx54ACdOnLjp65R6nXXk81L6NbZq1Srcdttt1kUXExMTsW3btlZf44zri4HIztLS0rBgwQK8/PLLyM7OxujRozFx4kTk5eU12z43NxeTJk3C6NGjkZ2djZdeegnPPvss0tPTnVy5NNr7eVmcOHEChYWF1ke/fv2cVLG0qqurERMTg+XLl7epvdKvL6D9n5mFUq+xPXv24Omnn8a3336LjIwMGI1GJCUlobq6usXXKPk668jnZaHUaywsLAxvvfUWDh06hEOHDuHuu+/G1KlTcfTo0WbbO+36Esmuhg8fLs6fP9/mWFRUlPjiiy822/5Pf/qTGBUVZXPst7/9rZiQkOCwGl1Jez+vXbt2iQDEsrIyJ1Tn2gCIGzdubLWN0q+vG7XlM+M1Zqu4uFgEIO7Zs6fFNrzOrmvL58VrrCl/f3/xX//6V7N/5qzriz1EdmQwGHD48GEkJSXZHE9KSsKBAweafU1WVlaT9hMmTMChQ4dQX1/vsFpdQUc+L4vY2FiEhIRg3Lhx2LVrlyPLlDUlX1+dxWusQUVFBQAgICCgxTa8zq5ry+dlwWsMMJlMWLduHaqrq5GYmNhsG2ddXwxEdlRSUgKTyYSgoCCb40FBQSgqKmr2NUVFRc22NxqNKCkpcVitrqAjn1dISAhWr16N9PR0bNiwAQMGDMC4ceOQmZnpjJJlR8nXV0fxGrtOFEUsXLgQd9xxBwYPHtxiO15nDdr6efEaA3766Sd07doVarUa8+fPx8aNGxEdHd1sW2ddX9zt3gEEQbB5Lopik2M3a9/ccXfVns9rwIABGDBggPV5YmIi8vPz8c4772DMmDEOrVOulH59tRevseueeeYZ/Pjjj9i3b99N2/I6a/vnxWus4TPIyclBeXk50tPTMXv2bOzZs6fFUOSM64s9RHbUvXt3eHh4NOndKC4ubpJuLYKDg5tt7+npicDAQIfV6go68nk1JyEhAadOnbJ3eW5BydeXPSnxGvv973+PLVu2YNeuXQgLC2u1La+z9n1ezVHaNebt7Y1bb70V8fHxWLx4MWJiYrBs2bJm2zrr+mIgsiNvb2/ExcUhIyPD5nhGRgZGjhzZ7GsSExObtN+xYwfi4+Ph5eXlsFpdQUc+r+ZkZ2cjJCTE3uW5BSVfX/akpGtMFEU888wz2LBhA7755htERkbe9DVKvs468nk1R0nXWHNEUURdXV2zf+a068uuQ7RJXLdunejl5SWuWbNGPHbsmLhgwQKxS5cu4rlz50RRFMUXX3xRnDVrlrX92bNnRV9fX/H5558Xjx07Jq5Zs0b08vIS169fL9VbcKr2fl5LliwRN27cKJ48eVL8+eefxRdffFEEIKanp0v1FpyqsrJSzM7OFrOzs0UA4nvvvSdmZ2eL58+fF0WR11dz2vuZKf0a+3//7/+JOp1O3L17t1hYWGh91NTUWNvwOruuI5+X0q+x5ORkMTMzU8zNzRV//PFH8aWXXhJVKpW4Y8cOURSlu74YiBxgxYoVYkREhOjt7S0OHTrUZvrl7NmzxTvvvNOm/e7du8XY2FjR29tb7N27t7hq1SonVyyt9nxeb7/9tti3b19Ro9GI/v7+4h133CF+/fXXElQtDct03Rsfs2fPFkWR11dz2vuZKf0aa+6zAiB+8skn1ja8zq7ryOel9GvsN7/5jfVnfo8ePcRx48ZZw5AoSnd9CaJ4bWQSERERkUJxDBEREREpHgMRERERKR4DERERESkeAxEREREpHgMRERERKR4DERERESkeAxEREREpHgMRERERKR4DERERESkeAxERUSNGo1HqEohIAgxERKRY586dgyAIWL9+PcaMGQO1Wo2NGzdKXRYRScBT6gKIiKSSk5MDAHj77bfx5ptvIjIyEj169JC2KCKSBAMRESnWkSNH0KVLF3z55Zfo3bu31OUQkYR4y4yIFCsnJwdTpkxhGCIiBiIiUq4jR47grrvukroMInIBDEREpEh6vR7nzp1DbGys1KUQkQtgICIiRTpy5AhUKhWGDBkidSlE5AIYiIhIkY4cOYKoqChoNBqpSyEiFyCIoihKXQQRERGRlNhDRERERIrHQERERESKx0BEREREisdARERERIrHQERERESKx0BEREREisdARERERIrHQERERESKx0BEREREisdARERERIrHQERERESKx0BEREREivf/AeX73L/pTVMBAAAAAElFTkSuQmCC",
      "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.12.4"
  },
  "vscode": {
   "interpreter": {
    "hash": "6834e2eb4702651dab439b72249722ae87da91a96342472ed2409ca0cf3240b1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}