{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "35354229",
   "metadata": {},
   "source": [
    "# Running a SepTop RBFE calculation\n",
    "\n",
    "This tutorial gives a step-by-step process to set up a relative binding free energy (RBFE) simulation campaign using the [Separated Topologies Protocol](https://pubs.acs.org/doi/10.1021/acs.jctc.3c00282) in OpenFE.\n",
    "In this tutorial we are performing a relative binding free energy calculation of two ligands binding to TYK2.\n",
    "\n",
    "The relative binding free energy is obtained through a thermodynamic cycle. The ligands are transformed into each other both in solvent, giving \n",
    "$\\Delta\\Delta G$ (solvent) and in the complex, giving $\\Delta\\Delta G$ (complex), which allows calculation of the relative binding free energy, $\\Delta\\Delta G$. \n",
    "Each ligand is represented with its own set of coordinates, meaning that the interactions of all atoms of one ligand are turned off while simultaneously turning on the interactions of all atoms of the other ligand. Therefore, restraints are required: \n",
    "Ligands are restrained to the protein in the complex states using orientational (Boresch-style) restraints; in the solvent states ligands are restrained to remain apart from each other using a single harmonic distance restraint between the ligands. Restraints are not depicted in the thermodynamic cycle below for simplicity.\n",
    "\n",
    "**Note:** In this `Protocol`, the coulombic interactions of the molecule are fully turned off (annihilated), while the Lennard-Jones interactions are decoupled, meaning the intermolecular interactions are turned off, while keeping the intramolecular Lennard-Jones interactions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23881797-0bd3-408b-9212-ac206d46cd39",
   "metadata": {},
   "source": [
    "<img src=\"septop_cycle.png\" width=\"600\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4df1661a",
   "metadata": {},
   "source": [
    "## 0. Setup for Google Colab\n",
    "\n",
    "If you are running this example in Google Colab, run the following cells to setup the environment. If you are running this notebook locally, skip down to `1. Loading the ligands`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e178d4d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NBVAL_SKIP\n",
    "# Only run this cell if on google colab\n",
    "import os\n",
    "if \"COLAB_RELEASE_TAG\" in os.environ:\n",
    "    # fix for colab's torchvision causing issues\n",
    "    !rm -r /usr/local/lib/python3.12/dist-packages/torchvision\n",
    "    \n",
    "    !pip install -q condacolab\n",
    "    import condacolab\n",
    "    condacolab.install_from_url(\"https://github.com/OpenFreeEnergy/openfe/releases/download/v1.7.0/OpenFEforge-1.7.0-Linux-x86_64.sh\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "111a9e15",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NBVAL_SKIP\n",
    "# Only run this cell if on google colab\n",
    "import os\n",
    "if \"COLAB_RELEASE_TAG\" in os.environ:\n",
    "    import condacolab\n",
    "    import locale\n",
    "    locale.getpreferredencoding = lambda: \"UTF-8\"\n",
    "    !mkdir inputs && cd inputs && openfe fetch rbfe-tutorial\n",
    "    # quick fix for https://github.com/conda-incubator/condacolab/issues/75\n",
    "    # if deprecation warnings persist, rerun this cell\n",
    "    import warnings\n",
    "    warnings.filterwarnings(action=\"ignore\", message=r\"datetime.datetime.utcnow\") \n",
    "    for _ in range(3):\n",
    "      # Sometimes we have to re-run the check\n",
    "      try:\n",
    "        condacolab.check()\n",
    "      except:\n",
    "        pass\n",
    "      else:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fea29c3",
   "metadata": {},
   "source": [
    "## 1. Loading the ligands\n",
    "\n",
    "First we must load the chemical models between which we wish to calculate free energies.\n",
    "In this example these are initially stored in a molfile (`.sdf`) containing multiple molecules.\n",
    "This can be loaded using the `SDMolSupplier` class from rdkit and passed to openfe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fc97de03",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import openfe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "41cf8be7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import Chem\n",
    "supp = Chem.SDMolSupplier(\"tyk2_ligands.sdf\", removeHs=False)\n",
    "ligands = [openfe.SmallMoleculeComponent.from_rdkit(mol) for mol in supp]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e5de19a",
   "metadata": {},
   "source": [
    "## 2. Charging the ligands\n",
    "\n",
    "It is recommended to use a single set of charges for each ligand to ensure reproducibility between repeats or consistent charges between different legs of a calculation involving the same ligand, like a relative binding affinity calculation for example. \n",
    "\n",
    "Here we will use some utility functions from OpenFE which can assign partial charges to a series of molecules with a variety of methods which can be configured via the `OpenFFPartialChargeSettings` class. In this example \n",
    "we will charge the ligands using the `am1bcc` method from `ambertools` which is the default charge scheme used by OpenFE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5219106c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Generating charges: 100%|███████████████████████| 10/10 [02:46<00:00, 16.66s/it]\n"
     ]
    }
   ],
   "source": [
    "from openfe.protocols.openmm_utils.omm_settings import OpenFFPartialChargeSettings\n",
    "from openfe.protocols.openmm_utils.charge_generation import bulk_assign_partial_charges\n",
    "\n",
    "charge_settings = OpenFFPartialChargeSettings(partial_charge_method=\"am1bcc\", off_toolkit_backend=\"ambertools\")\n",
    "\n",
    "charged_ligands = bulk_assign_partial_charges(\n",
    "    molecules=ligands,\n",
    "    overwrite=False,  \n",
    "    method=charge_settings.partial_charge_method,\n",
    "    toolkit_backend=charge_settings.off_toolkit_backend,\n",
    "    generate_n_conformers=charge_settings.number_of_conformers,\n",
    "    nagl_model=charge_settings.nagl_model,\n",
    "    processors=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6963be83",
   "metadata": {},
   "source": [
    "## 3. Creating the `LigandNetwork`\n",
    "\n",
    "The first step is to create a `LigandNetwork`. Here, we will be using the same process as in the relative hybrid topology `Protocol`, including the use of a mapper which is required by the scorer. **The mappings will not be used in the `Protocol`.** This is a temporary solution until we have developed a scorer specifically for the `SepTopProtocol`. Alternatively, the user can also manually define the edges they want to run when creating the transformations below, without creating a `LigandNetwork` first.\n",
    "\n",
    "The pipeline for creating a `LigandNetwork` can involve three components:\n",
    "\n",
    "1. **Atom Mapper**: Proposes potential atom mappings (descriptions of the alchemical change) for pairs of ligands. We will use the `LomapAtomMapper`. *The atom mapping will only be used to score the potential edges, the atom mapping is not used outside of the scorer*.\n",
    "2. **Scorer**: Given an atom mapping, provides an estimate of the quality of that mapping (higher scores are better). We will use `default_lomap_scorer`.\n",
    "3. **Network Planner**: Creates the actual `LigandNetwork`; different network planners provide different strategies. We will create a minimal spanning network with the `generate_minimal_spanning_network` method.\n",
    "\n",
    "Each of these components could be replaced by other options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5a3cf244",
   "metadata": {},
   "outputs": [],
   "source": [
    "mapper = openfe.LomapAtomMapper(max3d=1.0, element_change=False)\n",
    "scorer = openfe.lomap_scorers.default_lomap_score\n",
    "network_planner = openfe.ligand_network_planning.generate_minimal_spanning_network"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acc13581",
   "metadata": {},
   "source": [
    "The exact call signature depends on the network planner: a minimal spanning network requires a score, whereas that is optional for a radial network (but a radial network needs the central ligand to be provided)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f6e7bce5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ligand_network = network_planner(\n",
    "    ligands=charged_ligands,\n",
    "    mappers=[mapper],\n",
    "    scorer=scorer\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7492637",
   "metadata": {},
   "source": [
    "Now we can look at the overall structure of the `LigandNetwork`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e6ca6131",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAIHCAYAAAAPeO1DAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZbpJREFUeJzt3QmczdX/x/G3fd/3nVB2Zc1WEhGFJEKFSPIjS7Sg0kKbENqQRItkTUpC9pIUKqGk7CpkSWSZ/+Nzvv9ry3Bnu+vr+Xjcx9y5M3Pvme9Y3nPO+XxOspiYmBgBAAAgaiQP9gAAAAAQWARAAACAKEMABAAAiDIEQAAAgChDAAQAAIgyBEAAAIAoQwAEAACIMgRAAACAKJPSn086deqUdu7cqUyZMilZsmRJPyoAAADEiZ3tcejQIeXPn1/JkydPeAC08FeoUKG4jQIAAAABt23bNhUsWDDhAdBm/nxPmDlz5sQZHQAAABLNwYMH3YSdL7clOAD6ln0t/BEAAQAAQpc/2/UoAgEAAIgyBEAAAIAoQwAEAACIMgRAAACAKEMABAAAiDIEQAAAgChDAAQAAIgyBEAAAIAoQwAEAACIMgRAAPFWt67Uq5d3v2hRacSIxLuY1sh+5szEez4AwBkEQACJYtUqqUuXxLuYu3ZJN96ooNi3T+rRQ7riCil9eqlwYen++6UDB879vE2bpGbNpJw57ahMqVYt6fPPgzNmAIgLAiCARJErlxeWEkvevFKaNAqKnTu929Ch0nffSRMmSHPnSp06nft5TZpIJ05ICxdKq1dLV14p3XSTtHt3cMYNAP4iAAJIFOcvAW/YINWuLaVNK5UpI82fH7dl3bM/99dfvfenTJHq1JHSpZOqVvVm4GzmsUoVKWNGqVEj6Y8/zn2e8eOlsmW9MJkvn9S9+6Vfu1w5ado06eabpeLFpXr1pMGDpdmzvcBn/vxT+vln6eGHpQoVpJIlpWeflY4ckX74we/LBgBBQQAEkOhOnZKaN/dmBFeulMaMkQYMSPjzPv64NHCg9M03UsqUUps20oMPSi+9JC1dKm3eLD322JnPf/VV6X//85ambSbvww+lEiXi99q2/GvLvPa6JkcOqXRpaeJE6e+/vWD4+utSnjxS5coJ/14BICn9/z9lAJB45s3zwtiiRd5SrrEZtAYNEva8fftKDRt693v29ALgggXe3jtjS7S2XOvz9NPSAw94n+tjM4dxtXev9NRT0r33nnnMZiQ/+8zbA5gpk5Q8uRf+bKk4a9Z4f4sAEBDMAAJIdBs3SoUKnQl/plq1hD+vLbX6WNgy5cuf+9jvv3v37a3t47v++oS95sGD3l4/W8a2GUifmBipWzcpd25v9vGrr7wwaHsArYAFAEIZARBAorNwZDNkiS1VqjP3fc9//mO2/Gxsn2BCHTrk7Su0/YUzZpz7Wlb48dFH0uTJ3gxkpUrSK694r/vWWwl/bQBISgRAAImuVClp61Zpz54zj1mxRiDZsqwVptgScXxn/m64QUqd2ts7aMUsZ7NiD2NLv2ez930hFABCFXsAASQ62+tn1bPt20vPP+/NpPmKQJJiZjA2gwZJXbt6y7TWU9DGsXy51+PvYuzzLPxZyHv7bS8M2s3X7iZFCqlGDSlbNu97tMITm/kbO1bassVbMgaAUMYMIIBEZwHJWrgcPuwVXXTu7FXvmvNn0pKShTNrTWNLs9YKxvbn/fTTpb/OevpZ9bJVDlvVsLWP8d22bfM+x5o/W8GHfY/WJsZa0SxbJs2aJVWsmOTfGgAkSLKYGNutc3EHDx5UlixZdODAAWW2PggAEEc282Z9Aa13ns0OXsyxY15QtCrb+vW51ADgj7jkNZaAASQJK5qw4glrkGyhz1qxWLHEpcKfLbVOn+7tpbO9hACAxEcABJAkbB+dNWm2JVNbLrWZvBdf9D42ZIh3uxBbMl6/XnruOalgwaQZ2zvvnNvT72xFinCSB4DIxxIwgIDbt8+7XYgVUxQokPTh9OwK5bNZqxcLgQAQblgCBhDSsmf3bsFiLWLsBgDRiipgAACAKEMABAAAiDIEQAAAgChDAARwUSdOnOAKAUCEIQACiJX1iW/cuLEqVaqkCRMmcKUAIELQBxBArObPn6/PPvtMGTJkULNmzbhSABAhmAEEEKuhQ4e6t506dVK2bNm4UgAQIQiAAC7ou+++07x585Q8eXL16tWLqwQAEYQACOCChg0b5t62aNFCxYoV4yoBQAQhAAL4j127dukdOzBXUt++fblCABBhCIAA/mP06NE6fvy4atWqperVq3OFACDCEAABnOPvv//Wq6++6u4/8MADXB0AiEAEQADnePPNN7V//34VL15cTZs25eoAQAQiAAI47eTJkxo+fLi737t3b6VIkYKrAwARiAAI4LRZs2bpl19+cT3/OnTowJUBgAhFAAxRdetKvtZrRYtKI0Yk3nMnSybNnJl4z4fI8eKLL7q33bp1c6d/AAAiEwEwDKxaJXXpknjPt2uXdOONCop9+6QePaQrrpDSp5cKF5buv186cODM5yxa5IXUC93sWiBpfPnll1qxYoVSp06t7t27c5kBIIJxFnAYyJUrcZ8vb14Fzc6d3s1OGCtTRvrtN6lrV++xqVO9z6lZ0wupZ3v0UTuXVqpSJSjDjqrZv3bt2ilvMP+QAACSHDOAYeD8JeANG6TataW0ab0QZcEoLsu6Z3/ur79670+ZItWpI6VLJ1WtKm3a5M22WeDKmFFq1Ej6449zn2f8eKlsWSlNGilfPsmfSaNy5aRp06Sbb5aKF5fq1ZMGD5Zmz5ZOnPA+J3VqL6T6bjlySB9+KN19tzdWJD7b9zd9+nR3v0+fPlxiAIhwzACGmVOnpObNvaXTlSulQ4esV1vCn/fxx72Qac9rQatNGylzZumll7yl2latpMcek/6/PZx7aznh2We95WRbwl2+PH6vbV9rr5Uylj+NFv7+/FOiJiHpjBgxQqdOnVLDhg1VzlI6ACCiEQDDzLx50ubN3j453yqdzaA1aJCw57XTvho29O737OkFwAULpFq1vMc6dZImTDjz+U8/7QVP+1wfmzmMq717paeeku69N/bPeeMNb2yFCsX9+XFp1vNvvE3n0vgZAKIGATDMbNzoBaGzt2hVq5bw561Q4cz9PHm8t+XLn/vY77979+2t7dm7/vqEvebBg1KTJt4yts1AXsj27dKnn3pL1Egar7/+ujv9o0KFCqpfvz6XGQCiAHsAw0xMTNLsg0uV6sx93/Of/5gtPxvbJ5hQtnRt+wptf+GMGee+1tnefNPbA8iBFEnj33//1ahRo04f+5aMTZYAEBUIgGGmVClp61Zpz54zjwW6NUqmTF5hii0Rx3fm74YbvGIP299nxSyxhV0LgHfdFXtARMJMnjxZO3fuVP78+XX77bdzOQEgSrAEHGZsr59Vz7ZvLz3/vDeTNmCA97FATt4MGuS1b8md2ysCsXFYEYj1+LsY+zwLf0eOSG+/7YVBu/na3Zx98tjChdKWLd7+QyS+mJiY061fevTo4fr/AQCiAwEwzFhAshYunTt7RReXXSa98ILXViW2mbSkYAH06FHJjo21ApKcOaWWLS/9datXe9XLpkSJcz9mYc9mFs8u/rCegKVLJ/Lg4cyfP1/r1q1zJ37ce7EqHABAxEkWY9MAl3Dw4EFlyZJFBw4cUGbr14GQYjNv1hfw55+92cGLOXbMC4qffSax3z+6NWrUSJ9++qmb/Rs5cmSwhwMASKC45DVmAMOQFU1Y8UTJkl7os1Ys1q7lUuHPllqt12/y5N5eQkSv77//3oW/5MmTq5fv0GkAQNQgAIYh20f34IPStm3e0qvN5P3/Vi4NGeLdLsSWjNevl557TipYMGnG9s47sff0K1JE+uGHpHldxM2wYcPc2xYtWugy20cAAIgqLAFHmH37vNuFWPuWAgWSPpyeXaF8NqvktRCI4Nq9e7eKFCniWsB88cUXuvrqq/mRAEAEYAk4imXP7t2CxVrE2A2ha/To0S781axZk/AHAFGKPoBAFLETP1555ZXTjZ8BANGJAAhEkQkTJrizf4sXL65mzZoFezgAgCAhAAJR4uTJkxpujRslV/mb4uyu2wCAqEIABKLEhx9+qM2bNytbtmzq2LFjsIcDAAgiAmAYmTp1qipVqqTKlStryZIlwR4Owozv2Lf77rvPnf4BAIheBMAw0rJlS9WoUUPffPONmjRp4pr5Av748ssvtXz5cnfeb/fu3bloABDlCIBhZsSIEbruuut0+PBhNW3aVH/++Wewh4Qwmv1r27at8uXLF+zhAACCjAAYZlKlSqUPPvjAVXFu2bJFt912m44fPx7sYSGE2Z+T6XYGoKQ+ffoEezgAgBBAAAxDOXLkcBv6M2XKpEWLFun+++8P9pAQ4rPGp06d0g033KDy5csHezgAgBBAAAxTZcqU0XvvvadkyZLptddeO93cFzib9fx744033H0aPwMAfAiAYcwKQZ599ll332YBFy5cGOwhIcSMGTPGnf5hM38NGjQI9nAAACGCABjm+vXrpzvvvNM1+bX9gNbnDTB23u/IkSNPz/7ZbDEAAIYAGObsP3Wb5alevbr27dunm2++WQcPHgz2sBACJk+erJ07d7qq3zZt2gR7OACAEEIAjABp06bVjBkzVKBAAf3444+u1YfNCCJ6xcTEnG790qNHD9f/DwAAHwJghLBZnpkzZ7owOGfOHPXv3z/YQ0IQLViwQOvWrVP69Ol177338rMAAJyDABhBqlSpogkTJrj7zz//vCZNmhTsISFIfLN/d999t7Jnz87PAQBwDgJghGndurUGDBjg7t9zzz1auXJlsIeEAPvhhx80d+5cJU+eXL169eL6AwD+gwAYgZ588kk1a9ZMx44dU/PmzbV9+/ZgDwkBNGzYMPf2lltucSfGAABwPgJgBLKZH1v+LVeunHbv3u1C4JEjR4I9LASA/bzffvttd5/GzwCA2BAAI5QdE2fHxeXMmVOrV69Wp06dXGUoItvo0aNd/78aNWq4GwAAF0IAjGDFihXT1KlTlTJlStcTbsiQIcEeEpKQnfjx6quvuvvM/gEALoYAGOGuvfZavfzyy+7+wIEDXasYRKa33nrLNQO/7LLL3LI/AACxIQBGgS5duqh79+7u/h133OH6wyGyWOPv4cOHu/tW+ZsiRYpgDwkAEMIIgFHCwsH111/vlgmbNm2qP/74I9hDQiKaPXu2fv75Z2XLlk0dO3bk2gIALooAGCVsH+CUKVNUokQJ/fbbb7r11ltdsQAiw9ChQ93brl27KmPGjMEeDgAgxBEAo4idCGGVwZkzZ9bSpUvdsjCVweHPmn0vX75cqVKlOr3UDwDAxRAAo0zp0qX13nvvKVmyZBo7dqxrG4LIOPatbdu2yp8/f7CHAwAIAwTAKNS4cWN3VrDp3bu35s+fH+whIZ62bNmiadOmuft9+vThOgIA/EIAjFLWJ+6uu+5y1aO33Xabfvrpp2APCfHw0ksv6dSpU2rQoIEqVKjANQQA+IUAGKVsCfj111/X1Vdfrb/++stVBh84cCDYw0Ic7N+/X+PGjXP3+/bty7UDAPiNABjF0qZNqxkzZqhgwYLasGGDbr/9djcjiPAwZswY19anfPnybgYQAAB/EQCjXN68eTVr1iylS5dOc+fO1UMPPRTsIcEP1sJn5MiRp/f+2YwuAAD+IgBClSpV0oQJE05XlNqRYght77//vnbu3OkCfJs2bYI9HABAmCEAwmnVqpUeffTR00fHrVixgisToqx3o6/1S48ePZQmTZpgDwkAEGYIgDht0KBBuuWWW9zyYosWLbRt2zauTghauHCh1q5dq/Tp07uTPwAAiCsCIM78YUieXBMnTnTtRPbs2aNmzZq5IgOEFt/s39133+1OdwEAIK4IgDiHnSNrx8XlypVL3377rTp27MhxcSHkhx9+0CeffOKKPnr16hXs4QAAwhQBEP9RpEgRTZ8+3Z0t+8EHH+ipp57iKoWIYcOGube2VF+8ePFgDwcAEKYIgLig2rVr69VXX3X3H3/88dPHjSF4du/erbfffvv0SS4AAMQXARCx6tSpk3r27Onu27Fxa9as4WoF0csvv+wKdOz0lpo1a/KzAADEGwEQFzV06FB3ysSRI0dcUcjvv//OFQsCu/6+GVlm/wAACUUAxEWlTJnSNR2+/PLLtXXrVtce5tixY1y1ALPm3Hv37lWxYsXc/j8AABKCAIhLypYtm6sMzpIli5YvX65u3bpRGRxAdj7z8OHD3f3evXsrRYoUgXx5AEAEIgDCL1dccYWbCbRegePHj9dLL73ElQuQ2bNn66efflLWrFldWx4AABKKAAi/NWzY0O0J9O1D+/TTT7l6AWz8bKd+WJ9GAAASigCIOLHmwzYLderUKbVu3VobN27kCiahr776SsuWLXM9Ge3cXwAAEgMBEHFiJ1BYNaq1ITlw4ICaNm2q/fv3cxWTePavTZs2yp8/P9cZAJAoCICIszRp0riTQgoVKqRNmzbp9ttv14kTJ7iSiezXX3/V1KlT3f0+ffpwfQEAiYYAiHjJkyePqwxOnz695s2bp379+nElE5kV2thSu/VhrFixItcXAJBoCICItyuvvFITJ05090eMGOGqg5E4/vrrL40bN87dp/EzACCxEQCRILfeeqsGDRp0ukrV+gQi4caMGaPDhw+rXLlyuuGGG7ikAIBERQBEgj366KNq2bKljh8/7k4KsRNDEH923u/IkSNP7/2zwhsAABITARAJ/0OUPLkmTJjgloTtrGCrDP7777+5svE0ZcoU7dixQ3nz5lXbtm25jgCAREcARKLIkCGDZs2apdy5c2vt2rVq3769K2BA3MTExJxu/dK9e3dXcQ0AQGIjACLRFC5cWDNmzFDq1Kk1bdo0Pfnkk1zdOPr888+1Zs0aV11teyoBAEgKBEAkKmsQ/frrr7v7TzzxhD744AOucBz4jtqz01Zy5MjBtQMAJAkCIBJdhw4dTjcutqXgb7/9lqvsh/Xr1+uTTz5xRR925B4AAEmFAIgk8fzzz6tRo0b6559/1KxZM+3Zs4crfQnDhg1zb5s3b64SJUpwvQAASYYAiCSRIkUKvffee7riiiu0bds23XLLLTp27BhXOxYWkCdNmuTu0/gZAJDUCIBIMlmzZtXs2bPd2y+++MIVNViVK/7r5Zdfdv3/qlev7vZRAgCQlAiASFIlS5Z0fe1sRtB6BQ4fPpwrfp4jR47olVdeOT37R+NnAEBSIwAiyTVo0OD0/rZ+/fq5Qgec8dZbb2nv3r0qVqyYWyoHACCpEQARED169FDnzp1dc+jbb79dGzZs4MpL7nr4ZkWt8jdlypRcFwBAkiMAIiBsWdP2udWpU0cHDx50x8Xt378/6q++7ZH86aef3D7Ju+++O+qvBwAgMAiACBjfCSFFihRxoadVq1Y6ceJEVP8EfMe+3XvvvcqYMWOwhwMAiBIEQARUrly59OGHH7qzg+fPnx/VLU9WrVqlpUuXumVfWyIHACBQCIAIuAoVKpzueTdy5EiNGzcuqmf/2rRpowIFCgR7OACAKEIARFBYtetTTz3l7nfr1s3NhEWTX3/99fQ5ydE8CwoACA4CIIJmwIABat26tY4fP64WLVq4UBQtXnrpJVcBXL9+fVWsWDHYwwEARBkCIIJaGTx+/HhVrlxZf/75pzsz+PDhwxH/E/nrr79OL3sz+wcACAYCIIIqffr0mjlzpvLmzat169bpzjvvdDNjkWzs2LEu6JYtW1YNGzYM9nAAAFGIAIigK1iwoGbMmKE0adK4MPj4448rUtlytxW+mD59+nDsGwAgKAiACAlXX321mxkzTz/9tN5//31FIjsXefv27cqTJ4/atWsX7OEAAKIUARAhw5Z/7axg06FDB61evVqRJCYmRkOHDnX3re+fzXgCABAMBECElGeeeUZNmjTR0aNHXVHIrl27FCk+//xzrVmzRunSpVPXrl2DPRwAQBQjACKkpEiRQu+++65Kly6tHTt2uH6BFgYjqfFzx44dlSNHjmAPBwAQxQiACDmZM2d2x8Vly5ZNK1euVJcuXdzyaTj78ccf9fHHH7uij969ewd7OACAKEcAREgqUaKEpk6d6mYE7dg43965cDVs2DD31pa17XsDACCYCIAIWfXq1XMnZpiHHnpIc+bMUTjas2fP6bOPafwMAAgFBECENDsn+N5773VLwG3atNH69esVbl5++WUdO3ZM1apVU61atYI9HAAACIAIbbZnbtSoUbr22mt16NAhNW3aVHv37lW4OHLkiF555RV3v2/fvjR+BgCEBGYAEfJSpUrl9gMWK1ZMmzdvVqtWrdyJGuFg4sSJLrAWLVrUVTQDABAKCIAICzlz5nSVwRkzZtTChQvDopLWzjQePny4u9+rVy+lTJky2EMCAMAhACJslCtXTu+8845bRrV9da+//rpC2UcffaRNmzYpS5Ysuvvuu4M9HAAATiMAIqzYHsDBgwe7+927d9eiRYsU6o2frYglU6ZMwR4OAACnEQARdh5++GG1bdtWJ06cUMuWLfXLL78o1KxatUpLlixxy7527i8AAKGEAIiwY0vA48aNU9WqVV2BhTVXtgrhUJz9s9Y1BQsWDPZwAAA4BwEQYSldunSaMWOG8uXLp++//1533HGHK7oIBb/99purWjY0fgYAhCICIMJWgQIFNHPmTKVJk8ZVCA8cOFChwE4vOXnypK6//npVrFgx2MMBAOA/CIAIa3a6xvjx4939Z555Ru+++25Qx3PgwAG3PG2Y/QMAhCoCIMKeFYRYYYjp1KmTK8AIlrFjx7r9iGXKlFGjRo2CNg4AAC6GAIiIYK1hbr75Zh09etQVhezcuTPgY7DTSWz51/Tp04dj3wAAIYsAiIiQPHly1yS6bNmy2rVrl5o3b65//vknoGOYMmWKtm/frty5c6tdu3YBfW0AAOKCAIiIYc2WrRgkR44cbhm4c+fOiomJCchr2+v4Wr9Y37+0adMG5HUBAIgPAiAiymWXXeZasFgDZisIee655wLyunYiybfffuva09x3330BeU0AAOKLAIiIU7duXY0aNcrd79+/v2bPnp3kr+mb/evQoYObgQQAIJQRABGRunbtqm7durmlWasStmbRSeXHH3/UnDlzXNFH7969k+x1AABILARARKwRI0bouuuu0+HDh9W0aVP9+eefSfI6w4cPd2/tNUqWLJkkrwEAQGIiACJipUqVSh988IGKFy+uLVu26LbbbnOtWhLTnj17NHHiRHefxs8AgHBBAEREs/14VhlsFcJWqHH//fcn6vO/8sorOnbsmDuRpHbt2on63AAAJBUCICKencphFcG2R++1115zoS0xWNPpZcuW6aqrrtJTTz1F42cAQNhIFuNHo7SDBw8qS5Ys7pzTzJkzB2ZkQCKzljB2ZFyKFCk0b9481atXj2sMAIgYcclrzAAiajz44IO64447dPLkSbcfcPPmzcEeEgAAQUEARNSwJeCxY8e6/Xr79u1zZwfbb0sAAEQbAiCiih3RNnPmTOXPn9/177MegTYjCABANCEAIurky5dPs2bNcmHQGjjbaSEAAEQTAiCiUpUqVTR+/Hh3//nnn9ekSZOCPSQAAAKGAIio1aZNm9Ozf/fcc4++++67i35+3bpSr17e/aJF7aQRBcygQdKVVwbu9QAAkY0AiKhm/fuaNWvmmjnH5SSPVaukLl0UMH37SgsWBOa1jh3zwmayZNKaNf/9+IQJUoUKtp9SyptX6t49MOMCACSelIn4XEDYSZ48uVv+rVmzZpzOCs6VSwGVMaN3C4QHH5Ty55fWrv3vx4YNk158UXrhBal6dWuGLf3yS2DGBQBIPMwAIurZMXF2XFy2bNn8vhbnLwFv2CDZSXA2K1amjDR/vjeDNnOmf8+3Y4fUurVkQ8iRQ2rWTPr119iXgDt0kJo3l4YMkfLkkbJmlZ54QjpxQurXT8qeXSpYUPr/bY5+++QTad48aejQ/35s/35p4EDJjj5u21YqXlwqW1a6+ea4vQYAIPgIgICkYsWKuWKQ+Dh1ygtj6dNLK1dKY8ZIAwb4//VHjkjXXefN8C1ZIi1b5t1v1Ej699/Yv27hQmnnTu9rbGbOQuJNN3kh0sbRtat327bNv3Hs2WN7ISWrh7Hv5XyffeZ9rxZWS5f2AmarVv4/PwAgdBAAgf9XuXLleF0LmzGzQ0VsZqxiRW8mcPBg/79+8mRbipbGjZPKl/fC1ZtvSlu3SosWxf51Nss3cqR0xRXS3Xd7by1MWl1LyZLSI49IqVNLy5dfegx2IKTNKlpgrFLlwp9jS70WAG3W0WY/p06V9u2TGjS4eFAFAIQe9gACCbRxo1SokFcQ4VOtmv9fv3q19PPPUqZM5z5u++sudlqdLb9acPSxpeBy5c68nyKFt5z8+++XHsOoUXaGpBcaY2Ph7/hxL3TecIP32Hvved/3559LDRte+nUAAKGBAAgkkM2e2X6/+LJgZZOP77wTt2KTVKnOfd/GcKHH7PkvxZaTv/xSSpPm3MdtNrBdO+mtt6yBtveY7XE8e3w5c3qzlQCA8EEABBKoVCkvANkeOpuF87WJ8VelStL770u5c0uZMwfnx2Gzek8/feZ921toM3o2Lqv2NbVqnZnxtP1/xpaArXi6SJEgDBoAEG/sAQQSyPbAWUVs+/bSunXenjtfEYg/M4M2w2azaFb5u3SptGWLtHix1LOntH17YH48hQt7y8e+2+WXe4/b9+ULe/aYjdHGtWKF9P333vdsAdiKWAAA4YMAiKh39OhRLV26VG+88Ua8roXttbN2L4cPS1WrSp07e+1SjLWFuRSruLVKXgthLVp4RSBW1PHPP8GbEYyNFbrYjGCTJtK113pLznPn/nfpGQAQ2pLFxNgOpos7ePCgsmTJogMHDihzqP2PBMTR4cOH9cUXX2jJkiXutnLlSncSyFVXXaVvvvkmUa6nzQJaNbAVd9gsWkJZcYbNDlqLGAAAEprX2AOIiLdv3z4tX778dOBbvXq1Tp48ec7n5MmTRzf4SlvjYcYMr3eftV+x0GfLpLZnLqHhz349s/YrdgzcVVcl7LkAAPAhACLi7Nq1yy3p+gLfd99995/PKVq0qK655prTtxIlSihZAkp5Dx3yjlCzpsi2n69+fe/INGN98+x2IXXqeKdvxObAAa/q1paWrb9ffCVkDACAyMMSMMKa7WD47bffToc9u/3000//+bxSpUqdDnt16tRRYdtwFyBWKWu3C0mXTipQIDrGAABIWiwBI6ID38aNG88JfNvOO4vMZvIqVqx4TuDLbT1WgsRO7LBbMIXCGAAAoYMlYIQ026tnS7hnB74//vjjnM9JmTKlqlatejrs1apVS1mzZg3amAEACHUEQISUf//91xVp+MKeFW9YNdPZ0qZNqxo1apye4atevboyZMgQtDEDABBuCIAIqiNHjrg2LL7AZ+1Z/rEGeGfJlCmTateufTrwVa5cWWnOP7MMAAD4jQCIgLLZvBUrVpwOfKtWrdLx48fP+ZwcOXKcU6Fr+/lSWLdlAACQKAiASFK2X2/ZsmWnA9+aNWt06tSpcz6nQIEC5wS+0qVLJ6glSyBZmO3Ro4crThkwYIBuvfXWYA8JAIBLIgAiUW3fvv102LNefOvXr//P51jPvbMDn/XkC5fAd76aNWuqZcuW6t+/v9q0aaNFixa5xwAACGX0AUS82azX5s2bz6nQ3bJly38+r1y5cue0ZMmfP3/EXYdWrVpp6tSpyps3rytiibTvEQAQ+ugDiCRhS7c2o3d24LNTN86WPHlyVapU6XTgs+IN29MXyWz28s0339SGDRv0/fffuxnBzz//nEIVAEDIYgkYsTpx4oTbs3f2kq6dq3u21KlTq1q1aqcDny1/WtVutMmYMaNmzpypKlWquErm+++/X6+//nqwhwUAwAURAHHa0aNHXVWu7xxd68F3+PDhc66Q9duzkOcLfBb+rC8fpOLFi+u9995T48aNNWbMGNeupkuXLlwaAEDIYQ9gFLNwZ7NVvhk+68d37Nixcz7HTtSwfXu+wHfVVVcpVapUQRtzOHjmmWdcUYhdJ4pCAAChuAeQABhF9u/ff05LFitWsKPWzpYnT55zKnStgMP29cF/FIUAAIKBIhA4u3fvPr2cazc7U9fCydmKFClyTuArWbJk2LZkCRUUhQAAQh17ACPIb7/9dk6F7qZNm/7zOaVKlTqnJUvhwoWDMtZIR1EIACCUEQDDlM3kbdy48ZzAt23btv/MRNkxame3ZLElXgQGRSFA4NWtK115pTRihFS0qNSrl3cLhEGDpJkzpTVrAvN6QEKwBzBM2F49W8I9O/DZMWtnS5kypWtD4gt8tWrVckUcCC6KQoDgBED7JzJDBil9+sC8tjVNsDq6QLQ+tdepXl1au1b69lvve/ZZtUp6+GFp9WqbCJCqVpWef/7cz0FkYg9gBDh+/Lh+/PFHffvtt/rmm29cP75Dhw6d/njBggXdDFP58uVd42W72f106dIFddz4r4cfftj9DO2kEDsrmJNCgMDIlSuwVzpjRu8WCA8+KNmBQxYAz2b/TTRsKDVrJr3yivVzlR5/3Hts+3aJJg44LcYPBw4csMoB9xZA3B06dCimXLly7u9RjRo1Yo4ePcplBJLAtdfGxPTs6d0vUiQmZvjwMx/78ceYmFq1YmLSpImJKV06Juazz6wqLiZmxgz/nnv79piYVq1iYrJmjYnJnj0mpmnTmJgtW858/PHHY2IqVjzzfvv2MTHNmsXEDB4cE5M7d0xMliwxMYMGxcQcPx4T07dvTEy2bDExBQrExLzxRty+x48/jokpVSom5ocfvPF/++2Zj61a5T22deuZx9at8x77+ee4vQ7CT1zyGv09gAAWhdiSvO+kEACBc+qU1Ly5txy8cqU0Zow0YID/X3/kiHTddd4M35Il0rJl3v1GjaR//4396xYulHbu9L5m2DBvn+BNN0nZsnnj6NrVu523hTtWe/ZI99wjTZp04aXtK66QcuaU3njDG9c//3j3y5a1rg/+f7+IfARAIMBFIVacYyeF2A1AYMybJ23eLE2cKFWsKNWuLQ0e7P/XT55sZ51L48ZJ5ctLpUtLb74pbd0qLVoU+9dlzy6NHOkFs7vv9t5amOzfXypZUnrkETtSU1q+/NJjsHm8Dh28wFilyoU/x07itPG8/bZkO4IspH76qfTxx7ZP3P/vF5GPAAgEUKNGjTT4///X6d69u1asWMH1BwJg40apUCEpb94zj1Wr5v/XW0HFzz97Acu318/C3dGjXrCMjc28nd1L3xoxWID0SZHCKxr5/fdLj2HUKNvk74XG2NiMnwXNWrWkL7/0gqWNoXFj72OAD78PAAFGUQgQeDZ7lpAe97aEXLmy9M47cSs2Ob/owsZwocfs+S/FlpMt1KVJc+7jNhvYrp301lvSu+9Kv/4qffHFmeBpj9mS86xZ0u23X/p1EB0IgECAcVIIEHilSnnLtbaHztcO1dql+KtSJen996XcuaVLHLGaZGwp+emnz7xvewututfGZS1hjC0vW/A7O+z63vcnZCJ6JA/1fk6+Bp7W0NP6OiUW+8tgDTuBYKAoBAisBg1sH67Uvr20bp23NOorAvFnZtBm2Ky4wtqrLF0qbdkiLV4s9ezptVcJBDu4qVy5M7fLL/cet++rYMEz3+f+/dL//if9+KP0ww9Sx47e/j8rYgHCIgCezX5T69Il8Z5v1y7pxhsVFPv2ST16eJuBrYrL/lJbUeiBA+d+XtOm3sfSppXy5ZPuvNP7jQ+RgaIQIHBsr5390m/Nmq0xcufO0sCB3sfs39hLsX+rrZLX/k1u0cIrArG9dravLlgzgrHNdM6e7YXcGjWkOnW8/zfmzvX+HwHC4iSQszu6R5Lvv/cac1o1V5kydoavV9VVoYI0deqZzxs+3PsLbH9pd+yQ+vb1HqduILJwUggQHDYLaNXAVtxhs2gJZcUZNjtoLWKAYIhLXgubGcDzl4A3bPD+4tpvbhai5s+P27Lu2Z9rG2bt/SlTvN+WrHTefkPctMmbebQNtr5+T+edvqbx470KK9uUa0Gte/dLv7ZN3U+bJt18s/ePTr16XjsC+63Nurb79O4tXX2117upZk3vaB/bAHz8uH/fI8KnKKRly5bu9Bc7KWQn07xAkpgxQ/rsM+/ffPs/w1aVrFo2oeHPplGsEnjBAu//AyAchE0ATMyGnrGxWTlbEvjmG2+/RJs23nE7L73k/VZnf8Efe+zM57/6qrfPwv4R+e476cMPpRIl4vfatvxrYT22Pk22bGzVZxYEOconMotCypUrp927d7sweMwO+gSQqOyYtG7dvGVSW4GxX/StMtYMGXKmvcv5t0ttF7J/v20iwvr5WX+/+ErIGICoqAL2NfS0Zpe+nk42g2abXxPCllitosrYxl4LgPYbnf2GaDp1kiZMOPP5Vo31wAPe5/rYPyhxtXev9NRT0r33/vdjDz0kjR7tVXbZbOBHH8X9+RE+RSFVqlQ5fVLI66+/HuxhASHn+++/d78sxcddd3m3C7FtOK1aXfhjlzpiPWtWKTF+Z0vIGIComAFMaEPP2NgePB9fm4CzG3baY75mnfbWVuquvz5hr2lNPZs08X57tBnI8/XrJ337rRd6bROz/eN16V2bCEcUhQAXN3XqVN1j56AlAWvqbCs4F7oVKBCYn0wojAHRI3k0NvSMzdlLq77nP/8xXx+lxPhtzJYjbF+hTe/b3pQLLe1a2wEr9bfZTTuKyI7zsX2AiEycFAL8l9UqDhkyRLfddhvbI4BoDoBnN/T0iUtDz8RgxwFZYYotEcd35u+GG7w9I7Z30J82BL6ZP7aHRTaKQoAzbD9sx44dNeD/N3q3bduWywNE6x7Asxt6Pv+8N5MWl4aeiWXQIG/PhnWGtw26Ng5rK2A9/i7GPs/Cn+3rswO7LQzazXekkC31fvWVd7NKZzvC55dfvAIU+76tNQwiFyeFAJ4///xTLVq00NKlS5UiRQqNGjVK9913H5cHiNYZwIQ29EwsFkCtNc0rr3il/zfdJP30k3+Hilv1slUO294Oax/ju23bdmaJefp0b4+hNYy2hqO279k6z59/DiQiDyeFINpt2LBBV199tQt/1s9szpw5hD8gWhpBJ1VDT1tCtaBo/aDq1w/UCIG4mzt3rho3buz2QFlVcJfEPA4HCFELFixw7ZD++usvFS1aVB999JHK0mAPiM5G0InV0NOWWt97zzsc2/YSAqGMohBEm7Fjx7o/9xb+atasqa+++orwBySB5NHW0NMO8rbees89d+bw7MRmDZtje31+iUVcURSCaHDy5En17dvXzXKfOHHCFXvYTGAu2xgNINFFzBLw+adm2O1CbG9dUvdTsnB6doXy2azVix3tBsTF4cOHVaNGDdcE195+/vnnSsNmUETQn+927drpQ2uJIOmJJ57Qo48+6gqiAPgvLnktIgMgEIk2b97sTgqxpTGbJeGkEESC7du36+abb9aaNWvcLzUTJkzQ7bffHuxhAWEpKvYAAtGGk0IQaVavXq1q1aq58Jc7d24tWrSI8AcECAEQCCMUhSBSTJ8+XXXq1NGuXbtckcfKlStd2xcAgUEABMIMRSEIZ7br6LnnntOtt96qf/75x/1Ss2LFCtfuBUDgEACBMD0ppFy5ctq9e7frl2bHZQGh7t9//1WnTp3cLzGme/fumj17NnvLgSAgAAJhiJNCEG727t2rG264wf3ykjx5cnesm91SpgzLE0mBsEcABMIURSEIF5s2bXLtixYvXqxMmTK5kz1s9g9A8BAAgTBGUQhCnfWstOKOn376SUWKFHH7/W688cZgDwuIegRAIMxRFIJQNX78eLfsu3//flWvXt1V+treVQDBRwAEwhxFIQg1p06d0kMPPeQKPuxYN2vsbDOBefLkCfbQAPw/AiAQASgKQaj4+++/XWX6888/795/7LHH9O677yqdncMJIGQQAIEIQVEIgm3nzp265pprNGPGDKVOnVpvv/22O9eXM32B0EMABCIIRSEIlm+++cYd62Zvc+bMqYULF6pdu3b8QIAQRQAEIgxFIQi0WbNmuWPdduzYoTJlyrhij1q1avGDAEIYARCIMBSFIJDHug0dOlS33HKLjhw54ip+rc3LZZddxg8BCHEEQCACURSCpHb8+HF16dJF/fr1c0Hwvvvu05w5c5QlSxYuPhAGCIBAhKIoBEnF+vrZftNx48a5Y91GjBihl19+mWPdgDBCAAQiGEUhSGw///yzO9nDijxspvnDDz9Uz549qfQFwgwBEIhwFIUgsSxZssSd6GFn+xYqVEjLly9XkyZNuMBAGCIAAhGOohAkhrfeekv169fXvn37XLuXr776ShUqVODiAmGKAAhEAYpCkJBj3fr3768OHTq4wo/bbrtNixYtUt68ebmoQBgjAAJRgqIQxJW1dmnVqpWeeeYZ9/7AgQM1efJkjnUDIgABEIgiFIXAX7t27dK1116radOmuWPdJk6cqKeeespV/QIIf/xNBqIMRSG4lLVr17p9fl9//bVy5Mih+fPn68477+TCARGEAAhEGYpCcDGzZ892x7ht375dpUqVcse62TFvACILARCIQhSF4Hx2msfw4cPVrFkz/f3337r++uvdsW62dxRA5CEAAlGKohD4WHVv165d1adPHxcE7Yi3Tz75RNmyZeMiARGKAAhEeVHIkCFD3P3u3bu7GR9El7/++kuNGzfWmDFj3PaAF198Ua+99ppSpUoV7KEBSEIEQCDKPfTQQ2rZsqWbBbr11lu1c+fOYA8JAbJ582bVqFHDFXlkyJBBM2fOdLOAFgQBRDYCIBDlKAqJTsuWLXPHum3YsEEFCxZ07zdt2jTYwwIQIARAABSFRJlJkya5Io+9e/eqSpUq7li3K6+8MtjDAhBABEAADkUh0XGsm53mcdddd+nff/9VixYttHjxYuXLly/YQwMQYARAAKdRFBK5/vnnH7Vp00aDBw8+3RD8gw8+UPr06YM9NABBQAAEcA6KQiLP7t27VbduXU2ZMsVV944fP96d78uxbkD0IgACOAdFIZFl3bp1rtjD9vllz55dn332mTp27BjsYQEIMgIggP/gpJDI8PHHH7tj3bZu3arLL79cX375pa699tpgDwtACCAAArggikLCl53mMXLkSN188806fPiwrrvuOn3xxRcqWbJksIcGIEQQAAHEiqKQ8HPixAl3qkvPnj1d1W+nTp00d+5ct/wLAD4EQAAXRVFI+Dhw4ICaNGmiV155xe3lfOGFFzR27FilTp062EMDEGIIgAAuiqKQ8LBlyxbVrFlT8+bNc61dpk+frr59+3KsG4ALIgACuCSKQkLbihUrXKXv+vXrlT9/fi1dulTNmzcP9rAAhDACIAC/UBQSmt59913Vq1dPf/zxh6666irX7qVSpUrBHhaAEEcABOA3ikJCq9J30KBBateunY4dO+Zm/Gzmr0CBAsEeGoAwQAAEECcUhQTf0aNH1bZtWz3xxBPu/QcffFDTpk1ThgwZgj00AGGCAAggTigKCa49e/a4vn6TJ09WypQpNW7cOD333HMc6wYgTgiAAOKMopDg+P77712xh53okS1bNlfxa33+ACCuCIAA4oWikMCyZs7W5uW3335TiRIlXAi0mUAAiA8CIIB4oygkMEaPHu0aPB86dMid5Wvhz872BYD4IgACSBCKQpL2WLcePXq4mx3r1qFDB7fsmyNHjiR8VQDRgAAIIEEoCkkaBw8eVNOmTd3sn3nmmWc0fvx4jnUDkCgIgAASjKKQxGX7/GrVqqVPPvlE6dKlcy1eHn74YY51A5BoCIAAEgVFIYnD9vdVq1bNVfzmy5dPS5YsUYsWLRLp2QHAQwAEkGgoCkkY6+1Xt25d/f7777ryyivdsW5VqlRJpJ8OAJxBAASQqCgKid+xbk8++aTatGnjjnW7+eab3bFuBQsW5E8ngCRBAASQqCgKifuxbnfeeacef/xx9/4DDzygGTNmuH2VAJBUCIAAEh1FIf75448/VL9+fb3zzjvuWLcxY8Zo6NChSpEiBX8qASQpAiCAJEFRyMWtX7/eHeu2fPlyZcmSxZ30cc899/CnEUBAEAABJBmKQi7ss88+U40aNbRlyxZddtllrvL3+uuv508igIAhAAJIUhSFnOu1117TjTfe6Bo9165dWytXrlSpUqX4UwggoAiAAJIURSGekydPqnfv3rrvvvvc/bvuukvz589Xzpw5+RMIIOAIgACSXLQXhRw6dEjNmjXTiBEj3PuDBw/WhAkTlCZNmmAPDUCUIgACCIhoLQrZunWrW+qdM2eO0qZNqylTpqh///4c6wYgqAiAAAIm2opC7CQPO9Zt3bp1ypMnjxYvXqzbbrst2MMCAAIggMCKlqKQqVOn6tprr9WePXtUoUKF02EQAEIBM4AAAirSi0LsWDeb5bSZPjvlo0mTJlq2bJkKFy4c7KEBwGkEQAABF6lFIRZkO3TooAEDBrj3e/XqpVmzZilTpkzBHhoAnIMACCAoIq0o5M8//1SDBg00ceJEd5Tbq6++quHDh3OsG4CQRAAEEDSRUhSyYcMGXX311Vq6dKkyZ86sjz/+WF27dg32sAAgVgRAAEEV7kUhCxYscMe6bd68WcWKFXNL2jfccEOwhwUAF0UABBBU4VwUYsvWDRs21F9//aVatWq5Y93KlCkT7GEBwCURAAEEXbgVhdhRbg888IDuvfded79du3buWLdcuXIFe2gA4BcCIICQEC5FIYcPH1aLFi00bNgw9/6TTz6pSZMmuVM+ACBcEAABhIxQLwrZvn276tSpow8//NCd4zt58mQ9+uijHOsGIOwQAAGElFAtCvn666/dSR5r1qxR7ty5tWjRIrVu3TrYwwKAeCEAAggpoVgUMn36dF1zzTXatWuXG5cVe1jbFwAIVwRAACEnVIpC7Fg3C6NPP/20SpUq5ZalbTxFixYNyngAILEki7F/4S7h4MGDypIliw4cOOCanAJAIMydO1eNGzd2Qez1119Xly5duPAAEIu45DVmAAGErFAvCgGAcEUABBDSQrUoBADCGQEQQEgLxaIQAAh3BEAAIS9UikIAIFIQAAGEhXA5KQQAwgEBEEDYFoWsXbv2op9ft67Uq5d33zq3jBihgBk0SLryysC9HgDEBQEQQNgWhfTr18/vr1u1SgpkF5m+faUFC5L2NZo2lQoXluwY4nz5pDvvlM6vkenZU6pcWUqThkAK4AwCIICwLQr5888//f66XLmk9OkVMBkzSjlyJO1rXHedNGWKtHGjNG2atHmz1LLluZ9jnV7vvlvi1DoAZyMAAgjbopBMmTL5/TXnLwFv2CDVru3NnpUpI82fb+FSmjnTv+fbscMLVdmyeUGvWTPp119jXwLu0EFq3lyyFew8eaSsWaUnnpBOnJBsIjN7dqlgQWn8eL+/JfXuLdmJdEWKSDVrSg8/LH35pXT8+JnPGTlS+t//pMsu8/95AUQ+AiCAsC0KeeaZZ+L1tadOeWHMZgRXrpSsnmTAAP+//sgRb/bNZvmWLJGWLfPuN2ok/ftv7F+3cKG3RGtfM2yYFxJvuskLkTaOrl2927Ztcf+e9u2T3nnHC4KpUsX96wFEFwIggLBV09JOPMyb5y2XTpwoVazozQQOHuz/10+eLCVPLo0bJ5UvL5UuLb35prR1q7RoUexfZ7N8NiN3xRXesqy9tTDZv79UsqT0yCNS6tTS8uX+j+Whh6QMGbxZSHv9WbP8/1oA0YsACCDq2J65QoWkvHnPPFatmv9fv3q19PPPkq1A28yf3SzcHT3qBcvYlC3rBUcfWwq2AOmTIoUX5H7/3f+x2PLxt996oda+/q67vH1/AHAxKS/6UQCIQBaQbL9ffNkSslXW2pLrhYpNYnP+0qyN4UKP2fP7K2dO73b55d5MpAVb2wdYo4b/zwEg+hAAAUSdUqW85dI9e7xZOF+bGH9VqiS9/76UO7eUObNChm/mj5PyAFwKS8AAok6DBlZEIrVvL61b5+258xWB+DMz2K6dN+tmlb9Ll0pbtkiLF3s997ZvV0B89ZU0erS0Zo3022/S559Lbdt639fZs3+2VG2fs3u39M8/3n27XaxYBUDkIwACiDq2V87avRw+LFWtKnXuLA0c6H3M2sJcilUPWyWvNWFu0cJberWiDgtYgZoRTJdOmj5duv76M0Ul5cp5QdSaPvvY93bVVdLrr0ubNnn37XZ+w2gA0SVZTMyltwsfPHhQWbJk0YEDB5Q5lNY7ACCR2CygVQPbjJnNoiWUVfTa7KC1iAGAQIhLXmMPIICoNGOGV71r7Vcs9Nnyba1aCQ9/9iv1L794x8DZTBsAhCKWgAFEpUOHpG7dvIIQO6XDloJ9PfTstA5fe5fzbzfeePHnPXDAO1nE+vlZf7/4SsgYAOBSWAIGELbmzJmjJk2aJPrz2qkadott712BAon+kiE5BgDhhSVgABHt1KlTGjhwoObOnZskAdCaOtstmEJhDAAiF0vAAMLK4cOHdeutt8b7HGAAAEUgAMLI1q1b1bRpU61du1apU6fW008/HewhAUBYYgYQQFj48ssvVa1aNRf+cufOrUWLFqlx48bBHhYAhCUCIICQ9+6776pu3bras2ePKlSooK+++ko1OOwWAOKNAAgg5Is92rVrp2PHjrnl3+XLl6tIkSKxfs0nn3yiSpUqqWvXrok2jiNHjmjq1Km67bbb3HP7bvfcc4/mz5+v48ePJ9prAUAgEAABhKS///7bBa7Bgwe79x966CHNmDFDGa0R3kXs3r1b3377rbZt25ZoY0mfPr1atmypKVOmaNiwYbrsssu0bt06jRs3Tg0aNFCxYsXcfkSboQSAcEAABBBytm/frjp16mj69Omu2GPChAl69tlnlTz5pf/JsplCY1+X2JIlS+aWom02cMuWLRowYIBy5cqlHTt26NFHH1WhQoV0xx13uP2KfpyyCQBBQwAEEFJsf1/VqlXdLJ6Fq4ULF6p9+/Z+f/2///7r3qZJkyYJRykX9mzWz2YaJ02apOrVq7ul4HfeecftT7TvwYLrP//8k6TjAID4IAACCBmTJ0/Wtdde65Zxy5Ur58JgLTugNw58M4BJHQB97HV8s342Xgur9tjq1avVsWNHFxQffvhh/fbbbwEZDwD4gwAIICSKPR5//HG1adNGR48e1U033aQVK1aoaNGicX6uQAfAs/lm/WxW0BpVFy5cWHv37tVzzz3n9g02b97cFY2wPAwg2AiAAILKKmxvv/12Pfnkk+79vn37aubMmcqUKVO8ni+YAdDHlq5t1m/z5s2ucKV+/fou5M6aNcsVjZQpU0ajR49253YCQDAQAAEEjRVPXHPNNfrggw+UKlUqvfHGG3rhhReUIkWKeD+nbw9gUhSBxFXKlCndrN9nn32m9evXq3v37q6KecOGDerRo4cKFCjgHvvxxx+DPVQAUYYACCAovv76a3eyh+2Vy5Ejh1savfvuuxP8vKEwA3ghpUuX1qhRo1zotdm/UqVKuXONX375ZTcjaLOENvN54sSJYA8VQBQgAAIIOOunZzN/O3fudOHHiifs/cQQqgHQJ3PmzPrf//7nZgQt9NoMobW3WbBggW655RYVL17ctbz5448/gj1UABGMAAggYKz4wfb6tW7d2rVHsbN8v/jiC1cgkVhCPQCe3VPw+uuvd3sEf/nlF7dn0GZCt27dqkceecRVD1tF8apVq4I9VAARiAAIICAs8LVt29ZV+5revXvrww8/dDNiiSkpG0EnFTvazqqGrQG2VRFXrlzZfR8TJ050y+TWY9B6Dfq+NwBIKAIggCS3a9cu19/P+vxZYcTYsWPdkWoJKfYIdiPopJA2bdrTs37WV9D6C1qQtSXyu+66y80K2ukjiXnMHYDoRAAEkKS++eYb1x/PQk327NldRWznzp2T7PXCZQn4UsvDvlk/WxK2E0esYtj2BQ4ZMsT1R7z11lv1+eef01MQQLwQAAEkmWnTpql27dqu8tWqXm0my87STUqREADPlidPHjfr9+uvv7oziO36WU9BOye5Xr167sSUV1991VUUA4C/CIAAkqTYY/DgwWrZsqXb+9ewYUO3pGkVrkktHPcA+sOWzn2zft999526du2qDBkyuGribt26uRnCnj17auPGjcEeKoAwQAAEkKjsKDfbuzZw4ED3voWSjz76SFmyZAnIlQ7nPYD+8s362czqSy+9pJIlS7pTRUaOHOlmWi1wz549WydPngz2UAGEKAIggESze/dut0T57rvvuhmr1157TSNGjHD3AyXSloAvxkL1/fff704W+fTTT3XzzTe7/YPz5s1T06ZNVaJECXeyip1HDABnIwACSBRr1qxxLUtWrlypbNmyuUBy7733BvzqRlMA9LFG0jfccINrq2PnD/fr18/9DGzf4IMPPqiCBQuqU6dOriAHANy/G1wGAAllR5jVqlXLtSe54oorXAi0AoVgiMYAeLZixYrp+eefdz0F7WzlK6+80i3Ljx8/3vUXtJ+TzdD6lsoBRCcCIIAEFXvYsWV2hNmRI0fUoEEDV+xhe9KCxRdsIq0IJK7Sp0/vzla2Wb/ly5erTZs2bil+xYoVateunQoXLqzHHnvM7SMEEH0IgADixWaVrGmxHVtmunfvro8//lhZs2YN6hWN9hnA89mewJo1a7pZP+sp+MQTTyhfvnzas2ePnnrqKXcKSatWrbRkyRJ6CgJRhAAIIM4sPNgSrzUqttM8Xn75ZY0aNSqgxR6xIQDGzoKfzfr99ttvev/991WnTh1XKfzBBx+4k1oqVqyoMWPG6O+//w7gTwxAMBAAAcTJunXrXLHHF1984Wb75s6d6/rQhQoC4KWlSpXq9Kzf2rVr1aVLF6VLl871F7TCHesp2KdPH/38888B+IkBCAYCIAC/WZWpLSfaUqLt87P9fvXr1w+pKxipjaCTSoUKFfT666+7vYB2PrM16z5w4ICGDx/ufsaNGzd2S/t2+giAyEEABOBXsYdVljZv3twtD15//fUu/FnFbyix5UxfUGEPYNxY25jevXtr06ZNLvBZ8DOffPKJmjRpossvv9wFxP379yfBTw5AoBEAAVxyRq1jx4566KGHXBC87777XCjInj17yF053+yfIQDGv6fgjTfeqDlz5uinn35yS8G21G/9BR944AG3PGxLxrZ0DCB8JYuxf9EvwY4Yso7ztiyQOXPmwIwMQND98ccfrsWLtRGxYGDHjlm1b6iy2SlfMLV2MLbXDQlns75WRWyFPrZP0Kd27druz0OLFi241kAIiEteIwACuKDvv//eHS1mp0nYPyhTpkxxp02E+lF0VulqbCnYWqAg8dh8wbJlyzR69GhNmzbt9FnDds2teMRmBn3XH0BoB0CWgAH8x0cffaQaNWq48GdFAbbfL9TD39lNoG35l/CX+OyaWusYayFjrWSspUyePHm0a9cuDRo0yDWXtobTNmPsx+ISgCAiAAI4zf7TfvHFF9W0aVMdPnxYdevWdce6lSpVKiyuEi1gAsf2AlpTaasIt+Vhqw4/ceKEJk+e7JaGK1Wq5I6isxNiAIQeAiCA07NnnTt3Vt++fV0QvOeeezRv3jzlyJEjbK4QATDwrN2Ob9Zv9erV7vi5tGnTas2aNe7PU8GCBdWvXz/98ssvQRgdgNgQAAHozz//dOf4jh8/3hV7jBgxwvWGC7ciCgJgcPlm/bZv3+7aBhUtWtQV5gwdOlQlSpRwe0o//fRTegoCIYAACES59evXu5M97FQI2zRs7T969uwZlnvofHsAaQIdXDZrbLN+dpKINQ+3/aM2q2x7Sxs1auS2FFhFuW1UBxAcBEAgilk/Pyv22LJliy677DJ3vJv9Bx2umAEMLXZOtG/Wb8OGDbr//vvdLxnWX7BXr15uH6H1lbSKcwCBRQAEopDNxthRXzfddJNrG3DNNde4Yo8yZcoonBEAQ5edGmOzfrY8/Morr6hs2bKuv+Brr72m8uXL67rrrnOtZayQBEDSIwACUcaWSa1nm53wYL3yOnXqpM8++0w5c+ZUuCMAhr5MmTK5WT9rKP3555/r1ltvdTOFixYtUsuWLVWsWDENHjxYv//+e7CHCkQ0AiAQRfbu3ev2Y40dO9bt8bOzXe1+pOyZ8wXASPl+Ipn9+bM2Q1OnTnVbEAYMGKBcuXK5GcKBAweqUKFCuvPOO93MND0FgcRHAASixI8//qjq1atr8eLFbhZm9uzZ6t27d1gWe/jTCBrhw8Le008/rW3btmnSpEnuz6n9LN9++21dffXVrkjprbfe0tGjR4M9VCBiEACBKGCb8O0/0s2bN7vWHCtWrFCTJk0UaVgCDm8W3O+44w538sxXX32l9u3bu8e+/vprdejQwfUUfOSRR9wpJAAShgAIRDBbOhs1apQaN27sij3shAb7j7VcuXKKRATAyFG1alVNmDDBzQo+88wz7pg528Lw7LPPuor15s2ba/78+SwPA/FEAAQi1PHjx91me2u9YcUeNoNi/2HaPqtIxR7AyGN/Xh9++GE3ez1jxgxdf/317s/zrFmzXPNyq1wfPXq0+wUHgP8IgEAE2rdvn+vnZ6d52B4/O5XBTvmI9L1x7AGMXClTpjw962fNy//3v/8pY8aMrr9gjx49XE/B7t27u72uAC6NAAhEmI0bN7pN9AsXLnT/QdpMiZ3KEEnFHrFhCTg6lC5d2s367dixw721k0UOHz6sl19+2c0I1q9fXzNnzqSnIHARBEAgglg/Pwt/dgSX7Zlavny5O4khWhAAo4udKmIzgTYjaDODNkNoZ1kvWLBAt9xyi4oXL+72DP7xxx/BHioQcgiAQISw2Y8bb7zRna9as2ZNrVq1ShUqVFA0IQBGJ5vdtr2Btkfwl19+cXsG7TzirVu3uqphazNje2Ctmjip1K0r9erl3S9aVBoxQgEzaJB05ZWBez1EBgIgEObs6CybBbH9TydPntRdd93lln9z586taOPbA0gj6OhVpEgRVzVsDaWtirhy5cruFwPrI2iVxdYOyfoL+n5ZSAqrVklduihg+vaVFiwIzGvZZbOwaTtK1qw58/jatVKbNtbTUUqXzpbppZdeCsyYED8EQCCM7d+/38362dmqNgtiy132n16kF3vEhhlA+KRNm9b1EbSZcOsraP0FU6VK5U4WsRNGbFbQTh+xNjOJzQrt06cP3M8iY0YpR47AvNaDD0r58//38dWrve/77belH36QBgyQHnlEGj06MONC3BEAgTC1adMmN5the58yZMig6dOn66GHHoqKYo/YEABxPvv7YPti7YQRC3t24ohVDNu+wCFDhrizh+08YjuXOLGOnDt/CXjDBql2bQulUpky0vz53gzazJn+Pd+OHVLr1lK2bF7Qa9ZM+vXX2JeAO3SQmjeXhgyR8uSRsmaVnnjCVgukfv2k7NmlggWl8ePj9n198ok0b540dOh/P3b33dLIkdK110qXXSbdcYfUsaM0fXrcXgOBQwAEwpBtcrf/1CwE2kyGFXvYBvhoRwDExeTJk8fN+v3666/uDGI7i9i2TdgvT/Xq1XMN0l999VVXUZxYTp3ywpjNCK5cKY0Z482O+evIEem667xZviVLpGXLvPuNGtmWh9i/buFCaedO72uGDfNC4k03eSHSxtG1q3fzdwJ0zx7pnnukSZP8n908cMALmwhNBEAgzLz22mtq2LCh/vrrLzcDaCd7VKxYMdjDCgk0goa/PQV9s37fffedunbtqvTp07tq4m7durkZwp49e7qWSgllM2abN0sTJ0r219RmAgcP9v/rJ0+WkieXxo2Typf39ta9+aa0dau0aFHsX2fBy2bkrrjCm52ztxYm+/eXSpb0lmdTp5aWL7/0GGxi1GYVLTBWqeLfuL/4QpoyRbr3Xv+/VwQWARAIo2IPO9XDTvewWYt27dq5/8Dy5s0b7KGFDBpBI658s37WU3DEiBEqWbKkO1Vk5MiRrr+g/bI1e/Zs93cuPixDWmHE2X9Nq1Xz/+ttb93PP0uZMnkzf3azcHf0qBcsY1O2rBccfWwp2AKkT4oU3nLy779fegyjRkl20IqFRn/YHkBbpn7sMalBA/++BoFHAATCgM32NWnSxJ3rawYPHuz2NNlGd5zBEjDiK2vWrG7Wz04WmTt3rm666Sa3f3DevHlq2rSpmlmiiQebPUvItlxbQq5c2au4Pfu2aZPUtm3sX5cq1bnv2xgu9Jg9/6XYcvKXX0pWW5YypVSihPe4zQa2b3/u565fL9Wr5y0XDxzo97eJIEgZjBcF4D9r6mzNnO0/JlumsuDXokULLuEFEACRUNZI2mb97LZlyxY3Ozhu3DjttA118VCqlLdca3vobBbO1ybGX5UqSe+/L1lXp8yZFRS2lPz002fet0vRsKE3rurVz535s/BnoTAuy9wIDmYAgRBmS7xW7GHhr2DBglq2bBnh7yLYA4jEZBXCdo629RR8/PHH4/UctgRavLgXitat8/bc+YpA/JkZbNdOypnTW1JdulTaskVavFjq2VPavl0BUbiwLZWfuV1+ufe4fV9WTewLf1asYt9vnz7S7t3ejUNYQhcBEAhRY8eO1Q033KB9+/apWrVqrtjjqquuCvawQhp7AJEUbOY9vkvAttfO2r1YYXHVqlLnzmeWRv3ZwWEVt1bJayHMJv6tCMSKOv75J3gzghfywQde2HvnHSlfvjM3+54RmlgCBkKw2KNfv35uQ7q5/fbbNX78eKWz9vq4KJaAESxnV+Se3aPPtwxs7Vt8fJW3vr10l2IFJG+9dfHTOaw4xGfChIuPL7ZxxqXP4fktE63NjN0QPgiAQAixc3zbtGmjT6zjqqQnn3xSAwcOjOrmznFBAEQomjHDC2jWfsUqem35tlYtbwk1ISyE/fKLdwwciwOIK5aAgRCxefNm1ahRw4U/m+374IMP9OijjxL+4oAAiFB06JDUrZs3E2j99GxZdNYs72N2Woevvcv5txtvvHSjZTtZxPr5WX+/+ErIGBC+ksX4cfaN9UTKkiWLm53IHEqbDoAIsWTJElfcsXfvXuXPn18ffvihO8QecT/p4ffff9e6detU/uymZ0CI2rfPu12I7fooUCA6xoDEEZe8xhIwEGRvvPGGa+58/PhxValSRbNmzXIhEHHHDCDCjTV1DvZxaaEwBgQeS8BAkNjJAg888IA6d+7swl+rVq20ePFiwl8CEAABwD8EQCBI0/R2usAwO6XdVdAN0uTJk127CcSP7WYhAAKAf1gCBgLMThewkz1++OEHd5TbhAkT1Lp1a34OidA+x7elObXtigcAxIoACASQneRxyy236M8//1S+fPncfr+qdEpN1CbQJo0dWgoAiBVLwECA2ExfvXr1XPirVKmSVq1aRfhLRL7lX0MABICLIwACASj2ePDBB9WxY0dX7NGyZUvX9qUAvRWSJAAmT55cKVOyuIHwcejQIY0ePVo1a9Z0vxza7f7779dPP/0U7KEhghEAgST+h92WfF944QX3vjV2fv/995UhQwauexIFQPb/IdxkypRJ3bt317Rp01S9enXXx3LUqFG64oordMcdd+gXO+4DSGQEQCCJ/Pbbb6pVq5Zmz57tliTfeecdd7SbzVAh6fYAsvyLcGX7gl999VX9+OOPrjDMiprs341SpUqpR48e2rNnT7CHiAjC/0RAElixYoWqVaum7777zp1OYf392rZty7VOQrSAQaQoWbKkawu1evVqNWzY0G0dsSXi4sWL67HHHnNtpICEIgACiWzSpEm67rrr3JFkV155pSv2sGUdJC0CICKN7QWcO3euFi5c6H6h/Pvvv/XUU0/psssu0/Dhw3X06NFgDxFhjAAIJJJTp07pkUce0V133eWWI23vn7V9KVSoENc4AAiAiFT2C+WXX37p9gjavkA7M7xPnz7uvnUXsEIzIK4IgEAiOHz4sG699VY9++yz7v3+/ftr6tSpFHsEYQ8gRSCIRMmSJVOLFi30/fffa9y4cSpYsKC2bt3qugtUqFDB9RT1NUIH/EEABBLI/hGuXbu2Zs6c6cKHLQEPHjyYYo8AYwYQ0cBaHHXq1EmbNm1y3QWyZcum9evXq3nz5q6NjO03BvxBAAQS4IsvvnB7c9auXavcuXNr0aJFrm0DAo8AiGiSLl069e3b17WIsRUHO0fclonr1q2rxo0ba82aNcEeIkIcARCIJ2vPYHtzrDWDLcF89dVXqlGjBtczSAiAiEZZs2Z1Kw4///yzunXr5mYIP/nkE1111VVq166dNm/eHOwhIkQRAIF4FHsMGDDAzfRZ6GjatKmWL1+uIkWKcC2DiEbQiPYegi+//LLrIdimTRv32Lvvvut6CFqT6d27dwd7iAgxBEAgDqwNw2233aYhQ4a49x966CHNmDFDGTNm5DoGGY2gAalEiRIu+H3zzTdq1KiRTpw44YKh9RC0k4gOHDjAZYJDAAT8tG3bNtWpU0fTp093xR7WfsGqfjnZIzSwBAycYUvAthT8+eefuz6kR44c0dNPP+16CL744ov0EAQBEPCH7e+zYo9vv/1WuXLlco1Z27dvz8ULIQRA4L+sKMSK1WylonTp0tq3b58rHrHTRsaPH+9mCBGdmAEELuHTTz9Vz5493R4bW/616jo74xehhT2AQOw9BK1NzLp161zos+b027dvd+1krIDNwiE9BKMPARC4BDuL036Dtj01U6ZMUf78+blmIYg9gMDFWYWwNY62HoK2DJw9e3ZXNGINpq2DgbWxQvQgAAKICCwBA/5JmzatO0rOeggOHDjQ9RBcuXKla2tlhSO21QWRjwAIICIQAIG4yZIli5566inXK/B///ufmyG0LS+VKlVyrWSstyAiFwEQQEQgAALxkzdvXo0ePVobNmxQ27Zt3WOTJ092RSPWXHrXrl1c2ghEAAQQUXsArUUPgLizXoF2wpEtAd94442uQvjVV191vQWt+T09BCMLARBARGAGEEgcV155pT7++GNXFHL11Ve7HoLW/N56CA4dOlT//PMPlzoCEAAR0erWlXr18u4XLSqNGBG41x40yP4hDdzrRTsCIJC4rr32Wq1YsUIzZ85UmTJlXA/Bfv366fLLL9cbb7xBD8EwRwBE1Fi1SurSJXCv17evtGBB0r6Ghdpkyc69PfzwuZ+zdat0881ShgxSzpzS/ffbcqkiDgEQSJoegs2aNXM9BN98883TPQQ7d+6s8uXLu5OR6CEYngiAiBq5cknp0wfu9ex44Bw5kv51nnxSsj3avtvAgWc+dvKk1KSJnWEsLVtmG7uladOkBx5QxKERNJB0UqRIoQ4dOrgegsOGDVOOHDlc0citt97qlontyDmEFwIgosb5S8AbNki1a1tPLKlMGWn+fG8GbeZM/55vxw6pdWspWzYv6DVrJv36a+xLwB06SM2bS0OGSHnySFmzSk88IdlJTP36SdmzSwULSuPHx+37ypTJqvjO3Cx4+sybJ61fL739tp0NKtWvL734ojR2rHTwoCIKjaCBwPQQ7N27t+sh+OijjypDhgzuqMx69eq5pvnWMB/hgQCIqHTqlBfGbEZw5UppzBhpwAD/v/7IEem667ywtWSJN7tm9xs1uvjy6sKF0s6d3tcMG+aFxJtu8kKkjaNrV++2bZv/Y3nuOS+AWtgcPPjc1//iC6lcOensw0saNrTZMmn1akUUloCBwMmcObOefPJJ10Owe/fuSpUqlebNm6fKlSvr9ttv108//cSPI8QRABGVbGZs82Zp4kSpYkVvJtDCk79sKTV5cmncOKl8eal0aenNN739dhc7Tclm+UaOlK64Qrr7bu+thcn+/aWSJaVHHrE2JtLy5f6No2dPbyy2+tK9uzfD2a3bmY/v3u3NNp7Nwqa9hn0skhAAgcDLkyePRo0a5ZaD27Vr5/YMvv/++65o5L777qOHYAgjACIqbdwoFSrkLZn6VKvm/9fb7Jk1ybflV5v5s5uFu6NHvWAZm7JlveDoY+HMAqRPihTebN7vv/s3jt69rVJPqlBB6txZeu016Y03pL17z3yOLWufLybmwo+HM/YAAsFjLWLefvtt10OwcePGrkL4tddec70F+/fvr7/++osfT4ghACIqJTQA2RJy5crSmjXn3jZtkv6/kf4FpUp17vs2hgs9Zs8fH1df7b31neBkAff8mb79+6Xjx/87Mxju2AMIBF/FihU1Z84cLV68WDVr1nQ9A5955hkXEJ9//nl6CIYQAiCiUqlS3nLtnj3ntonxV6VKkm1xyZ1bKlHi3FuWLAoa3xnu+fJ5b2vUkL7/3qsOPnv5O00aL8BGEpaAgdBxzTXXaNmyZZo1a5bKli2r/fv366GHHlLJkiU1btw4egiGAAIgolKDBnbskdS+vbRunbfnzlcE4s/MYLt2Xk89q/xdulTaskVavNjbk7d9uwLCCjyGD/dmHu31p0yR7r1XatpUKlzY+5wbbvAqnO+80wuH1pfQ+hPec49t4lZEIQACocX2AzZt2lRr167VhAkTVLhwYe3YsUP33HOPypUrp6lTp9JDMIgIgIhKttfO2r0cPixVrertn/P1z7O2MJdi1cNWyWtBq0ULrwjEijrshKRABSubxXv/fe+0Ewt5jz3mBbv33jv3+5wzx/ueatWSWrXyqp+HDlXEIQACodtDsH379q6H4IgRI5QzZ05t3LhRt912m6pVq6YFSd0xHxeULMaPFt4HDx5UlixZ3EHQVvoNRCKbBbRqYNs/Z7ODCWUVvTY7aC1ikPTs3yj7t8r+k7FlJgChyf6eWjPpF198UYftt3C3KtPA7RW0NjKIv7jkNWYAEbVmzJA++8xr3mxNoO2YOJslS2j4s1+prBLYfqm1ql8EBjOAQHiwYDJo0CDXQ/D+++93PQQ/++wzValSRa1bt3a/xCHpEQARtQ4d8nrmWUGIndJhS8GzZnkfs9M6fO1dzr/deOPFn/fAAW9J1nrtWX+/+ErIGKKNLWQQAIHwkjt3br300ktuOfjOO+90ewanTJniegjee++92mld85FkWAIGLmDfPu92IenSSQUKRMcYwqkFTBrbFOna3OxXVjtnD0BYWbdunQYMGKCPPvrIvZ8uXTo3Q2jVw9msgz0SdQmYAAgg7Nk+okzWlVvS33//rfRWpQMgLFn7mIcffljL//9IJPuFzt7v0aMHf7cvgT2AAKKKb/nX+GYCAYSn2rVra+nSpfrwww9duxg7RcQCoBV3jRkzRsetkz0SjD2AACImAFq7CbsBCG+2H/Dmm2/WmjVrNHHiRBUpUsTtCbS9gRYKP/jgA3oIJhABEEDYowAEiEz2C50ViFihiBWM5MqVy1UJt2rVSlWrVtV8a+GAeCEAAoiYAJjaSq8BRBzb2mEFIdY6xlrIZMyYUatXr3b9A+329ddfB3uIYYcACCAiqoAN+/+AyGbFXo8//rgLgj179nS/9NksoM0G2skiNlMI/xAAAYQ9loCB6OshaMfKWeC766673J5BO1u4bNmy6tKliztzGBdHAAQQ9giAQHQqWrSo3nrrLddD0IpGTp48qbFjx6pEiRKuf+C+2JqpggAIIPwRAIHoZpXB1jbGeghaG5mjR4/q+eefV/HixfXss8/qyJEjwR5iyGEGEEDE7AGkCASIbrVq1dKSJUvcaSLly5d3PQQfeeQRNyP4+uuvR3wPweTJ/Y91BEAAYY8ZQAA+th+wSZMm+vbbbzVp0iS3TLxr1y517drV7RG084ZPnToVkRfMqqP9RQAEEPYIgAAu1EPwjjvu0IYNGzRy5EjXQ/Cnn35S69atXdXwvHnzorqZNAEQQNgjAAKIjbWHsnOErXXME0884VrJfPPNN2rYsKHq16+vVatWReXFIwACCHvsAQRwKRb8HnvsMRcEe/Xq5fYML1y4UNWqVVPLli3dTGE0IQACCHvMAALwly0FDx8+3B0p16FDB1c4MW3aNLc/8J577tH27duj4mISAAGEPQIggLgqUqSI3nzzTddDsFmzZq4wZNy4cSpZsqQefPDBiO8hSAAEEPYIgADiq2zZspo5c6aWL1+uOnXquB6CL7zwgi677DI988wz+vvvvxN8cevWlXr18u4XLSqNGBH8nxcBEEDEBED6AAKIr5o1a2rx4sWaM2eOKlSooAMHDqh///6uh+Crr76aaD0EreakS5fg/5wIgAAipgjEqv0AICE9BBs3bux6CL799tsqVqyYdu/erW7duql06dKaPHlygnsI5solpU8f/J8RARBA2GMJGEBiSp48udq1a+cqg0eNGqXcuXO76uE2bdqoSpUq+vTTT+PdQ/D8JWArPq5dW0qbVipTRpo/34KoNHOmf8+3Y4fUurWULZuUI4f/4yAAAgh7BEAASSF16tTq3r27C39PPvmkayVjs4ONGjVSvXr1tHLlygQ9v00mNm/uzQjaU40ZIw0Y4P/X2xHH111nJ4BIS5ZIy5b5/7UEQABhjwAIIKmPWHv00Uf1yy+/qE+fPi4YLlq0SFdffbVatGihLVu2xOt5582TNm+WJk6UKlb0ZgIHD/b/6ydPttlKadw4qXx5qXRp/7+WAAgg7NEIGkAg5MyZUy+++KI7Uq5jx45uqXjGjBm67bbb4vV8GzdKhQpJefOeeaxaNf+/fvVq6eefrcm1NwsYh6OACYAAwh8zgAACqXDhwho/frzrIdi8efN4F4bYNkLb7xdf9rKVK0tr1py5+Stl/F8WAEIDARBAsHoIzpgxwwXB+ChVStq6VdqzR8qTx3ssLkcTV6okvf++lDu3lDlz3F6bJWAAYY8ACCCYKlSoEK+va9BAKl5cat9esgy5fPmZIhB/ZgbbtbNlaalZM2npUikuWxEJgADCHnsAAYSjFCm8di+HD0tVq0qdO0sDB3ofs7Ywl2LVw1b9W7iw1KJF3IpAWAIGEPaYAQQQyhYtOnP/11//uwx8dvsWmwU0JUr499xWQPLWW2fe93dPIQEQQNgjAAIIVzNmeNW7JUt6Fb09e0q1anlLw0mJJWAAYY8ACCBcHTokdevmzQR26OAtBc+a5X1syJAz7V3Ov914Y8JelxlAABETAK05KwCEk7vu8m4X0rWr1KrVhT+WLl3CXpcACCBiikDSpEkT7KEAQKLJnt27JQWWgAGEPZaAASBuCIAAwh4BEADihgAIIOwRAAEgbtgDCCDslSlTRsePH1eGDBmCPRQACAsEQABhb+XKlcEeAgAE3WE7UsRPLAEDAABEgFOnTvn9uQRAAACAKEMABAAAiDIEQAAAgChDAAQAAIgyBEAAEaNuXalXL+9+0aLSiBGBe+1Bg6Qrrwzc6wFAQhAAAUSkVaukLl0C93p9+0oLFgTmtY4d88JmsmTSmjVnHp8wwXvsQrfffw/M2ACEB/oAAohIuXIF9vUyZvRugfDgg1L+/NLatec+3rq11KjRuY916CAdPSrlzh2YsQEID8wAAohI5y8Bb9gg1a4tpU1rJ4dI8+d7M2MzZ/r3fDt2eAErWzYpRw6pWTPp119jXwK24NW8uTRkiJQnj5Q1q/TEE9KJE1K/flL27FLBgtL48XH7vj75RJo3Txo69L8fS5dOypv3zC1FCmnhQqlTp7i9BoDIRwAEEPGsN6qFsfTp7dQQacwYacAA/7/+yBHpuuu8Gb4lS6Rly7z7Ntv277+xf52Fr507va8ZNswLiTfd5IVIG0fXrt5t2zb/xrFnj3TPPdKkSd73cikTJ3qf17Kl/98rgOhAAAQQ8WzGbPNmLxBVrOjNBA4e7P/XT54sJU8ujRsnlS8vlS4tvfmmtHWrtGhR7F9ns3wjR0pXXCHdfbf31sJk//5SyZLSI49IqVNLy5dfegwxMd6sogXGKlX8G7fNLrZt680MAsDZ2AMIIOJt3CgVKuQti/pUq+b/169eLf38s5Qp07mP2946C5axKVvWC44+thRcrtyZ922J1paT/SnQGDVKOnjQC43++OILaf16L/QCwPkIgAAins2e2X6/hCwhV64svfNO3IpNUqU6930bw4Ue8+f4TltO/vJLKU2acx+32cB27aS33jr3cZuttD2JNm4AOB8BEEDEK1XKW661PXQ2C+drE+OvSpWk99/3KmkzZ1ZQ2FLy00+fed/2FjZs6I2revVzP/fwYWnKFOmZZwI+TABhgj2AACJegwZS8eJS+/bSunXenjtfEYg/M4M2w5Yzp1f5u3SptGWLtHix1LOntH27AqJwYW/52He7/HLvcfu+rJr4bBYKrdrYxg0AF0IABBDxbK+dtXuxmbGqVaXOnaWBA72PWVuYS7FKWqvktRDWooVXBGJFHf/8E7wZwYt54w1vnFZtDAAXkiwmxnbHXNzBgweVJUsWHThwQJlD8V87AIgjmwW0amAr7rBZtISy4gybHbQWMQAQDHHJa+wBBBAVZszwevdZ+xULfbZ8W6tWwsOf/Qr9yy/eMXBXXZVYowWApMUSMICocOiQ1K2bVxBi/fRsKXjWLO9jdlqH7yi382833njx5z1wwDtZxPr5WX+/+ErIGAAgrlgCBhD19u2Tu12INVEuUCA6xgAgvLEEDABxYCd22C3axwAgerAEDAAAEGUIgAAAAFGGAAgAABBlCIAAAABRhgAIAAAQZQiAAAAAUYYACAAAEGUIgAAAAFGGAAgAABBlCIAAAABRhgAIAAAQZQiAAAAAUYYACAAAEGUIgAAAAFGGAAgAABBlCIAAAABRJqU/nxQTE+PeHjx4MKnHAwAAgHjw5TRfbktwADx06JB7W6hQofiMBwAAAAFiuS1LliwX/ZxkMX7ExFOnTmnnzp3KlCmTkiVLlphjBAAAQCKwSGfhL3/+/EqePHnCAyAAAAAiB0UgAAAAUYYACAAAEGUIgAAAAFGGAAgAABBlCIAAAABRhgAIAAAQZQiAAAAAii7/B5CuBx+BZWSAAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from openfe.utils.atommapping_network_plotting import plot_atommapping_network\n",
    "plot_atommapping_network(ligand_network)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30b276b6",
   "metadata": {},
   "source": [
    "You can output the ligand network to the same `graphml` format as we saw in the CLI tutorial with the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2263838f",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"ligand_network.graphml\", mode='w') as f:\n",
    "    f.write(ligand_network.to_graphml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "056924a3",
   "metadata": {},
   "source": [
    "## 4. Creating a single `Transformation`\n",
    "\n",
    "The `LigandNetwork` only knows about the small molecules and the alchemical connections between them. It doesn't know anything about environment (e.g., solvent) or about the `Protocol` that will be used during the simulation.\n",
    "\n",
    "That information in included in a `Transformation`. Each of these transformations corresponds to a thermodynamic cycle of one ligand transformation, including the complex and the solvent leg.\n",
    "\n",
    "In practice, this will be done for each edge of the `LigandNetwork` in a loop, but for illustrative purposes we'll dive into the details of creating a single transformation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0cb1329",
   "metadata": {},
   "source": [
    "### Creating `ChemicalSystem`s\n",
    "\n",
    "OpenFE describes complex molecular systems as being composed of `Component`s. For example, we have `SmallMoleculeComponent` for each small molecule in the `LigandNetwork`. We'll create a `SolventComponent` to describe the solvent, and a `ProteinComponent` for the protein.\n",
    "\n",
    "The `Component`s are joined in a `ChemicalSystem` to define the entire system.\n",
    "\n",
    "In state A of the `SepTopProtocol`, the ligand A is fully interacting in the complex and the ChemicalSystem contains the ligand A, the protein, and the solvent. Ligand B is fully decoupled in this state. In the other endstate, state B, ligand B is fully interacting in the complex while ligand A is decoupled. Therefore, the ChemicalSystem in state B only contains the ligand B, protein and the solvent.\n",
    "\n",
    "Note that for SepTop simulations, we are not separately defining the end states of the solvent leg, but the `Protocol` creates that based on the complex states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9d2fbc22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# defaults are water with NaCl at 0.15 M\n",
    "solvent = openfe.SolventComponent()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3f1706ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "protein = openfe.ProteinComponent.from_pdb_file(\"./tyk2_protein.pdb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "710285ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "systemA = openfe.ChemicalSystem({\n",
    "    'ligand': charged_ligands[0],\n",
    "    'solvent': solvent,\n",
    "    'protein': protein\n",
    "})\n",
    "systemB = openfe.ChemicalSystem({\n",
    "    'ligand': charged_ligands[1],\n",
    "    'solvent': solvent,\n",
    "    'protein': protein    \n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "340d1a6e",
   "metadata": {},
   "source": [
    "### Creating a `Protocol`\n",
    "\n",
    "The actual simulation is performed by a `Protocol`. We'll use an OpenMM-based separated topologies relative free energy `Protocol`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3f394a0d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from openfe.protocols.openmm_septop import SepTopProtocol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56658a3a",
   "metadata": {},
   "source": [
    "There are various different parameters which can be set to determine how the SepTop simulation will take place.\n",
    "\n",
    "The easiest way to customize protocol settings is to start with the default settings, and modify them. Many settings carry units with them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c2b3de6d-1f11-432d-82f3-9deff16c78dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from openff.units import unit\n",
    "settings = SepTopProtocol.default_settings()\n",
    "\n",
    "# Run only a single repeat\n",
    "settings.protocol_repeats = 1\n",
    "# Change the min and max distance between protein and ligand atoms for Boresch restraints to avoid periodicity issues\n",
    "settings.complex_restraint_settings.host_min_distance = 0.5 * unit.nanometer\n",
    "settings.complex_restraint_settings.host_max_distance = 1.5 * unit.nanometer\n",
    "# Set the equilibration time to 2 ns (which is also the default)\n",
    "settings.solvent_simulation_settings.equilibration_length = 2000 * unit.picosecond\n",
    "settings.complex_simulation_settings.equilibration_length = 2000 * unit.picosecond"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7adf42d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "protocol = SepTopProtocol(settings)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "318ff872",
   "metadata": {},
   "source": [
    "### Creating the `Transformation`\n",
    "\n",
    "Once we have the mapping, the two `ChemicalSystem`s, and the `Protocol`, creating the `Transformation` is easy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "44ba94ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "transformation = openfe.Transformation(\n",
    "    systemA,\n",
    "    systemB,\n",
    "    protocol,\n",
    "    mapping=None,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4283dfe4",
   "metadata": {},
   "source": [
    "To summarize, this `Transformation` contains:\n",
    "- chemical models of both sides of the alchemical transformation in `systemA` and `systemB`\n",
    "- a description of the exact computational algorithm to use to perform the estimate in `Protocol`\n",
    "\n",
    "**The `mapping` is set to `None` since no atoms are mapped in the SepTop protocol.**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e29d1c8",
   "metadata": {},
   "source": [
    "## 5. Creating the `AlchemicalNetwork`\n",
    "\n",
    "The `AlchemicalNetwork` contains all the information needed to run the entire campaign. It consists of a `Transformation` for each edge of the campaign. We'll loop over all the edges to make each transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "66666a80",
   "metadata": {},
   "outputs": [],
   "source": [
    "transformations = []\n",
    "for edge in ligand_network.edges:\n",
    "    # use the solvent and protein created above\n",
    "    sysA_dict = {'ligand': edge.componentA,\n",
    "                 'protein': protein,\n",
    "                 'solvent': solvent}\n",
    "    sysB_dict = {'ligand': edge.componentB,\n",
    "                 'protein': protein,\n",
    "                 'solvent': solvent}\n",
    "    \n",
    "    # we don't have to name objects, but it can make things (like filenames) more convenient\n",
    "    sysA = openfe.ChemicalSystem(sysA_dict, name=f\"{edge.componentA.name}\")\n",
    "    sysB = openfe.ChemicalSystem(sysB_dict, name=f\"{edge.componentB.name}\")\n",
    "    \n",
    "    prefix = \"rbfe_\"  # prefix is only to exactly reproduce CLI\n",
    "    \n",
    "    transformation = openfe.Transformation(\n",
    "        stateA=sysA,\n",
    "        stateB=sysB,\n",
    "        mapping=None,\n",
    "        protocol=protocol,  # use protocol created above\n",
    "        name=f\"{prefix}{sysA.name}_{sysB.name}\"\n",
    "    )\n",
    "    transformations.append(transformation)\n",
    "\n",
    "network = openfe.AlchemicalNetwork(transformations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c61fe36",
   "metadata": {},
   "source": [
    "## 6. Running the SepTop simulations using the OpenFE CLI\n",
    "\n",
    "We'll write out each transformation to disk, so that they can be run independently using the `openfe quickrun` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d6cebd9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "# first we create the directory\n",
    "transformation_dir = pathlib.Path(\"transformations\")\n",
    "transformation_dir.mkdir(exist_ok=True)\n",
    "\n",
    "# then we write out each transformation\n",
    "for transformation in network.edges:\n",
    "    transformation.to_json(transformation_dir / f\"{transformation.name}.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b96b57a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rbfe_lig_ejm_31_lig_ejm_42.json rbfe_lig_ejm_42_lig_ejm_43.json\n",
      "rbfe_lig_ejm_31_lig_ejm_46.json rbfe_lig_ejm_46_lig_jmc_23.json\n",
      "rbfe_lig_ejm_31_lig_ejm_47.json rbfe_lig_ejm_46_lig_jmc_27.json\n",
      "rbfe_lig_ejm_31_lig_ejm_48.json rbfe_lig_ejm_46_lig_jmc_28.json\n",
      "rbfe_lig_ejm_31_lig_ejm_50.json\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/atravitz/micromamba/envs/openfe-notebooks/lib/python3.13/pty.py:95: DeprecationWarning: This process (pid=80258) is multi-threaded, use of forkpty() may lead to deadlocks in the child.\n",
      "  pid, fd = os.forkpty()\n"
     ]
    }
   ],
   "source": [
    "!ls transformations/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c30e8ae2",
   "metadata": {},
   "source": [
    "Each of these individual `.json` files contains a `Transformation`, which contains all the information to run the calculation.  These could be farmed out as individual jobs on a HPC cluster."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32178871-7539-4042-9861-4d37b5ba497b",
   "metadata": {},
   "source": [
    "You can run the SepTop simulation from the CLI by using the `openfe quickrun` command. It takes a transformation JSON as input, and the flags -o to give the final output JSON file and -d for the directory where simulation results should be stored. For example,\n",
    "\n",
    "`openfe quickrun path/to/transformation.json -o results.json -d working-directory`\n",
    "\n",
    "where path/to/transformation.json is the path to one of the files created above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0968657-7b60-4b10-ab3b-5303cf08caa1",
   "metadata": {},
   "source": [
    "## 7. Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c7cd90f-c0e5-4191-b959-70918f50297e",
   "metadata": {},
   "source": [
    "Finally now that we've run our simulations, let's go ahead and gather the free energies for both phases.\n",
    "If you ran the simulation using the CLI (i.e. by calling openfe quickrun ) you will end up with a JSON output file for each transformation. To get your simulation results you can load them back into Python in the following manner:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "de7d7920-b171-42fe-9b82-5ff5694d26e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment the following lines if you've run the transformation\n",
    "\n",
    "# import json\n",
    "\n",
    "# outfile = \"results/rbfe_lig_ejm_31_lig_ejm_42.json\"\n",
    "# with open(outfile) as stream:\n",
    "#     results = json.load(stream)\n",
    "#     estimate = results['estimate']\n",
    "#     uncertainty = results['uncertainty']\n",
    "# estimate"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}