01-annihilate-Elec-vacuum.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c8347bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from pymbar import MBAR, timeseries\n",
    "\n",
    "import simtk.unit as simtk_unit\n",
    "from simtk.openmm.app import (\n",
    "    AmberPrmtopFile,\n",
    "    AmberInpcrdFile,\n",
    "    Simulation,\n",
    "    HBonds,\n",
    "    NoCutoff,\n",
    "    DCDReporter,\n",
    "    StateDataReporter,\n",
    "    PDBFile,\n",
    ")\n",
    "from simtk.openmm import (\n",
    "    XmlSerializer,\n",
    "    LangevinIntegrator,\n",
    "    MonteCarloBarostat,\n",
    "    CustomNonbondedForce,\n",
    "    CustomBondForce,\n",
    "    NonbondedForce,\n",
    "    Platform,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b41cb42c",
   "metadata": {},
   "source": [
    "### Functions to add `lambda_electrostatics`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e58dab1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_lambda_elec(system, molecule, lambda_value=1.0):\n",
    "    \"\"\"\n",
    "    * Adds lambda scaling to both inter- and intramolecular electrostatic\n",
    "      interactions in `NonbondedForce` for the molecule.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    system: openmm.System\n",
    "        The OpenMM system containing the nonbonded force.\n",
    "    molecule: list\n",
    "        A list atom indices (starts from 0) for the molecule of interest.\n",
    "    lambda_value: float\n",
    "        The starting lambda value for the Elec scaling.\n",
    "    \"\"\"\n",
    "    nonbonded = [\n",
    "        force for force in system.getForces() if isinstance(force, NonbondedForce)\n",
    "    ][0]\n",
    "\n",
    "    nonbonded.addGlobalParameter(\"lambda_electrostatics\", lambda_value)\n",
    "\n",
    "    for atom_i in molecule:\n",
    "        charge, sigma, epsilon = nonbonded.getParticleParameters(atom_i)\n",
    "        nonbonded.addParticleParameterOffset(\n",
    "            \"lambda_electrostatics\", atom_i, charge, 0.0, 0.0\n",
    "        )\n",
    "\n",
    "    for exception_i in range(nonbonded.getNumExceptions()):\n",
    "        atom_i, atom_j, charge, sigma, epsilon = nonbonded.getExceptionParameters(\n",
    "            exception_i\n",
    "        )\n",
    "        if atom_i in molecule and atom_j in molecule:\n",
    "            nonbonded.addExceptionParameterOffset(\n",
    "                \"lambda_electrostatics\", exception_i, charge, 0.0, 0.0\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d0d682e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_custom_LJ(system, molecule, force_group=5):\n",
    "    \"\"\"\n",
    "    * move the Intramolecular LJ to `CustomBondForce`\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    system: openmm.System\n",
    "        The OpenMM system containing the nonbonded force.\n",
    "    molecule: list\n",
    "        A list of atom indices (starts from 0) for the molecule of interest.\n",
    "    force_group: int\n",
    "        Force group number for the intermolecular LJ\n",
    "    \"\"\"\n",
    "    nonbonded = [\n",
    "        force for force in system.getForces() if isinstance(force, NonbondedForce)\n",
    "    ][0]\n",
    "\n",
    "    # Intramolecular LJ\n",
    "    intramol_LJ = CustomBondForce(\"4*epsilon*((sigma/r)^12 - (sigma/r)^6);\")\n",
    "    intramol_LJ.addPerBondParameter(\"sigma\")\n",
    "    intramol_LJ.addPerBondParameter(\"epsilon\")\n",
    "    intramol_LJ.setForceGroup(force_group)\n",
    "\n",
    "    for exception_i in range(nonbonded.getNumExceptions()):\n",
    "        atom_i, atom_j, charge, sigma, epsilon = nonbonded.getExceptionParameters(\n",
    "            exception_i\n",
    "        )\n",
    "        if atom_i in molecule and atom_j in molecule:\n",
    "            intramol_LJ.addBond(atom_i, atom_j, [sigma, epsilon])\n",
    "\n",
    "    system.addForce(intramol_LJ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "67fb56e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def turn_off_nonbonded(system, molecule):\n",
    "    \"\"\"\n",
    "    * Turn off the LJ parameters in `NonbondedForce` for the molecule\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    system: openmm.System\n",
    "        The OpenMM system containing the nonbonded force.\n",
    "    molecule: list\n",
    "        A list of atom indices (starts from 0) for the molecule of interest.\n",
    "    \"\"\"\n",
    "    nonbonded = [force for force in system.getForces() if isinstance(force, NonbondedForce)][0]\n",
    "    \n",
    "    # Turn off Elec and LJ parameters\n",
    "    for atom_i in molecule:\n",
    "        charge, sigma, epsilon = nonbonded.getParticleParameters(atom_i)\n",
    "        nonbonded.setParticleParameters(atom_i, 0.0, sigma, 0.0)\n",
    "        \n",
    "    # Turn off Elec and LJ Exceptions\n",
    "    for exception_i in range(nonbonded.getNumExceptions()):\n",
    "        atom_i, atom_j, charge, sigma, epsilon = nonbonded.getExceptionParameters(exception_i)\n",
    "        if atom_i in molecule and atom_j in molecule:\n",
    "            nonbonded.setExceptionParameters(exception_i, atom_i, atom_j, 0.0, sigma, 0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b84f3141",
   "metadata": {},
   "source": [
    "# Initial Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2386fadc",
   "metadata": {},
   "outputs": [],
   "source": [
    "build_folder = \"build_files\"\n",
    "work_folder = \"elec-vac\"\n",
    "os.makedirs(work_folder, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "32f199c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "platform = Platform.getPlatformByName(\"CPU\")\n",
    "#platform.setPropertyDefaultValue(\"DeviceIndex\", \"0\")\n",
    "#platform.setPropertyDefaultValue(\"Precision\", \"mixed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cc04039f",
   "metadata": {},
   "outputs": [],
   "source": [
    "prmtop = AmberPrmtopFile(f\"{build_folder}/ethanol_vacuum.prmtop\")\n",
    "inpcrd = AmberInpcrdFile(f\"{build_folder}/ethanol_vacuum.rst7\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "58f537de",
   "metadata": {},
   "outputs": [],
   "source": [
    "ethanol = [atom.index for atom in prmtop.topology.atoms() if atom.residue.name == \"MOL\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "88efc5fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = prmtop.createSystem(\n",
    "    nonbondedMethod=NoCutoff,\n",
    "    constraints=HBonds,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "160cdf96",
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 298.15 * simtk_unit.kelvin\n",
    "friction = 1.0 / simtk_unit.picosecond\n",
    "timestep = 2.0 * simtk_unit.femtosecond\n",
    "kT = temperature * simtk_unit.BOLTZMANN_CONSTANT_kB * simtk_unit.AVOGADRO_CONSTANT_NA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "76d2879d",
   "metadata": {},
   "outputs": [],
   "source": [
    "add_lambda_elec(system, ethanol, lambda_value=1.0)\n",
    "create_custom_LJ(system, ethanol)\n",
    "turn_off_nonbonded(system, ethanol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e1f5d122",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(f\"{work_folder}/system.xml\", \"w\") as f:\n",
    "    f.write(XmlSerializer.serialize(system))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aeff92c",
   "metadata": {},
   "source": [
    "# MD Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "649101a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrator = LangevinIntegrator(temperature, friction, timestep)\n",
    "\n",
    "simulation = Simulation(\n",
    "    prmtop.topology,\n",
    "    system,\n",
    "    integrator,\n",
    "    platform,\n",
    ")\n",
    "simulation.context.setPositions(inpcrd.positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92dd5c4f",
   "metadata": {},
   "source": [
    "## 01. Minimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "892b85e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation.minimizeEnergy(\n",
    "    tolerance=1.0 * simtk_unit.kilojoules_per_mole, maxIterations=5000\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66834eaf",
   "metadata": {},
   "source": [
    "## 02. Equilibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7ce2f524",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_freq = 5000\n",
    "equil_steps = 5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "4bd9c6f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "dcd_reporter = DCDReporter(f\"{work_folder}/equilibration.dcd\", output_freq)\n",
    "state_reporter = StateDataReporter(\n",
    "    f\"{work_folder}/equilibration.log\",\n",
    "    output_freq,\n",
    "    step=True,\n",
    "    kineticEnergy=True,\n",
    "    potentialEnergy=True,\n",
    "    totalEnergy=True,\n",
    "    temperature=True,\n",
    "    volume=True,\n",
    "    speed=True,\n",
    "    separator=\",\",\n",
    ")\n",
    "simulation.reporters.append(dcd_reporter)\n",
    "simulation.reporters.append(state_reporter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2ced4d96",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation.step(equil_steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "41d1e692",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation.saveState(f\"{work_folder}/equilibration.xml\")\n",
    "with open(f\"{work_folder}/equilibration.pdb\", \"w\") as f:\n",
    "    PDBFile.writeModel(\n",
    "        prmtop.topology,\n",
    "        simulation.context.getState(getPositions=True).getPositions(),\n",
    "        f,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e0568f7",
   "metadata": {},
   "source": [
    "## 03. Alchemical Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d2283f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_values = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]\n",
    "n_lambda = len(lambda_values)\n",
    "\n",
    "output_freq = 500\n",
    "equil_steps = 5000\n",
    "prod_steps = 50000\n",
    "\n",
    "total_iterations = int(prod_steps/output_freq)\n",
    "\n",
    "u_kln = np.zeros((n_lambda, n_lambda, total_iterations))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "83fc039f",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrator = LangevinIntegrator(temperature, friction, timestep)\n",
    "\n",
    "simulation = Simulation(\n",
    "    prmtop.topology,\n",
    "    system,\n",
    "    integrator,\n",
    "    platform,\n",
    ")\n",
    "simulation.loadState(f\"{work_folder}/equilibration.xml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9a27db18",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running simulation for lambda: 0.0 ...\n",
      "Running simulation for lambda: 0.1 ...\n",
      "Running simulation for lambda: 0.2 ...\n",
      "Running simulation for lambda: 0.3 ...\n",
      "Running simulation for lambda: 0.4 ...\n",
      "Running simulation for lambda: 0.5 ...\n",
      "Running simulation for lambda: 0.6 ...\n",
      "Running simulation for lambda: 0.7 ...\n",
      "Running simulation for lambda: 0.8 ...\n",
      "Running simulation for lambda: 0.9 ...\n",
      "Running simulation for lambda: 1.0 ...\n"
     ]
    }
   ],
   "source": [
    "for k, lb1 in enumerate(lambda_values):\n",
    "    print(f\"Running simulation for lambda: {lb1} ...\")\n",
    "    \n",
    "    simulation.context.setParameter(\"lambda_electrostatics\", lb1)\n",
    "    \n",
    "    # Short equilibration\n",
    "    simulation.step(equil_steps)\n",
    "    \n",
    "    # Add reporters\n",
    "    dcd_reporter = DCDReporter(f\"{work_folder}/production-{k}.dcd\", output_freq)\n",
    "    state_reporter = StateDataReporter(\n",
    "        f\"{work_folder}/production-{k}.log\",\n",
    "        output_freq,\n",
    "        step=True,\n",
    "        kineticEnergy=True,\n",
    "        potentialEnergy=True,\n",
    "        totalEnergy=True,\n",
    "        temperature=True,\n",
    "        volume=True,\n",
    "        speed=True,\n",
    "        separator=\",\",\n",
    "    )\n",
    "    simulation.reporters.append(dcd_reporter)\n",
    "    simulation.reporters.append(state_reporter)\n",
    "    \n",
    "    # Run production\n",
    "    for n in range(total_iterations):\n",
    "        simulation.context.setParameter(\"lambda_electrostatics\", lb1)\n",
    "        simulation.step(output_freq)\n",
    "        \n",
    "        for l, lb2 in enumerate(lambda_values):\n",
    "            simulation.context.setParameter(\"lambda_electrostatics\", lb2)\n",
    "            u_kln[k,l,n] = simulation.context.getState(getEnergy=True).getPotentialEnergy() / kT\n",
    "    \n",
    "    simulation.saveState(f\"{work_folder}/production-{k}.xml\")\n",
    "    \n",
    "    simulation.reporters.pop(-1)\n",
    "    simulation.reporters.pop(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "97172608",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(f\"{work_folder}/reduced_elec_vacuum_potential.npy\", u_kln)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b44fce49",
   "metadata": {},
   "source": [
    "# Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f2969768",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subsample data to extract uncorrelated equilibrium timeseries\n",
    "N_k = np.zeros([n_lambda], np.int32)\n",
    "\n",
    "for k in range(n_lambda):\n",
    "    [nequil, g, Neff_max] = timeseries.detectEquilibration(u_kln[k,k,:])\n",
    "    indices = timeseries.subsampleCorrelatedData(u_kln[k,k,:], g=g)\n",
    "    N_k[k] = len(indices)\n",
    "    u_kln[k,:,0:N_k[k]] = u_kln[k,:,indices].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "210cf203",
   "metadata": {},
   "outputs": [],
   "source": [
    "mbar = MBAR(u_kln, N_k)\n",
    "[DeltaF_ij, dDeltaF_ij, Theta_ij] = mbar.getFreeEnergyDifferences(return_theta=True)\n",
    "ODeltaF_ij = mbar.computeOverlap()['matrix']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "d6c9dcf1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 = 0.00 +- 0.00 kcal/mol\n",
      "0.1 = -0.69 +- 0.00 kcal/mol\n",
      "0.2 = -1.29 +- 0.01 kcal/mol\n",
      "0.3 = -1.80 +- 0.01 kcal/mol\n",
      "0.4 = -2.20 +- 0.01 kcal/mol\n",
      "0.5 = -2.49 +- 0.01 kcal/mol\n",
      "0.6 = -2.68 +- 0.01 kcal/mol\n",
      "0.7 = -2.76 +- 0.01 kcal/mol\n",
      "0.8 = -2.73 +- 0.01 kcal/mol\n",
      "0.9 = -2.59 +- 0.01 kcal/mol\n",
      "1.0 = -2.35 +- 0.01 kcal/mol\n"
     ]
    }
   ],
   "source": [
    "convert_kcalmol = kT.value_in_unit(simtk_unit.kilocalorie_per_mole)\n",
    "for k, lb in enumerate(lambda_values):\n",
    "    print(f\"{lb} = {DeltaF_ij[0,k]*convert_kcalmol:.2f} +- {dDeltaF_ij[0,k]*convert_kcalmol:.2f} kcal/mol\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "db50b1f6",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'overlap fraction')"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdQAAAGMCAYAAABj+lP2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwxklEQVR4nO3debwcVZ3+8c+TIJEl4ACySIwi44LbAMbfKLvjRJnMiBvIIDNjgBkU3BAYNcpA0NG4AgIDqEQQt4FRR2UUBQUdBGQJsjlRtrDvIiQsSci9398fpzp0Ol19q+pW9+3b/bx51au7q05VnXM7nG+fU6dOKSIwMzOz8Zky0RkwMzMbBA6oZmZmNXBANTMzq4EDqpmZWQ0cUM3MzGrggGpmZlYDB1QzM7MaOKCamZnVwAHVzMysBg6oNpQkzZUUks6c6LwMA0nzs7/3/InOi1m3OKBaKZLeJOlsSXdIWi7pYUmLJH1S0hYTnb9BImmPLAg1lleOkf63TWnPrDEP8yXtUcfxzAaZA6oVImljST8FfgS8A3gmcD3wILADcBRwk6R3TFwuB94/5G2Q9DJg+y6ccw/gmOx1PB4C/pC9mg0kB1Qbk6R1gQuANwK3AXOALSLi1RHxYuB5wH8C04HvSHr7ROV1QN0DLAPeKSnv/9l/zF7/0JsslRMRJ0fESyLi5InOi1m3OKBaEccCrwbuBXaNiPOi6TFFEXFnROwHnEH6N3W6u39r9STwfWBr4HWtGyUJeCfwOPDfvc2amTU4oFpHkp4FvC/7eGRE3NUh+QdJXXrN+yDpu9l1vSM7nOdNWZqr22ybIelESTdKelLSI5IukrR3zrF+mR1rD0nbZ+e/X9KopLkFyjxb0smSrs2uES+XdIukUyXNzNnnzOyccyW9KLvO/ECW399KOnCs847hm9nrP7bZtgfwXFIwfbyuckkKUncvwDEt13PPbEp3W7bu+ZJeJ+k8SQ81voMszVqDkiRtJunebP1a5ZK0nqTF2faPdf7zmE08B1QbyxxgQ1Kg/K9OCSNiGfCt7GPztdRvZ6/7ddi9se07zSsl7Q7cALwfmAHcBCwlBZH/kvSFDsfcDfgNqav6TmBJp/w3OQ84FNgSuD075xbAe4CrJb20w74vBK4A3pyd837Stc2Fkk4seP52LgTuBt4maf2WbY1rq98Y4xhly3UJqQxkr5c0LTe2Of5+wM+BvwRuBTr9+CIiHgIOyj6eLOl5LUk+D7wEuBT47BhlM5t4EeHFS+4CnAwE8MOC6d+WpQ9g02zdNOCRbN2L2uyzPvAYMAo8t2n9c4A/ZuvnAdOatu1EqrAD+LuW4/0yW78K+DKwftO29bLXuVmaM9vk52DgOS3r1gM+lu1zUZt9zsy2PUUKfs9u2rY3sDLb/rcl/vZ7ZPvcnH3+XPZ5v6Y0zwQeJV1nnUoaHFZnueZn2+Z3yOdtTX/vo4F1svVqfGedjgOclm37FTAlW/eG7HtfBrxgov8/8OKlyOIWqo1l6+z1loLpm9NtDRARK0jXAKF9K3UvYAPg1xFxZ9P6I4BNgBMiYkF2HLJjXkpqWQF8KCcvNwCHRMQTTfs9OVYBIuIrEXFPy7onI+LTwK+BPSRt3X5vRoB3RsSDTft+Fzgp+/jhsc7fQaMF2tw9+mZgI+A7ETHSaedxlquIn0TEJyJiVXbsaP7OOjiC1FreDThS0iak6/ECDouIW8eRJ7OeWWeiM2B9b3r2mnttrkVzuulN778NHEAKqMe27LNfU5pmb8teT885109JLb+dJK3TqMibfDMiRgvluoWkWaSW5UuBjUmtP0hdugCvJHXBtvp+RNzXZv0pwOHAzpI2iIiif8/VIuJ6SdcBsyVtHhEPULy7FxhXuYo4q8pOEfF4dg31EuCTpMsMzwF+FBELK+bFrOccUG0sy7LXDQqmb063rOn9hcB9wIsl7RARv4XVg572JHUXfreRWNKGwPOzj19JA1lzPRPYlHS9stnignleLRsxezLpWmMnm+SszzvnrcAKUvf3tsB1ZfOW+Sap63c/Sd8i/e1+FxHXdNqphnIVUfrv3RARl0v6NPBvwO7AA8C/jCMvZj3nLl8bS6O1sm3B9M3pVrd0spbi2dnH5m7ftwPrAudHGqTSsHHT+507LOtmadZrk5fSrUBSd+qh2b6Hklpu60eEIkI8PejqGTn7P9BuZUQEaRIMWLPlXta3SNcW/wH4e9KP4iKt0/GWq4gqf+9mFza9/5+sBW42aTig2lguzV53klSkR2O37PWmiPhjy7bGCN6/19NNzraje0mDlBrWbVT8HZbbihSmgP2z1yMi4tSIuLnluutzx9j/2e1WZuVtbFvWLk0R2TXQC4FZwL+Sguu3Ou6UjLdcXZX1SDS6d0eBuZJ2ncAsmZXmgGpj+Qmp5bEZsE+nhJKm83TFfXbr9oi4nDRo6bnALpK2JI1kfRL4QUvaxshVgJdVzn15z89eL23dIOkZwHZj7J+3fRtSd+8oxQd45WnckzoT+FV0vje44fnZa9lyRc76un0JeAFpRq7DSHXTWdm/KbNJwQHVOoqIR4D/yD5+UdKMDsm/RAq8jzbt06rREt0P2Jc0KObciHisTdrGyODDSmR5vBqttnYzPR1ATgu0ydtzZolqXLu8pMqApBbfA84HfgEUvbe1arka+7XrUq+FpL2AA0m3Vh1IutZ7PulHwJe6dV6zujmgWhFHA1cDWwH/K2nPpi7bxkxGjVG8ARycM9IVnu6e3IenR6i2ju5t+CzwMPAuScdlA5hWk7SJpAMlHVWlUDl+nb3+u6TVQUbSnqSJBpaPsf9U4FuSNmva962kiSnIjjEuEfFYRLwxIv46In5QcLeq5WrcslK0y78USZsDX80+HhoRd2XXmw8gffcHSHpL3ec16wYHVBtTdi/h60mz4GxDmnHnfklXSvo9cAepxfkYsH9EnNPhWL8HriG1ZGeRWiXn5aS9i3SP6kOke00fkHSdpN9IuiVbvxB4eQ3FbPgcqSL/S+D2bNrAJVkeF5Fah518nlSuOyVdle37fdLgqVMi4twa81pG1XKdD/wJ2AW4Q9KvlaZ2/GhN+foqsDlwTkSsvo6eXSs+JPv4lSzwmvU1B1QrJCIeiYjZwFtIt7esIN2zuCVwLfBp4IXNlWIHzS3S70XEyg7nvYR0z+SngP8jBfRXkq5F/pTUlfrBsuXpcL47gNeSguBK0tR3y0lz2jZu7+nkRuD/AeeSrnFuRfr7HEzT/Ma9VrVcEbGUNGvReaRrwK8l3dbykvHmSdJBpB9M9/J08Gw+9zmkfyvPJv9eZLO+odS7YmbjkU0W/y7ggIg4c2JzY2YTwS1UMzOzGjigmpmZ1cAB1czMrAZ9GVAlHSppSfYA5EVjzZgiafcs3XJJt0p6T6f0ZmbWPyTtJulcSfdkD5R/S4F9+q7e77uAKmlf4ATSqM4dgIuB8yTNzEm/DWk2n4uz9J8GTpT09p5k2AyIiLnZFIhnTnRezCahDUij4QuNhO/Xer/vRvlKuhy4OiIOaVq3GPhBRMxrk/6zwF4RsV3TutOAv4iI1/Yiz2ZmVg9JAby106Ql/Vrv99Xj2yStC7wK+EzLpvOBnXJ2e222vdnPgIMkPSMinmpznmmke+qabUK68d3MbLKYDtwTXWoZSXomTz/RqQ4rCj50fiyl6/1e6KuASpo9ZyprP9fyftIEAu1smZN+nex497bZZx7phnYzs8luBtUfCp9L0jO33Hzqk/c9MFL1EI8BG7asOxaYP558ZarU+13XbwG1ofXXltqsGyt9u/UNC4Djmj5PB+7a/QWHsM6U1oZrFzxS+eld1Tz7z3p2Kj3a27LF9KLPPR8/LRvvnPblxAbr9+5kPf7e2Ki1nu2iPz3as1NpvWf27FyrRlfyy3sWwjgeBziGde97YIQli57HRtPLDbdZumyUbV51+4akYN+cvzpapw1l6/2u67eA+hAwwtqt0c1Z+9dIw3056VcBrc/jBFbPTbv6i23M877OlGmsM7UHAXVKnf+mCuhFmTKakjuLYFdET8s21qyD9epl2Qb53yRT6uyx7Ey9+EHeYxtsmJYyRp4Oacuy6SvrVrre74W+GuWbzem6CJjdsmk2bZ7jmLmsTfo3AFdNVD+6mdmgGCUqLV3Wl/V+XwXUzHHAP2eP5dpO0vGkScZPA5C0QNJZTelPA56XPd5rO0kHAgcBX+h5zs3MBsxoxf/KkLShpO0lbZ+t2ib7PDPbPinq/X7r8iUizpa0KekZnFsBNwBzIuL2LMlWpADbSL9E0hzgeOC9wD3AByJirMdsmZlZf5gFXNT0uTHG5evAXCZJvd93ARUgIk4BTsnZNrfNul8BO3Y5W2ZmQ2ckgpGSd+WUTR8Rv+TpQUXtts9ts67v6v2+DKhmZtYfqlwT7cE11L7kgGpmZrlGCUYcUAtxQDUzs1xuoRbngGpmZrl6cQ11UPTjbTNmZmaTjluoZmaWazRbyu4zjBxQzcws10iFQUll0w8KB1QzM8s1EmvMzVt4n2HkgGpmZrnc5VucA6qZmeUaRYzkT2KUu88wckA1M7Nco5GWsvsMI982Y2ZmVgO3UM3MLNdIhS7fsukHhQOqmZnlckAtzgHVzMxyjYYYjZKDkkqmHxQOqGZmlsst1OIcUM3MLNcIUxgpOX51pEt56XcOqGZmlisqdPnGkHb5+rYZMzOzGriFamZmuXwNtTgH1GYRaRk0g1imiaDhrCS6wn/LSWMkpjASJa+hDmmV44BqZma5RhGjJa8OjvrxbWZmZmtyl29xDqhmZparWpfvcLZQPcrXzMysBm6hmplZrnQN1c9DLcIB1czMco1WmCnJg5LMzMxa+BpqcQ6oZmaWa5Qpvm2mIAdUMzPLNRJipOTcvGXTDwoHVDMzy1XtaTPD2UL1bTNmZmY1cAvVzMxyjcYURksOShr1oCQzM7M1ucu3OAdUMzPLNUr5QUaj3clK33NANTOzXNVumxnO4TkOqGZmlqvaxA7DGVCHs9RmZmY1cwvVzMxyeXL84hxQzcwsl7t8i3NANTOzXNVum3FANTMzW8NoiNGyt814Ll8zM7M1VXseqluoZmZma6g29eBwBtThLLWZmVnN3EI1M7NcI4iRkrfBlE0/KBxQzcwsl7t8i3NANTOzXCOUb3GOdCcrfc8BtVlEWrp+nmF9FoOVMmU4u81qpx62lnr6nfXmXG6hFteXpZZ0qKQlkpZLWiRp1w5p3ybpAkkPSloq6TJJb+xlfs3MBlVjpqSyS1ll6v0s/f6SrpX0hKR7JZ0hadPKBa1B3wVUSfsCJwCfAnYALgbOkzQzZ5fdgAuAOcCrgIuAcyXt0P3cmpnZeJWt9yXtApwFLAReBuwDvBo4vRf5zdN3ARU4HFgYEadHxOKIOAy4EzikXeKIOCwiPhcRV0bETRHxMeAm4E29y7KZ2WCKbHL8MkuU744uVe8DrwFui4gTI2JJRPwa+DIwq2o569BXAVXSuqRW5vktm84Hdip4jCnAdODhDmmmSdqosWTpzcysxTi7fKc317WSprUev2K9fykwQ9IcJVsAewM/rqfU1fRVQAU2A6YC97esvx/YsuAxjgA2AM7pkGYe8GjTcle5bJqZDYfGXL5ll8xdrFnXzmtzitL1fkRcCuwPnA2sBO4DHgHeP77Sjk+/BdSG1qG2arNuLZL2A+YD+0bEAx2SLgA2blpmVMummdlgazxtpuySmcGade2CDqcqXO9LeilwIvAJUut2T2Ab4LTKBa1Bv9028xDpFqbWXyWbs/avlzVkF7UXAvtExM87pY2IFcCKpn0rZdbMbNCN82kzyyJi6RjJq9T784BLIuLz2efrJD0OXCzpqIi4t1SGa9JXLdSIWAksAma3bJpN6jNvK2uZngm8MyImtA/dzGyQjDKl0lJUxXp/faD1hv7GfBIT1kLqtxYqwHHANyRdBVwGHAzMJGvKS1oAbB0R/5R93o80fPqDwG8kNX7lPBkRj/Y682ZmVlqpeh84F/iqpEOAnwFbkW67uSIi7ulx3lfru4AaEWdnN+ceTfoj3QDMiYjbsyRbkf7QDe8mleM/sqXh68DcrmfYzGyAjYQYKdnlWzZ92Xo/Is6UNB14H/BF0oCkC4GPlDpxzfouoAJExCnAKTnb5rZ83qMHWTIzG0rjvIZaWJl6P1t3EnBS6RN1UV8GVDMz6w9RYS7fGNK5fB1Qzcwsl5+HWpwDqpmZ5RqN8l24oz14aFc/Gs52uZmZWc3cQjUzs1x+HmpxDqhmZpar8QSZsvsMIwdUMzPL1Yv7UAeFA6qZmeVyl29xDqhmZpZrlAoTO7jL18zMbE1R4RpqDGlAHc52uZmZWc3cQjUzs1y9mst3EDigmplZLg9KKs4B1czMcrmFWpwDqpmZ5fLEDsU5oDZRBOrBrM49nzc6hnSmaitM6m0FGFOGs8KdjNxCLc4B1czMcjmgFjecV47NzMxq5haqmZnlcgu1OAdUMzPL5YBanAOqmZnlCsqP2h3WYZAOqGZmlsst1OIcUM3MLJcDanEe5WtmZlYDt1DNzCyXW6jFOaCamVkuB9TiHFDNzCxXhIiSAbJs+kHhgGpmZrk8OX5xDqhmZpbLXb7FOaCamVkud/kW59tmzMzMauAWqpmZ5RrULl9JU4A/BzanpXEZEf9b5ZgOqGZmlmsQu3wlvQb4NvA8WGsEVQBTqxzXAdXMzHJFhRZqvwdU4DTgKuBvgXupaT5/B1QzM8sVQJQMN5PgaTMvBPaOiJvrPKgHJZmZWa7Gfahllz53Oen6aa3cQjUzs2FzEvBFSVsC1wNPNW+MiOuqHNQB1czMcg3ioCTge9nr15rWBWmAkgclmZlZ/UZDaPBum9mmGwd1QDUzs1wRFQYl9fmopIi4vRvHdUA1M7NcA9rli6RtgcOA7UjdvIuBL0XELVWP6VG+ZmaWqxFQyy79TNIbgf8D/h9wHXAD8JfA7yTNrnpct1DNzCzXgF5D/QxwfER8tHmlpM8AnwUuqHJQB9RmVS4W2JpGRyc6B2ZmY9kOeEeb9V8jdQNX4i5fMzPL1WhnlF363IPA9m3Wbw88UPWgbqGamVmuFCDLDkrqUmbq81XgK5JeAFxKGpS0C/AR4ItVD+qAamZmuQZ0lO8ngWXAEcCCbN09wHzgxKoHdUA1M7NcQfnJ7vu9gRoRARwPHC9perZu2XiP64BqZma5BrSFulodgbTBAdXMzAaepKuB10fEnyT9lg4N6YjYsco5+nKUr6RDJS2RtFzSIkm7FtxvZ0mrJF3T5SyamQ2HqLiUVLbelzRN0qck3S5phaRbJB3YYZcfAiua3ndaKum7FqqkfYETgEOBS4B3A+dJemlE3NFhv42Bs4BfAFv0IKtmZoOvysxHJdNXrPfPIdX1BwE3A5vTIaZFxLFN7+eXymBB/dhCPRxYGBGnR8TiiDgMuBM4ZIz9vgx8G7isy/kzMxsaPboPtVS9L2lPYHdgTkT8PCJui4grIuLSIieTdKukTdusf5akW0vnPtNXAVXSusCrgPNbNp0P7NRhvwOAbYFj89K0pJ8maaPGAkyvmGUzs4E2zrl8pzfXtZKmtR6/Yr2/F3AV8GFJd0u6UdIXJK1XsFjPp/0zT6cBMwoeYy391uW7GamQ97esvx/Yst0Okl5Impdx14hYJRXqapgHHDOOfJqZDYdQ6S7cpvR3tWw5lnSvZ7PS9T7wAtJEDMuBt2bHOAXYBMi9jippr6aPb5T0aNPnqcDrgSV5+4+l3wJqQ2uHgdqsQ9JUUjfvMRFxY4njLwCOa/o8nbW/eDOzoTfO56HOIE2g0LBircRNu7V8blvvZ6Zk2/aPiEcBJB0OfFfSeyPiyZz9ftB0rq+3bHsKuI002UMl/RZQHwJGWPtXyeas/esFUiCcBewg6eRs3RRAklYBb4iIC1t3iogVNH2xBVu1ZmZWzrKIWDpGmrL1PsC9wN2NYJpZTArCM4Cb2u0UEVMAJC0BXh0RD42Rt1L66hpqRKwEFgGtz6ObTZpvsdVS4BWkCY0by2nAH7L3l3clo2Zmw6LLt81UqPchjQR+jqQNm9a9CBilQG9jRGxTdzCF/muhQuqK/Yakq0gjdg8GZpICJZIWAFtHxD9FxCjpwbCrSXoAWB4RN2BmZuPSo5mSCtf7WfpvA/8GnCHpGNI11M8DX+vQ3buapBOBmyPixJb17wP+PBtlXFrfBdSIODsbznw0sBUpYM6JiNuzJFuR/tBmZtYLXZ6ct2y9HxGPSZoNnEQa7ftH0n2pRxU85dtJI4VbXQp8lIrPRO27gAoQEaeQRmy12zZ3jH3ns/YoMjMzq6BXc/mWrfcj4ves3U1c1KbAo23WLyW1divpq2uoZmbWZ3o09WCP3Qzs2Wb93wCVJ3boyxaqmZn1C2VL2X362nHAyZKeDTTuBHk96ZaZw6oe1AHVzMyGSkR8LZu16eOkwU2Q7kE9JCLOqnpcB1QzM8s3iE8YByLiVODUrJX6ZEQ8Nt5jOqCamVm+AQ2oDRHxYF3HckA1M7N845vLt29J2ht4B+l2nHWbtw3UA8bNzKw/9OjxbT0l6QPAGcADwA7AFaR7WV8AnFf1uA6oZmaWbzBvmzkUODgi3gesBD4XEbOBE4GNqx7UAdXMzIbNTJ6eJ/hJnn4m9jeA/aoe1NdQzYrq9VOJ/BQk6weDeQ31PtJsSbdny2uAa4FtGMdNtG6hmplZLkW1pc9dCLwpe78QOF7SBcDZwH9XPahbqGZmlm8wb5s5mKxBGRGnSXoY2AU4l+wJN1U4oJqZWb4B6/KVtA5phqSvAXcCRMQ5pKfVjIu7fM3MLN+AjfKNiFXAvwJT6z62A6qZmeUbsICa+TmwR90HdZevmZkNm/OABZJeDiwCHm/eGBE/qnJQB1QzM8s3mIOSTs1eD2+zLajYHeyAamZm+QZsUBJARHTlcqevoZqZWa5BuQ9V0sOSNsvef03S9LH2KcsB1czM8g3OoKR1gY2y9+8Cnln3CbrS5ZtF/pcDr4yIL3fjHGZmZiVcBvxA0iLS9IInSnqyXcKIOLDKCcYVUCVNBV4MvKJpeSVp4uEGB1QzM5to/wB8CNiW1IbemJpbqYUDqqRNgVmkgNkIni/h6QezNq5CrwSuJ000fG1tOTUzs54T5a+J9uOQpIi4H/gogKQlwD9GxB/rPEehgCrpEOBLpKHEzX+rAP4XuJIUPK8DFmczUZiZ2WQ3mKN8t+nGcYu2UD+WvZ5FerL5Hdm61wKbAd+NiCvqz56ZmU2owbwPtSuKjvJ9DnBqRBwQEadGxI8jYmfgA6TrpZdI+qKk9bqWUzMz673BGeXbdUUD6ieBi1pXRsTJwMuAn5Iu9l4v6a/qy56ZmU2kQbkPtRcKBdSImB8RP8jZdldEvAnYD9gQuEDSVyVtXF82zcxsQriFWlhtEztExNmkUb9nAQcBv5P05rqOb2ZmVhdJUyXtLenfJB2VvR/XraS1TuwQEY8AB0j6Jun+0+/ThWfOmZlZjwzgoKTsKTM/BLYE/pCtfhHwoKS9IuL6KsftytSDEfEL0kxJX+jG8c3MrDcG9Brq6cDvgBkRsWNE7Ag8l3Tr51eqHrRrT5uJiOXAR7p1fDMz64EBvA8V+AtgVkT8qbEiIv4k6eOkeRUq8eT4ZmaWbzAHJf0B2KLN+s2Bm6se1M9DNTOzXFW6cCdBl+/HSJPjzwd+k617DXA08BFJjafSEBFLix7UAdXMzIbN/2Sv5/B0e7rRT31u0+egxMBaB1QzM8s3gKN8gdd146AOqGZmlq/KqN0+D6gR8atuHNcBtVlEWnpxHjOzyWAwW6gASFqfNB/9us3rI+K6KsdzQDUzs3wDGFAlPRs4A/ibnCSVJiTybTNmZpZrQCd2OAH4M9LI3ieBPYF3ATcBe1U9qFuoZmY2bP4KeHNEXClpFLg9Ii6QtBSYB/y4ykHdQjUzs2GzAfBA9v5h4NnZ++uBHase1AHVzMzyDe5MSS/O3l8DvFvS1sB7gHurHtRdvmZmlmtAZ0o6Adgqe38s8DNgf2AlMLfqQR1Qzcyss/4PkKVExLea3v9W0vNJz/O+IyIeqnpcB1QzM8s3gLfNtIqIJ4Crx3scB1QzM8s1KF2+ko4rmjYiDq9yDgdUMzPLNzgt1B0KpqucewdUMzMbeBHRlQnxm/XlbTOSDpW0RNJySYsk7TpG+mmSPiXpdkkrJN0i6cBe5dfMbFD1aqaksvV+0347S1ol6ZqC6dfJ0r+8fC4767uAKmlf0pDmT5Ga6BcD50ma2WG3c4DXAweR7i3aD/h9d3NqZjYEenAfasV6H0kbA2cBvyh6rohYBdxOxfl6O+m7gAocDiyMiNMjYnFEHAbcCRzSLrGkPYHdgTkR8fOIuC0iroiIS3uXZTOzAdWbiR1K1ftNvgx8G7is5Pn+HVggaZPSOe2grwKqpHWBVwHnt2w6H9gpZ7e9gKuAD0u6W9KNkr4gab0O55kmaaPGAkyvI/9mZoNmnF2+05vrWknT1jp+tXofSQcA25ImZijrA8CuwD2S/iDp6ualwvGA/huUtBmpGX5/y/r7gS1z9nkBsAuwHHhrdoxTgE2AvOuo84BjxptZM7OBN75Rvne1bDkWmN+yrnS9L+mFwGeAXSNilaSSGeQHZXcoot8CakPr16c26xqmZNv2j4hHASQdDnxX0nsj4sk2+ywAmu9Jms7aX7yZmY3PDGBZ0+cVHdIWqvclTSV18x4TETdWyVREVGnVjqnfAupDwAhr/yrZnLV/vTTcC9zdCKaZxaQvYwbp+XZriIgVNH2xFX7dmJkNh/G1UJdFxNIxUpet96cDs4AdJJ2crZsCSNIq4A0RceFYWZT0LGBvUrfx5yPiYUk7AvdHxN1j7d9OX11DjYiVwCJgdsum2UDeIKNLgOdI2rBp3YuAUdzqNDMbl27fNlOh3l8KvALYvmk5jfQEme2By8csk/RK4EbgI8CRwLOyTW8l9WBW0lcBNXMc8M+SDpS0naTjgZmkPxiSFkg6qyn9t4E/AmdIeqmk3YDPA1/L6e41M7OiejPKt3C9HxGjEXFD80J6tuny7PPjBc93ZkS8kDT+puE8YLfSuc/0W5cvEXG2pE2Bo0mP17mBdEvM7VmSrUh/6Eb6xyTNBk4ijfb9I+m+1KN6mnEzswHUi7l8y9b7NXg18O426+8mfwDsmPouoAJExCmkkbrtts1ts+73rN1dYGZm49WjuXzL1vst2+ez9ujhTpYDG7VZ/2LgwRLHWUM/dvmamVm/6E2Xb6/9EDha0jOyz5HNyvQZ4HtVD+qAamZmw+ZI4Nmka6/rAb8Cbibd4vPxqgftyy5fMzPrD8qWsvv0s+xWnl0k/RWwI6lxeXVE/Hw8x3VANTOzfIPzPNTVJD0/m/f9QmDMe1aLcpevmZnl6tXj23rsVkm/lvTuOifId0A1M7N8gzkoaRbpCTVHkSbI/6GkfdpN3l+GA6qZmXU2WMGUiLg6Iv6VdG/r35AGJ30ZeEDS16oe1wHVzMyGUiQXRcS/AH8N3Aq8q+rxHFDNzCzXgF5DBUDScyV9WNI1wJXA48D7qh7Po3zNzCzfYI7yPRjYH9iZNKn+t4C3RMRt4zmuA6qZmeXqxVy+E+DfgP8EPhgR19R1UAdUMzPLN4AtVGBmRNSeSwdUMzPLNYgt1IiI7AHjBwHbkX4CLAYWRsSjVY/rQUlmZpZvAO9DlTQLuAX4ELAJsFn2/hZJO1Y9rluoZmY2bI4HfgT8S0SsApC0DnA6cAIVHzLugGpmZvkG8xrqLJqCKUBErJL0OeCqqgd1l6+ZmeUa0PtQl5JmSWr1XNIj3CpxQDUzs3wDeA0VOBtYKGnfbHKHGZL+ntTl+52qB3WXr5mZ5VIEKnmHSdn0E+BIUtg/i6fj4FPAqcBHqx7UAdXMzPIN4DXUiFgJfFDSPGBb0jPRb46IJ8ZzXAdUMzMbSlkAvb6u4zmgmplZrkGc2KFbHFDNzCzfAHb5dosDqpmZ5XILtTgHVDMzy+cWamEOqGZmlsst1OIcUM3MLJ9bqIV5piQzM7MauIVqZmYdDWsXblkOqGZmli8iLWX3GUIOqGZmlsuDkopzQDUzs3welFSYA6qZmeXSaFrK7jOMPMrXzMysBm6hmplZPnf5FuaAamZmuTwoqTgHVDMzy+fbZgpzQDUzs1xuoRbngGpmZvl8DbUwB1QzM8vlFmpxvm3GzMysBm6hmplZPg9KKswB1czMcrnLtzgHVDMzy+dBSYU5oJqZWS63UItzQDUzs3yjkZay+wwhB1QzM8vnLt/CfNuMmZlZDdxCNTOzXKLCNdSu5KT/9WULVdKhkpZIWi5pkaRdx0i/v6RrJT0h6V5JZ0jatFf5NTMbWI37UMsuJZWp9yW9TdIFkh6UtFTSZZLeOK5y1qDvAqqkfYETgE8BOwAXA+dJmpmTfhfgLGAh8DJgH+DVwOm9yK+Z2SBrjPItu5Q6R8l6H9gNuACYA7wKuAg4V9IO1UpZj74LqMDhwMKIOD0iFkfEYcCdwCE56V8D3BYRJ0bEkoj4NfBlYFZvsmtmNsCi4lJOqXo/Ig6LiM9FxJURcVNEfAy4CXhT6TPXqK8CqqR1Sb82zm/ZdD6wU85ulwIzJM1RsgWwN/DjDueZJmmjxgJMryH7ZmYDRxGVlsz05rpW0rS1jl+t3m89xhRSPf5w9ZKOX18FVGAzYCpwf8v6+4Et2+0QEZcC+wNnAyuB+4BHgPd3OM884NGm5a7xZNrMzNq6izXr2nlt0pSu99s4AtgAOKdaNuvRbwG1obXDQG3WpQ3SS4ETgU+QfuXsCWwDnNbh+AuAjZuWGePMr5nZYBqtuCQzWLOuXdDhTIXr/TUSSfsB84F9I+KBsQvUPf1228xDwAhr/yrZnLV/vTTMAy6JiM9nn6+T9DhwsaSjIuLe1h0iYgWwovFZGtZB3mZmnbV04RbeJ7MsIpaOkbxKvZ/OkwYzLQT2iYifl8pkF/RVCzUiVgKLgNktm2aTrpW2sz7Nv4eSkezVkdLMbDy6PCipYr3faJmeCbwzInLHzPRSv7VQAY4DviHpKuAy4GBgJlkXrqQFwNYR8U9Z+nOBr0o6BPgZsBVp+PUVEXFPj/NuZjZYevM81FL1fhZMzwI+CPxGUqN1+2REPFr25HXpu4AaEWdnkzIcTQqONwBzIuL2LMlWpD90I/2ZkqYD7wO+SBqQdCHwkV7m28xsEPXiaTNl633g3aT49R/Z0vB1YG65s9en7wIqQEScApySs21um3UnASeN+8RSWszaqTD7y6Q6n1k7vWmhlqr3I2KP0ifogb66hmpmZjZZ9WUL1czM+oNG01J2n2HkgGpmZvl61OU7CBxQzcwsnx8wXpgDqpmZ5RrnxA5DxQHVzMzyucu3MI/yNTMzq4FbqGZmli9Ye3LXIvsMIQdUMzPL5WuoxTmgmplZvqDCNdSu5KTvOaCamVk+D0oqzAHVzMzyjVL+QZieKcnMzGxNvoZanG+bMTMzq4FbqGZmls/XUAtzQDUzs3wOqIU5oJqZWT4H1MIcUM3MLJ9H+RbmgGpmZrk8yrc4j/I1MzOrgVuoZmaWz9dQC3NANTOzfKMBKhkgRx1QzczM1uQWamEOqGZm1kGFgDqkj5txQDUzs3xuoRbmgGpmZvlGg9ItTl9DtZ5R2bukzcys3zmgmplZvhhNS9l9hpADqpmZ5fM11MIcUM3MLJ+voRbmgGpmZvncQi3MAdXMzPIFFQJqV3LS9xxQzcwsn1uohflpM2ZmZjVwC9XMzPKNjlL6ieGjvm3GzMxsTe7yLcwB1czM8jmgFuaAamZm+XwfamEOqGZmlitilCg5lWDZ9IPCo3zNzMxq4BaqmZnliyjfhetrqGZmZi2iwjVUB1QzM7MWo6MgP76tCAdUMzPL5xZqYQ6oZmaWK0ZHiZIt1GEd5euAamZm+dxCLcy3zZiZmdXALVQzM8s3GiC3UItwQDUzs3wRlH7ajAOqmZnZmmI0iJIt1BjSgNp311Al7SbpXEn3SApJbymwz+6SFklaLulWSe/pQVbNzAZfjFZbSpJ0qKQlWT2+SNKuY6Tvu3q/7wIqsAFwLfC+IoklbQP8BLgY2AH4NHCipLd3LYdmZkMiRqPSUoakfYETgE+R6vGLgfMkzcxJ35f1ft91+UbEecB5AJKK7PIe4I6IOCz7vFjSLOBI4HvdyKOZmdXqcGBhRJyefT5M0huBQ4B5bdL3Zb3fdwG1gtcC57es+xlwkKRnRMRTrTtImgZMa1o1HWDV6MquZbJZ9Og8q42s6Nmp1OOyRU/L1rtzAcTI1J6da5C/N3pYNo0WagTUolf11apYUboLdxWrq93pLQ2jFRGxxpcvaV3gVcBnWg5zPrBTzilK1/u9MAgBdUvg/pZ195PKthlwb5t95gHHtK785R1fqT1zfeGPE52BLrpvojNglbT+H2vjsQmwtAvHXQnc92t+smXF/R8D7mpZdywwv2XdZsBU2tfjeeeuUu933SAEVFh7Gg/lrG9YABzX9Hk66YufASyrN2t9a9jKPGzlheEr87CW9+FuHDwilmfXKtet8bCduiba1eOdLsaWrfe7bhAC6n2s/Stmc2AVOW2zrMth9Rfb1CWxLCK68Uuv7wxbmYetvDB8ZR7i8nZNRCwHlnf5NA8BI7Svx/P6MkrX+73Qj6N8y7oMmN2y7g3AVRPVj25mZsVExEpgEWvX47OBS3N268t6v+8CqqQNJW0vafts1TbZ55nZ9gWSzmra5TTgeZKOk7SdpAOBg4Av9DbnZmZW0XHAP0s6MKvHjwdmkur3SVPv92OX7yzgoqbPjWudXwfmAluR/tAARMQSSXOA44H3AvcAH4iIMkOnV5Aulvd2GOfEGrYyD1t5YfjK7PJOUhFxtqRNgaNJdfwNwJyIuD1L0o16v3Ya1imizMzM6tR3Xb5mZmaTkQOqmZlZDRxQzczMauCAamZmVoOhCaiD8GigssqUWdLbJF0g6UFJSyVdlk1OPWmU/Y6b9ttZ0ipJ13Q5i7Wq8G96mqRPSbpd0gpJt2S3G0waFcq8v6RrJT0h6V5JZ2SjSfueH2U5+QxFQB2URwOVUbbMwG7ABcAc0kTVFwHnStqh+7kdvwrlbey3MXAW8Itu57FOFct7DvB60v16Lwb2A37f3ZzWp8L/x7uQvtuFwMuAfYBXA6e3S9+H/CjLySYiBn4BLgdObVm3GFiQk/6zwOKWdacBl010WbpV5pxj/A44eqLL0s3yAv8JfJI0Yfc1E12ObpUX2BN4BNhkovPewzIfCdzSsu79wJ0TXZYKZQ/gLWOkmfT11mRfBr6F2vRooNZH/VR5NNAsSc+oN4f1q1jm1mNMIU2+3ZWJt+tUtbySDgC2Jd0cP2lULO9ewFXAhyXdLelGSV+QtF4Xs1qbimW+FJghaY6SLYC9gR93L6cTalLXW4OgH2dKqtvAPBqohCplbnUEqcvpnBrz1S2lyyvphaTnL+4aEat6MdF4jap8vy8AdiFNdP7W7BinkB79NRmuo5Yuc0RcKml/4GzgmaT/f39EaqUOosleb016A99CbTLpHw1UQdkyp0TSfqQu0H0j4oEu5KtbCpVX0lTg28AxEXFjLzLWJWW+3ynZtv0j4oqI+AlwODB3srRSM4XLLOmlwInAJ0it2z2Bbcjmhx1Qg1BvTVrD0EIdmEcDlVClzMDqgR8LgX0i4ufdyV7typZ3OmnO6B0knZytmwJI0irgDRFxYbcyW4Mq3++9wN0R8WjTusWkCncGcFPdmaxZlTLPAy6JiM9nn6+T9DhwsaSjImLQWmyTvd6a9Aa+hRoD9GigoiqWudEyPRN4Z0RMmutMFcq7FHgFsH3Tchrwh+z95V3JaE0qfr+XAM+RtGHTuhcBo6SHVPe1imVen1S+ZiPZ66Tq4y9oUtdbA2GiR0X1YgH2BVaSrhVtR3pCwWPA87LtC4CzmtJvAzxOetLNdtl+K4G3T3RZuljm/YCngENJv3Iby8YTXZZulLfN/vOZXKN8y36/GwJ3Av8FvJR0m9SNwFcnuixdLPPc7N/0IaRryDsDVwKXT3RZCpZ3Q57+wRfAh7L3M3PKO+nrrcm+THgGelbQFChuIz3qaBGwW9O2M4FftqTfHbg6S78EeM9El6GbZQZ+mf1P27qcOdHl6NZ33LLvpAqoVcoLvIR0r/ETWXD9IrDeRJejy2V+P+n2rydIj/j6JrD1RJejYFn36PT/5KDWW5N58ePbzMzMajDw11DNzMx6wQHVzMysBg6oZmZmNXBANTMzq4EDqpmZWQ0cUM3MzGrggGpmZlYDB1QzM7MaOKCadYmkYyWNSrpL0qETnR8z6y7PlGTWJZLeBrwOeA9pMvatIuLBic2VmXWLW6hmXRIR34+I9wPfIj0ce9YEZ8nMusgB1az7rsheXzmhuTCzrnJANeu+O7LXV0xoLsysqxxQzbpIkoB52UcHVLMB5oBq1l3vBXbK3r9E0jMmMjNm1j0e5WvWJZJmAP9Hepj3g6SHP78iIm6Y0IyZWVe4hWrWPacCGwAHAVdm69ztazagHFDNukDSvsDfASdFxG+Aa7NNDqhmA8oB1axmkv4M+BJwG/DxbPV12asDqtmAWmeiM2A2gI4DtgDeEBGPZ+sWAyvxvahmA8stVLMaSXo9MBc4IyIuaKyPiKeA3wMzJW08Qdkzsy5yQDWriaT1gC8D9wFHtEnSuI768p5lysx6xgHVrD7HAtsC742IP7XZ3gio7vY1G0C+D9XMzKwGbqGamZnVwAHVzMysBg6oZmZmNXBANTMzq4EDqpmZWQ0cUM3MzGrggGpmZlYDB1QzM7MaOKCamZnVwAHVzMysBg6oZmZmNXBANTMzq8H/B7JCDcaqBFZxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x400 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(dpi=100)\n",
    "plt.imshow(ODeltaF_ij, extent=[0, 1, 1, 0])\n",
    "plt.xlabel(r\"$\\lambda$\", fontsize=14)\n",
    "plt.ylabel(r\"$\\lambda$\", fontsize=14)\n",
    "plt.title(\"Overlap Matrix\", fontsize=16)\n",
    "plt.clim(0, 1)\n",
    "cbar = plt.colorbar()\n",
    "cbar.ax.set_ylabel(\"overlap fraction\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f05a823e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}