{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "i8YiintZl8vd" }, "source": [ "# Tutorial\n", "\n", "\n", " \"Open\n", "\n", "\n", "In this notebook we will perform a multimap targeted free energy perturbation (TFEP) on a toy system using the `tfep` library and compare its performance with standard free energy perturbation.\n", "\n", "To run this you need to install the `tfep` library first together with the OpenMM optional package (see the [installation guide](https://tfep.readthedocs.io/en/latest/installation.html)). On Google Colab, the following cell will install the packages." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "DLIrGhOJl8vh", "outputId": "5174d235-82e3-4011-9ea7-c4dc8ec82c51" }, "outputs": [], "source": [ "import os\n", "\n", "if os.getenv(\"COLAB_RELEASE_TAG\"):\n", " try:\n", " import tfep\n", " except ModuleNotFoundError:\n", " !pip install git+https://github.com/andrrizzi/tfep.git\n", " !pip install openmm" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "UiYCiQe3up4p", "outputId": "fb83f03a-a730-4b31-a872-4e25c853344d" }, "outputs": [], "source": [ "# -------------- #\n", "# Global imports #\n", "# -------------- #\n", "import functools\n", "import os\n", "\n", "import lightning\n", "from matplotlib import pyplot as plt\n", "import numpy as np\n", "import MDAnalysis\n", "import openmm\n", "import pint\n", "import seaborn as sns\n", "import torch\n", "\n", "import tfep.analysis" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "id": "gglRfJo6rB3f" }, "outputs": [], "source": [ "# -------------------- #\n", "# Global configuration #\n", "# -------------------- #\n", "\n", "# TODO: what is the default context?\n", "sns.set_context('notebook')\n", "\n", "# Create the subfolder where to save all the generated data.\n", "MAIN_DIR = 'tfep-tutorial'\n", "os.makedirs(MAIN_DIR, exist_ok=True)\n", "\n", "# Create a registry for unit and define the thermodynamic constants.\n", "UNITS = pint.UnitRegistry()\n", "TEMPERATURE = 298.15 * UNITS.kelvin # Temperature of the system.\n", "RT = (TEMPERATURE*UNITS.molar_gas_constant).to(UNITS.kJ/UNITS.mol).magnitude # RT in units of kJ/mol" ] }, { "cell_type": "markdown", "metadata": { "id": "1jEMqeNBl8vk" }, "source": [ "## MD simulation of a triatomic molecule\n", "\n", "To perform a TFEP calculation we need samples from the perturbed state. We will consider a toy triatomic molecule in vacuum modeled using two harmonic bonds and a harmonic angle (think ozone) but no nonbonded interactions. We will collect samples for this system by running molecular dynamics using OpenMM." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "id": "_ILIdAWSl8vk" }, "outputs": [], "source": [ "def create_triatomic_molecule(\n", " r0=1.278*openmm.unit.angstroms, # equilibrium bond length.\n", " K_r0=290.1*openmm.unit.kilocalories_per_mole/openmm.unit.angstrom**2, # bond harmonic constant.\n", " alpha0=2.038*openmm.unit.radian, # equilibrium angle.\n", " K_alpha0=900.0*openmm.unit.kilocalories_per_mole/openmm.unit.radian**2, # angle harmonic constant.\n", " mass=15.999*openmm.unit.amu, # atoms mass.\n", "):\n", " \"\"\"Create a triatomic molecule OpenMM System object.\n", "\n", " Parameters\n", " ----------\n", " r0 : openmm.unit.Quantity\n", " The equilibrium bond length betwen atoms.\n", " K_r0 : openmm.unit.Quantity\n", " The force constant for the harmonic bond.\n", " alpha0 : openmm.unit.Quantity\n", " The equilibrium angle between atoms 0-1-2.\n", " K_alpha0 : openmm.unit.Quantity\n", " The force constant for the harmonic angle.\n", " mass : openmm.unit.Quantity\n", " The mass of the three atoms.\n", "\n", " Returns\n", " -------\n", " system : openmm.System\n", " The OpenMM System object to simulate.\n", " topology : openmm.app.Topology\n", " The OpenMM Topology object.\n", " starting_positions : openmm.unit.Quantity\n", " The starting positions for the simulations (the equilibrium\n", " configuration).\n", "\n", " \"\"\"\n", " from openmm.app import Topology, Element\n", "\n", " # Create an empty system object.\n", " system = openmm.System()\n", "\n", " # Add three particles to the system.\n", " n_atoms = 3\n", " for i in range(n_atoms):\n", " system.addParticle(mass)\n", "\n", " # Add a harmonic bond between atoms 0-1 and 1-2.\n", " force = openmm.HarmonicBondForce()\n", " force.addBond(0, 1, r0, K_r0)\n", " force.addBond(1, 2, r0, K_r0)\n", " system.addForce(force)\n", "\n", " # Add a harmonic angle between atoms 0-1-2.\n", " force = openmm.HarmonicAngleForce()\n", " force.addAngle(0, 1, 2, alpha0, K_alpha0)\n", " system.addForce(force)\n", "\n", " # Set the initial positions for a molecule in equilibrium.\n", " starting_positions = np.zeros([n_atoms, 3]) * openmm.unit.angstroms\n", " starting_positions[0, 0] = r0\n", " starting_positions[2, :2] = np.array([np.cos(alpha0._value), np.sin(alpha0._value)]) * r0\n", "\n", " # Create topology.\n", " topology = Topology()\n", " element = Element.getBySymbol('O')\n", " chain = topology.addChain()\n", " residue = topology.addResidue('O3', chain)\n", " atom0 = topology.addAtom('O', element, residue)\n", " atom1 = topology.addAtom('O', element, residue)\n", " atom2 = topology.addAtom('O', element, residue)\n", " topology.addBond(atom0, atom1)\n", " topology.addBond(atom1, atom2)\n", "\n", " return system, topology, starting_positions\n", "\n", "\n", "def run_md():\n", " \"\"\"Run an MD simulation of the toy system.\n", "\n", " This runs Langevin dynamics at temperature `TEMPERATURE` for 1ns using\n", " a timestep of 2fs and saving energies and configurations every 100fs.\n", "\n", " The results are saved in:\n", " - MAIN_DIR/md/traj.dcd : The trajectory holding the sampled configurations.\n", " - MAIN_DIR/md/energies.csv : The energies of the sampled configurations in units of kJ/mol.\n", "\n", " The function also saves the OpenMM Topology as a pdb file in MAIN_DIR/md/triatomic.pdb\n", "\n", " \"\"\"\n", " from openmm.app import Simulation, DCDReporter, StateDataReporter, PDBFile\n", "\n", " # Load the system from files.\n", " system, topology, positions = create_triatomic_molecule()\n", "\n", " # Create a simulation.\n", " simulation = openmm.app.Simulation(\n", " topology, system,\n", " integrator=openmm.openmm.LangevinMiddleIntegrator(\n", " TEMPERATURE.magnitude*openmm.unit.kelvin, # from pint unit to OpenMM unit framework.\n", " 1/openmm.unit.picosecond, # frictionCoeff\n", " 2.0*openmm.unit.femtosecond, # stepSize\n", " ),\n", " )\n", " simulation.context.setPositions(positions)\n", "\n", " # Configure output trajectory file. Output every 0.1 picosecond.\n", " os.makedirs(f'{MAIN_DIR}/md', exist_ok=True)\n", " simulation.reporters.append(openmm.app.DCDReporter(f'{MAIN_DIR}/md/traj.dcd', reportInterval=50))\n", " simulation.reporters.append(openmm.app.StateDataReporter(f'{MAIN_DIR}/md/energies.csv',\n", " reportInterval=50, step=True, potentialEnergy=True))\n", "\n", " # Write also a pdb file for readibility.\n", " PDBFile.writeFile(topology, positions, f'{MAIN_DIR}/md/triatomic.pdb')\n", "\n", " # Run for 1ns.\n", " simulation.step(500000)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "id": "gKd7cARGl8vl" }, "outputs": [], "source": [ "# Run the molecular dynamics simulation to sample configurations and potential energies.\n", "run_md()" ] }, { "cell_type": "markdown", "metadata": { "id": "NTdH3Wc_l8vl" }, "source": [ "## The free energy of the triatomic system\n", "\n", "Let's analyze the system. Because the system is in vacuum and has three atoms, the potential energy is invariant to rigid translations (3 degrees of freedom) and rigid rotations (3 DOFs). So the system has only 3 DOFs:\n", "- The distance between atoms 0 and 1.\n", "- The distance between atoms 1 and 2.\n", "- The angle between the vectors 0-1 and 1-2\n", "we can plot their sampled distributions." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 391 }, "id": "tXfJRYRcl8vm", "outputId": "639839be-67bd-40c8-90d6-5e66b2da9394" }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/andrea/miniforge3/envs/tfep/lib/python3.11/site-packages/MDAnalysis/coordinates/DCD.py:165: DeprecationWarning: DCDReader currently makes independent timesteps by copying self.ts while other readers update self.ts inplace. This behavior will be changed in 3.0 to be the same as other readers. Read more at https://github.com/MDAnalysis/mdanalysis/issues/3889 to learn if this change in behavior might affect you.\n", " warnings.warn(\"DCDReader currently makes independent timesteps\"\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAGACAYAAABbdrqFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABfq0lEQVR4nO3dfVxUdd7/8fcAAt6goCgiKChGoOEGZqXtz7tuVMxsa8vMx5qZtl6ppRZp6HqXmptY6WZl6GZtm+m666UmWm3qXnuZrd2oUVBTSEmamiKKJijD+f3hNRMDw/0Mwwyv5+PBQ+fM95z5Hmb4zPmc753JMAxDAAAAAADAJXzcXQEAAAAAALwZiTcAAAAAAC5E4g0AAAAAgAuReAMAAAAA4EIk3gAAAAAAuBCJNwAAAAAALkTiDQAAAACAC5F4AwAAAADgQn7urkBjFxwcrOLiYoWHh7u7KgCc5Mcff1RAQIAKCgrcXRWPRWwEvBPxsf6Ij4D3cUZsJPGuRnFxsUpKStxdDQBOxN90/REbAe/E33X9ER8B7+OMv2kS72pY71YePnzYzTUB4CzdunVzdxU8HrER8E7Ex/ojPgLexxmxkTHeAAAAAAC4EIk3AAAAAAAuROINAAAAAIALkXgDAAAAAOBCJN4AAAAAALgQiTcAAAAAAC7k9sR73bp1MplMFX5mzZplVy4jI0OJiYkKDAxU9+7d9dJLLzk8XlpamqKjoxUYGKg+ffpoz549DXAWAAAAAAA41mjW8d65c6fatGljexwREWH7/759+zRy5EiNHTtWzz33nPbu3aupU6fK399fEyZMsJVLS0tTamqqlixZoqSkJKWnp2vYsGHav3+/EhISGvR8AAAAAACQGlHi3bt3b4WGhjp8buHChUpKStLatWslSYMGDdKRI0c0d+5cjR8/Xj4+PiouLtaiRYs0bdo0PfHEE5KkAQMGKCEhQYsXL9bbb7/dYOcCAAAAAICV27uaV6e4uFi7du3SfffdZ7d9zJgx+vHHH3XgwAFJ0ocffqizZ89q9OjRtjK+vr4aNWqUMjIyZBhGg9YbAAAAAACpESXePXv2lK+vr7p166ZnnnlGFotFkpSTk6NLly4pPj7ernyPHj0kSdnZ2Xb/xsXFVShXWFioo0ePuvoUAAAAAACowO1dzcPDw7VgwQLdcMMNMplM2rp1q+bMmaOjR4/qxRdf1JkzZyRJwcHBdvuFhIRIkvLz8yVJZ86cUUBAgJo3b15pucjISId16NatW6X1y8vLU+fOnet0bgAAAAAAuD3xHjJkiIYMGWJ7fNttt6l58+Z6/vnnNXv2bNt2k8nkcP+y2x2VsXYxr2x/oKSkxNZjwio+Pl5+fm7/8wAAlyL+AUDtWSwW5eTk2B7HxMTI19fXjTWCJ2iU36z33nuv0tLSdPDgQUVFRUmSreXbyvrY2qIdEhKioqIiFRUVKTAw0FauoKDArpwjhw8frvS5qlrD4R2ys7M1adU2BYV1kSQVnjiiVyaLmfABeD3iHwDUXk5Ojiau2q5WoZ10/tQxpU8ertjYWHdXC41coxnjXVbZidBiYmLk7+9f4Y58VlaWJNnGflv/dVQuKCjIbnkyoLygsC4KjohRcESM7QIUcKd169bJZDJV+Jk1a5ZduYyMDCUmJiowMFDdu3fXSy+95PB4aWlpio6OVmBgoPr06aM9e/Y0wFnAExD/AKD2WoV2UlBYF7UK7eTuqsBDNMoW7w0bNsjX11eJiYkKCAjQ4MGDtXHjRk2fPt1WZv369QoPD1diYqIkqV+/fmrTpo02bNhg22axWLRx40YlJyfT1Rw25btWms1mZr1Ho7Vz5061adPG9rjsTcR9+/Zp5MiRGjt2rJ577jnt3btXU6dOlb+/vyZMmGArl5aWptTUVC1ZskRJSUlKT0/XsGHDtH//flo2AQAAGoDbE+8hQ4bo5ptv1jXXXCNJ2rp1q1599VU99thj6tixoyRp7ty56t+/vyZOnKgxY8Zo7969Sk9P1+rVq+Xjc6XRPiAgQHPmzFFqaqrat2+vpKQkrVmzRocPH2YNb9gp37XyeNZ+tYmuPPlgDCTcqXfv3goNDXX43MKFC5WUlKS1a9dKkgYNGqQjR45o7ty5Gj9+vHx8fFRcXKxFixZp2rRpeuKJJyRJAwYMUEJCghYvXkx8BAAAaABuzxzi4uK0Zs0a/fDDDyotLVVsbKxeeOEFTZ061Vamb9++2rJli1JTU/XGG28oMjJSK1eutGvRkaTHH39chmFo5cqVOnHihBISEpSRkUGLDiqwdq2UpMITeVWWZQwkGqPi4mLt2rVLS5cutds+ZswYpaen68CBA+rdu7c+/PBDnT17VqNHj7aV8fX11ahRo7R8+XIZhkGPIAAAABdze+K9YsUKrVixotpyycnJSk5OrrKMyWRSSkqKUlJSnFU9QJJ9og40pJ49e+rUqVOKiorSxIkT9eSTT8rX11c5OTm6dOmSbX4Lqx49eki6csOod+/ett4acXFxFcoVFhbq6NGjlS61CAAAAOdwe+INAKgoPDxcCxYs0A033CCTyaStW7dqzpw5Onr0qF588UXbyg7BwcF2+1lXcMjPz5d0ZQWIgIAANW/evNJylSXeVa3qkJeXp86dO9fp3OA5yg+1YZgNAAB1w7cnADRCQ4YM0ZAhQ2yPb7vtNjVv3lzPP/+8Zs+ebdteWTfxstsdlbFOKEg3c1Sl7FAbhtkAAFB3JN4A4CHuvfdepaWl6eDBg4qKipIkW8u3lfWxtUU7JCRERUVFKioqUmBgoK1cQUGBXTlHDh8+XOlzVbWGw7sw1AYAgPprlOt4AwAqKrvsXUxMjPz9/SvMuJ+VlSVJtrHf1n8dlQsKCrJbngwAGrt169bJZDJV+Jk1a5ZduYyMDCUmJiowMFDdu3fXSy+95PB4aWlpio6OVmBgoPr06aM9e/Y0wFkAaIpo8QYAD7Fhwwb5+voqMTFRAQEBGjx4sDZu3Kjp06fbyqxfv17h4eFKTEyUJPXr109t2rTRhg0bbNssFos2btyo5ORkupo3MeXHbJvNZrsbOoCn2Llzp9q0aWN7XPYm4r59+zRy5EiNHTtWzz33nPbu3aupU6fK39/fbkWctLQ0paamasmSJUpKSlJ6erqGDRum/fv3M6QCgNOReANAIzRkyBDdfPPNuuaaayRJW7du1auvvqrHHntMHTt2lCTNnTtX/fv318SJEzVmzBjt3btX6enpWr16tXx8rnRoCggI0Jw5c5Samqr27dsrKSlJa9as0eHDh1nDuwkqvzzi8az9ahNNggHP07t3b4WGhjp8buHChUpKStLatWslSYMGDdKRI0c0d+5cjR8/Xj4+PiouLtaiRYs0bdo0PfHEE5KkAQMGKCEhQYsXLyY+AnA6Em8AaITi4uK0Zs0a/fDDDyotLVVsbKxeeOEFTZ061Vamb9++2rJli1JTU/XGG28oMjJSK1eutGvRkaTHH39chmFo5cqVOnHihBISEpSRkUGLThNVdsx24Yk8u+dKSy0ym822x7SIw9MUFxdr165dWrp0qd32MWPGKD09XQcOHFDv3r314Ycf6uzZsxo9erStjK+vr0aNGqXly5fLMAx6BAFwKhJvAGiEVqxYoRUrVlRbLjk5WcnJyVWWMZlMSklJUUpKirOqBy914adjemZbkUKjLkqiRRyNV8+ePXXq1ClFRUVp4sSJevLJJ+Xr66ucnBxdunTJNr+FVY8ePSRd6fXRu3dv25CLuLi4CuUKCwt19OjRSpdaBIC6IPEGAAA2LdtHVtoiDrhbeHi4FixYoBtuuEEmk0lbt27VnDlzdPToUb344ou2lR2Cg4Pt9rOu4JCfny/pygoQAQEBat68eaXlqkq8q1rZIS8vT507d671uQHwbiTeaBLKTihE10kAADzTkCFDNGTIENvj2267Tc2bN9fzzz+v2bNn27ZX1k287HZHZazXB3Qzb5osFotycnLstsXExMjX19dNNYI3IfFGk1B2QiG6TgIA4D3uvfdepaWl6eDBg4qKipIkW8u3lfWxtUU7JCRERUVFKioqUmBgoK1cQUGBXbnKHD58uNLnqmoNR+OWk5Ojiau2q1VoJ0nS+VPHlD55uGJjY91cM3gD1vFGk2GdUKhlu3B3VwUAADhJ2V5sMTEx8vf3t1s2T5KysrIkyTb22/qvo3JBQUF2y5OhaWkV2klBYV0UFNbFloADzkDiDQCAlyopKVFmZqbth6E28EYbNmyQr6+vEhMTFRAQoMGDB2vjxo12ZdavX6/w8HAlJiZKkvr166c2bdpow4YNtjIWi0UbN25UcnIyXc0BOB1dzQEA8FKs2w1vM2TIEN1888265pprJElbt27Vq6++qscee0wdO3aUJM2dO1f9+/fXxIkTNWbMGO3du1fp6elavXq1fHyutDkFBARozpw5Sk1NVfv27ZWUlKQ1a9bo8OHDrOENwCVIvAEA8GJVrdsNeJq4uDitWbNGP/zwg0pLSxUbG6sXXnhBU6dOtZXp27evtmzZotTUVL3xxhuKjIzUypUrNWHCBLtjPf744zIMQytXrtSJEyeUkJCgjIwMJSRwcwqA85F4wyuVncVcqt1M5qWlFpnN5jrtCwAAXGfFihVasWJFteWSk5OVnJxcZRmTyaSUlBSlpKQ4q3oAUCkSb3il+nSvvPDTMT2zrUihURdrvS8AAAC8T9mlxnJzc0WbDGqLxBteqz7dK1u2j6x03/It4vHx8fLz408JAADAW5VdauzkNwfVOjJOkmSUlio3N9euLGt/wxGyBaCWyraIF544olcmi/FgAAAAXs661Nj5U8ds2y7kH9f8zUfUNqJAEmt/o3Ik3kAdlG0RBwAAQNPVom1H2/BGoDIk3gAAoFrlh9lIDLUBAKCm+LYEAADVKj/xJENtAACoORJvAABQIwyzAQCgbki8gXqg6yUAAACA6pAdAPVA10sAAAAA1SHxBuqJrpcAAADep+wa3bm5uTIMN1cIHo3EGwAAAADKKbtG98lvDqp1ZFy1+5RN1iUpJiZGvr6+rqwmPASJNwAAAAA4YF2j+/ypYzUqXzZZP3/qmNInD1dsbKyLawlPQOINAAAAAE5iTdaBsnzcXQEAAAAAALwZiTcAAAAAAC5EV3PAiVjXG4C7lZSUKDs7W5JkNptlMA0vAABuRzYAOBHregNwt+zsbE1atU1BYV10PGu/2kQTfwAAcDcSb8DJWNcbgLsFhXVRcESMCk/kubsqAABAjPEGAAAAAMClSLwBAAAAAHAhupoDAODByk6mJjGhGgAAjRGJNwAAHqzsZGqSmFANAIBGiMQbXoEWHwBNmXUyNUlMqAYAQCNE4g2vQIsPAAAAasJisSgnJ8f2OCYmRr6+vm6sEZoCEm94DVp8AAAAUJ2cnBxNXLVdrUI76fypY0qfPFyxsbHurha8HIk3AAAAgCalVWgnW09JoCGwnBgAAAAAAC5E4g0AAAAAgAvR1RwAAA/CKg4AAHgeEm8AADwIqzgAgGcwSkuVm5tre8zs6U1bo+pqfv78eUVGRspkMumTTz6xey4jI0OJiYkKDAxU9+7d9dJLLzk8RlpamqKjoxUYGKg+ffpoz549DVBzAAAajnUVh+CIGLVsF+7u6gAAHLiQf1zzNx/U9A0HNHHVdrslzND0NKrE++mnn1ZJSUmF7fv27dPIkSOVlJSkHTt2aNy4cZo6darWrFljVy4tLU2pqamaMmWKMjIy1L17dw0bNkyZmZkNdQoAAAAAPIS1VdpsNis3N1fOHrnTom1HBYV1UavQTs49MDxOo+lq/tVXX2nVqlVavny5Jk2aZPfcwoULlZSUpLVr10qSBg0apCNHjmju3LkaP368fHx8VFxcrEWLFmnatGl64oknJEkDBgxQQkKCFi9erLfffrvBzwkAAABA43WlVfqI2kYU6OQ3B9U6Ms7dVYKXajQt3o8++qgmTZqkq6++2m57cXGxdu3apfvuu89u+5gxY/Tjjz/qwIEDkqQPP/xQZ8+e1ejRo21lfH19NWrUKGVkZDDxDNyitNQis9mszMxM24+jXh0AAABwD2urdIuQDu6uCrxYo2jx3rRpkw4dOqRNmzbps88+s3suJydHly5dUnx8vN32Hj16SLoyyUzv3r1tM7zGxcVVKFdYWKijR48qMjLShWcBVHThp2N6ZluRQqMuSpIKTxzRK5OlhAQmQgIAAACaCrcn3j///LNmzJihZ555Rq1bt67w/JkzZyRJwcHBdttDQkIkSfn5+bZyAQEBat68eaXlKku8u3XrVmn98vLy1Llz55qdDOBAy/aRCo6IcXc1AMCprD16rKy9efz8frm0iI+Pt3sMAEBT5fZvw0WLFiksLEzjxo2rspzJZKp2u6My1i7mle0Pz8Q6tmhqzp8/r7i4OB09elQff/yxrrvuOttzGRkZmj17trKzsxUZGakZM2bokUceqXCMtLQ0vfjiizp+/LgSEhK0bNkyDRw4sAHPAt6kfI+e41n75dcyRKFRV0mihw8AAGW5NfH+/vvvtXz5cm3evFnnzp2TdOXi0vrv+fPnbS3W1pZvK+tj6/MhISEqKipSUVGRAgMDbeUKCgrsyjly+PDhSp+rqjUc7sM6tmhqqlv1YezYsXruuee0d+9eTZ06Vf7+/powYYKtnHXVhyVLligpKUnp6ekaNmyY9u/fT2KEOivbo6fwRJ78gkLp4QMAgANunVwtNzdXly5d0vDhwxUSEqKQkBCNGDFC0pWZy2+55RbFxMTI39/frnVTkrKysiTJNvbb+q+jckFBQYqIiHD16aCBsY4tmgrrqg8LFiyo8FzZVR8GDRqkOXPm6KGHHtLcuXNVWloqSRVWfRg8eLDefPNNde3aVYsXL27o0wEAAGhy3Jp4X3vttdq9e7fdz/PPPy9JeuWVV/TSSy8pICBAgwcP1saNG+32Xb9+vcLDw5WYmChJ6tevn9q0aaMNGzbYylgsFm3cuFHJycl0NQfgsVj1AQAA72GxXJkjo+yPxWJxd7XgYm7tah4cHFzp+MLevXsrKSlJkjR37lz1799fEydO1JgxY7R3716lp6dr9erV8vG5cu8gICBAc+bMUWpqqtq3b6+kpCStWbNGhw8fZg1vAB6LVR8AAPAuOTk5mrhqu1qFdpIknT91TOmThys2NtbNNYMruX1ytZro27evtmzZotTUVL3xxhuKjIzUypUr7cYvStLjjz8uwzC0cuVKnThxQgkJCcrIyGD8IgCP5O5VH1jxAQAA12gV2sk2VxGahkaXeA8cONBht8fk5GQlJydXua/JZFJKSopSUlJcVT0AaDCs+gAAAOAd3DrGGwDgmHXVhwULFujcuXMqKChwyqoPZVW36sPhw4cr/aG1G0BjcP78eUVGRspkMumTTz6xey4jI0OJiYkKDAxU9+7d9dJLLzk8RlpamqKjoxUYGKg+ffpoz549DVBzAE0NiTcANEKs+gAA1atuqcWkpCTt2LFD48aN09SpU7VmzRq7ctalFqdMmaKMjAx1795dw4YNU2ZmZkOdAoAmotF1NQcA/LLqQ1kHDx7U9OnT9corr6hPnz52qz5Mnz7dVq6qVR+s21j1AYCnsy61uHz5ck2aNMnuubJLLUpXblgeOXJEc+fO1fjx4+Xj41NhqUVJGjBggBISErR48WIm5/UyFotFOTk5kq7c3GZBDzQ0Em8AaIRY9QEAqlbdUotLly612z5mzBilp6frwIED6t27d5VLLS5fvlyGYXBj0ouUnUn85DcH1ToyrvqdACci8QYAD8aqDwCaIpZaRF1YZxI/f+qYu6uCJojEGwA8BKs+AID7l1qUWG4RQO0xuRoAAAA8BkstAvBEtHgDAADAI1iXWty8ebPOnTsnSU5ZajEwMNBWrrqlFqUryy1WpqrWcABNFy3eAAAA8AgstQjAU9HiDQBAI1ZSUmKXGJjNZodj/Rub0lKLzGaz3bb4+Hj5+XHpgbpjqUUAnopvPwAAGrHs7GxNWrVNQWFdJEnHs/arTXTjn43+wk/H9My2IoVGXZQkFZ44olcmi5n0US8stQhPZZSWKjc3VxLriDdVJN4AADRyQWFdFBwRI0kqPJHn5trUXMv2kbZ6Aw2JpRbR2FzIP675m4+obUQB64g3USTeAAAA8FgstQhP0aJtR9YRb8KYXA0AAAAAABci8QYAAAAAwIVIvAEAAAAAcCESbwAAAAAAXIjEGwAAAAAAFyLxBgAAAADAhUi8AQAAAABwIRJvAAAAAABcyM/dFQAAAN6vtNQis9lsexwfHy8/Py5DAABNA994AADA5S78dEzPbCtSaNRFFZ44olcmSwkJCe6uFgAADYLEGwAANIiW7SMVHBHj7moAANDgGOMNAAAAAIAL0eINj1BSUqLs7GzbY7PZLMMw3FgjAAAAAKgZEm94hOzsbE1atU1BYV0kScez9qtNNGMDAQAAADR+JN7wGEFhXWxjAwtP5Lm5NgAAAABQMyTeAAA0IgytAQDA+5B4AwDQiDC0BgAA70PiDQBAI8PQGgAAvAuJN9CASkstMpvNdtvi4+Pl58efIgAAAOCtuNoHGtCFn47pmW1FCo26KEkqPHFEr0yWEhLoRgoAAOAsFotFOTk5tse5ubliugy4E4k30MBato+0dSEFAACA8+Xk5Gjiqu1qFdpJknTym4NqHRnn5lqhKSPxBgAAAOB1WoV2sk1Uef7UMTfXpnJGaalyc3Ntj2NiYuTr6+vGGsEVSLwBAAAAwE0u5B/X/M1H1DaiQOdPHVP65OGKjY11d7XgZCTeAAAAAOBGLdp2tLXOwzv5uLsCAAAAAAB4MxJvAAAAAABciMQbAAAAAAAXIvEGAAAAAMCFSLwBAAAAAHAhEm8AAAAAAFyIxBsAAAAAABci8QYAAAAAwIVIvAEAAAAAcCG3J97vvvuuBgwYoPbt2ysgIEDdunXTjBkzdPbsWbtyGRkZSkxMVGBgoLp3766XXnrJ4fHS0tIUHR2twMBA9enTR3v27GmAswAAAAAAwDG3J975+fnq16+fXn31Vb377ruaMWOG3njjDd1zzz22Mvv27dPIkSOVlJSkHTt2aNy4cZo6darWrFljd6y0tDSlpqZqypQpysjIUPfu3TVs2DBlZmY29GnBCUpKSpSZmanMzEyZzWYZhuHuKgEAAABArfm5uwKjR4/W6NGjbY8HDhyogIAAPfzwwzp27Jg6deqkhQsXKikpSWvXrpUkDRo0SEeOHNHcuXM1fvx4+fj4qLi4WIsWLdK0adP0xBNPSJIGDBighIQELV68WG+//bZbzg91l52drUmrtikorIuOZ+1Xm+gEd1cJAAAAAGrN7S3ejrRr106SdPnyZRUXF2vXrl2677777MqMGTNGP/74ow4cOCBJ+vDDD3X27Fm7JN7X11ejRo1SRkYGraUeKiisi4IjYtSyXbi7qwIAAAC4lFFaqtzcXJnNZpnNZlksFndXCU7SaBJvi8WioqIiffbZZ1q4cKFGjBihqKgo5eTk6NKlS4qPj7cr36NHD0lXWkXL/hsXF1ehXGFhoY4ePdoAZwEAAAAAdXMh/7jmbz6o6RsOaOKq7crJyXF3leAkjSbxjoqKUvPmzdW7d2+Fh4dr/fr1kqQzZ85IkoKDg+3Kh4SESLoyRtxaLiAgQM2bN6+ynCPdunWr9CcvL88p5wcAtcHEk/BmpaUWmc1m2zwemZmZKikpcXe1AKBRaNG2o4LCuqhVaCd3VwVOVKfE29fXV/v373f43KeffipfX99aHzMjI0N79+7Vq6++qi+//FIjRoyw61phMpkc7ld2u6My1i7mle0PAM7izNjIxJPwZhd+OqZntn2ulE2HlLLpkCat2mbruQbv5IprR8Db0e3cu9RpcrWqxkuXlpbWKcnt1auXJKlfv35KSkrSddddp82bN9u6lFtbvq2sj60t2iEhISoqKlJRUZECAwNt5QoKCuzKOXL48OFKn+vWrVutzwVA0+TM2MjEk/B2LdtHKjgixt3VQANxxbUjUJ7FYrF1zc7NzZWnT/F0pdv5EbWNKND5U8eUPnm4YmNj3V0t1FGdu5pXFiA//fRTtWnTps4VkqRrr71Wvr6++vbbbxUTEyN/f/8Kd8KzsrIkyTb22/qvo3JBQUGKiIioV50AoCZcGRuZeBKAJ3NlfAQkKScnRxNXbdf0DQc05697VFxU5O4q1Rvdzr1HjRPvFStW2MY9m0wm3XnnnRXGQ4eHh2vy5Mm65ZZb6lWpffv2yWKxqFu3bgoICNDgwYO1ceNGuzLr169XeHi4EhMTJV1pKW/Tpo02bNhgK2OxWLRx40YlJydzJxWAS7g6NjLxJABP5Yr4yPwXqE6r0E4KCuuiFiEd3F0VwE6Nu5p36NBBPXv2lCR999136tatW4UJzwICApSQkKDHHnusxhW46667dN1116lXr15q3ry5Dh06pGeffVa9evXSnXfeKUmaO3eu+vfvr4kTJ2rMmDHau3ev0tPTtXr1avn4+Nhee86cOUpNTVX79u2VlJSkNWvW6PDhw3SlBOAyroqNVlFRUbbkeOjQoS6ZeDIyMtLha1c11CYvL0+dO3eu5dkAaEpcER+t819MmzZNISEh+uKLLzR//nx98cUXeu+99yT9Mv/F2LFj9dxzz2nv3r2aOnWq/P39NWHCBNuxrPNfLFmyRElJSUpPT9ewYcO0f/9+JSQkOOeXAAD/p8aJd9nxhoMGDdLLL79coQWlLq6//npt2LBBS5cuVWlpqaKjo/Xwww/riSeekL+/vySpb9++2rJli1JTU/XGG28oMjJSK1eutAuekvT444/LMAytXLlSJ06cUEJCgjIyMgieAFzGVbHRKiMjQ+fPn9eXX36pp59+WiNGjND7779ve56JJwE0Vq6Ij8x/AcBT1Wlytd27dzutArNmzdKsWbOqLZecnKzk5OQqy5hMJqWkpCglJcVZ1QNcyrqkjlV8fLz8/Or0Z4lGwJmx0YqJJwF4A1fERytH818sXbrUrsyYMWOUnp6uAwcOqHfv3lXOf7F8+XIZhsGNSQBOVecrfMMw9PHHH+v777/XxYsXKzw/duzYelUMaAquLKlTpNCoiyo8cUSvTBY9NDycK2Nj2YknR4wYYZt4cujQobYyVU08aZ0Tw1qOiScBNCRnxkeLxaLLly8rKyvLbv6LrKysaue/6N27d43mv6hsGA4A1EWdEm+z2aw77rhD33zzjcMZcU0mE4k3UEMsqeM9XB0bK5t4cvr06bYyVU08ad3GxJMAGpqz46M757+QmAMDQO3VKfGePHmyioqKtGHDBvXq1UsBAQHOrhcAeBxnxkYmngTgTZx97cj8FwA8TZ0S7/379ys9PV2//e1vnV0fAPBYzoyNTDzZtJSUlNi6vprNZtZXh9dx9rWjO+e/kJgDA0Dt1SnxbtWqlVq3bu3sugCAR3NmbGTiyaYlOztbk1ZtU1BYFx3P2q820dwUgXdx5bUj818A8AQ+ddnpwQcf1FtvveXsugCARyM2oj6CwrooOCJGLduFu7sqgNO5Mj5WNv9FWVXNf2HF/BcAXKlOLd7XXHON1q9frzvuuEMjRoywLeNQ1l133VXvygGAJyE2AoBjzoqPzH8BwFPVKfG+//77JUm5ubl65513KjxvMplksVjqVzM0OWXHOEpNb5xj+TW9Jdb19jTERgBwzFnxkfkvAHiqOl3R796929n1AOzGOEpqcuMcy67pLYl1vT0QsREAHHNWfGT+CwCeqk6J94ABA5xdD0DSL2McJanwRJ6ba9PwWNPbsxEbAcAx4iOApq5Ok6sBAAAAAICaqVOL9+DBg6t83mQy6YMPPqhThQDAUxEbAcAx4iOApq5OiXdpaWmFZRZOnTqlr7/+Wh06dFBsbKxTKgcAnoTYCACOER8BNHV1Srz37NnjcLvZbNbIkSM1b968+tQJADwSsREAHCM+AmjqnDrGOzY2VikpKXryySedeVgA8GjERgBwjPgIoKlw+uRq0dHR+uKLL5x9WADwaMRGAHCM+AigKahTV/Oq/P3vf1enTp2cfVigySkttchsNttti4+Pl5+f0/9s0QCIjQDgGPER9WGxWJSTkyNJys3NlWG4uUJAJep0BT9+/PgK24qLi/X5558rKytLzz77bL0rBjR1F346pme2FSk06qIkqfDEEb0yWUpISHBzzVAZYiMAOEZ8hKvk5ORo4qrtahXaSSe/OajWkXHurhLgUJ0S7127dlWYmTIwMFDR0dF66qmndP/99zulckBT17J9pIIjYtxdDdQQsREAHCM+wpVahXZSUFgXnT91zN1VASpVp8T7u+++c3I1AMDzERtRUyUlJcrOzrY9NpvNMugfCS9GfATQ1DFYFACABpadna1Jq7YpKKyLJOl41n61iWYYCQAA3qrOiXd+fr6ef/55ffDBBzp9+rRCQ0N1yy23aNq0aQoJCXFmHQHAYxAbUVNBYV1sQ0kKT+S5uTaA6xEfATRldVpO7OjRo0pKStLixYt19uxZdenSRQUFBXr66aeVlJSkY8cYXwGg6SE2AoBjxEcATV2dEu/U1FRdvHhR//nPf/Tll1/q/fff15dffqn//Oc/unjxolJTU51dTwBo9IiNAOAY8RFAU1enxHvnzp1atGiR+vTpY7e9T58+WrhwoXbs2OGUygGAJyE2AoBjxEc4k8VikdlsltlsZu1ueIw6jfE+e/asoqOjHT7XtWtXnT17tj51AgCPRGwEAMeIj3Am1u6GJ6pTi3fXrl21fft2h8/t2LFDXbt2rVelAMATERsBwDHiI5zNunZ3i5AO7q4KUCN1avF+8MEHNWvWLJWWluqBBx5QeHi4fvzxR7355pv605/+pKVLlzq7ngDQ6BEbAcAx4iOApq5OiXdKSopycnL04osvatWqVbbthmHo4Ycf1hNPPOG0CgKApyA2AoBjxEcATV2dEm+TyaTVq1drxowZ2r17t06fPq127dpp8ODBio2NdXYdAcAjEBsBwDHiI4CmrsZjvM+cOaO7775b77zzjm3b1VdfrUmTJmn27NmaNGmSzGaz7r77bp0+fdollQWAxobYCACOER8B4Bc1TrzXrFmjQ4cOaejQoZWWGTp0qDIzM+26EAGANyM2AoBjxEfAeYzSUuXm5tqWUTObzbJYLO6uFmqhxon322+/rYkTJ8rPr/Le6X5+fpo4caK2bt3qlMoBQGNHbAQAx4iPgPNcyD+u+ZsPavqGA5q+4YAmrtqunJwcd1cLtVDjxNtsNuu6666rtlxSUpLMZnO9KgUAnoLYCACOER8B52rRtqOCwrooKKyLWoV2cnd1UEs1TrxLSkrUrFmzass1a9ZMly9frlelAMBTEBsBwDHiIwD8osaJd3h4uLKysqot9+WXX6pjx471qhSahpKSEmVmZtp+zGazDMNwd7WAWiE2AoBjxEcA+EWNE+8BAwbopZdeqvKO5OXLl/Xyyy9r0KBBTqkcvFt2drYmrdqmlE2HlLLpkJ5++18quljk7moBtUJsBADHiI8A8IsaJ97Tp0/XV199pd/85jc6duxYheePHTumO++8U19//bWmT5/u1ErCOzhq4W7VobOCI2IUHBGjlu3C3V1FoNaIjQDgGPERAH5R+TST5fTq1UurVq3SI488oq5du6p3797q2rWrJCk3N1effvqpSktL9fLLLyshIcFlFYbnsrZwB4V1kSQdz9qvNtF8VuDZiI0A4BjxEQB+UePEW5ImTpyoa665RkuWLNHu3bv10UcfSZJatGihoUOH6qmnntKNN97okorCOwSFdVFwRIwkqfBEnptrAzgHsRHllZSUKDs7225bfHx8lcsqAd6I+AgAV9T6CqBv377atm2bSktLderUKUlSaGiofHxq3GsdALwOsRFlle/hU3jiiF6ZLFr1KlFaaqmwnBQ3KrwH8REA6pB4W/n4+KhDhw7OrAuAKnBh6hmIjbAq28MHVbvw0zE9s61IoVEXJUnnfszV40PMio2NtZUh3nk+4iOApoxvMMBDlL8wpQUNgDdp2T7SbijSM9s+J94BALwGiTfgQcpemAKANyPeAQC8CYNrAAAAAABwIRJvAAAAAABcyO2J99/+9jfdeeed6ty5s1q2bKlevXrp5ZdfVmlpqV25jIwMJSYmKjAwUN27d9dLL73k8HhpaWmKjo5WYGCg+vTpoz179jTAWQAAUDnr5IiZmZnKzMyU2WyWYRjurhYAAGggbh/jvXz5ckVFRWnZsmUKCwvT7t279eijj+rw4cNatmyZJGnfvn0aOXKkxo4dq+eee0579+7V1KlT5e/vrwkTJtiOlZaWptTUVC1ZskRJSUlKT0/XsGHDtH//fiZkAQC4TfnJEY9n7VebaL6XAABoKtyeeG/btk3t27e3PR40aJDOnz+vF198UYsWLVJAQIAWLlyopKQkrV271lbmyJEjmjt3rsaPHy8fHx8VFxdr0aJFmjZtmp544glJ0oABA5SQkKDFixfr7bffdsv5AQAgVZy1GwAANB1u72peNum2SkxMVFFRkfLz81VcXKxdu3bpvvvusyszZswY/fjjjzpw4IAk6cMPP9TZs2c1evRoWxlfX1+NGjVKGRkZdOkD4FEYhgMAAOA93J54O/Lvf/9bbdu2VYcOHZSTk6NLly4pPj7erkyPHj0kSdnZ2Xb/xsXFVShXWFioo0ePNkDNAcA5li9froCAAC1btkzvvPOO7rzzTj366KOaOXOmrYx1GE5SUpJ27NihcePGaerUqVqzZo3dsazDcKZMmaKMjAx1795dw4YNU2ZmZkOfFgAAQJPk9q7m5X3yySd67bXXNG/ePPn6+urMmTOSpODgYLtyISEhkqT8/HxJ0pkzZxQQEKDmzZtXWi4yMtLha3br1q3S+uTl5alz5851OhcAqCuG4QAA8AuLxaKcnBxJUm5urujMCk/TqFq8jx8/rrvvvlvXX3+9XauOJJlMJof7lN3uqIy1i3ll+wNAY8QwHACoiGE4TVdOTo4mrtqu6RsOaM5f96i4qMjdVQJqpdEk3mfPntWwYcPUokULbd26Vc2aNZP0S4u1teXbyvrY+nxISIiKiopUVO6PsKCgwK6cI4cPH670h9ZuAI0Fw3AANHUMw2naWoV2UlBYF7UI6eDuqgC11ii6mhcVFemOO+7QiRMntG/fPrVr1872XExMjPz9/ZWdna2hQ4fatmdlZUmS7aLT+m92drYSExPtygUFBSkiIqIhTgUAXIJhOGjKrOuglxUfHy8/v0ZxGYMGxDAcAJ7K7S3eJSUluvfee3Xo0CHt3LlTUVFRds8HBARo8ODB2rhxo9329evXKzw83JZk9+vXT23atNGGDRtsZSwWizZu3Kjk5GS6mgPwWAzDQVN3ZR30z5Wy6ZBSNh3SpFXbbL050LQwDAeAp3L7reLJkydr27ZtevbZZ/Xzzz/ro48+sj3Xo0cPtW7dWnPnzlX//v01ceJEjRkzRnv37lV6erpWr14tH58r9w4CAgI0Z84cpaamqn379kpKStKaNWt0+PBh7lq6UUlJie3iyGw280XmRLQANQ3OHIYTGBhoK1fTYTiVqao1HHCFsuugA2WVHYbz9ddfVzsMp3fv3jUahlNZbyAAqAu3X6G/++67kqQnn3yywnO7d+/WwIED1bdvX23ZskWpqal64403FBkZqZUrV2rChAl25R9//HEZhqGVK1fqxIkTSkhIUEZGhhISEhrkXFBRdna2Jq3apqCwLjqetV9tonkvnOVKC1CRQqMuSpIKTxzRK5PF592LMAwHAKrmjmE4EkNxANSe2xPv7777rkblkpOTlZycXGUZk8mklJQUpaSkOKFmcJagsC4KjohR4Yk8d1fF69AC5L3KDsP5n//5nyqH4UyfPt22vaphONZtDMMB4A0YhgPAk7g98QYAVMQwHAConDuH4UgMxQFQeyTeANAIMQwHqFz5OS6Y36JpYRgOAE/EtxQANEIMwwEqV3aOC+a3aFoYhgPAU5F4AwAAj8McF00Tw3AAeCoSbwAAAHgEhuF4N4vFopycHNvjmJgY+fr6urFGgPOQeAMAAMAjMAzHu+Xk5Gjiqu1qFdpJhSd/0Ozbr1HXrl0lSbm5ufq/SecBj0TiDQAAAKBRaBXaSUFhXXT+1DHN33xQbSMKJEknvzmo1pFx7q0cUA8k3gAAAAAanRZtOyoorIsk6fypY26uDVA/JN4AAAAA4EGM0lLl5ubaHjMevvEj8QYAAAAAD3Ih/7jmbz6ithEFOn/qmNInD1dsbKy7q4UqkHjDqUpKSpSdnW17bDabZTATBgAAAOBUZbvio/Ej8YZTZWdna9KqbbYgcDxrv9pEsywHAAAAgKaLxBtOFxTWRcERMZKkwhN5bq4NAAAAALiXj7srAAAAAACANyPxBgAAAADAhUi8AQAAAABwIRJvAAAAAABciMQbAAAAAAAXIvEGAAAAAMCFWE4MAIB6KikpUXZ2tu2x2WyWYRhurBEAAGhMSLwBAKin7OxsTVq1TUFhXSRJx7P2q010gptrBQAAGgsSbwAAnCAorIuCI2IkSYUn8txcGwBAU2GUlio3N9f2OCYmRr6+vm6sERwh8QYAAAAAD3Uh/7jmbz6ithEFOn/qmNInD1dsbKy7q4VySLwBAAAAwIO1aNvRNtwJjROzmgMAAAAA4EIk3gAAAAAAuBCJNwAAAAAALkTiDQAAAACAC5F4AwAAAADgQiTeAAAAAAC4EIk3AAAAAAAuROINAAAAAIAL+bm7AgAAAACA+jNKS5Wbm2t7HBMTI19fXzfWCFYk3gAAAADgBS7kH9f8zUfUNqJA508dU/rk4YqNjXV3tSASbwAA6qSkpETZ2dmSJLPZLMMw3FwjAACkFm07Kiisi7urgXJIvFEvZS88JS4+ATQd2dnZmrRqm4LCuuh41n61iU5wd5WapNJSi8xms922+Ph4+flxiQMAaDz4VkK9lL3wlMTFpxuVv/jkwhNwvaCwLgqOiFHhiTx3V6XJuvDTMT2zrUihURclSYUnjuiVyVJCAt9FAIDGg6ty1Jv1wlMSF59uVPbikwtPAE1Jy/aRtu8hAAAaIxJvwItYLz7pegkAAAA0HlyFA16IrpcAAABA40HiDXgpul4CAAAAjYOPuysAAAAAAIA3o8UbAAAAgFtYLBbl5ORIknJzc8WqtPBWJN4AAAAA3CInJ0cTV21Xq9BOOvnNQbWOjHN3lQCXoKs5AAAAALdpFdpJQWFd1CKkg7urArgMiTcAAAAAAC7k9sT722+/1aRJk3TttdfKz89P11xzjcNyGRkZSkxMVGBgoLp3766XXnrJYbm0tDRFR0crMDBQffr00Z49e1xYewAA0JiUllpkNpuVmZlp+ykpKXF3tQAATZzbx3h/+eWX2r59u2644QaVlpaqtLS0Qpl9+/Zp5MiRGjt2rJ577jnt3btXU6dOlb+/vyZMmGArl5aWptTUVC1ZskRJSUlKT0/XsGHDtH//ftYvBgCgCbjw0zE9s61IoVEXJUmFJ47olcniOgAA4FZub/EeMWKE8vLytGnTJiUlJTkss3DhQiUlJWnt2rUaNGiQ5syZo4ceekhz5861JerFxcVatGiRpk2bpieeeEKDBw/Wm2++qa5du2rx4sUNeUoA4BT0CALqpmX7SAVHxCg4IkZBYV3cXR0AANyfePv4VF2F4uJi7dq1S/fdd5/d9jFjxujHH3/UgQMHJEkffvihzp49q9GjR9vK+Pr6atSoUcrIyJDB2gQAPIy1R1D37t3Vo0cPh2WsPYKSkpK0Y8cOjRs3TlOnTtWaNWvsyll7BE2ZMkUZGRnq3r27hg0bpszMzIY4FQBwGm5KAvBEbk+8q5OTk6NLly4pPj7ebrv1IjQ7O9vu37i4uArlCgsLdfTo0QaoLQA4Dz2CAKAibkoCNWOUlio3N1dms9n2Y7FY3F2tJsvtY7yrc+bMGUlScHCw3faQkBBJUn5+vq1cQECAmjdvXmm5yMhIh6/RrVu3Sl8/Ly9PnTt3rlPdAaA+atojaOnSpXbbx4wZo/T0dB04cEC9e/euskfQ8uXLZRiGTCaTS84BAJxtxIgRGjlypCRp3Lhx+uSTTyqUKXtTUpIGDRqkI0eOaO7cuRo/frx8fHwq3JSUpAEDBighIUGLFy/W22+/3XAnBbjAhfzjmr/5iNpGFEiSzp86pvTJwxUbG+veijVRjb7F26qyi8Ky2x2VsXYx56ISgLehR1DDKikpsZsp22w2M4wJcAOGKQI116JtRwWFdVFQWBe1Cu3k7uo0aY2+xdvaYm1t+bayPrY+HxISoqKiIhUVFSkwMNBWrqCgwK6cI4cPH670uapawwHAnVzdI4jeQPays7M1adU222Rdx7P2q000M2UDjU1Nbkr27t27RjclK+stCXgia9dzq5iYGPn6+rqxRk1Lo0+8Y2Ji5O/vr+zsbA0dOtS2PSsrS5JsQdX6b3Z2thITE+3KBQUFKSIiogFrDQANhx5BDScorIuCI2IkSYUn8txcGwCOMEwRcKxs13O6nTe8Rp94BwQEaPDgwdq4caOmT59u275+/XqFh4fbkux+/fqpTZs22rBhg22bxWLRxo0blZyczIUlAK/j6h5B9AYC4Mm4KQlUZO16jobn9sT7559/VkZGhiTp+++/17lz57Rp0yZJVya4aN++vebOnav+/ftr4sSJGjNmjPbu3av09HStXr3aNs4nICBAc+bMUWpqqtq3b6+kpCStWbNGhw8fZnIMJyopKbF1zZLEGEfAjegRBAAVMUyx8bFYLMrJybF7LF0ZU5+bmysuJdEUuD3xPnnypO655x67bdbHu3fv1sCBA9W3b19t2bJFqampeuONNxQZGamVK1dqwoQJdvs9/vjjMgxDK1eu1IkTJ5SQkKCMjAwlJDAGz1kY4wg0HvQIAoCKuCnZ+OTk5Gjiqu22yb1OfnNQfi2C1TYiWie/OajWkXHVHAHwfG5PvKOjo2vUYpqcnKzk5OQqy5hMJqWkpCglJcVZ1YMDjHH0PKWlFpnNZrtt8fHx8vNzewhAFegRBAC1x03JxqlVaCdbw835U8fk17KtgsK66PypY26uGdAwuOoGmoALPx3TM9uKFBp1UZJUeOKIXpkseoM0cvQIAoCKuCkJwBOReANNRMv2kbaeCvAM9AgCgIq4KQnAE5F4AwAAwGNwUxKoP9b0bngk3qgSs5gDAAAA3oU1vRseiTeqxCzmAABPxuSSAOAYa3o3LL51UC1mMfc+XIgCaCqYXBIA0BhwlQ00QVyIAmhKmFwSAOBuJN5AE8WFKFA15rgAAADOQuINAIADzHEBAACchcQbAIBKMMcFAABwBh93VwAAAAAAAG9G4g0AAAAAgAuReAMAAAAA4EKM8QbAut4AAACAC3FVjQrKLqHD8jlNA+t6AywfBgAAXIfEGw4vNp977ysFdYxi+ZwmhHW90dSxfBgAAHAVEm9UerEZHBHD8jkAmhSWDwMAAK5A4g1JXGwCAJqm8r2+mN8CAOAKfLMAAIAmo/xkkmWHVzG/BeA8FotFOTk5kqTc3FwxZQaaOhJvAADQZJSfTLLs8CoAzpOTk6OJq7arVWgnnfzmoFpHxrm7SoBbkXgDAIAmpexkkgyvAlynVWgnBYV10flTx9xdFVTBKC1Vbm6u7XFMTIx8fX3dWCPvROINAAAAAE3Uhfzjmr/5iNpGFOj8qWNKnzxcsbGx7q6W1yHxBlCl8hMPSUw+BAAA4E1atO1oW+EIrsGVM4AqlV9ujsmH4C3K31Qym80ymP2nSSs/8ZrEjUYAgHPwTQKgWmWXmwO8RfmbStZJttB0lZ94jRuNAABnIfFugmjlAYAryt5UYpItSPYTrwEA4Cwk3k0QrTyoTtnultyYAQAANcHa3UDlSLybKFp5UJWy3S25MQMAAGqCtbuBypF4A3DI2t2SGzMAAKCmWLvbs5Vf01tiXW9nIfEGAAAAANit6S2Jdb2diMQbANBklJ1ckvkLAACoiDW9XYPE2wuVn7VcYh1SAJDsJ5dk/gIAqD8mVANqhkzMC5WftZx1SAE0VY6WT2zVoTPzFwCAkzChGlAzJN5equys5WWXhpLoXgmg6WD5RABwPSZUA6pH4t0ElF0aSuLCE0DTwvKJqKvyN64lhm4BAOqGb44mwro0lMSFJwDv5ahrOT18UFflb1wzdAtAU1N+eTHr0mJlx/aX3Y7KkXgDALwGXcvhbGVvXANAU1N2ebGyS4uVHdvPkmM1Q+INAPAqdC2Hq9D1HFCFlk5mMvd+lS0vZh3bj5rhmwJArVR34clydgC8FV3PAftZzCUxkzlQQ1wJA6iV6i48Wc4ODYkx3WhodD0H7Fs6mckcqBkSbwC1Vt2FZ9muvoArMaYbAAB4AhJvD0V3XgC4gjHdAACgsSNL81B054WnKD8mnBtEqK+yNx7pWg53YrI1eJvKlogqu53J1JquskuL8TmoPb4ZPISjcYytOnS2tfKU/fLnQhSNSdkx4dwggjOUvfFI13K4U20mW6OnGjxBZUtEld3OZGpNV9mlxfgc1B7R3kNUN46x7Jc/F6JoSOVbfBzd+LGOCad1CM5i7V5O13K4W00nW6OnGjyFdeK08q2bLdtd2c5kak2bdWmxsp+Dsp8VK2tvCfyCq10PUt04RuuXPxeiaEjlW3yquvHDUjwAmjImnoS7VNaFvCq0bqKmyn5WJNn1lsAvvC7xNpvNevTRR/Xvf/9bLVu21OjRo7V06VI1b97c3VUDvFbZFp/qbvyULVu+BbykpESS7FrAaRF3Dk+OjSwZBk/BnBaeyZPjY02V7SpeePIHzb79GnXt2lVS1Um4o9ZNwBHrZ0Wq2AJO6/cVXvVtUFBQoMGDBysqKkp///vfdfLkSc2YMUOnT5/Wm2++6e7qASjHUWu5X8sQhUZdJYkWcWdp7LGxurGvLBkGT8GcFp6nscdHZ7J2IT9/6pjmbz6othEFtEzCJcq2gPMZ+4VXJd6rV6/WmTNndPDgQYWGhkq60nI2ZswYzZ49W/Hx8W6uYeXKX3iWb/mjhQfeqnxruV9QKF0xnayxx8aajH1lyTB4ipqO+Ubj0Njjo6uUbZ0EXMH6GaP1+xdelXhnZGTolltusQVOSbr77rs1fvx4ZWRkNOrg6ahFp2zLHy08gD1mCK65xhAby79f5d+rsol1TSbsAxo7PseeoTHER3cqnxSxRBScrWzrd9lhDhaLRZLsknBvT8q96go1Oztb48ePt9sWEBCgmJiYChforlDbRKD8WrRllwcr3/JHCw+aoqrGgJvNZj333lcK6hgliW7pVXF3bLTWwXpz8dyPuXp8iNnW7ax8QlKbCfuAxqq6z3F9Vnmo7nqjuhtd+IW742PZSc/KJiJVJSU13afs48oS6vKTYjGJGlyh7FwB1mEOJ785KL8WwWobES1JlSblNf27qKm6TDToLF4Vhc+cOaPg4OAK20NCQpSfn1/pft26dav0udzcXPn5+VVZxurSpUv66cw5mf7vzTMsFgW3aq5mzZo5LH/58mUVnL8ok6+vLJcvy8evmXz9rpS1XC6WTD51elyffRvrsRpLPThWwx9rT6khH7//u4i4fFkmk498/Cr+zZRaSjTsnZfk7++v6uTl5TWpi1B3x0bpSnw8de5n+fj6yXK5WKPftH9f6xr/PPmzzbE4Vtn4Vt01Q1llrx8c7Vv2ecNiUfuQ1jWKjRLx0aqh4uPly5d1sqBQPj6+spRclslkuhIny/xfunKjpm2r5vLz81NJSYnyz1+sdp/yz/n4NpOP35UYbDL5Vvi/pEqf88R9PKWe7PPLPne8ZlT4LNf076Kmyv79lJZa1CE4qEZx1xmx0esiq8lkqrDNMAyH22t6vJr+kv39/RURFlp9wTLlW7ZsWcmzrerxuD77Vn2svLwzkqSIzp0buF6uOydnHCsv70qPhM6230vjqJe7j1Xx91LfejmHn5+fAgICXHLsxsqdsVG6Eu86hVov/Kt7XxtnXKk6/jWOv7mGPJb7vg8a57HqH+9qrurrh+qfrwrx8YqGio/NmjVTRPu2tTp+s2bNFFHJjOt5eXkqlaPrkbJaVvL/qp5r+H3y8q7c+OhkO5favk7jObeK8aHx1K22++Tl5cmi8p8xZ7+O65T9+8nLy9Px4xeq+Xu5whmx0asS75CQEJ05c6bC9oKCgirH6Bw+fNiV1fIq1ru3/M7s8XtxjN9L40BsdA4+z/b4fdjj9+GZvC0+etPnkHNpnDiXuvNpkFdpIPHx8RXG4xQXFysnJ8frJ8cAgMoQGwHAMeIjgIbiVYl3cnKyPvjgA50+fdq2bfPmzSouLlZycrIbawYA7kNsBADHiI8AGopXJd6///3vFRwcrJEjR+rdd9/VX/7yF02dOlVjxozhriWAJovYCACOER8BNBSvSryDg4O1a9cutWzZUnfddZdmzJih0aNHKz093d1VAwC3ITYCgGPERwANxasmV5Ok2NhYvfvuu+6uBgA0KsRGAHCM+AigIXhVizcAAAAAAI2NyTAMw92VAAAAAADAW9HiDQAAAACAC5F4AwAAAADgQiTeAAAAAAC4EIk3AAAAAAAuROINSdK3336rSZMm6dprr5Wfn5+uueaaavc5d+6c5s+frxtuuEHBwcFq3769hg4dqs8++6wBatww6vJ7kaSZM2eqZ8+eCgoKUuvWrdWnTx+9/fbbLq5tw6nr76WszZs3y2Qy1WlfwJmIf/aIe/aId3CHun7uLl26pJkzZ6pTp05q3ry5rr/+en3wwQcVyplMpgo/HTt2dPZp6G9/+5vuvPNOde7cWS1btlSvXr308ssvq7S0tNp9X3/9dcXFxSkwMFDXXHON/va3v1Uoc/nyZT311FMKDw9XixYtNGjQIH3++edOPw/J9efSUO+JVPdz2bBhg+6++25FRETIZDIpLS3NYTlPeF9qei7OfF9IvCFJ+vLLL7V9+3Z1795dPXr0qNE+R44c0erVq3XLLbdow4YNeu2112SxWNSvXz+vuPiU6vZ7kaQLFy5o0qRJ+vvf/66//e1vuvbaazV69Gi99dZbLqxtw6nr78Xq4sWLmjFjhsLCwlxQO6B2iH/2iHv2iHdwh7p+7qZNm6ZVq1Zp5syZ+u///m9169ZNycnJDuPS1KlTtW/fPttPRkaGM09BkrR8+XIFBARo2bJleuedd3TnnXfq0Ucf1cyZM6vcb9OmTRo3bpx+85vfaMeOHbr55ps1atQovffee3blpk+frlWrVmnhwoXasmWL/Pz8dPPNN+v48eMedy5Sw7wn9T2Xw4cPa8SIEVWW85T3pSbnIjnxfTEAwzAsFovt/w888IDRs2fPavc5f/68ceHCBbttFy9eNMLDw41x48Y5vY7uUJffS2X69etn3Hrrrc6oltvV9/fyhz/8wejfv3+9f6eAMxD/7BH37BHv4A51+dz98MMPhq+vr7Fy5UrbttLSUiMhIcG444477MpKMpYtW+a8Clfi5MmTFbZNnz7dCAwMNIqKiirdLy4uzrjnnnvstt12223GDTfcYHtsPd9Vq1bZtp07d85o166dMXPmTCfU3p4rz8UwGu49MYy6n0vZz2Vl9fWU96Um51Ldc7VFizckST4+tf8otGzZUi1atLDbFhgYqPj4eB07dsxZVXOruvxeKtOuXTtdvnzZacdzp/r8XnJycrR8+XKtXLnSiTUC6o74Z4+4Z494B3eoy+fu888/l8Vi0ZAhQ2zbTCaTbrvtNr377ru6dOmSM6tYI+3bt6+wLTExUUVFRcrPz3e4T25urr766iuNHj3abvv999+v/fv369SpU5Kk9957TxaLRffdd5+tTFBQkEaMGKHt27c78SyucOW5NLS6nItUs8+lJ7wvknO/62qKxBtOdeHCBR04cEDx8fHurorbGYahkpISFRQU6C9/+Yvee+89TZ482d3VcrvHHntMY8eO1a9+9St3VwVwKuIfca884h0aUlFRkSTJ39/fbntAQICKi4uVm5trt33p0qVq1qyZgoODNWrUKB05cqRB6vnvf/9bbdu2VYcOHRw+n52dLUkVYmmPHj1kGIa++uorW7mwsDC1bdu2Qrmvv/66RmOv68tZ52LlrvdEqv5casoT3pfactb74ueU2gD/Z86cOfr55581ZcoUd1fF7T744APdeuutkiQ/Pz+9+OKL+u1vf+vmWrnXtm3b9OGHH8psNru7KoDTEf+Ie2UR79DQYmNjJUn79+9XdHS0bftHH30kSXatf2PHjtXtt9+usLAwffHFF3r66af161//WocOHVJISIjL6vjJJ5/otdde07x58+Tr6+uwzJkzZyRJwcHBdtut9bKex5kzZyqUsZa7fPmyzp8/r9atWzuv8uU481wk970nUs3OpaY84X2pDWe+LyTecJq33npLL7zwglatWqXu3bu7uzpud8MNN+jjjz/W2bNntWPHDk2ZMkV+fn566KGH3F01tygqKtK0adO0YMEChYaGurs6gFMR/64g7l1BvIM79OzZUwMHDtTMmTMVGRmpq6++Wq+99pr+9a9/SbLvWvv666/b/t+/f3/9+te/VlJSktLT0/Xkk0+6pH7Hjx/X3Xffreuvv77aia+kK93ky7oy3NZ+e/kylZVzNlecizveE6n251ITnvK+1IQz3xcSbzjF+++/rwcffFApKSl65JFH3F2dRiEoKEjXXXedJOnmm29WcXGxZsyYoXHjxjnlDpyneeGFF+Tj46PRo0eroKBA0pVlT0pLS1VQUKAWLVpU6B4HeALi3y+Ie1cQ7+Au69at0z333KObbrpJkhQVFaW5c+dq3rx5VS6B1KtXL1199dX69NNPXVKvs2fPatiwYWrRooW2bt2qZs2aVVrW2op45swZu9UArH9L1udDQkJsLcplFRQUqFmzZmrZsqUTz+AXrjgXR1z9nki1O5ea8oT3pT7q874wxhv1tn//ft11112655579Mc//tHd1Wm0evfurXPnzumnn35yd1Xc4quvvtK3336r9u3bKyQkRCEhIVq/fr2ys7MVEhKiP//5z+6uIlBrxL+qNdW4R7yDu0RFRWn//v3Kzc3Vl19+qZycHDVv3lzh4eGKioqqcl9ri6SzFRUV6Y477tCJEye0c+dOtWvXrsry1vHQ1vHRVllZWTKZTIqLi7OVO3nyZIUJtLKysnT11Ve7ZPIsV51LZVz1nki1P5ea8oT3pb7q+r6QeKNesrOzlZycrJtuukmvvfaaS7uPeLr//d//VevWrZtst8NZs2Zp9+7ddj9DhgxRdHS0du/erTvuuMPdVQRqhfhXvaYa94h3cLfo6Gj16NFDly5d0tq1azVhwoQqyx88eFBms1l9+vRxaj1KSkp077336tChQ9q5c2e1yb8kde3aVXFxcdqwYYPd9vXr1+v666+3xZPbbrtNPj4+2rhxo63M+fPntW3bNg0fPtyp5yG59lwccdV7ItXtXGrKE96X+qjP+0JXc0iSfv75Z9ti8N9//73OnTunTZs2SZIGDBig9u3b66GHHtLrr7+ukpISSdLJkyc1ZMgQNWvWTCkpKXZdLgICApSYmNjwJ+Jkdfm9fP7555o5c6buueceRUdH24LN2rVrtXTpUvn5ef6fXV1+L3FxcRXu7K5bt04//PCDBg4c2KD1B8oi/tkj7tkj3sEd6vK5k6QXX3xRbdq0UefOnfXdd9/pueeeU2BgoN1417S0NB0+fFgDBgxQhw4d9MUXX2jx4sXq3LlztQl6bU2ePFnbtm3Ts88+q59//tk20Zt0ZZbr1q1bOzyPhQsXatSoUYqJidGtt96qLVu26L333tPOnTttZSIiIjRp0iTNnDlTfn5+ioqKUlpamiRp2rRpTj0PV59LQ74n9TmXrKwsZWVl2R5nZmZq06ZNatmypYYNGybJc96XmpyL098Xp6wGDo+Xm5trSHL4s3v3bsMwDOOBBx4wyn5kdu/eXek+UVFR7jkRJ6vL7+X48ePGfffdZ0RFRRkBAQFGhw4djP79+xv//d//7aazcL66/F4ceeCBB4yePXs2QI2ByhH/7BH37BHv4A51/dylpaUZ3bp1M/z9/Y3w8HBj8uTJRn5+vl2ZrVu3GjfeeKMREhJi+Pn5GeHh4cb48eONY8eOOf08oqKi6vz3s27dOiM2Ntbw9/c3evToYWzcuLFCmeLiYmPmzJlGWFiYERgYaAwYMMA4ePCg08/D1efSkO9Jfc5l3rx5Nfre84T3pSbn4uz3xWQYLhw8AAAAAABAE8cYbwAAAAAAXIjEGwAAAAAAFyLxBgAAAADAhUi8AQAAAABwIRJvAAAAAABciMQbAAAAAAAXIvEGAAAAAMCFSLwBAAAAAHAhEm8AAAAADeK7776TyWTSunXrnHrcf/7zn+rbt69atGih0NBQjRs3TidPnqzx/i+88ILuuusude3aVSaTSQMHDqzV6+/atUvjx49XXFycWrZsqYiICI0cOVKffvppjfYfOHCgTCaTTCaTbr/99lq9dn2YTCbNnz/f9njdunUymUz67rvvXPq61157rVvO151IvAEAAAB4rH/9618aNmyYwsLCtGXLFq1YsUL//Oc/dfPNN6u4uLhGx3jllVf0/fffa/DgwWrfvn2t6/Dyyy/ru+++02OPPaaMjAytWLFCJ0+e1I033qhdu3bV6BiJiYnat2+fli9fXuvXd5bhw4dr3759Cg8Pd+nr/OUvf9G+ffvUsWNHl75OY+Ln7goAsLdp0yY9+uijunjxop566ik9+eST7q4SADQKxEcAjqSkpCg2NlabNm2Sn9+V9KZr16666aab9Oc//1n/9V//Ve0xsrKy5ONzpU3ymmuuqXUdVq1apQ4dOthtGzp0qLp3764lS5Zo8ODB1R6jdevWuvHGG6stZ7FYVFJSooCAgFrXszrt27ev042H2kpISJAkl5xDY0WLN9CIFBYW6uGHH1ZaWpreeOMNLVu2TJmZme6uFgC4HfERcJ5vv/1WDz74oK666iq1aNFCERERGjFiRIW/qT179shkMmn9+vWaPXu2OnXqpNatW+uWW27R119/bVfWMAwtWbJEUVFRCgwM1HXXXaf3339fAwcOrFG37W+++Ub333+/OnTooICAAMXHx2vVqlXV7nf06FF9/PHH+t3vfmdLuiWpX79+io2N1ebNm2v0O7Em3XVVPumWpFatWqlHjx7Ky8ur83GtXfOfffZZLVq0SF27dlVAQIB2796toqIiPf7447r22mvVpk0btW3bVn379tWWLVsqHOfcuXOaOHGi2rVrp1atWmno0KEym80Vyjnqav7+++9r5MiRioyMVGBgoLp3767f//73OnXqlN2+8+fPl8lk0pdffqnRo0erTZs2CgsL0/jx43X27Nk6/w68BS3eQCNy8eJFGYahXr16qVWrVgoODlZJSYm7qwUAbkd8BJzn2LFjateunZYuXar27dsrPz9fr7/+um644QYdOHBAV199tV351NRU3XTTTVqzZo3OnTunmTNnasSIEcrOzpavr68kafbs2XrmmWf08MMP66677lJeXp4mTJigy5cvKzY2tsr6ZGVlqV+/furSpYuWL1+ujh076t1339Wjjz6qU6dOad68eZXu+8UXX0iSevXqVeG5Xr16ae/evbX99TjN2bNn9dlnn9Wotbs6K1euVGxsrNLS0tS6dWtdddVVKi4uVn5+vp544glFRETo0qVL+uc//6m77rpLr732msaOHSvpyk2RO++8Ux9++KHmzp2rPn36aO/evRo2bFiNXjsnJ0d9+/bVhAkT1KZNG3333Xd67rnn9Otf/1qZmZlq1qyZXfm7775bo0aN0kMPPaTMzEw99dRTkqQ///nP9f49eDQDaGDz5s0zyn70XnvtNUOSkZubW+Nj7N2715g3b55x5syZGu9TWFhoPPbYY0Z4eLgREBBg/OpXvzLWr19f4/3PnTtnpKSkGLfeeqsRGhpqSDLmzZtX4/0feOABQ5IhyejZs2el5caNG2crN3v27ArPb9682fa8JOPjjz+ucR0ANG7uiI/1jW2GYRgffPCB8eCDDxpXX3210aJFC6NTp07GHXfcYXzyySc12p/4CLhXSUmJcenSJeOqq64ypk+fbtu+e/duQ5KRnJxsV37jxo2GJGPfvn2GYRhGfn6+ERAQYIwaNcqu3L59+wxJxoABA2zbcnNzDUnGa6+9Zts2ZMgQIzIy0jh79qzd/lOmTDECAwON/Pz8Suv+17/+1a4uZT388MOGv79/tedfXs+ePe3qXFdjxowx/Pz8ahQLBwwY4PA1rb+vmJgY49KlS1Ueo6SkxLh8+bLx0EMPGYmJibbtO3bsMCQZK1assCu/ePHiCjG/uu+d0tJS4/Lly8b3339vSDK2bNlie876Hfbss8/a7fPII48YgYGBRmlpaYXjRUVFGcOHD6/yvLwFXc3hdnWZxOHDDz/UggULVFBQUON97rrrLr3++uuaN2+eduzYoT59+mj06NF66623arT/6dOn9eqrr6q4uFh33nlnjV+3rI4dO2rfvn1VvmZpaant/7/+9a8rPD9gwADt27dPc+bMqVMdAHiOhoiPzohtzphUiPgINJySkhItWbJEPXr0kL+/v/z8/OTv769vvvlG2dnZFcrfcccddo+trcvff/+9JOmjjz5ScXGx7r33XrtyN954o6Kjo6usS1FRkT744AP95je/UYsWLVRSUmL7SU5OVlFRkT766KNqz8lkMlW7veyxS0pKZBhGtcctyzCMCseozB/+8Af99a9/1fPPP6/evXvX6nUcueOOOyq0LEvS3/72N910001q1aqV/Pz81KxZM61du9bufdy9e7ckacyYMXb73n///TV67ZMnT2rSpEnq3Lmz7TWioqIkqcafl6KiolrNMu+N6GoOt2uISRwyMjL0/vvv66233tLo0aMlSYMGDdL333+vlJQUjRo1ytZVqjJRUVE6c+aMTCaTTp06pTVr1tS6HgEBAVVOmrF371795S9/0fDhw7V9+3YdPHhQQ4cOtSsTEhKiG2+8UV999VWtXx+AZ2mI+OiM2OaMSYWIj0DDmTFjhlatWqWZM2dqwIABCgkJkY+PjyZMmKCLFy9WKN+uXTu7x9YJsaxlT58+LUkKCwursK+jbWWdPn1aJSUl+tOf/qQ//elPDsuUH0vsqG7WOpSVn5+vtm3b2h6XT1xfe+01jRs3rsr6lfWvf/1LgwYNstuWm5tb4ebCggULtGjRIi1evFhTpkyp8fGr4ugG7D/+8Q/de++9uueee5SSkqKOHTvKz89PL7/8sl237tOnT8vPz6/C+1iTGcVLS0t122236dixY/rDH/6ghIQEtWzZUqWlpbrxxhvr9Hlpqki84VLbt2/X7NmzlZ2drU6dOmny5MkVyqxbt04PPvigXeD66aefNHv2bO3YsUMnT55U69atFRsbqwULFuh///d/tWDBAklXZqy02r17d6WTd2zevFmtWrXSPffcY7f9wQcf1P3336///Oc/6tevX5XnUtmdVGexWCyaMmWKIiMj9frrryssLEwHDx506WsCcJ/GEh+dEdtcNamQFfERcK4333xTY8eO1ZIlS+y2nzp1SsHBwbU+njXROnHiRIXnjh8/XmWrd0hIiHx9ffW73/3OYRyU7ONZedYZyDMzM5WcnGz3XGZmpt0M5R9//HGNj+tI7969KxyjU6dOdo8XLFig+fPna/78+UpNTa3V8aviKFa/+eab6tq1qzZs2GD3fPkl1Nq1a6eSkhKdPn3aLik+fvx4ta/7xRdf6NChQ1q3bp0eeOAB2/Zvv/22LqfRpJF4w2U++OADjRw5Un379tXbb78ti8WiZ5991mFQLu93v/udPvvsMy1evFixsbEqKCjQZ599ptOnT2vChAnKz8/Xn/70J/3jH/+w3QHs0aNHpcf74osvFB8fbzfbpfRLV6kvvvii2sTb1V555RUdPHhQb7/9ttq1a6errrqKC0vASzWm+OgqzpxUiPgIOJfJZKqwjNP27dt19OhRde/evdbHu+GGGxQQEKANGzborrvusm3/6KOP9P3331eZeLdo0UKDBg3SgQMH1KtXL/n7+9fqtSMiInT99dfrzTff1BNPPGHrwfjRRx/p66+/1rRp02xlr7vuulodu7ygoKAqj/H0009r/vz5mjNnTpUTwjmLyWSSv7+/XdJ9/PjxCrOaDxo0SM8++6z++te/6tFHH7Vtr8lwS+uxy39eVq9eXZ+qN0kk3nCZ2bNnKywsTO+//74CAwMlSUOGDKl2rI90pUvhhAkTNHHiRNu2kSNH2v7fpUsXSVJiYmKNjnf69Gl169atwnZr9yNH3ZMa0k8//aQ//OEPGjBggEaNGiXpyk2BTZs26eeff1aLFi3cWj8AztWY4qOrTJ48WRcuXNDs2bPrdRziI+B8t99+u9atW6e4uDj16tVLn376qZYtW6bIyMg6Ha9t27aaMWOGnnnmGYWEhOg3v/mNfvjhBy1YsEDh4eHVLtW1YsUK/frXv9b/+3//T//1X/+l6OhoFRYW6ttvv9W2bduqnSvij3/8o2699Vbdc889euSRR3Ty5EnNmjVL11xzjR588MEancMnn3xiW0Lr3LlzMgxDmzZtkiT16dPHNqa5MsuXL9fcuXM1dOhQDR8+vMK49Jqsz11bt99+u/7xj3/okUce0W9/+1vl5eXp6aefVnh4uL755htbudtuu039+/fXk08+qQsXLui6666zDd+pTlxcnGJiYjRr1iwZhqG2bdtq27Ztev/9951+Pt6OxBsuceHCBX388cd65JFHbBeV0pU7hSNGjNDrr79e5f7XX3+91q1bp3bt2umWW25R7969HU4oURtVdacsP/FGWb6+vi7vZj5r1iydO3dOK1eutG3r1auXNm7cqM8//9wlwRqAezTG+FhTNY2P1kmF/vSnP9V7UiHiI+B8K1asULNmzfTMM8/o/PnzSkpK0j/+8Y96TUy4ePFitWzZUq+88opee+01xcXF6eWXX9bs2bOr7b7eo0cPffbZZ3r66ac1Z84cnTx5UsHBwbrqqqsqdB93ZODAgcrIyNDcuXM1YsQItWjRQrfffruWLVtWoaW2Mi+++GKF+GsdoliTseDbtm2TJO3cuVM7d+6s8HxtJ3KriQcffFAnT57UK6+8oj//+c/q1q2bZs2aZbvpYeXj46OtW7dqxowZevbZZ3Xp0iXddNNNysjIUFxcXJWv0axZM23btk2PPfaYfv/738vPz0+33HKL/vnPf9pu9KKG3DqnOrxWXl6eIclYtGhRhedmzpxZ7XI5P/30k/HYY48ZUVFRhiSjVatWxu9+9zvjxx9/NAzDMJYtW1arJXZuvPFGo0+fPhW2f/HFF4YkY/Xq1YZh/LJkQ9mf3bt3V9jvp59+qtNyYlFRURW2/+c//zFMJpPxwAMPGGfOnLH9vPXWW4Yk4+WXX3Z4POvvjeVyAM/S2OJjWVXFtprGx/nz5xuSjMWLF9f4dYmPgHc6fPiw4e/vX6t40FQNGDDA6N+/v3H58mXDYrG4uzouZ136rCktJ0aLN1wiJCREJpPJ4aQNNZnIITQ0VC+88IJeeOEFHTlyRFu3btWsWbN08uRJh3cRq5OQkKD169erpKTEbpx3ZmampF8m5ujUqVOFSTOuvvrqWr9eTZWWlmry5MkyDEOvv/66w5YuxjEC3qWxxceaqkl8dOakQsRHwLMcOnRI69evV79+/dS6dWt9/fXXevbZZ9W6dWs99NBD7q6eR/if//kfNWvWTMOHD9c777zj7uq4VO/evXXo0CFJspsAz5uReMMlWrZsqeuvv17/+Mc/tGzZMlt3ysLCQltXnJrq0qWLpkyZog8++EB79+6VVPtlCX7zm98oPT1df//7321jBCXp9ddfV6dOnXTDDTdIkvz9/es98UZtpKen65NPPtGCBQvUv3//Cs/ffffdXFgCXqaxxceaqi4+OntSIeIj4FlatmypTz75RGvXrlVBQYHatGmjgQMHavHixdUuKYYrk5UVFhZKUp1mlvc0b731ln7++WdJTeN8JRJvuNDTTz+toUOH6tZbb9Xjjz8ui8WiP/7xj2rZsqXy8/Mr3e/s2bMaNGiQ7r//fsXFxSkoKEgff/yxdu7caZspMyEhQdKVMUoPPPCAmjVrpquvvlpBQUEOjzls2DDdeuut+q//+i+dO3dO3bt31/r167Vz5069+eab1a7hbbVjxw5duHDBFhizsrJsE28kJyfXapKf/Px8zZ49W/369dMf/vAHh+Mkf/WrX+k///mPSktLq52YBIDnaEzxUap/bHP2pELER8DzdO/eXf/85z/dXQ2P5coelo2RO1bbcDt393WHd9u6davRq1cvw9/f3+jSpYuxdOlSY968eVWOYSwqKjImTZpk9OrVy2jdurXRvHlz4+qrrzbmzZtnXLhwwbbfU089ZXTq1Mnw8fGpdKxhWYWFhcajjz5qdOzY0fD39zd69eplrF+/vlbnYx1T6einuvGU5ccw/v73vzf8/PyMzMzMSveZNm2aIcnIzs6u8BxjGAHP1pjiY31im2FcGZtY2f41udQgPgIAvJ3JMFwwxR6ACsaNG6c9e/bo22+/lclkqnEre3mGYchiseiNN97QQw89pI8//rhBu8cDgLMRHwEA3o6+WUAD+v7779WsWTP96le/qvMxtmzZombNmjFRCQCvQnwEAHgzWryBBvLdd9/p1KlTkqTmzZurZ8+edTpOQUGBvv32W9vjHj161GpsOQA0NsRHAIC3I/EGAAAAAMCF6GoOAAAAAIALkXgDAAAAAOBCJN4AAAAAALgQiTcAAAAAAC5E4g0AAAAAgAuReAMAAAAA4EIk3gAAAAAAuBCJNwAAAAAALkTiDQAAAACAC5F4AwAAAADgQv8fLGw5/5uN0XIAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def analyze_configurations():\n", " # Compute the three degrees of freedom.\n", " dist_10 = np.empty(10000)\n", " dist_12 = np.empty(10000)\n", " angles = np.empty(10000)\n", " with MDAnalysis.coordinates.DCD.DCDReader(f'{MAIN_DIR}/md/traj.dcd') as trajectory:\n", " for idx, timestep in enumerate(trajectory):\n", " positions = timestep.positions\n", " r10 = positions[0] - positions[1]\n", " r12 = positions[2] - positions[1]\n", " dist_10[idx] = np.linalg.norm(r10)\n", " dist_12[idx] = np.linalg.norm(r12)\n", " angles[idx] = np.arccos(np.dot(r10/dist_10[idx], r12/dist_12[idx]))\n", "\n", " # Plot the histograms of the three degrees of freedom.\n", " fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 4))\n", " for ax_idx, data in enumerate([dist_10, dist_12, angles]):\n", " sns.histplot(data, ax=axes[ax_idx])\n", " axes[0].set_xlabel('dist 0-1 [$\\AA$]')\n", " axes[1].set_xlabel('dist 1-2 [$\\AA$]')\n", " axes[2].set_xlabel('angle 0-1-2 [radian]')\n", " fig.tight_layout()\n", " plt.show()\n", "\n", "analyze_configurations()" ] }, { "cell_type": "markdown", "metadata": { "id": "N9njgFNkl8vm" }, "source": [ "The three DOFs have quite narrow distributions, and if we perturb the equilibrium bond length $r_0$ we can expect the overlap between the reference and pertrubed distributions to decrease quite fast.\n", "\n", "In this simple case, it is possible to predict the free energy difference as a function of $r_0$. We can remove the rigid rototranslational contributions by placing atom 1 in the origin, atom 0 on the $x$ axis, and atom 2 on the $x-y$ plane. Then, after converting the coordinates of atom 2 from cartesian to polar, it is easy to show that the free energy difference (in reduced units) of perturbing $r_0 \\to r'_0$ is\n", "\n", "$$\n", "\\Delta f = - \\log \\frac{\\int e^{-K_{r_0}(x - r'_0)^2} \\mathrm{d}x \\int e^{-K_{r_0}(r - r'_0)^2} r \\mathrm{d}r \\int e^{-K_{\\alpha_0}(\\alpha - \\alpha_0)^2} cos(\\alpha) \\mathrm{d}\\alpha}{\\int e^{-K_{r_0}(x - r_0)^2} \\mathrm{d}x \\int e^{-K_{r_0}(r - r_0)^2} r \\mathrm{d}r \\int e^{-K_{\\alpha_0}(\\alpha - \\alpha_0)^2} cos(\\alpha) \\mathrm{d}\\alpha} = -\\log \\frac{r'_0}{r_0}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "id": "exIRyhT4l8vm" }, "source": [ "## Standard free energy perturbation\n", "\n", "We can now check the performance of standard FEP as the perturbation in $r_0$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "u06QnB8Wl8vm" }, "outputs": [], "source": [ "def compute_fep_estimate(target_r0):\n", " \"\"\"Standard free energy perturbation.\n", "\n", " Compute the free energy difference between the reference distribution and\n", " the distribution with equilibrium bond length target_r0.\n", "\n", " Parameters\n", " ----------\n", " target_r0 : openmm.unit.Quantity\n", " The equilibrium bond length of the target potential energy function.\n", "\n", " \"\"\"\n", " # First, we need to re-evaluate the potential energy of the configurations\n", " # sampled with the reference MD simulation using target_r0.\n", " openmm_context = openmm.Context(\n", " create_triatomic_molecule(r0=target_r0)[0],\n", " openmm.openmm.VerletIntegrator(2.0*openmm.unit.femtosecond),\n", " openmm.Platform.getPlatformByName('CPU'),\n", " )\n", "\n", " # Compute the target potentials for all configurations.\n", " target_potentials = []\n", " with MDAnalysis.coordinates.DCD.DCDReader(f'{MAIN_DIR}/md/traj.dcd') as trajectory:\n", " for timestep in trajectory:\n", " # MDAnalysis reads positions in Angstrom but OpenMM accepts nanometers.\n", " openmm_context.setPositions(timestep.positions / 10)\n", " state = openmm_context.getState(getEnergy=True)\n", "\n", " # OpenMM return the energies in kJ/mol.\n", " target_potentials.append(state.getPotentialEnergy()._value)\n", "\n", " # Convert from list to numpy array.\n", " target_potentials = np.array(target_potentials)\n", "\n", " # Now load the energies sampled with the reference MD.\n", " steps, reference_potentials = np.loadtxt(f'{MAIN_DIR}/md/energies.csv', delimiter=',').T\n", "\n", " # Compute the free energy difference.\n", " work = torch.tensor(target_potentials - reference_potentials)\n", " df = tfep.analysis.fep_estimator(work, kT=RT)\n", "\n", " # The TFEP library has functions to perform (bayesian) bootstrap\n", " # analysis to get an approximate estimate of the uncertainty.\n", " fep_estimator = functools.partial(tfep.analysis.fep_estimator, kT=RT)\n", " bootstrap_analysis = tfep.analysis.bootstrap(work, fep_estimator, n_resamples=2000, bayesian=True)\n", "\n", " df_mean = bootstrap_analysis['mean']\n", " ci_low = bootstrap_analysis['confidence_interval']['low']\n", " ci_high = bootstrap_analysis['confidence_interval']['high']\n", " print(f'Df = {df} ({df_mean}) [{ci_low}, {ci_high}] kJ/mol')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "ALaYy8JJl8vn", "outputId": "0c4976a1-239e-439f-b120-0e366d1b5ee1" }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.10/dist-packages/MDAnalysis/coordinates/DCD.py:165: DeprecationWarning: DCDReader currently makes independent timesteps by copying self.ts while other readers update self.ts inplace. This behavior will be changed in 3.0 to be the same as other readers. Read more at https://github.com/MDAnalysis/mdanalysis/issues/3889 to learn if this change in behavior might affect you.\n", " warnings.warn(\"DCDReader currently makes independent timesteps\"\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Df = -0.2860998560439967 (-0.2855375363935697) [-0.3484750311933955, -0.2271823693093556] kJ/mol\n", "Correct Df = -0.06130654404637665 kJ/mol\n" ] } ], "source": [ "compute_fep_estimate(target_r0=1.31*openmm.unit.angstroms)\n", "print(f'Correct Df = {-RT*np.log(1.31/1.278)} kJ/mol')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "IgKXnhV_l8vn", "outputId": "dab4115b-8754-40d9-b470-7eee5819c658" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Df = 109.42345178743534 (110.74018250759458) [106.15834910440587, 117.06674738063921] kJ/mol\n", "Correct Df = -0.7073255071448409 kJ/mol\n" ] } ], "source": [ "compute_fep_estimate(target_r0=1.7*openmm.unit.angstroms)\n", "print(f'Correct Df = {-RT*np.log(1.7/1.278)} kJ/mol')" ] }, { "cell_type": "markdown", "metadata": { "id": "fPsvuTJul8vn" }, "source": [ "As expected, the performance of standard FEP decreases very rapidly with the size of the perturbation." ] }, { "cell_type": "markdown", "metadata": { "id": "gtN1mYUFl8vo" }, "source": [ "## Targeted free energy perturbation with the `tfep` library\n", "\n", "Now we can test TFEP instead. The following two functions implement a single-epoch multimap TFEP calculation as described in [this paper](https://doi.org/10.1073/pnas.2304308120)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "4szLMpuVl8vo" }, "outputs": [], "source": [ "def run_tfep_training(output_dir_path, target_r0=1.278*openmm.unit.angstroms, n_epochs=1):\n", " \"\"\"Trains a TFEP map and saves the potential energies evaluated with r0=target_r0.\n", "\n", " Parameters\n", " ----------\n", " output_dir_path : str\n", " Where to save the output files.\n", " target_r0 : openmm.unit.Quantity\n", " The equilibrium bond length of the target potential energy function.\n", " n_epochs : int\n", " The number of training epochs.\n", "\n", " \"\"\"\n", " import tfep.app\n", " import tfep.potentials.openmm\n", "\n", " # Create TFEP map.\n", " tfep_map = tfep.app.CartesianMAFMap(\n", " # Target potential energy function.\n", " potential_energy_func=tfep.potentials.openmm.OpenMMPotential(\n", " system=create_triatomic_molecule(r0=target_r0)[0],\n", " platform='CPU', # Change this to 'CUDA' to run on the GPU.\n", " positions_unit=UNITS.angstrom, # Units of the input positions.\n", " energy_unit=UNITS.kJ/UNITS.mole, # Units of the returned energies.\n", " # Useful for caching the OpenMM Context and avoid re-initializing it at every batch.\n", " system_name='triatomic',\n", " precompute_gradient=True, # Necessary to compute gradients for training.\n", " ),\n", " # Data.\n", " topology_file_path=f'{MAIN_DIR}/md/triatomic.pdb',\n", " coordinates_file_path=f'{MAIN_DIR}/md/traj.dcd',\n", " # We don't really need to map the rototranslational degrees of freedom.\n", " # origin_atom and axes_atoms force the map to place atom 1 in the origin\n", " # and atoms 0, 2 on the x-y plane and maps only the unconstrained DOFs.\n", " origin_atom=[1],\n", " axes_atoms=[0, 2],\n", " # We need to flag atom 1 as an atom whose coordinates are not changed by\n", " # the map but affect the transformation of the other DOFs.\n", " conditioning_atoms=[1],\n", " # Other options.\n", " temperature=TEMPERATURE, # Temperature.\n", " batch_size=16, # Batch size.\n", " tfep_logger_dir_path=f'{output_dir_path}/tfep_logs', # Where to save the target potential energies.\n", " )\n", "\n", " # Optimize the map and save the computed target potential energies.\n", " trainer = lightning.Trainer(\n", " max_epochs=n_epochs, # Run for n_epochs.\n", " default_root_dir=output_dir_path, # Where to save the Lightning logs.\n", " )\n", " trainer.fit(tfep_map)\n", "\n", "\n", "def compute_single_epoch_mtfep_estimate(dir_path):\n", " \"\"\"This function reads the potentials saved during training and compute the Df.\n", "\n", " It performs a bootstrap analysis of the data as a function of the number of\n", " samples to evaluate the convergence.\n", "\n", " \"\"\"\n", " # Read the data from the Read the data from the TFEP logger (DFT potentials).\n", " tfep_logger_dir_path = f'{dir_path}/tfep_logs'\n", " tfep_logger = tfep.io.TFEPLogger(tfep_logger_dir_path)\n", "\n", " # This reads the potential energies evaluated during the first epoch of training\n", " tfep_data = tfep_logger.read_train_tensors(epoch_idx=0)\n", "\n", " # Read the vacuum potentials (saved in kJ/mol).\n", " steps, reference_potentials = np.loadtxt(f'{MAIN_DIR}/md/energies.csv', delimiter=',').T\n", " reference_potentials = torch.tensor(reference_potentials)\n", "\n", " # The samples where shuffled during training so we need to match them.\n", " reference_potentials = reference_potentials[tfep_data['trajectory_sample_index']]\n", "\n", " # Potentials and log_det_Js in the tfep logger are stored in in units of kT\n", " # while OpenMM saves energies in kJ/mol.\n", " target_potentials = tfep_data['potential'] * RT\n", " log_det_J = tfep_data['log_det_J'] * RT\n", "\n", " # Compute the generalized work (in units of kT).\n", " work = target_potentials - log_det_J - reference_potentials\n", "\n", " # Compute the free energy and work statistics with the whole dataset.\n", " result = {'df': tfep.analysis.fep_estimator(work, kT=RT)}\n", "\n", " # Run the bootstrap analysis as a function of the number of samples.\n", " fep_estimator = functools.partial(tfep.analysis.fep_estimator, kT=RT)\n", " result['bootstrap_sample_size'] = torch.arange(1000, 11000, 1000)\n", " result['df_bootstrap'] = tfep.analysis.bootstrap(\n", " data=work,\n", " statistic=fep_estimator,\n", " n_resamples=2000,\n", " bootstrap_sample_size=result['bootstrap_sample_size'],\n", " take_first_only=True,\n", " batch=1000,\n", " bayesian=True,\n", " )\n", "\n", " return result\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 712, "referenced_widgets": [ "2f034474d13a49b6b8065d3caaadbbf3", "1d97d8e663044c0c8ad072aead4afe0f", "6a250960606e43ff8849ca4e7694b888", "d1ba10853d6942eeb638882da3fc2dfe", "affadb1606fd4833a61dbcaff8c550c8", "e7494aac57b54e1abf008651d19dd24c", "75cd8e35ace64d398411ae1317940e92", "cd0faa7c5a5e42fbafa26f8496d67d3e", "8eff48866de6405ab2a3923d8d832caa", "a52f60082adb49d0a32c9c27f6593bc6", "a2a4c9492ab74f2c880e4d28aff68338" ] }, "id": "M_y50_4Bl8vo", "outputId": "5e84393f-3892-4621-884f-37e60d5d9467" }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO: GPU available: False, used: False\n", "INFO:lightning.pytorch.utilities.rank_zero:GPU available: False, used: False\n", "INFO: TPU available: False, using: 0 TPU cores\n", "INFO:lightning.pytorch.utilities.rank_zero:TPU available: False, using: 0 TPU cores\n", "INFO: HPU available: False, using: 0 HPUs\n", "INFO:lightning.pytorch.utilities.rank_zero:HPU available: False, using: 0 HPUs\n", "WARNING: Missing logger folder: tfep-tutorial/mtfep/lightning_logs\n", "WARNING:lightning.pytorch.loggers.tensorboard:Missing logger folder: tfep-tutorial/mtfep/lightning_logs\n", "/usr/local/lib/python3.10/dist-packages/tfep/nn/flows/oriented.py:139: UserWarning: Using torch.cross without specifying the dim arg is deprecated.\n", "Please either pass the dim explicitly or simply use torch.linalg.cross.\n", "The default value of dim will change to agree with that of linalg.cross in a future release. (Triggered internally at ../aten/src/ATen/native/Cross.cpp:62.)\n", " plane_normal_vector = torch.cross(axis_vector, plane_axis_vector)\n", "INFO: \n", " | Name | Type | Params | Mode \n", "------------------------------------------------------------------------\n", "0 | _potential_energy_func | OpenMMPotential | 0 | train\n", "1 | _loss_func | BoltzmannKLDivLoss | 0 | train\n", "2 | _flow | CenteredCentroidFlow | 252 | train\n", "------------------------------------------------------------------------\n", "252 Trainable params\n", "0 Non-trainable params\n", "252 Total params\n", "0.001 Total estimated model params size (MB)\n", "INFO:lightning.pytorch.callbacks.model_summary:\n", " | Name | Type | Params | Mode \n", "------------------------------------------------------------------------\n", "0 | _potential_energy_func | OpenMMPotential | 0 | train\n", "1 | _loss_func | BoltzmannKLDivLoss | 0 | train\n", "2 | _flow | CenteredCentroidFlow | 252 | train\n", "------------------------------------------------------------------------\n", "252 Trainable params\n", "0 Non-trainable params\n", "252 Total params\n", "0.001 Total estimated model params size (MB)\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2f034474d13a49b6b8065d3caaadbbf3", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Training: | | 0/? [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def plot_free_energy_trajectory(analysis, ref_r0, target_r0):\n", " # Plot the the free energy trajectory.\n", " fig, ax = plt.subplots()\n", " bootstrap_sample_size = analysis['bootstrap_sample_size']\n", "\n", " # Plot true free energy.\n", " ref_df = - RT* np.log(target_r0 / ref_r0)\n", " ax.hlines(ref_df, bootstrap_sample_size[0], bootstrap_sample_size[-1], label='true', color='black', linestyle=':')\n", "\n", " # Plot mean bootstrap statistics.\n", " mean_bootstrap_df = [x['mean'].tolist() for x in analysis['df_bootstrap']]\n", " ax.plot(bootstrap_sample_size, mean_bootstrap_df, label='bootstrap mean', marker='.')\n", "\n", " # Plot bootstrap confidence interval.\n", " low_ci = [x['confidence_interval']['low'].tolist() for x in analysis['df_bootstrap']]\n", " high_ci = [x['confidence_interval']['high'].tolist() for x in analysis['df_bootstrap']]\n", " ax.fill_between(bootstrap_sample_size, low_ci, high_ci, alpha=0.25)\n", "\n", " ax.set_ylabel('$\\Delta f$ [kJ/mol]')\n", " ax.set_xlabel('# samples')\n", " ax.legend()\n", "\n", " print(f'Df = {analysis[\"df\"]} ({mean_bootstrap_df[-1]}) [{low_ci[-1]}, {high_ci[-1]}] kJ/mol')\n", "\n", "plot_free_energy_trajectory(mtfep_analysis, ref_r0=1.278, target_r0=1.71)" ] } ], "metadata": { "colab": { "provenance": [] }, "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.11.9" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "1d97d8e663044c0c8ad072aead4afe0f": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "HTMLModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "HTMLModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "HTMLView", "description": "", "description_tooltip": null, "layout": "IPY_MODEL_e7494aac57b54e1abf008651d19dd24c", "placeholder": "​", "style": "IPY_MODEL_75cd8e35ace64d398411ae1317940e92", "value": "Epoch 0: 100%" } }, "2f034474d13a49b6b8065d3caaadbbf3": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "HBoxModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "HBoxModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "HBoxView", "box_style": "", "children": [ "IPY_MODEL_1d97d8e663044c0c8ad072aead4afe0f", "IPY_MODEL_6a250960606e43ff8849ca4e7694b888", "IPY_MODEL_d1ba10853d6942eeb638882da3fc2dfe" ], "layout": "IPY_MODEL_affadb1606fd4833a61dbcaff8c550c8" } }, "6a250960606e43ff8849ca4e7694b888": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "FloatProgressModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "FloatProgressModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "ProgressView", "bar_style": "success", "description": "", "description_tooltip": null, "layout": "IPY_MODEL_cd0faa7c5a5e42fbafa26f8496d67d3e", "max": 625, "min": 0, "orientation": "horizontal", "style": "IPY_MODEL_8eff48866de6405ab2a3923d8d832caa", "value": 625 } }, "75cd8e35ace64d398411ae1317940e92": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "DescriptionStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "DescriptionStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "" } }, "8eff48866de6405ab2a3923d8d832caa": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "ProgressStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "ProgressStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "bar_color": null, "description_width": "" } }, "a2a4c9492ab74f2c880e4d28aff68338": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "DescriptionStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "DescriptionStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "" } }, "a52f60082adb49d0a32c9c27f6593bc6": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "affadb1606fd4833a61dbcaff8c550c8": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": "inline-flex", "flex": null, "flex_flow": "row wrap", "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": "100%" } }, "cd0faa7c5a5e42fbafa26f8496d67d3e": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": "2", "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "d1ba10853d6942eeb638882da3fc2dfe": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "HTMLModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "HTMLModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "HTMLView", "description": "", "description_tooltip": null, "layout": "IPY_MODEL_a52f60082adb49d0a32c9c27f6593bc6", "placeholder": "​", "style": "IPY_MODEL_a2a4c9492ab74f2c880e4d28aff68338", "value": " 625/625 [00:32<00:00, 18.98it/s, v_num=0]" } }, "e7494aac57b54e1abf008651d19dd24c": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } } } } }, "nbformat": 4, "nbformat_minor": 4 }