ljmartin/calculate solvation free energy .ipynb

## calculate solvation free energy .ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1f6a9d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import Chem\n",
    "from rdkit.Chem import Draw\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "from simtk.openmm import app\n",
    "from simtk import openmm\n",
    "from simtk import unit\n",
    "\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "\n",
    "from openff.toolkit.topology import Molecule\n",
    "from openmmforcefields.generators import GAFFTemplateGenerator\n",
    "\n",
    "TEMPERATURE = 298*unit.kelvin\n",
    "FRICTION = 1/unit.picosecond\n",
    "TIMESTEP = 2*unit.femtosecond"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2bbb30c2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0                                       mobley_1017962\n",
       "1                                         CCCCCC(=O)OC\n",
       "2                                     methyl hexanoate\n",
       "3                                                -2.49\n",
       "4                                                  0.6\n",
       "5                                                 -3.3\n",
       "6                                                 0.03\n",
       "7                                    10.1021/ct050097l\n",
       "8                             10.1021/acs.jced.7b00104\n",
       "9     Experimental uncertainty not presently availa...\n",
       "Name: 0, dtype: object"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.read_csv('/Users/ljmartin/Documents/projects/code_projects/freesolv.txt', sep=';', header=None).iloc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f7b9af7e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAACWCAIAAADCEh9HAAAABmJLR0QA/wD/AP+gvaeTAAAWl0lEQVR4nO3de3CU1f0G8GdzvxEEEwiXgoAI5RLCRQjhVlHB0UBtlTqtrqjFxdGZHYuMS23H1bHYjdUh1M5vTNVqwrRTg1INQULRekkIyD0XoKVJSGqAJBBCEkhIwu75/XGWJSQhu8nuvpfN8xn/WJOzu1+yOc++Ofue72sQQoCIiPorSO0CiIj0jTFKROQVxigRkVcYo0REXmGMEhF5hTFKROQVxigRkVcYo0REXmGMEhF5hTFKROQVxigRkVdC1C6AiPro/HnU1aG1FYMHY/RoRESoXdBAx6NRIp2or8dLL2HyZMTHY+pUzJmDiRMxeDCWLsXf/gaHQ+36Bi4DOzwR6UBeHn7+c1y86Pzf4GBER6Op6fqA+fPx6acYNkyV6gY4Ho0Sad6uXVi5EhcvIiICL76IkhJcvYrGRrS2YudO/OhHALB3LxYvviFYSSk8GiXStvp6TJ2K2loMHoxduzBvXtcBQuCFF7BpEwA8+ST+8hflaxzgeDRKpG3vvIPaWgDYtKmHDAVgMODNN7FgAQBkZeHUKUXLI8Yokda9/z4AjBqFxx+/6ZigIPzmNwBgt+PDD5Wpi1wYo0Qadvq08+hyxQoEB/c2ctkyxMQAQEGBEoVRJ4xRIg0rKnLemDHDzcjgYCQm3nAXUgpjlEjD6uudN4YPdz9YjmlogN3ux5KoG8YokYa1tDhvREa6HxwdDQAOx/V7kSIYo0QaJpc7AY+S8dIl4NqZ+aQgxiiRht16q/OGPOepd3LM0KEI4rxWFH/cRBrm+mTp6FE3I+12FBcDQFKSf0uibhijRBo2YgQmTACAHTvcfHC0ezcuXwaAhQuVKIw6YYwSaduTTwLA6dP4619vOkYIpKUBQEgInnhCmbrIhTFKpG1r1yI+HgDWrbvpOaGvvYavvwaAJ57AmDGKlUYSY5RI2+Li8N57CA5GfT0WL8Ybb1z/uEkIHD6Mhx+G1QoAt9+Ot95SsdIBix2eiPTgk09gNKK1FQCCgpCQgOho1NWhsdE5YNYs5OZixAgVaxywGKNEOlFdjddew7ZtOH/+hq8nJuKZZ/D00wjhNYHUwRgl0hWHAydOoK4OTU2Ii8O4cRg5Uu2aBjrGKJGGlZdjyxZMn46HHlK7FLopfsREpGFFRXj11d5OdXJ56y2kpV2/WBMpiIspRBp24QLQaUtoL/7wB9TWYvVqf1dE3fFolEjDZKM8T2K0oQEAhg71bz3UE8YokYbJGHUbjs3NaG/HoEEIC1OgKOqCMUqkYR4ejXp+0Ep+wBgl0rA+xSj/olcJY5RIw3g0qgeMUSINY4zqAWOUSMMYo3rAGCXSMHka05AhboYxRlXFGCXSqqYmdHQgNtb9aUyMUVUxRom0yvNwZIyqijFKpFWen8bk+Z5R8gPGKJFW8WhUJxijRFrV1xjl6fcqYYwSaRWPRnWCMUqkVR6G49WraGpCSAgGD1agKOqOMUqkVR7G6IULEAJDhsBgUKAo6o4xSqRV3MKkE4xRIq3y8DQmxqjaGKNEWsWjUZ1gjBJplYenMTFG1cYYJdIqHo3qBGOUSJM6OtDcjJAQxMa6GcmdoGpjjBJpkjyNaehQ96cxcQuT2hijRJrELUz6wRgl0iTGqH4wRom0KL+p6bnExA+nTXM/lDGqNsYokRadrKn5v+Li/PZ2tyPXxsZaZs1qZIyqJ0TtAoioBxcuXABwqwfhmHnoUFtb26v8iEk9PBol0qL6+noAQ92F46VLl9ra2qKjoyMiIhSpi3rAGCXSIhmjbo9GPRxGfsUYJdIiD/PR87/9yX8Yo0RaxKNRHWGMEmlRn2LU7RIq+RVjlEiLeDSqI4xRIi2Si55uDzMZo1owsM4bra2t/eabb0aOHDl79uzIyEi1yyHqWXNzc3t7e0xMTHh4eO8jGaNaMFBitLa2dtOmTW+//XZ7e3tYWFhCQsKHH364aNEitesi6oHn4cgY1YLA/6O+vLzcZDKNGTMmLS2ttbV18eLFo0aNqqioWLp06W9/+9t2DzbbESnM89OYGKNaEMgxWl5evnbt2smTJ7/77rtXr15NTU09ePDgl19+WVxcbLFYhBAbN26cPXv24cOH1a6U6AY8GtWXwIzR4uLixx9/fNKkSX/+858NBoPRaDx27Nj27dtnzZoFICIiwmazFRQU3HHHHaWlpcnJyRs2bOjo6FC7aiInz09j4un3WhBoMVpYWLhixYqkpKQtW7YEBwcbjcbjx49nZWVNnjy5y8jk5OSjR49aLBa73Z6WlrZw4cJ///vfqtRM1EVfj0Z53qi6AidGCwoKVqxYsWDBgtzc3KioKLPZXFFRkZWVdfvtt9/sLpGRkTabbdeuXWPGjNm/f//MmTPT0tIcDoeSZRN152GM2u32xsbGoKCgW265RZG6qGeBEKNffPFFSkrKokWLcnNzBw0aZDaby8rKNm/ePGrUKE/ufs8995SUlJhMpitXrmzYsGHRokVlZWX+rpmoFx7GaENDg8PhGDJkSFBQIExk/dLxT9/hcGzfvn3u3Ln33nvv3r174+LirFZrVVXV5s2bExIS+vRQsbGxGRkZO3fuHDVqVGFh4YwZMzZv3iyE8FPlRL3jFiZ90WWMOhyOrVu3Tp8+feXKlQcOHBg2bJjVai0vL3/llVeGDBnS74e97777SktLjUZjS0vL888/f99991VXV/uw7ADW0dGRlZV15513/vrXv/7Tn/7EdyAvMUZ1RuhKW1tbZmbmHXfcIYsfO3Zsenp6S0uLb58lOzs7Li4OwODBgzMyMnz74AGmpaXlj3/84w9+8AP5isi/LpcuXVpZWal2aTo2d+5cAN99913vw3JycgCkpqYqUxXdjG5i9MqVKxkZGa7pOn78+PT09CtXrvjp6Wpqan784x/L53rggQfOnDnjpyfSr+bm5vT09JEjR8qf0rRp0zIzM3NycuRX5DqJw+FQu0xdmjBhAoCysrLeh33wwQcAVq9erUhRdFM6iNEep2tHR4cCT52dnS1XCeLj4z/55BMFnlEXGhsbbTab6ySbGTNmZGZm2u12+d26urqHHnpIfmv58uXV1dXqVqtH8pP3Cxcu9D7szTffBLBu3TplqqKb0XSMdpmuSUlJnaerMqqqqu6++25ZwKpVq86fP6/ks2vNuXPnrFar6/SaBQsW5OTk9HjImZ2dLdfsbrnlFi6M9ElHR4fBYAgODnb7q/7SSy8B2LhxozKF0c1oNEbr6uo8nK4KcDgcGRkZMTExABISEnJyclQpQ101NTUWiyUqKqrzK9L7Xc6ePbty5UrXO9C5c+eUKVXv6urqAMTFxbkduXbtWgDvvPOOAlVRLzQXo92n6+7du9UuSgghKioqFi9eLKsyGo1NTU1qV6SQU6dOmc1meeFJg8GQmpq6b98+z++emZk5aNAgAMOGDfvHP/7hvzoDxokTJwBMmjTJ7ciHH34YQHZ2tgJVUS80FKNeTlcF2O329PR02QLytttu+9e//qV2Rf5VVlZmMplCQkLkR/Cyt0s/HqeysvKuu+5yHZbW19f7vNRAUlBQACAlJcXtSPlT/fLLLxWoinqhiRj11XRVRmlp6Zw5c2TWm0ymS5cuqV2R7xUXFxuNxuDgYAChoaFGo/HEiRPePKBcGImOjgYwYsSI3NxcX5UaePLy8iIiIrqcxrR///66urouHzolJiYCOHr0qLIFUlcqx6jPp6syOjo6bDZbWFgYgAkTJuTn56tdkc8cPnx41apVBoMBQFhYmNFoPHnypK8evLy8XLbKlu9Azc3NvnrkwNPe3i5v5Ofnp6amAkhOTh4+fPhnn33mGjNt2jQA//vf/1SqkZxUi9Hu0/W///2vWsX0T1FR0YwZMwCEhIRYLBb/ncSqDNd0BRAdHW02m7///nufP8vVq1dtNptcGBk3btzXX3/t86cIGDt37ly4cKF8RQYNGjR+/Hj5DrRmzRrX0rzdbufJuapTIUa7T1f9nlrY2tpqsVjk0fS0adMOHz6sdkX9kZ+fv3TpUtd0NZvN/t5uUFJSInu/ysPSy5cv+/Xp9MXhcOTk5MiNTABuvfVWq9VaX18vF0bkp69jx47lkqh2KBqj3afr2bNnlSzATwoLC+X+1NDQUIvF4vpzTON6nK5uT/n2lfb2dpvNFhoaCmDKlCkHDhxQ5nm1zG635+TkyDcYuenDarVevHix85hjx44F/NK87igRo3K63nnnnapMV2W0tLRYLBa5o3zevHkaX+HtMl1lb5cu01UZ3333neyoLRdG2tralK9BC9rb2zMzMydNmiRfkTFjxqSnp9/sIL3z0vz48eO//fZbhaulLvwbo3a7PTs7e8qUKapPV2Xs3r1b7vqX1ylReMOVJ3qcrj7v7dIncmFEvgMlJiYeOXJExWKUJ7vtuJqLjxs3Lj09vbW11e0di4uLk5KSAAQHBwfA0ryu+StG5XR1tWLSwnRVRmNjo8lkcu0d0M7nZnK6yp4Xrumqnbm3Z8+eiRMnyoURq9V69epVtSvyu0uXLqWnp7uai0+dOrWvzSLa29utVqtcmp86deqhQ4f8Vy31wvcx2qUVk9amqzI+//xz2UslKioqPT1d3c9Su/R26cd0Vcbly5fNZrM8eSM5Ofk///mP2hX5S1NTU3p6uqu5eGJiYmZmZr/fOfbu3Sv/vJALI3pZmg8kvoxRFVsxaVBDQ8Njjz0mfxTLly/3x8lDbjU1NdlsNldbX9mKSeMHert27Ro9ejSuXSlLgwsj3pC9XVzNxVNSUnzSLKLz0vzcuXOPHz/uk2oHgjNnzqxbt87LHee+idHeO6cNZK5GRwp3gPa8FZMGXbx40bUwcs8991RVValdkQ/U1tZardbY2NjOr4hvnyI/P18u2sileY2/X6qusrLy2WeflbvPFy1a5M1DeRujup6uylC4A3Q/WjFp044dO0aMGIFrHaDVLqf/KisrzWZzZGSk641h7969fnouuTQvF0ZSUlJ8uAMtkJSXl5vNZrkHRO4+379/vzcP6FWMbtiwwTVdly1b9s0333jzaIHN1QF62LBh27Zt88dTKDldlVFXV/fTn/5U/nPkpbHUrqhvuk9XZU6PlRdn1MjSvKaUlJQYjUZX+45Vq1YdO3bM+4f1KkafffbZwJiuyqisrHTtPvBto6Py8nKTySRPZVdyuiqjcwforKwstcvxiJyu8jN0OV0VXq9saGhwLYzce++93Hd/9OhRo9Eol4/l7nMffobpVYxWV1cXFRX5qpSBoEsH6O3bt3v5gJ17u6gyXZVx9uzZFStWuN6BtNwB+siRI12aRah4ygEvziiu7T6Xr0h4eLjJZPL5m4omGuUNNOXl5bIDtNzP178O0N2na8AvhLk6QA8fPlyDHaC7T1dVTs/oora29sEHH5TvQPfff//p06fVrkg5ndt3xMTEmM1mP/3zGaPq6NIB+quvvvL8vp1/ObQzXZVx6tSpzh2gNbKlOD8/33XBLjldtXYp2ezs7KFDh0aEhDQsXCg+/ljtcvxL7j6fN2+efEViY2MtFotfm4UzRtVUWlo6e/ZseNzoSPvTVQGdO0CPGTPmiy++ULEShaerN6qqqvJNJgEIQDz6qNDGO5BvyWYRck7hWm+XhoYGfz8vY1RlndtM/PCHP+zxxIsurZg0Pl2VUVZWJntxqtIBusfpqoNmEQ6HyMgQMTECEAkJolMHaL2Tu89lmxu57GOz2RRrwMgY1YQuHaBdjY486Zw2YMl3IFcHaGXOt1N3uvpGRYVYssR5WGo0Cp1fnLFLb5fbbrvNw94uPsQY1YrOHaCnT5++f/9+zzunDWQlJSUzZ86UJyqYzWb/dW/QwnT1GXlYGhUlADF2rNDnxRm79HaZMGFCRkaGKrvPGaPa8u2338prRcg9avKX47333huwjTg90bkD9NSpU31+zqx2pquPHTsm5swRgDAYhMkkdNQB+uJF8bvfLZ4zR74iSUlJW7duVXH3OWNUc5qbm5955pmXX35Zs62YtGnfvn0+7wDt21ZMWtTRIWw2ERYmADF+vNB+B+jz54XVKoYMEcCWBQtmzpyZnZ2t+jYtxqh2qf7LoTtdOkB7c+VhP7Vi0qiiIpGUJAARHCwsFqHNtpZnzogXXnB+PgaIu+5yaOZqVIxRCjQFBQVyBbN/jY4UaMWkRe3twmoVwcECENOmCU11gK6qEmaziIx0BuiCBUIzASoxRikAde4APX/+fA+3Y1ZVVXXp7VJYWOjvUrVl714xaZIAREiIsFiE6h2gKyqE2SzCw50LuKmpwrtWTH7CGKWAlZeX52EH6IqKCt92TtOxlhZhsYigIAGIuXOFWhdnLC0VRqMICRGACAoSq1YJX7Ri8hPGKAWyLh2gu/ek8FPnNN3LzxcTJghAREQIm00o+ana0aPCaHQuL4SGCqNRaP5yMoxRCnxbt26Nj4/HjR2gO3dOCw0NVbcVkxY1NgrX5tGUFKFA45uCApGaKgwGAYjwcGEyCZ309zMIIUAU6Gpra9euXfvZZ58BWLJkSVBQ0FdffQUgMjLy6aefXr9+vesijHSDvDysWYPTpxEVhddfh9kMg6HnkU1NOHEC9fUAEBeHKVMQE+PpsxQUIC0NubkAEBODp56CxYJrV3XTAbVznEg5stGRvGiYXzunBZSGhuuHpcuWie7txHJzxZIlznVM139hYeLuu8U//+nmwXfvFsnJzrvExgqLReiwWQSPRmlgqa6u/v777w8ePPjoo4+6LsJI7n30EZ57DvX1GDcOJ08iJAQA2tqwejU++sg5JiICo0bB4cCZM2hrc35x9Wq8+y5CQ294NIcDO3bg1Vdx6BAAxMXhuefw/PO4dlU3nVE7x4lIJ2pqxIMPii1brn/lJz9xHkjOni0+/1y4do61topPPxVTpji/+4tfXL9Le7vIzBSTJzu/NXy4sNmEzptF8GiUiPrl/fexZg0ArFiBjz9GWFjXAZcu4f77kZ8PAH//Ox55BADmzsWBAwAwfjwsFqxejfBwJav2B8YoEfWdw4FJk1BWhmHDcOIEbrY8Ul2NKVPQ3IzERBQVAcDGjfjgA7z4Ip56yrkyoH9BahdARDq0Zw/KygDgl7+8aYYCGD0ajz0GAMXFOHwYANavx8mTMJkCJkPBGCWi/igsdN544AE3I1eudN7YswcAwsMRFGixE2j/HiJSQmmp80ZiopuRM2Y4bxw75sd6VMUYJaK+u3ABAMLDMWiQm5Hx8TfcJRAxRomo71pbAXj0IXtIiPOk0cuX/VuSehijRNR3gwcDniXjlSvo6ACAaz2wAw9jlIj6Li4OAOx21Na6GXnmzA13CUSMUSLqu5kznTeOHHEzUp7nBODadcIDD2OUiPpu8WLnjW3b3Iz8+GMAMBiu3yXgcBcTEfXL/PnYtw9RUSgtxbhxPY8pLcWsWejowPLlyMtTtj7l8GiUiPrl1VdhMKClBT/7GRoaehhQU4NHHkFHB4KD8corSpenIMYoEfXLsmX41a8A4OBBJCYiIwM1Nc5vnT6Nt99GUhKOHwcAqxXJyarV6X/8o56I+ksIvPwyfv972O3Or0RGQghcueL839BQvP461q9Xq0BlMEaJyDulpXjjDeTl4dy5619MSEBqKl58ERMnqleZQhijROQLDgfOnsX58wAQH48RI2561aaAwxglIvIKP2IiIvIKY5SIyCuMUSIirzBGiYi8whglIvIKY5SIyCuMUSIirzBGiYi8whglIvIKY5SIyCuMUSIir/w/RHXVsazbAKAAAADFelRYdHJka2l0UEtMIHJka2l0IDIwMjEuMDMuMwAAeJx7v2/tPQYg4AFiRgYI4ARiDiBuYGRjSACJM0NoJibctAaQZmbhgNBMMBqmnxtoNiMTAxMzUA0DCysDK5sGEyu7AjsHgwiDOMxaBs5JiTMPcM7O2QfiHAyIOTBz/XU7EPv0e54D621vgMXnfo3cf7cuGyxeKXPc7o2Y7H4Q+9nqe/bKVXK2YL3n3thn5R2xB7FfbrNzYOL3A6tZNqXL4YyjDNgcMQASyi0FKSV7SAAAANx6VFh0TU9MIHJka2l0IDIwMjEuMDMuMwAAeJydkk1qAzEMhfdzCl0gQn8eS+tMVyEtdNE7dJ/7E409MbNIoVgY8x7y+5CNF9jre7v9PmCUbMsCEAAOQG9XRMCPEFGeg4thDZVdEYoy7wzC7BJc4S/EeTWKoplYz5qwzlEEhYT7LCYySSH0wscsopmdoWRC3aPzVtZ1jsJYNV5Z9TpL8ULcVVE7z/L1f4pg1Dhu5GQxRzFUpbVTuMjUu+RVpO27SaXDpLJhUpVhUq3DpKoteRgfHW9/v5k7wMfntjwBhZZ8SJcm3fYAAACYelRYdFNNSUxFUyByZGtpdCAyMDIxLjAzLjMAAHicFY5LDsUwCAOv8patlCIwnwRFXfUAvVAP/8DL0djwPJ3jfs/3+X3HZTRTJcbFBBWOsS8lM5gOJoMgUAgEhlhbBqi1xrRcjK08IHMWY9JVyfZCNFTGFpqaNcLEuiJ5NlrOdUrIa6r2QTkzpWuLLeudbaTK4Y3E4enj/P4uXCYLVP570QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<rdkit.Chem.rdchem.Mol at 0x1a8177880>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smi = 'CCCCCC(=O)OC'\n",
    "Chem.MolFromSmiles(smi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "add550c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an OpenFF Molecule object for benzene from SMILES\n",
    "molecule = Molecule.from_smiles(smi)\n",
    "# Create the GAFF template generator\n",
    "gaff = GAFFTemplateGenerator(molecules=molecule)\n",
    "# Create an OpenMM ForceField object with AMBER ff14SB and TIP3P with compatible ions\n",
    "from simtk.openmm.app import ForceField\n",
    "forcefield = ForceField('amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml')\n",
    "# Register the GAFF template generator\n",
    "forcefield.registerTemplateGenerator(gaff.generator)\n",
    "# You can now parameterize an OpenMM Topology object that contains the specified molecule.\n",
    "# forcefield will load the appropriate GAFF parameters when needed, and antechamber\n",
    "# will be used to generate small molecule parameters on the fly.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "003a8a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = molecule.to_rdkit()\n",
    "AllChem.EmbedMolecule(m)\n",
    "Chem.MolToPDBFile(m, 'mol.pdb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "bca7d4d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#instantiate a Modeller object using the topology and xyz coordinates,\n",
    "pdb = app.PDBFile('mol.pdb')\n",
    "pdb.topology.setPeriodicBoxVectors(np.eye(3)*2.4)\n",
    "\n",
    "\n",
    "#and add salt water:\n",
    "modeller = app.Modeller(pdb.topology, pdb.positions)\n",
    "modeller.addSolvent(forcefield, \n",
    "                    ionicStrength = 0.15*unit.molar,\n",
    "                   padding=1.2*unit.nanometer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "03942ab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_cnb(original_nbforce, system, solute_idx):\n",
    "    \"\"\"\n",
    "    original_nbforce: the nonbonded force that openmm sets up along with the system automatically\n",
    "    system: the system object\n",
    "    longRangeDispersion: bool, i.e. whether or not to apply long range dispersion correction\n",
    "    \"\"\"\n",
    "    \n",
    "\n",
    "    ##Implement Coulomb force:\n",
    "    \n",
    "    ONE_4PI_EPS0 = 138.935456\n",
    "    epsilon_solvent = original_nbforce.getReactionFieldDielectric()\n",
    "    r_cutoff = original_nbforce.getCutoffDistance()\n",
    "    k_rf = r_cutoff**(-3) * ((epsilon_solvent - 1) / (2*epsilon_solvent + 1))\n",
    "    c_rf = r_cutoff**(-1) * ((3*epsilon_solvent) / (2*epsilon_solvent + 1))\n",
    "    \n",
    "    \n",
    "    energy_expression  = \"select(condition, 1, electro_scalingFactor)*electrostatics;\"\n",
    "    energy_expression += \"condition = delta(soluteFlag1-soluteFlag2);\" #solute must have flag int(1)\n",
    "    \n",
    "    energy_expression += \"electrostatics=ONE_4PI_EPS0*chargeprod*(r^(-1) + k_rf*r^2-c_rf);\"\n",
    "    energy_expression += \"chargeprod = charge1*charge2;\"\n",
    "    energy_expression += \"k_rf = %f;\" % (k_rf.value_in_unit_system(unit.md_unit_system))\n",
    "    energy_expression += \"c_rf = %f;\" % (c_rf.value_in_unit_system(unit.md_unit_system))\n",
    "    energy_expression += \"ONE_4PI_EPS0 = %f;\" % ONE_4PI_EPS0 \n",
    "    \n",
    "    electrostatics_force = openmm.CustomNonbondedForce(energy_expression)\n",
    "    electrostatics_force.addGlobalParameter('electro_scalingFactor', 1)\n",
    "    electrostatics_force.addPerParticleParameter('charge')\n",
    "    electrostatics_force.addPerParticleParameter('soluteFlag')\n",
    "\n",
    "    electrostatics_force.setNonbondedMethod(original_nbforce.getNonbondedMethod())\n",
    "    electrostatics_force.setCutoffDistance(original_nbforce.getCutoffDistance())\n",
    "    electrostatics_force.setUseLongRangeCorrection(False) #aka DispersionCorrection in the NonbondedForce\n",
    "\n",
    "    for index in range(system.getNumParticles()):\n",
    "        soluteFlag = 1 if index in solute_idx else 0\n",
    "        [charge, _, __] = original_nbforce.getParticleParameters(index)\n",
    "        electrostatics_force.addParticle([charge, soluteFlag])\n",
    "        \n",
    "    ##Implement LJ forces:\n",
    "    #softcore\n",
    "    energy_expression = \"sterics;\"\n",
    "    energy_expression += \"sterics=scaling*4*epsilon*x*(x-1.0);\"\n",
    "    energy_expression += \"x = (sigma/reff_sterics)^6;\"\n",
    "    energy_expression += \"reff_sterics = sigma*(0.5*(1.0-scaling) + (r/sigma)^6)^(1/6);\"\n",
    "    energy_expression += \"epsilon = sqrt(epsilon1*epsilon2);\"\n",
    "    energy_expression += \"sigma = 0.5*(sigma1+sigma2);\"\n",
    "    energy_expression += \"scaling = select(condition, 1, lj_scalingFactor);\"\n",
    "    energy_expression += \"condition = delta(soluteFlag1-soluteFlag2);\" #solute must have flag int(1)\n",
    "    \n",
    "    \n",
    "    lj_force = openmm.CustomNonbondedForce(energy_expression)\n",
    "    lj_force.addGlobalParameter('lj_scalingFactor', 1)\n",
    "    lj_force.addPerParticleParameter('sigma')\n",
    "    lj_force.addPerParticleParameter('epsilon')\n",
    "    lj_force.addPerParticleParameter('soluteFlag')\n",
    "    \n",
    "    lj_force.setNonbondedMethod(original_nbforce.getNonbondedMethod())\n",
    "    lj_force.setCutoffDistance(original_nbforce.getCutoffDistance())\n",
    "    lj_force.setUseLongRangeCorrection(True) #aka DispersionCorrection in the NonbondedForce\n",
    "\n",
    "\n",
    "    \n",
    "    for index in range(system.getNumParticles()):\n",
    "        soluteFlag = 1 if index in solute_idx else 0\n",
    "        [_, sigma, epsilon] = original_nbforce.getParticleParameters(index)\n",
    "        lj_force.addParticle([sigma, epsilon, soluteFlag])\n",
    "    \n",
    "    \n",
    "    ## Exceptions:\n",
    "    \n",
    "    energy_expression = \"4*epsilon*((sigma/r)^12 - (sigma/r)^6) + ONE_4PI_EPS0*chargeprod/r;\"\n",
    "    energy_expression += \"ONE_4PI_EPS0 = {:f};\".format(ONE_4PI_EPS0)  # already in OpenMM units\n",
    "    custom_bond_force = openmm.CustomBondForce(energy_expression)\n",
    "    custom_bond_force.addPerBondParameter('chargeprod')\n",
    "    custom_bond_force.addPerBondParameter('sigma')\n",
    "    custom_bond_force.addPerBondParameter('epsilon')\n",
    "\n",
    "    for index in range(original_nbforce.getNumExceptions()):\n",
    "        idx, jdx, c, s, eps = original_nbforce.getExceptionParameters(index) #c,s,eps = charge, sigma, epsilon\n",
    "        lj_force.addExclusion(idx, jdx)\n",
    "        electrostatics_force.addExclusion(idx, jdx)\n",
    "        c_value = c/unit.elementary_charge**2\n",
    "        eps_value = eps/(unit.kilojoule/unit.mole)\n",
    "        if c_value != 0 or eps_value!=0:\n",
    "            custom_bond_force.addBond(idx, jdx, [c, s, eps])\n",
    "        \n",
    "    return electrostatics_force, lj_force, custom_bond_force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "50adbc9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Quantity(value=-25090.217760919884, unit=kilojoule/mole)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "\n",
    "system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.CutoffPeriodic,\n",
    "         nonbondedCutoff=0.9*unit.nanometer, constraints=app.HBonds)\n",
    "\n",
    "solute = [c for c, i in enumerate(modeller.topology.atoms()) if i.residue.name=='UNL']\n",
    "for c, i in enumerate(system.getForces()):\n",
    "    if isinstance(i, openmm.NonbondedForce):\n",
    "        electro, lj, cbf = create_cnb(i, system, solute) #<---- note LRD is TRUE now !\n",
    "        system.removeForce(c)\n",
    "        \n",
    "system.addForce(electro)\n",
    "system.addForce(lj)\n",
    "system.addForce(cbf)\n",
    "\n",
    "\n",
    "integrator = openmm.LangevinIntegrator(TEMPERATURE, FRICTION, TIMESTEP)\n",
    "platform = openmm.Platform.getPlatformByName('OpenCL')\n",
    "simulation = app.Simulation(modeller.topology, system, integrator, platform)\n",
    "simulation.context.setPositions(modeller.positions)\n",
    "\n",
    "simulation.context.getState(getEnergy=True).getPotentialEnergy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "fefc5494",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def scaling_function(simulation, level):\n",
    "\n",
    "    #Free Energy calculation control parameters\n",
    "    ele_schedule = [1, 0.75, 0.5, 0.25, 0, 0,0,0,0,0,0,0,0,0,0,0]\n",
    "    vdw_schedule = [1,1,1,1,1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0]\n",
    "\n",
    "    simulation.context.setParameter('lj_scalingFactor', vdw_schedule[level])\n",
    "    simulation.context.setParameter('electro_scalingFactor', ele_schedule[level])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c702355",
   "metadata": {},
   "source": [
    "# Estimate free energy of solvation with MBAR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "03b83d9c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-36-d51f46f6e760>:12: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0\n",
      "Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`\n",
      "  for k in tqdm.tqdm_notebook(range(nstates)):\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f5d53e1be313477e8dbe944f159978c0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/16 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-36-d51f46f6e760>:13: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0\n",
      "Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`\n",
      "  for iteration in tqdm.tqdm_notebook(range(niterations)):\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5869d800dea44faa99eb83e479408f40",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4c7174cbe15b4cba96ca93c65ec71432",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8646c8c98d7e4255ab5424b60828d9d2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c90f113dd6624e30b0050d4bcc8b227f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a786cac3b1854e569c9de5527856d6a9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fab0cf49cb644f3abc6ccfc15882075b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "088817d2ce7c46d59a5b415991a1ad48",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "529ba471ed294bfc8c09bd6c53b25adb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ea73f3b30df34f2e87a03db71bf5a42d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7eeb6816b18a43dcbf57f07595df6cad",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b5ccac5976ca4acb9c39cd9b43cc56c6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "78494362371c42cd942e7e51839f44c3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "807fbd14b16d48e89ef983357dee78a5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "33dfad9877734bae96e455e52fd2b8f2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d2cac6e2182477eb8007fc2d99cb826",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "383f4d3af9dd40858c740f658d1ac1dd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import tqdm\n",
    "\n",
    "nsteps = 500 # number of steps inbetween samples\n",
    "niterations = 100 # number of samples to collect at each level. \n",
    "\n",
    "nstates = 16\n",
    "\n",
    "u_kln = np.zeros([nstates,nstates,niterations], np.float64)\n",
    "kT = unit.AVOGADRO_CONSTANT_NA * unit.BOLTZMANN_CONSTANT_kB * simulation.integrator.getTemperature()\n",
    "\n",
    "\n",
    "for k in tqdm.tqdm_notebook(range(nstates)):\n",
    "    for iteration in tqdm.tqdm_notebook(range(niterations)):\n",
    "        # Set the separation distance:\n",
    "        scaling_function(simulation, k)\n",
    "        # Sample in MD-space for a bit:\n",
    "        simulation.step(nsteps)\n",
    "        # Compute what the energy of the current configuration would be at ALL separations,\n",
    "        # this is used to calculate energetic overlap between the states. \n",
    "        for l in range(nstates):\n",
    "            scaling_function(simulation, l)\n",
    "            u_kln[k,l,iteration] = simulation.context.getState(getEnergy=True).getPotentialEnergy() / kT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "ca3a23e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate free energy of Lennard-Jones particle insertion\n",
    "from pymbar import MBAR, timeseries\n",
    "# Subsample data to extract uncorrelated equilibrium timeseries\n",
    "N_k = np.zeros([nstates], np.int32) # number of uncorrelated samples\n",
    "for k in range(nstates):\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\n",
    "# Compute free energy differences and statistical uncertainties\n",
    "mbar = MBAR(u_kln, N_k)\n",
    "[DeltaF_ij, dDeltaF_ij,] = mbar.getFreeEnergyDifferences()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "d63a86f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(16, 16)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DeltaF_ij.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "b59c2ca4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "b594d756",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1ad5be7f0>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD4CAYAAAAjDTByAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAARM0lEQVR4nO3de4yc5XXH8e+Zy67tNWAcE7CxWwNCRDRqC7IQSSoalRIRinAq9Q9Q07pNpChSaaFKlDhCavJn07TpNUpEAyltEUhNoEERtCBCFFUKNOCaizEJl1IwNhhCZLPrtXdn5vSPeWnGy6x3n/NebPf5faTVzs68Z9+zz8zZd25njrk7IpKf1vFOQESODxW/SKZU/CKZUvGLZErFL5KpTpM7a6+e8s7atelxk73kmBWd9BiArvWTY1oWe8VkwmI5RuIif9cwLhRGm/TAViAGwElf/9ngq1xvDVaE4vqefpz1wHoc2HuIQz89sqzARou/s3YtGz51Y3Lcaef9NDnmPe/anxwDsH7FgeSYyVasiH9u8iehuM0TbyTHbOik/10AZ7Vj/zROa00kx0xaN7Svvg+SY3bNz4X29d2Z94TiDvRWJcfMezs55hvXPbTsbXW3XyRTKn6RTJUqfjO70sx+ZGbPmdn2qpISkfqFi9/M2sBXgA8DFwLXmdmFVSUmIvUqc+S/BHjO3V9w9zngTmBrNWmJSN3KFP/ZwMsjP+8pzjuKmX3CzB41s0f70zMldiciVSpT/ONeS3zHi6fufrO7b3H3Le3VUyV2JyJVKlP8e4BNIz9vBPaWS0dEmlKm+H8InG9m55jZBHAtcE81aYlI3cLv8HP3npldD/w70AZudfddlWUmIrUq9fZed78XuLeiXESkQXqHn0immu3qOwynPpv+/+YApyfH/Odb6Y0UACtXHUmO6bTSG0sA1k4dCsWdsXI6OeasFQdD+1rXTd8XwGnt2eSYVa30tQfoB45hTx/aENrXjjc2Lb3RGLPz6U1L/UGgq2/uB8veVkd+kUyp+EUypeIXyZSKXyRTKn6RTKn4RTKl4hfJlIpfJFMqfpFMqfhFMqXiF8mUil8kU4029rR6MPVa+gSYudPSJ5ccZjI5BmB6ZWBqTDs2+ml6Kjb66c2p9KalV1eeGtrXmsn0Bh2A1d30Jp2V7fnQvgae3gDz/MF1oX3tfWNNKK5/JP02TKBfrD+//P3oyC+SKRW/SKZU/CKZKjOxZ5OZPWRmu81sl5ndUGViIlKvMk/49YBPufsOMzsFeMzMHnD3pyvKTURqFD7yu/s+d99RnH4L2M2YiT0icmKq5DG/mW0GLgIeGXPZ/43rmj8S+zw4Eale6eI3s9XAt4Ab3f0dnxI5Oq6rO7m67O5EpCKlit/MugwL/3Z3v6ualESkCWWe7TfgFmC3u3+5upREpAlljvwfAH4H+DUz21l8XVVRXiJSszKz+v6D8WO6ReQkoHf4iWSq0a4+6zvd6fSuvu5MekdU71DsTkkv8P/Qg119vVZs+Wcs1rEY0RvEjg+zvfTuyMlOL7SvSFffgdlYR2V/Jnad2VxgHSNT4BJGfOnIL5IpFb9IplT8IplS8YtkSsUvkikVv0imVPwimVLxi2RKxS+SKRW/SKZU/CKZUvGLZKrRxh4cWr30JhhL7wUKxQBYr7kuZe/H9jUINNvM9wPjoog39swN0vdn/ViDlAcae6J/V0rjzFEC17VFGnsS6MgvkikVv0imVPwimario7vbZvZfZvadKhISkWZUceS/geG0HhE5iZT93P6NwG8AX68mHRFpStkj/18BnyH2aWMichyVGdpxNbDf3R9bYrufzeqbn4nuTkQqVnZoxzVm9iJwJ8PhHf+8cKOjZvV1p0rsTkSqVGZE9+fcfaO7bwauBb7r7h+tLDMRqZVe5xfJVCXv7Xf37wHfq+J3iUgzdOQXyVTz47oOziXHdacnkmPmV8e6ryKdVIEGNgB6xAIjzW+zwRdj+/3guK5u+riuiU6sFXMQWI9D07GRZ+3p2Hq05prp6kvpZtWRXyRTKn6RTKn4RTKl4hfJlIpfJFMqfpFMqfhFMqXiF8mUil8kUyp+kUyp+EUypeIXyZSKXyRTDXf1DWgdOJQcNzGzKjlmfjrWMWeBDrFBOzi/LRgW6QYcBObZARwJdvXNT6TneLgdaz30yHU2k951CDA5E+3qS48JdfUlxOjIL5IpFb9IplT8IpkqO7FnjZl908yeMbPdZva+qhITkXqVfcLvr4F/c/ffMrMJIP2ZORE5LsLFb2anApcBvwfg7nNA4DlNETkeytztPxd4HfhGMaL762b2jpE8o+O65vrpL/OJSD3KFH8HuBj4qrtfBMwA2xduNDqua6KtRwUiJ4oyxb8H2OPujxQ/f5PhPwMROQmUmdX3KvCymV1QnHU58HQlWYlI7co+2/+HwO3FM/0vAL9fPiURaUKp4nf3ncCWalIRkSY12tjDoI+9NZMc1p1emxzTOSX6iCa9AWYQXUVrriGoF3yEFxmFBdDvRdYxtrNIY09rNrYendlQGO0jgaBIn5Mae0RkKSp+kUyp+EUypeIXyZSKXyRTKn6RTKn4RTKl4hfJlIpfJFMqfpFMqfhFMqXiF8mUil8kU8129fUHDKbTu/o6s73kmO5sbFyXB/4dtmK7woNjviJxkb9rKBbog0COvWALYSCsfTi29u1oV99cepKhcV0Ju9GRXyRTKn6RTKn4RTJVdlzXH5vZLjN7yszuMLMVVSUmIvUKF7+ZnQ38EbDF3d8LtIFrq0pMROpV9m5/B1hpZh2Gc/r2lk9JRJpQ5nP7XwH+HHgJ2AcccPf7F2531LguPxzPVEQqVeZu/+nAVuAcYAMwZWYfXbjdUeO69JSAyAmjzN3+Xwf+291fd/d54C7g/dWkJSJ1K1P8LwGXmtkqMzOG47p2V5OWiNStzGP+RxgO59wBPFn8rpsryktEalZ2XNfngc9XlIuINEjv8BPJVKNdfe4DfDa9Lap1pJ8cE+miAmh3A91owa6+wVwwbiI9pjUX7CAMjhOMjJkjuI6RJKPrEb1dtQLXdairT7P6RGQpKn6RTKn4RTKl4hfJlIpfJFMqfpFMqfhFMqXiF8mUil8kUyp+kUyp+EUypeIXyVSz47ocvJc+esvm0xt7WvPBBozAyCgPTpmy9KUo4gKNLMF9DbqxuEiOUZEGmOh6NBkX+btS6MgvkikVv0imVPwimVqy+M3sVjPbb2ZPjZy31sweMLNni++n15umiFRtOUf+fwCuXHDeduBBdz8feLD4WUROIksWv7t/H3hzwdlbgduK07cBH6k2LRGpW/Qx/5nuvg+g+P7uxTYcHdc1z5Hg7kSkarU/4Tc6rqvLZN27E5Flihb/a2a2HqD4vr+6lESkCdHivwfYVpzeBny7mnREpCnLeanvDuAHwAVmtsfMPg78KXCFmT0LXFH8LCInkSXf2+/u1y1y0eUV5yIiDdI7/EQy1WxXH4AFur0CMR79txZIzyN/U3Bf0bjo2K14joFWx+C+Itd1dD2avF2F136ZdOQXyZSKXyRTKn6RTKn4RTKl4hfJlIpfJFMqfpFMqfhFMqXiF8mUil8kUyp+kUyp+EUy1Whjj3U7dM44Mznu0LqVyTGza9vJMQC9VYEmotiumF8di+tNpTfN9FbGZooNVgZnRnXT46wTyzEyLm2O2Byy9lzseNmOfHxlYOkHCbdFHflFMqXiF8mUil8kU9FxXV8ys2fM7Akzu9vM1tSapYhULjqu6wHgve7+i8CPgc9VnJeI1Cw0rsvd73f3XvHjw8DGGnITkRpV8Zj/Y8B9i104Oq5rbjBbwe5EpAqlit/MbgJ6wO2LbTM6rmuilf56vYjUI/wmHzPbBlwNXO4eeZuFiBxPoeI3syuBzwK/6u6Hqk1JRJoQHdf1d8ApwANmttPMvlZzniJSsei4rltqyEVEGqR3+IlkqtGuPp/oMr85vatvekN6mofOjM066q1Kf+4y2tUX6c4D8Kne0hst0F6VHgMwtWI+FDfZTY+b6PRD++oP0o9hb05OhfY1y4pQXOtI+u0xMvHME5oVdeQXyZSKXyRTKn6RTKn4RTKl4hfJlIpfJFMqfpFMqfhFMqXiF8mUil8kUyp+kUyp+EUypeIXyVSjXX20jMFkegtcfyJ9V4NATDTO28HuvMnYHDybSI+bmIh19UW68wBWdtP3N9mJ5Tjw9I656YnJ0L4OTwRnHgZiLNDkmLIUOvKLZErFL5Kp0Liukcs+bWZuZuvqSU9E6hId14WZbQKuAF6qOCcRaUBoXFfhL4HPAPrMfpGTUOgxv5ldA7zi7o8vY9ufjeuam4nsTkRqkPxSn5mtAm4CPrSc7d39ZuBmgFNP3ah7CSIniMiR/zzgHOBxM3uR4YTeHWZ2VpWJiUi9ko/87v4k8O63fy7+AWxx9zcqzEtEahYd1yUiJ7nouK7RyzdXlo2INEbv8BPJVLPjugBvpTdhxGKSQ8Jx0X2F5jEB1kqPawViANrBuE4rvZWlY7FGp4Gl3z6i6+HBOAK34RA19ojIUlT8IplS8YtkSsUvkikVv0imVPwimVLxi2RKxS+SKRW/SKZU/CKZUvGLZErFL5IpFb9Ipsy9uY/VM7PXgf9Z5OJ1wInwaUDK42jK42gneh4/7+5nLOcXNFr8x2Jmj7r7FuWhPJRHM3nobr9IplT8Ipk6kYr/5uOdQEF5HE15HO3/TR4nzGN+EWnWiXTkF5EGqfhFMtVo8ZvZlWb2IzN7zsy2j7nczOxvisufMLOLa8hhk5k9ZGa7zWyXmd0wZpsPmtkBM9tZfP1J1XmM7OtFM3uy2M+jYy6vdU3M7IKRv3OnmR00sxsXbFPbepjZrWa238yeGjlvrZk9YGbPFt9PXyT2mLenCvL4kpk9U6z73Wa2ZpHYY16HFeTxBTN7ZWT9r1okNm093L2RL6ANPA+cC0wAjwMXLtjmKuA+hh9AfCnwSA15rAcuLk6fAvx4TB4fBL7T0Lq8CKw7xuW1r8mC6+hVhm8UaWQ9gMuAi4GnRs77M2B7cXo78MXI7amCPD4EdIrTXxyXx3Kuwwry+ALw6WVcd0nr0eSR/xLgOXd/wd3ngDuBrQu22Qr8ow89DKwxs/VVJuHu+9x9R3H6LWA3cHaV+6hY7Wsy4nLgeXdf7F2YlXP37wNvLjh7K3Bbcfo24CNjQpdzeyqVh7vf7+694seHGQ6lrdUi67EcyevRZPGfDbw88vMe3ll0y9mmMma2GbgIeGTMxe8zs8fN7D4z+4W6cmA4y+R+M3vMzD4x5vIm1+Ra4I5FLmtqPQDOdPd9MPxnzchg2BGN3laAjzG8BzbOUtdhFa4vHn7cusjDoOT1aLL4x80SWfg643K2qYSZrQa+Bdzo7gcXXLyD4V3fXwL+FvjXOnIofMDdLwY+DPyBmV22MNUxMZWviZlNANcA/zLm4ibXY7mavK3cBPSA2xfZZKnrsKyvAucBvwzsA/5iXJpjzjvmejRZ/HuATSM/bwT2BrYpzcy6DAv/dne/a+Hl7n7Q3aeL0/cCXTNbV3Uexe/fW3zfD9zN8O7bqEbWhOENd4e7vzYmx8bWo/Da2w9tiu/7x2zT1G1lG3A18NtePLheaBnXYSnu/pq79919APz9Ir8/eT2aLP4fAueb2TnFUeZa4J4F29wD/G7xDPelwIG37/5VxcwMuAXY7e5fXmSbs4rtMLNLGK7TT6rMo/jdU2Z2ytunGT7B9NSCzWpfk8J1LHKXv6n1GHEPsK04vQ349phtlnN7KsXMrgQ+C1zj7ocW2WY512HZPEaf4/nNRX5/+npU8QxlwjOZVzF8dv154KbivE8CnyxOG/CV4vIngS015PArDO8OPQHsLL6uWpDH9cAuhs+YPgy8v6b1OLfYx+PF/o7XmqxiWMynjZzXyHow/IezD5hnePT6OPAu4EHg2eL72mLbDcC9x7o9VZzHcwwfR799O/nawjwWuw4rzuOfiuv+CYYFvb6K9dDbe0UypXf4iWRKxS+SKRW/SKZU/CKZUvGLZErFL5IpFb9Ipv4XiH8c4EYwuLkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(DeltaF_ij)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "377301ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1acb5e610>]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlDklEQVR4nO3dd3zV5d3/8dcni0xIIAk7hL1EBAIi4AAcCKi1rXXRuhAntdNqbW/v4a+3Vm21raM4igtsFWrFAYqiqCAQwLA3AUKABEJYIWRdvz8SeyMGEnJO8j3j/Xw88kjyPSc57yTkzTfX+V7XZc45REQk+ER4HUBERBpGBS4iEqRU4CIiQUoFLiISpFTgIiJBKqopHyw1NdVlZmY25UOKiAS9pUuX7nXOpZ14vEkLPDMzk+zs7KZ8SBGRoGdm22o7riEUEZEgVWeBm9mLZlZgZqtOOD7ZzNab2Woz+33jRRQRkdrU5wx8KjDm+ANmNhK4AjjTOdcXeMz/0URE5FTqLHDn3Hyg6ITDdwAPO+eO1dynoBGyiYjIKTR0DLwHcK6ZLTKzT81s8MnuaGaTzCzbzLILCwsb+HAiInKihhZ4FJACDAV+CfzDzKy2OzrnpjjnspxzWWlp37oKRkREGqihBZ4HzHTVFgNVQKr/YomISF0aeh34W8Ao4BMz6wHEAHv9FUpETt/ew8f4ZH0hew6WkhIfQ8uEb74kx0UTEVHrH8oSpOoscDObDlwApJpZHvAg8CLwYs2lhWXADU4Li4s0Kecca3Yd5OO1BXy0roCcvGJO9VsYYZASH0PK16UeH0PLxJrXCbW/xEZHNt0XJKetzgJ3zl17kpsm+DmLiNThaFklX2zay0frCpi3roDdB0sB6N+hBT8Z3YPRvdPpmpZI8dEy9h0uY39JGUVH/u/tfUfK2H+k+vXmwsMsya0+XnWS4o+PiaRbeiJ//eEg2raIa8KvVOqjSafSi8jpy9tfwrx1BXy8roAFm/dxrKKKhJhIzu2exqje6VzQM430pNhvfExcTFy9C7eqynHgaDlFNWV//Mu+w2X8fcl2Jk9bzvRJQ4mO1OTtQKICFwkwlVWO5dv383FNaa/bfQiATq3iue7sDEb3as3gzik0i/LP8EZEhJGSUD200rWWC8XOykjmx9OX8+ic9fx6bG+/PKb4hwpcJAAcOFrO/A2FfLyugE/WF7C/pJyoCCMrM4UHxvZmVO90uqQmcJKrdRvV5f3bsWRrEVPmbyGrUwoX923T5BmkdipwEY8cKi3nzaV5zF61m+xt+6mscrRMiGFkz3RG9U7n3O5ptIiL9jomAL8Z35uvdhTzizdyeLdtczq2jPc6kgDWlBePZGVlOS0nK+Gu4FApU7/I5ZUvt3GotIJebZIY3TudUb1ac1bHZCID9FK/HUUljP3TZ2S2SuDNO87x2xCO1M3Mljrnsk48rjNwkSayde8RpszfwoxleZRXVnHpGW247byu9O+Y7HW0eunYMp7Hr+rPpFeW8tA7a/mf75zhdaSwpwIXaWQ5O4p59tPNzF69m+jICL4/qAOTzu1CZmqC19FO28V92zDpvC5Mmb+FwZ1bcnn/dl5HCmsqcJFG4Jzj0w2FPPvpZr7cUkTz2CjuvKArNw7rTFpSM6/j+eSXl/Rk2bb93D9jBX3bNadrWqLXkcKWxsBF/Ki8sop3V+zi2U83s273Ido0j2XiuZ25ZkgGic1C53xp14GjjPvT56QlNuOtu4YTF6Px8MakMXCRRlRSVsHfl+zg+c+2srP4KN3TE3nsqv5c3r8dMVGhN/mlbYs4nrj6LG7422J++69VPHZVf68jhSUVuIgPio6UMXVBLi8vzKW4pJzBmSn81+V9GdUrPeQXjjqvRxqTR3XnTx9tZEjnlvwgq6PXkcKOClykAXYUlfDcZ1v4R/YOSsuruLB3a+64oAuDOrX0OlqTumd0d7Jzi/jtW6vo174Fvds29zpSWNEYuMhpWJ1/gL9+uoV3V+4iwuA7Z7XntvO70C09yetonik8dIxxf/qMxGZR/Ovu4STFBsbko1CiMXARHxw4Ws7/vreW15fsILFZFLeM6MzNwzvTpkVs3R8c4tKSmvHnawdw7XNfcv/Mlfz52gGeTPkPRypwkTrMWb2b3761ir2Hj3HbeV24c2S3gJniHijO7tKKX1zSk9/PXs+Qzi350TmZXkcKCypwkZMoPHSM/3x7Ne+u3EWvNkk8f0MWZ3ZI9jpWwLr9vK5k5+7nf95ZQ/8OyUEzwzSY1Xl9k5m9aGYFNbvvnHjbL8zMmZn2w5SQ4ZxjxtI8LvzDp3y4Zg+/uLgHsyaPUHnXISLCePyq/qQnxXLXtGUcKCn3OlLIq88FqlOBMSceNLOOwEXAdj9nEvFM3v4SbvjbEn7+Rg7d0hN5754R3D2quzYyqKeUhBj+ct0A9hws5edvfIV2Wmxcdf6rdM7NB4pquemPwL2AfkIS9KqqHC8vzOWSP84nO7eI/7q8L2/cdk5YX13SUAMyUvj12N7MXVvAlPlbvI4T0ho0Bm5mlwM7nXM5erZZgt2mgsPcN2MF2dv2c16PNH535Rl0SNF61764cVgmS3KL+P2c9QzslMLgzPC6Pr6pnPbfhWYWDzwA/Ec97z/JzLLNLLuwsPB0H06k0ZRXVvHUvE2MffIzNhYc5vGr+vPSTYNV3n5gZjz8vTPpmBLH3dOWsffwMa8jhaSGDOx1BToDOWaWC3QAlplZrfssOeemOOeynHNZaWm1bLgn4oFVOw9wxV++4NE567moT2vm/ux8vjeog65f9qPmsdE8ff0g9peU85PXv6KySqOt/nbaBe6cW+mcS3fOZTrnMoE8YKBzbrff04n4WWl5JY/MXscVT31B4eFjPDthEE9dPzDol3gNVH3aNee/L+/L55v28uePN3odJ+TUOQZuZtOBC4BUM8sDHnTOvdDYwUT8bfHWIu6bsYIte49wdVZHfj22Ny3iNSGnsV09uCOLtxbx5EcbyerUkhHdddWxv9RZ4M65a+u4PdNvaUQawaHSch6ZvY5Xv9xOx5ZxvDbxbIZ3U4k0FTPjoSvPYOXOA9zz+nLeu+dcWjfXEgT+oItbJaTNW1fAxX+cz7RF25k4ojNzfnKeytsD8TFRPDNhIEfLK5k8bTkVlVVeRwoJKnAJSccqKrl/5kpumrqEpNgoZtwxjN+M70N8jFaP8Eq39CT+97v9WJxbxAufb/U6TkhQgUvIyS8+yg/++iXTF2/n9vO7MmvyCAZkpHgdS4ArzmrPyJ5pPP3JZg4c1VR7X6nAJaQs3LyPy/78OZsLDvPshEHcd2kvmkVpv8ZA8stLenGwtJy/frrZ6yhBTwUuIcE5x/OfbWHCC4tIjo/mrbuGM+aMWqcmiMf6tGvOFf3b8eIXW9lzsNTrOEFNBS5Br6Ssgnte/4qH3l3Lhb3Teeuu4XRLT/Q6lpzCzy7qSWWV408f6dpwX6jAJaht23eE7z69gFkr8vnlJT155vpB2tIrCGS0iue6IRm8vmQHW/ce8TpO0FKBS9Cat76Ay/78ObsOlDL1piHcNbJbyO8EH0ruHtWdZlERPP7Beq+jBC0VuASdqpo/vW+euoT2KfHMunsE5/fQOjvBJi2pGRNHdOadFbtYmXfA6zhBSQUuQeVgaTm3vbqUP3y4gSv6t2PmHcPIaKXVA4PVred1ISU+mt/PWed1lKCkApegsXHPIb7zly/4eF0BD17Whz9efRZxMbpEMJglxUZz18hufLZxLws27fU6TtBRgUtQeH/lLr7z1BccLC1n2sSzuWl4Zy39GiImDO1EuxaxPDJ7nbZgO00qcAlolVWOh99fxx2vLaNHmyTemXwuZ3dp5XUs8aPY6Eh+elEPcvIOMHuVVqU+HSpwCVj7j5Rxw4uLefbTzVx3dgavTxpKmxZaxS4UfXdgB7qnJ/LoB+u10NVpUIFLQFq18wDj//w5i7cW8cj3+vG7K/tpSnwIi4wwfnlJT7YUHuHNpXlexwkaKnAJODOW5vG9ZxZQ5Rxv3H4OVw/O8DqSNIGL+rRmQEYyT8zdSGl5pddxgoIKXAJGWUUVD/5rFT9/I4cBGcnMmjyC/h2TvY4lTcTM+NWYXuw+WMpLC3K9jhMU6ixwM3vRzArMbNVxxx41s3VmtsLM/mlmyY2aUkLe/iNlTHhhES8t3Mat53bm1VvOJjVR+1SGm6FdWnGBlputt/qcgU8Fxpxw7EPgDOfcmcAG4H4/55Iwsn1fCd97ZgFfbS/myWvO4oFxfYiK1B+H4eqXl/TkwFEtN1sfdf6WOOfmA0UnHPvAOVdR8+6XQIdGyCZh4KsdxVz59BfsO1LGqxPP5oqz2nsdSTzWt10LrjirernZAi03e0r+OM25GXj/ZDea2SQzyzaz7MLCQj88nISKD1bv5popC4lvFsnMO4cxpHNLryNJgPjZRT2oqHQ8qeVmT8mnAjezB4AK4LWT3cc5N8U5l+Wcy0pL04JDUm3qF1u57dWl9GzTnJl3DKdrmtbvlv/TqVVC9bX/Wm72lBpc4GZ2AzAeuN5p/qvUU1WV46F31vCfs9ZwYe/WvH7rUNKS9GSlfNvdo7oRE6nlZk+lQQVuZmOAXwGXO+dK/BtJQlVpeSV3TVvG859v5cZhmTw7YZAWo5KTSk+KZeK51cvNrtqp5WZrU5/LCKcDC4GeZpZnZrcAfwGSgA/N7Csze7aRc0qQ23f4GNc99yWzV+/mN+N68+BlfYjU5gtSh6+Xm31ktpabrU1UXXdwzl1by+EXGiGLhKjcvUe48W+L2XWglKevG8il/dp6HUmCRPOa5WYfenctCzbtZVi3VK8jBRRdbCuNaum2Iq58+gsOHC1n2q1DVd5y2iYM7UTbFrE8Mme9lps9gQpcGs37K3dx7XOLaBEXzcw7hzOoU4rXkSQIxUZH8tMLe5Czo5g5q7Xc7PFU4OJ3zjme/2wLd05bRt92zZlxxzA6pyZ4HUuC2HcHtqdbeiKPztFys8dTgYtfVVY5/mvWGh56dy1j+rZh+q1DaaU1TcRHUZER/OLinmwuPMKMZVpu9msqcPGbo2WV3P7qUqYuyGXiiM48dd1AYqN1maD4xyV9W3NWRy03ezwVuPjF3sPHuOa5L5m7dg8PXtaH34zvQ4QuExQ/+nq52V0HSnl5Ya7XcQKCClx8trnwMFc+/QXrdx/k2QmDuGl4Z68jSYg6p2srzu+RxlPztNwsqMDFR4u3FvHdpxdQcqyS6bcO5ZK+bbyOJCHu6+Vmp8zXcrMqcGmwWTn5THh+Ea0SYvjnncMZkKHLBKXxndG+BZf3b8eLn+eG/XKzKnBpkBc/38rk6cvp37EFM+4YRkareK8jSRj52UU9KK+s4k8fh/dysypwOS3OOR6bs57/fmcNF/dpzSu3nE1KQozXsSTMZKYmcO2QDF5fvIPcMF5uVgUu9VZZ5XjgrVX8Zd4mrs7qyNPX6zJB8c7k0d2Ijozg8Q83eB3FMypwqZdjFZVMnr6MaYu2c/v5XXn4e/20b6V4Kj0plltGdGZWTn7YLjer30Cp0+FjFdw8dQnvrdzNA2N7c9+lvTDTNd7ivUnndyE5PpqJL2XzypfbOFYRXhN8VOBySl+v4/3lliIeu6o/t57XxetIIv/WPDaaF28cTPuUOH771ipGPvoJr365jbKK8FgvxZpyecasrCyXnZ3dZI8nvtlZfJQfvrCInfuP8tR1A7mwT2uvI4nUyjnHZxv38se5G1i+vZj2yXHcObIrVw3qSExU8J+nmtlS51zWicfrsyPPi2ZWYGarjjvW0sw+NLONNa91AXCI2VRwiO8/s4DCg8d4+eYhKm8JaGbGeT3SmHnHMF66eQhpSc144J+rGPnYJ0xbtD1kz8jr81/TVGDMCcfuAz5yznUHPqp5X0LE8u37+f6zCymvdPz9tnM4u0srryOJ1IuZcX6PNP555zCm3jSYtKRm/PqfKxn52CdMX7yd8hBbirZeQyhmlgm845w7o+b99cAFzrldZtYW+MQ517Ouz6MhlMA3f0Mht7+6lNTEZrxyyxA6tdI63hK8nHN8sqGQJ+ZuJGdHMR1S4rh7ZDe+N6gD0UF0FdXJhlAaWuDFzrnk427f75yrdRjFzCYBkwAyMjIGbdu2rUFfgDS+d1bk89O/f0XXtERevnkI6c1jvY4k4hfOOT5ZX8gTczeQk3eADilxTB7Vje8ODI4i96zAj6cz8MD1ysJc/uPt1WR1SuH5GwbTIi7a60gifuecY976Ap6Yu5EVeQfo2DKOySO7c+XA9gFd5A1+EvMk9tQMnVDzusCXcOId5xxPzt3Ib/+1mlE903n55rNV3hKyzIxRvVrzr7uG88INWSTHxXDvjBWMfvxT/pG9I+i2a2togb8N3FDz9g3Av/wTR5pSVZXjP99ezR/nbuC7A9vz7A8HERejqfES+syM0b1b8/bdw3n+R1k0j4vi3jdXMPoPn/JGEBV5nUMoZjYduABIBfYADwJvAf8AMoDtwFXOuaK6HkxDKIGjrKKKX7yRw9s5+Uwc0Zlfj+2tHXQkbDnnmLu2gCfmbmB1/kFSE5vRu20SXdMS6ZqeSNe0BLqlJZKW1MyTWcg+jYH7iwo8MJSUVXD7q8uYv6GQe8f05I7zu2pqvAjVRf7hmj28t3IXmwuPsLnwMCVl/zc9P6lZFF1qCr1rWiJd0xLplp5ARsuERp0wpAIXAIpLyrhp6hJydhTzuyv7cc2QDK8jiQQs5xy7D5ayuaC6zP/9UnCE3cdtJhEZYXRqGU+XtES6fV3w6Yl0TU2kRbzvzymdrMCjfP7MEjR2HyjlRy8uIndvCU9fP5AxZ7T1OpJIQDMz2raIo22LOEZ0T/3GbYePVbDluEL/utznbyik7Lgx9NTEZnRNS+DeMb0Y1Mm/k9ZV4CHOOUfuvhKWbC3iyY82cuBoOVNvHsywrql1f7CInFRisyjO7JDMmR2Sv3G8orKKvP1Hv3G2vrnwMM0aYYhFBR5iyiurWJN/kCW5RWTn7id7WxF7D5cB0LZFLNNvHUq/Di08TikSuqIiI8hMTSAzNYHRvRt3DSEVeJA7fKyC5dv3syR3P9m5RSzfXszR8uonXTq2jOO87mlkZbZkcGYKXdMSdaWJSAhRgQeZgoOlLMndX32Gva2INfkHqXIQYdC7bXOuHtyRrMwUsjq1pE0LTYUXCWUq8ADmnGNz4RGyc4tYXDMksr2oBIDY6AgGdEzh7pHdyMpsyYCMZJJiNYNSJJyowANQeWUVj81Zzz+yd7C/pByAVgkxZGWm8KNzOpGV2ZK+7ZoH9NoNItL4VOABprikjLumLeOLTfsY168t5/dIIyszhc6pCZpsIyLfoAIPIJsKDjHxpWzyi0t59PtnclVWR68jiUgAU4EHiHnrCpg8fTmx0ZFMn3Q2gzq19DqSiAQ4FbjHnHNMmb+Fh2evo0/b5jz3oyzaJcd5HUtEgoAK3EOl5ZX8euZKZi7fybh+bXn0qjOJj9GPRETqR23hkYKDpUx6ZSlf7Sjmpxf24Meju+lJShE5LSpwD6zIK2bSy0s5cLScZ64fyKX9tKiUiJw+FXgTezsnn1++kUNqYjNm3DGMPu2aex1JRIKUTwVuZj8FJgIOWAnc5JwrPfVHhaeqKsfjH67nqXmbGZyZwjMTBpGa2MzrWCISxBo8lc/M2gM/BrJqdquPBK7xV7BQcvhYBbe9upSn5m3m6qyOvDZxqMpbRHzm6xBKFBBnZuVAPJDve6TQsqOohFtfzmbDnkM8eFkfbhyWqScrRcQvGlzgzrmdZvYY1ZsaHwU+cM59cOL9zGwSMAkgIyO8tu/6css+7nxtGRWVVbx08xDO7Z7mdSQRCSG+DKGkAFcAnYF2QIKZTTjxfs65Kc65LOdcVlpa+BTYtEXbmfD8IpLjo3nrruEqbxHxO1+GUC4EtjrnCgHMbCYwDHjVH8GCVXllFQ+9s4aXFm7j/B5p/OnaAbSI0zKvIuJ/vhT4dmComcVTPYQyGgjrLeeLS8q487VlLNi8j4kjOnP/2N5EagccEWkkvoyBLzKzN4FlQAWwHJjir2DBZuOeQ0x8OZtdWklQRJqIT1ehOOceBB70U5agte/wMa7660KiIiK0kqCINBnNxPSDJ+Zu5FBpBe/fcy49Wid5HUdEwoT25PLRpoJDTFu8nevPzlB5i0iTUoH76HfvrSM+OpJ7Rnf3OoqIhBkVuA8+37iXj9cVcPeobrTS1HgRaWIq8AaqrHI89O4aOraM44ZhmV7HEZEwpAJvoDeX7mDd7kP8akwvYqMjvY4jImFIBd4Ah49V8NgHGxiYkcw4bcYgIh5RgTfAXz/dTOGhY/xmfB+tLCginlGBn6b84qM899kWLu/fjoEZKV7HEZEwpgI/TY/NWU+Vg3vH9PQ6ioiEORX4aViRV8zM5Tu5ZURnOqTEex1HRMKcCryenHM89O5aWiXEcOcFXb2OIyKiAq+vOav3sHhrET+7uAdJsVrfW0S8pwKvh7KKKv73/bV0T0/kai0TKyIBQgVeDy8vzGXbvhIeGNebqEh9y0QkMKiN6lBcUsafP97EeT3SuKBnutdxRET+zacCN7NkM3vTzNaZ2VozO8dfwQLFkx9t5FBpOQ+M7e11FBGRb/B1Q4cngdnOue+bWQwQUtfWbSk8zCsLt3H14Ax6ttFa3yISWBpc4GbWHDgPuBHAOVcGlPknVmB4+P11NIuK4GcX9fA6iojIt/gyhNIFKAT+ZmbLzex5M0s48U5mNsnMss0su7Cw0IeHa1oLN+/jgzV7uHNkN9KStNa3iAQeXwo8ChgIPOOcGwAcAe478U7OuSnOuSznXFZaWpoPD9d0qmrW+m6fHMctIzp7HUdEpFa+FHgekOecW1Tz/ptUF3rQm7l8J6vzD3LvmJ5a61tEAlaDC9w5txvYYWZfr+o0Gljjl1QeKimr4LE56+nfMZnLzmzndRwRkZPy9SqUycBrNVegbAFu8j2St56bv5XdB0v5y3UDiIjQWt8iErh8KnDn3FdAln+ieG/PwVKe/XQz4/q1JSuzpddxREROSTMxj/P4B+uprHL8akwvr6OIiNRJBV5jdf4B3liax43DM8loFVLzkUQkRKnAqV7r+/+9u5bkuGjuGtnN6zgiIvWiAgc+WlvAgs37+MmFPWgRp7W+RSQ4hH2Bl1dW8bv319IlLYHrzs7wOo6ISL2FfYFPW7SdLYVHeGBsb6K11reIBJGwbqwDR8t5Yu4GhndrxaheWutbRIJLWBf4U/M2UXy0nAfG9sFMk3ZEJLiEbYFv23eEqV/kctWgDvRp19zrOCIipy1sC/yR2euIijR+fnHPuu8sIhKAwrLAV+QV897K3dx+fldaN4/1Oo6ISIOEZYHPWJpHs6gIbhqe6XUUEZEGC7sCr6xyvLtyN6N7p5MUq0k7IhK8wq7AF23Zx97Dx7TWt4gEvbAr8Fkr8kmIiWSkrvsWkSAXVgVeVlHF+6t2c3HfNtoqTUSCns8FbmaRNbvSv+OPQI3pi017KS4pZ/yZbb2OIiLiM3+cgd8DrPXD52l0s3LyaR4bxbnd07yOIiLiM58K3Mw6AOOA5/0Tp/GUllfywZo9XHpGW2KiwmrkSERClK9N9gRwL1Dle5TG9cn6Ag4fq+Cy/rr6RERCQ4ML3MzGAwXOuaV13G+SmWWbWXZhYWFDH85ns3J2kZoYw9Au2qxYREKDL2fgw4HLzSwXeB0YZWavnngn59wU51yWcy4rLc2bsecjxyr4aN0exvZrS5TW/BaRENHgNnPO3e+c6+CcywSuAT52zk3wWzI/mrt2D6XlVRo+EZGQEhano7NydtG2RSyDMlK8jiIi4jd+KXDn3CfOufH++Fz+dqCknE83FDCuX1siIrRpg4iEjpA/A5+zZjfllU7DJyISckK+wGfl5JPRMp4zO7TwOoqIiF+FdIHvPXyMBZv3cVn/ttrzUkRCTkgX+PurdlNZpeETEQlNIV3gs3Ly6Z6eSM/WSV5HERHxu5At8N0HSlmSW8Rl/dtp+EREQlLIFvi7K3fhHFo6VkRCVsgW+KycfPq2a06XtESvo4iINIqQLPAdRSV8taNYT16KSEgLyQKftSIfgHH9NHwiIqErNAs8ZxcDM5Lp2DLe6ygiIo0m5Ap8U8Eh1u46qOETEQl5IVfgs3J2YabhExEJfSFV4M45Zq3IZ2jnVqQ3j/U6johIowqpAl+76xBbCo9o+EREwkJIFfisFflERRhjzmjjdRQRkUYXMgXunGNWTj7Du6XSMiHG6zgiIo3Ol13pO5rZPDNba2arzewefwY7XV/tKCZv/1ENn4hI2Ijy4WMrgJ8755aZWRKw1Mw+dM6t8VO20zIrZxcxkRFc3Le1Fw8vItLkfNmVfpdzblnN24eAtUB7fwU7HZVVjndW5HNBzzSax0Z7EUFEpMn5ZQzczDKBAcCiWm6bZGbZZpZdWFjoj4f7liW5RRQcOqbhExEJKz4XuJklAjOAnzjnDp54u3NuinMuyzmXlZaW5uvD1WpWTj5x0ZGM7p3eKJ9fRCQQ+VTgZhZNdXm/5pyb6Z9Ip6eisor3V+3mwj6tiY/xZUhfRCS4+HIVigEvAGudc3/wX6TTs2DzPoqOlHGZNm4QkTDjyxn4cOCHwCgz+6rmZayfctXbrJx8kppFcX7PxhmeEREJVA0ec3DOfQ54utnksYpKZq/ezcV929AsKtLLKCIiTS6oZ2LO37CXQ6UVXNZfwyciEn6CusBn5eSTEh/N8G6pXkcREWlyQVvgJWUVfLhmD5f2a0t0ZNB+GSIiDRa0zffxugKOlldy2ZmavCMi4SloC3xWTj7pSc0Y0rml11FERDwRlAV+qLSceesLGXdmWyIjPL0QRkTEM0FZ4B+u2UNZRRXjNXwiImEsKAt8Vk4+7ZPjGJiR7HUUERHPBF2B7z9Sxmcb9zK+f1uqZ/OLiISnoCvw2at3U1HldPWJiIS9oCvwWTn5dElNoG+75l5HERHxVFAVeMHBUhZu2cf4/u00fCIiYS+oCvy9lbtwDi0dKyJCkBX4rBW76NUmie6tk7yOIiLiuaAp8J3FR1m6bb/2vRQRqRE0Bf7uinwAxmv4REQE8H1PzDFmtt7MNpnZff4KVZtZObvo36EFnVolNObDiIgEDV/2xIwEngIuBfoA15pZH38FO97WvUdYufOAhk9ERI7jyxn4EGCTc26Lc64MeB24wj+xvumdnOrhk3EaPhER+TdfCrw9sOO49/Nqjn2DmU0ys2wzyy4sLGzQA7VuHssPsjrQtkVcw5KKiIQgXwq8tpk07lsHnJvinMtyzmWlpTVs5/gfDO7I77/fv0EfKyISqnwp8Dyg43HvdwDyfYsjIiL15UuBLwG6m1lnM4sBrgHe9k8sERGpS1RDP9A5V2FmdwNzgEjgRefcar8lExGRU2pwgQM4594D3vNTFhEROQ1BMxNTRES+SQUuIhKkVOAiIkFKBS4iEqTMuW/NvWm8BzMrBLY18MNTgb1+jNMYAj1joOeDwM8Y6PlAGf0h0PJ1cs59ayZkkxa4L8ws2zmX5XWOUwn0jIGeDwI/Y6DnA2X0h0DP9zUNoYiIBCkVuIhIkAqmAp/idYB6CPSMgZ4PAj9joOcDZfSHQM8HBNEYuIiIfFMwnYGLiMhxVOAiIkEqKAq8KTdPPl1m1tHM5pnZWjNbbWb3eJ3pZMws0syWm9k7Xmc5kZklm9mbZrau5nt5jteZTmRmP635Ga8ys+lmFhsAmV40swIzW3XcsZZm9qGZbax5nRJg+R6t+TmvMLN/mlmyV/lq8nwr43G3/cLMnJmlepGtLgFf4E25eXIDVQA/d871BoYCdwVYvuPdA6z1OsRJPAnMds71AvoTYDnNrD3wYyDLOXcG1UsoX+NtKgCmAmNOOHYf8JFzrjvwUc37XpnKt/N9CJzhnDsT2ADc39ShTjCVb2fEzDoCFwHbmzpQfQV8gdOEmyc3hHNul3NuWc3bh6gunm/tDeo1M+sAjAOe9zrLicysOXAe8AKAc67MOVfsaajaRQFxZhYFxBMAO1A55+YDRSccvgJ4qebtl4DvNGWm49WWzzn3gXOuoubdL6nezcszJ/keAvwRuJdatooMFMFQ4PXaPDkQmFkmMABY5HGU2jxB9T/GKo9z1KYLUAj8rWaI53kzS/A61PGcczuBx6g+G9sFHHDOfeBtqpNq7ZzbBdUnGEC6x3lO5Wbgfa9DnMjMLgd2OudyvM5yKsFQ4PXaPNlrZpYIzAB+4pw76HWe45nZeKDAObfU6ywnEQUMBJ5xzg0AjuDtn/3fUjOOfAXQGWgHJJjZBG9TBTcze4DqIcjXvM5yPDOLBx4A/sPrLHUJhgIP+M2TzSya6vJ+zTk30+s8tRgOXG5muVQPQY0ys1e9jfQNeUCec+7rv1zepLrQA8mFwFbnXKFzrhyYCQzzONPJ7DGztgA1rws8zvMtZnYDMB643gXeZJSuVP9HnVPzO9MBWGZmbTxNVYtgKPCA3jzZzIzqsdu1zrk/eJ2nNs65+51zHZxzmVR//z52zgXM2aNzbjeww8x61hwaDazxMFJttgNDzSy+5mc+mgB7ovU4bwM31Lx9A/AvD7N8i5mNAX4FXO6cK/E6z4mccyudc+nOucya35k8YGDNv9OAEvAFXvNkx9ebJ68F/hFgmycPB35I9VntVzUvY70OFYQmA6+Z2QrgLOB33sb5ppq/Dt4ElgErqf7d8Xy6tZlNBxYCPc0sz8xuAR4GLjKzjVRfRfFwgOX7C5AEfFjz+/KsV/lOkTEoaCq9iEiQCvgzcBERqZ0KXEQkSKnARUSClApcRCRIqcBFRIKUClxEJEipwEVEgtT/B2qqdo68ugTAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(DeltaF_ij[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "585f5ad2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.236805457595942"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DeltaF_ij[0][-1] * 0.239006"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "3ed68e58",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(16, 16, 100)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_kln.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "3a03b024",
   "metadata": {},
   "outputs": [],
   "source": [
    "def estfrnrg(dat, nit):\n",
    "    u_kln= dat[:,:,:nit].copy()\n",
    "    \n",
    "    # Subsample data to extract uncorrelated equilibrium timeseries\n",
    "    N_k = np.zeros([nstates], np.int32) # number of uncorrelated samples\n",
    "    for k in range(nstates):\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\n",
    "    # Compute free energy differences and statistical uncertainties\n",
    "    mbar = MBAR(u_kln, N_k)\n",
    "    [DeltaF_ij, dDeltaF_ij,] = mbar.getFreeEnergyDifferences()\n",
    "    return DeltaF_ij[0][-1]* 0.239006"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "a3f643c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "nrgz = list()\n",
    "x = np.arange(2,100)\n",
    "for i in x:\n",
    "    \n",
    "    nrgz.append(estfrnrg(u_kln, i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "c36b6c64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1adcedcd0>]"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABHBklEQVR4nO29d5RkZ33n/X0qx+6u7q7O3dOTgyYoDCMUEEogIViSMQZszOsXXh0MxmIPNpiXs2vv2uu1F1uwXpK1wEvSK9ZIsggrwAIJBRQnaHLumekcq7srx/vsH/c+t25V3UrdVVPp9zlnznRX3a56bnfVt373+/wC45yDIAiCaHwMtV4AQRAEURlI0AmCIJoEEnSCIIgmgQSdIAiiSSBBJwiCaBJMtXri7u5uPjo6WqunJwiCaEgOHTq0yDn36t1XM0EfHR3FwYMHa/X0BEEQDQlj7Eq++8hyIQiCaBJI0AmCIJoEEnSCIIgmgQSdIAiiSSBBJwiCaBJI0AmCIJoEEnSCIIgmgQSdqCueODKF1Uii1ssgiIaEBJ2oG6ZXIvj0/3odjx6arPVSCKIhIUEn6oaFQAwAML4UqvFKCKIxIUEn6gZfKA4AmFiO1HglBNGYkKATdcNiUI7QJ3zhGq+EIBqTooLOGBtmjD3DGDvNGDvJGHugwLFvYIylGGPvq+wyiVZAROiTyxHQrFuCKJ9SIvQkgM9wzncCeCOATzLGdmUfxBgzAvh7AL+s7BKJVmFJEfRIIoXFYLzGqyGIxqOooHPOZzjnh5WvAwBOAxjUOfRTAB4DMF/RFRItw5JGxCeWyXYhiHIpy0NnjI0CuA7AK1m3DwJ4D4BvFPn5+xljBxljBxcWFspcKtHs+EIxOCxGAOSjE/XLkfFl/OrUXK2XoUvJgs4Yc0GOwD/NOfdn3f1lAJ/jnKcKPQbn/CHO+X7O+X6vV3fgBtHCLIXi2D3YDoAEnahf/u7nZ/Affnyi1svQpaSJRYwxM2Qxf5hz/rjOIfsB/JAxBgDdAO5jjCU5509UaqFE87MUjOPGTZ0YW7Bgwkepi0T9kUhJODq5glhSQiyZgtVkrPWSMigq6ExW6W8BOM05f1DvGM75Rs3x3wHwMxJzolyWQjF0OS0Y8jjIQyfqktMzfkQTEgBgZiWK0W5njVeUSSkR+i0APgzgOGPsdeW2/xfACABwzgv65gRRCuF4EtGEhE6nFcOdDrw+sVzrJRFEDoeupF+Xk8uRxhN0zvkLAFipD8g5/7/WsyCiNREZLl0uC0Y67Xjy+AySKQkmI9W+EfXD4fEV2M1GRBKpuryKpHcLUReIHPQupwXDHgdSEsfMarTGqyKITA5fWcbt270wGRgmSdAJQh9fSC7773LJlgtAmS5EfTGzGsHUSgQHNnaiv8NWlxv3JOhEXSAqQ0WEDlBxEVFfHL6yAgC4YYMHwx4HRegEkQ/Rx6XTaUF/hw0GhrqMgIjW5dCVZdjMBuzsb8OQx16XXUFJ0Im6wBeKw2Y2wGExwmw0oL/dThE6UVccGl/G3qEOmI0GDHkcWAjEEE0UrKW86pCgE3XBYjCGLqcVSnEahjvt5KETdUM0kcLJqVXcsMEDQH59AsDUSn1F6SToRF3gC8XR5bKo3490OurykpZoHf7ltQk8c0buNXhschVJieOGEVnQhzz1uXFfUuk/QVQbXyiOTmda0IeVS9pIPAW7pb7Kq4nm58ysH5997BgA4O6dPRjokCPy60WErgj6ZJ0FHRShE3XBUjBL0DvFG6a+IiCiNfj6by7CaTHiM2/ZhhcvLuF7L13Bpm6n+hrtcVthNrK62+chQSdqDuccS6EYul1W9TbhUdbbG6bRWArG8I1nL9IEqDIYXwrjp0en8aEbR/Cpu7bi1595M35v/zA++ia1ZRUMBobBDnvdRehkuRA1JxxPKX1cMi0XAHj45XE8/PI4Do0v44MHRvC5e3fUapkNycOvjOPBp87hLbt6sdnrqvVyGoKHnr8Ik8GAj71pEwCgv92Ov3/f3pzjhjsddSfoFKETNcenKfsXeN1WeBxm/PrMPC4thsAAvD6+UpsFNjCvXFoCkDkNisjPfCCKfzk4id+5YRC9bbaCxw557JikTVGCyETt46LJcmGM4ecP3AajgcHrtuL/+d7BussoqHfiSUntDrgUjNV4NY3Bt1+4jGRKwv23bS567JDHgaVQHOF4Eg5LfUgpRehEzRFi0+m0Ztze126D1y3f1uW0qMJPlMbxqRW1d3ex391/ffI0Xrq4dDWWVbdEEyk8/PIVvG1PPzaW0BZ3yCPv89ST7UKCTtScJR3LJZtOpwXLoTht7pXBy2M+9etClksyJeGfnxvDE0emcu777ouX8cVfnqnK+uqNyeUwArEk3rKzt6Tjhzz1l4lFgk7UHG0v9Hx0Oi1IShz+aPJqLavheeWSD9t6XWizmdRulnqsRhIA9DOKnnh9Co+8OlG1NdYTUytyu+ZBJfIuhpqJVUc9h0jQiZrjC8VgNxsL+pAiA8ZHtktJJFMSDl324caNXehyWbFY4Pe2HJbv07MOJnxh+ELxlvDgp5UyflFEVAyvywqryUAROkFoWcqqEtXDQ4JeFiem/QjFU7hxUye6nBb4Clguy2E5Qp9eiSAlpS2tUCyptjW+MB+s7oLXwbHJFcxWYBjK9EoERgNDr9ta/GDIG/dDnsK56BcXgvj848ex+y9/qbYRqCYk6ETNWQrGC9otQNpfJ0EvjVfG5A3OAxs70eWyYKmA5SJ+p0mJY2Y1LU5aC+Z8HQv6R797EF/+1bl1P87UcgR9bbayxh7mG2geS6bwiYcP4e4Hn8VjhycRT0n45cnZda+xGCToRM3xheIFN0QBreXS/Jf+leCVSz5s8jrR47ah02kt+EG4Ek7fp402x5fSQlXpCD0STyEcL7wfMrkcxicePlTQ7oklU1gIxCoyrnBqJYKBjsK559kMeey6HvpLF5fw5PFZfOSmUbz4F3fitq3dePWST+cRKgsJOlFzloKxnJTFbNKCnrgaS2poUhLHa5dk/xwAul0W+EJxSJJ+hpD2d6rN9R9Xvh7pdFRc0P/80aP45MOHCx7zr4en8OTxWXzj2Yt5j1kIxDL+Xw/Tq5GS/XPBQIcdq5FETl90sZ4/umUU3S4rDmzsxNhiCPOB6s7JJUEnaorcxyWO7iKWi8Nigs1syInQJ3xh/Pr0XDWX2HCcmvYjEEvijZs6AcgfhhIHViL6H4Yr4TgsRoM8JUoToU/4wnDbTNg/6sH5+UBF13hhPogLC4U/JJ49twAA+N5LVzDv1xfCOb/8elhc56ZtSuKYXY2WLejidZv9gSL2HkR/ogPKh+trl5bXtc5ikKATNSUcTyGWlIpuigJAl9OaE6F/8/kx3P/9Q0Uv31sJUe4vIvQuRVTyWReidXFfmy2jlH3cF8ZIpwNbe9yY88fgj5Z/dbQUjCGRknJunw/EMO+P5a0r8EcTODKxgnfs7UdS4vh6nih9QYl4l0LxjA3dclkMxpBIcQyWKeii8C37A2UxKGduOa1y5tY1A21wWIx49VJ1i7dI0ImaInLQSxF0j9OcE6FPKZkZxydXq7K+RuT8XBDdLgv62mU/WOxP5KsWXQ4n0OEwY8iT2WwqLehyU69ybZdYMoU7/uE3+O6LlzNujycl+EJxxJISAjH9D+IXLywiJXH84U2jeO91g3j4lXHdTBYRoackrqZfrgVx3uUKuojAcyP0GLrd6de02WjADRs8eKXKPjoJOlFTRPZFsSwXALqbe9NKMciRiZWKr61RueILYUTpJw+kf7f5qkWXw3KEPtSZnuMqSRwTyxEMdzqwtVcR9LnyBP3ifAj+aDLng0Abzc779a8anj23ALfVhOtGOvCnd22FJHF8/TcXco7TetLrsV3KzUEXpCP0zN/tUjCe0Q4aAA6MduLsXCBjE7rSkKATFeVXp+bw2wuLJR+vVokW2RQFgE6HGb6sN4NIs6NOjGnGl8LY0JXuRVIsQ2g5HIfHYcGwx4FZfxSxZArzgRjiSQnDnQ4MeRywmAxFPe9szs75ASAnA2VeE83qbWZyzvHcuUXcvKULZqMBw50O/O7+YTzy6oQqvII5f+HHKpW0oJeX5SJet3oRevZr+sDGTnAOHLxcPR+dBJ2oKF/61Tl86anSc4J9Op0W89HptGYUyETiKbUo5shEdTebGoVYMoUZfzQjQu90yL/b7ChSsByKw+M0Y7jTAc6BmZVoRoaL0cCw2evC+bnyNkbPzMrHz2VtaGo3OBd0ouqLC0FMrURw2zavetv/fcso4ikJL2QFC/OBGNw2k3J+6xP0NpsJbpu5rJ+zmAzocJh1PXSvO/M1vW+4AxaTAa9erp7tQoJOVJRwPFVW97l0Y67iEXqXy4JQPKWmiE0r0fm+4Q7M+WMZRTGtyoQvAs6BDV1pQTcZZdHRKy6SJI7VSAIeh0XtHjixHM4QdADY0uMqu7jozIws6IUidL3slWfPyaJ929a0oI92O2E0sJwWyvP+KK4ZaANQeoTujybwnJJBI5Bz0MuzWwTdLmvGc6ckDl8o13KxmY24drhDLfqqBiToREUJxZKYC8iX7aWwpGQDlDII2qNEmmLza0bxz9++pw8AcIRsF4z7QgAyBR2QN0b1iov80QQkLv9uxRzXCV8E474wGEtvEm7tcWFqJVJWNtFZJUJfjSQQiadfD/OBGBgDLEaDboT+3LkFbPI61fUA8qbiQIdN/aDRPtYmrwtWk6EkQb+yFMJ7vvpb/OG3X8Xh8fRV3dRKtOwNUYHXZc2I0H2hOCSOHEEHgBs3duLEtB/BPJvB64UEnagokXgKnKc3K4vhCxUv+xcIL1j47iJCv2tnLywmA16njVFcWRKRdWY/7y6nVddyESLvcZrlsncDw+RyGBO+MAba7bCYZInY2uMC58DYQqikdayGE5j1R7G91w0AmNXaLIEoupwWeN1WLGRtikYTKbw8tpQRnQtGOh0Zgi6yZXrdct/8fJaS4OWxJbzrq79Vrwp/ez5t30yvRErusphNt9ua8cEkxF1P0A9s7ERK4jh8pToWIQk6UTE45wgpEVyp04UWSyj7FwhBFxH69EoEjMnzR3cPtOHIOPnoV5bCcFiMOYVaXS79CF3sQXgcFhgNDAMddkwsyxG6aA8LyJYLgJILjM7Myhuib94uC7M25XDeH4NXEeHsCP3VSz7EkhLevE1f0LWvK/GzvW1W+bEKROjPnJ3Hh7/1CrqcFjzxiVuwq78NLyoDPYKxJFYjiTVbLl6XFYsBPUHPfV1fP+KB0cCq1gaABJ3I4MULi/jUI0dwYqr8vO5YUoKo7cj20X90cAL3fvm5nEISXyimFr4UI7uF7sxKFN0uKywmA64d9uD41KpuEUsrIXLHGWMZt3c6LbqFRcsiQlfsrOFOOyZ8YfVxBBu6nDAZWMm56GeVDdTbFWGe9adfD/OBGHrcVvS4rTlpi8IGuVGpctUy3OnAYjCOkGJXiM3W3jZbjo+tJZmS8Nc/PYUNXU48/olbMNrtxC1bunBofBnRRGrNKYuCbre8tyPsKFXQdbo2Oq0m/M27d+Pe3X1req5ikKATKvOBKP7kkSP46dFpvON/vIBP//BIWb2ewxqfNLsD3fPnF3FmNqAOUxAsBYu3zhVkd1ycXo1gQCmeuW6kA9GEpPq2tSaaSCFZgw+XK0uhHP8ckKtFVyKJnDWJqx1V0D0OXFwIYiEQyxB0i8mA0W4nzpeYi35mNoB2uxn7hjsAALOrmo3QQBQ9bqtuhH5lKYyBdptub3yxHvHaEhuqXuWx8mW5PH54CmOLIXz2nu1ot8tZLDdv7kY8KeHg5WVMrYiiovJSFgVeJSBZDMQz/tezXADggwdGsHuwfU3PVQwSdAKAbJd87tFjCMaSeOyPb8Ynbt+Mn5+Yxb1ffr7k4QbaDbPsCF1kSGjzhkUfl1I99Ha7GQamEXRNZsK1inDUS4HRO7/yAv72yas7uk0UA2lz0AVdTgs4T1ssAlXQnbLQDXc6EFCmQmk3JQFgi9eFk9P+nA9lPc7OBrC9zw2n1QS3zYRZZb8jJXEsBuPoUWwSXyiecVV1eSmku34gLeiiC6TIlhERui8cz7lCiyVT+O+/Po99wx14y670aLkDGzthMjD89uKiGqEPduR+EJaCiMQXgvIHzGIwBovRgDbb1R8cXVTQGWPDjLFnGGOnGWMnGWMP6Bzz+4yxY8q/Fxlj+6qzXKJaPPzKOJ45u4DPv20HbtjgwWfv3YGvfuh6BGNJXCxxIywjQtd4nSmJ46JSlKLdHAvFU4gnpZI9dIOBweOQh0VzzjGzGkV/uyzoQx47ul3WuvDRfaE4zs0F8eTxmas6A3XWH0U8KWVE1gLxoZntoy+HEzAZGFxKz5EhzcZg9uPcvasXUysR3Pp3T+PBfzubt+KRc46zswHs6JM3RPvbberf3af0XOlx29DjliNibWR9ZSmM0TwDmlVBV15bc/4ojAambrBynnt+j7wyjqmVCP7srdsybCin1YRrhzvw4gVZ0E0GplZ9lotXLf+Xn3shGEO3y5Jje10NSonQkwA+wznfCeCNAD7JGNuVdcwlAG/mnO8F8NcAHqrsMolqcnEhiL/536dw2zYvPnLTqHr7kLIpVmp+r/A2e9usOT1B4kk5ctIWmYjIv1jrXC1iWLQ/kkQ4nlIr+xhjuG6ko+xMl5nVCFbDlW3Je3Ja3n+Y9UfV4pqrgchw0bNc0hlCmX9LuagoLT5i8DGQK+jvu2EI//tPb8WtW7vxT09fwFu/9JxuGuPkcgTBWBLbFUHvbbOpm6KiVF9YLkD69bUaScAXimNUZ/2AfIXmtpnUYGHeH4PXZYXBwDSimj6/cDyJrzxzETdu7MStW7pzHu/mzV04PrWKMzMB9LXbYDSsTYDV8wiKzo9xXf/8alBU0DnnM5zzw8rXAQCnAQxmHfMi51yERi8DGKr0Qonq8dWnL8BsNOCL79sLg+ZFrXqDJVouItd4W68bi8GY+r22wlBbSLJURpWowOOUI3SRsigidEC2XcYWQiV3BUykJLz7q7/FX/30ZMnPXwqnpv3q1785u1DgyMqi5qB35ka4ws/NbtAll/2nqyNFZovTYtTd27hmoB1f/4Mb8N/etxfzgRguLeZevYl9DL0IXdgkPW3ypiiQ7ucyrn4g6UfojLGM1MW5QAw9bfJjiKpMrSf/g5evYDEYw5/fs103Wr55SzckLveNWeuGKCB/WDIGNdNlMRDL659Xm7I8dMbYKIDrALxS4LCPAvh5np+/nzF2kDF2cGHh6r3QicK8dsWHW7d0o7ctc1NIpLKVHKErAi5yj6dW5Dee8M/tZmOGh57u41K6oHcpEbrwPfs1G1kioix1vuQzZ+Yx54/haIV995PTfgx22LGzvw2/OVv9OZKCK0thmAxMtx9J/gg9oW6IAunBx8M6mTJadvbJ1Zl6VcEiw2Wb8jroa7NhIRBDMiWpeec9StoikBbhS0vyh8Nod34vWyvo8/6oatt4XYp9o3mtPn1mHrsH27B/NDdjBpA30m1mA5ISx9A6BN1sNMDjsKjnsRSKFe3vXy1KFnTGmAvAYwA+zTn35znmDsiC/jm9+znnD3HO93PO93u9uXmmRGUoZ+7mfCCKCV8E1494cu4zKP5ktqBLEsff/OwULmTlJIvLb3GpLUZzXZgPor/dhg1djgwP3ad2Wiw9mvEoFY/Timhrq/uyL+GL8S8HJwDIQlLJfuonp1exs78Nd2z34uCV5bxXDOF4ck3pofm44gtj0GPXnYnpcchRZK6HHs8QdMYYtvW61b9hPtQ2ATr1BmdmAxjy2NW+KH3tdkhcFm5huXjd1pzWs1cW819hCEY6HZhYjkCSOOYDMfQqEXp3VoQuSRwnpvzqZrkeVpMRb1DEfj0ROpDORZckrttp8WpRkqAzxsyQxfxhzvnjeY7ZC+CbAN7FOa9uF/c6ZyEQw8EqNuApxIQvjDf8l1+VPGH88JUVAMD1G3IFHYBuatmMP4pvvnAJT53KfA6xKSrEQKQ8np8PYEuPC71ttgzLZXGtEXo4jqlleSNL+8bpLsMimvdH8czZBWzrlSsgK5XuGI4nMbYYwjUDbbh9ew9SEs+oSNTy0HNjeO/XXswZX7ZWxpfCuhuiAGBUNpQX9QQ96/f/rY/sx39+5+6Cz9XhMMNlNelH6LN+1W4BgL52+e8yuxpVm2nZzEZYTAZ4HGZV5C8vhdHXZivYBmK404F4UsLUSgS+UFyN0B0WE5wWo/rhcGkphGAsib1DHQXP4+bNsre+XkHvdssR+mokgaTEywpSKkkpWS4MwLcAnOacP5jnmBEAjwP4MOd8/eO3G5z/8fR5/N5DL+PKUmnZIZXk5LQfKYnjpRIbAB0ZX4bFaMDuwTbd+/Uq8GYUuyMYy4w8haCPdDpgNRnUSOrCfBBbe9zobbNmWC6+UBwOixE2c/E+LgIxTu3MrB+9bZkbWeVE6I8enkRK4vjC2+X9/dMzlRH0M7MBcC5PqLl+pANumwnP5LFdXp9YQTwllZQGWAr5ctAFXU5LRrdKzjlWwokMDx0AetpsaHcU7jrIGMOQx54j6LFkCmMLoYwIv69NFsvZ1Sjm/THVOwcyX1/F1g+kbbVDSum8iNDFY4kg4djkCgBg71DhfO+7dvbAbGRqg6+1Ivq5FKoSvRqUEqHfAuDDAO5kjL2u/LuPMfZxxtjHlWP+I4AuAF9T7j9YrQU3AkcnV5GSeMHhttVibFH2q8ULuhiHx5dxzWAbrCZ9UfXqVOCJ7nnBaKZNEVayXBwWk/JmD2NqJYJoQsLWXjlCXwjG1FFh5fRxEQgv+MSUP8crbrOZYDHpN3zSwjnHjw5O4sBoJ960pRsuqwmnZ3RdxLI5qWyI7hpog8lowG1bvXj23EJO+iLnXLVbViqQZbMSjsMfTRa0KzqdloyOi4FYEkmJl1zYlY34G2u5tBhCUuKqfw5AnZw0649iIRhTo2pA9tLFRunlpTBG82yICoSgv6ZcAfe0ZV6hiZF0xyZXYTcbscXrKvh423rdOP5X96gFUGtFVKqK1563XiN0zvkLnHPGOd/LOb9W+fck5/wbnPNvKMd8jHPu0dy/v/pLr08SKQmnZ/ywmAx49NBkTkP+aiOaJ52Y8ued8i6IJyUcm1zV9c8FXrcVS6FYxmOJTcdAlqCH4ilYjAZYTAYMeRyY8EXU3h9be1zoabMhJXF1Y24xGCsrZRFIC/piMJaR4QLIUaPeB1A2r17y4dJiCO9/wzAMBoYdfe6KCfqpaT/a7WbV23/zdi/m/LGcK4A5f0yNJtczOk2gNuUqEOF2u6wZU4tE2X+HY62CLo+s035YiUpSraB7HGZYTAbFcolmiLCI0IOxJBaDMWwosCEKyNaIgaUjdO2HgzbaPza5imuUD9VilHOFmA+v24poQlL/DnWbtkiUx4X5IOJJCZ++eys4B/75KkfpIo0sGEuq0Xo+Ts34EUtKRQU9keIZtoBIGcyeBxmJJ1X/c1gZZybe4Ft6XOhTsmiE7eILxdFdZnSo3cDT8z27izRpAoB/OTgJl9WE+5S2u7sG2nBmNlD0A7AUTk2vYld/m5ohInqZZNsu2s3QSkToV3z5c9AFXS5LRtpiujFXeUMdBEMeu9rYSnBhPggDAzZqioMYY+hrs2FGx3LpcVsxH4jhsvK6LRahW0wG9Lfb1Uya7A+HxWAcyZSEk9OrRf3zSiL2b84ogUFdb4oSpXNceaPec00ffuf6ITzy2kTG3MNqM7YQxH5lg/NY1uDk1UhCLfABoLbwvH5DR97HUzMRNDbGbB7LJRRPwakI+pDHgZVwAkfGV+B1W9HhsKh+p8h0KaePi0Br0eil58leZv6IdzWSwJPHZ/Dv9vWr/UJ29rchGEtm9J958eIivv6b8j6MkykJZ2YDGX5sT5sN1wzkpi8e1wj6aiRzvfP+KG79+6fLmhA0ruzX5NsUBeSrm9VIQi2PVxtzrdlykZ9L66NfmA9ipNORE/X2tdtwfj6IWFLKiarjSUn9fRQTdEA+R86hVIlmWi6rkQROzfgRTUhF/fNKIvZvTs8GYDQwdNjX9iG5XkjQK8zJqVU4LUZs7HLij2/fjGRKwjefv3RVnns5FMdyOIG7d/XCYTFmCHpK4nj7Pz2PTzx8SL1EPjy+jIF2W451oUVvo1H10HMi9FQ6Qlfe7C9cWFSnxveqEXoUnMtTXTrL9NC1Ebreur3u3DRLLY8emkQkkcLv37hBvW1nvyzAWtvlvz55Bn//izNlTbofWwwhlpRwTdYG8x3be3B4fCWjIvXk9KraWCy7v8rZuQAmlyM5H8iFuLIURrfLqtvUSiAyL4SQZzfmKheRuqj10S/MB9VWu1r62mw4lyeqBtKeeLFNUSD9oeV1WXU3xZ9WMrz2XEVB10boXU5LRoHe1YQEvcKcmPbjmoF2GAwMo91OvHPfAL7/0pWK5Dk/emgSn3/8eN77x5TL1q09LuweaM/YGH3tsg+TyxH86vQ8fnZsBoA84ee6POmKAn1BF1ku2RF6Es6sniDBWFIV9C6nBQYmR6CBWBLxlITuMj10m9moXgX0t+tH6L5QeuNViyRxfP+ly7h+pCOj2932XjcMDDil+Nzn5gJqxPjwK1dKXpso+d/Vnykkd+yQ0xefO58upjs+tYoDGzthMRpyLBeRK15sc1fLufmg+nvOR5e6/xDPeJ7ONQr6cFaEnkxJGFsMYktPbg57f7tN/Zt43fqC7nVb1ddPIcQ+gfaDAUhvRD59Zh5uqwkbS4j2K4U4D380WTO7BSBBrygpiePUtD8jQrt1qxeRRCqn73O5XFwI4gv/ehyPvDqed6N1TGmAtcnrwp6hdpyc9quX1z87Ng2b2YBrBtrwn356Emdm/Zha0S8o0pIt6ImUpGYlBLIKZsKxFByqh56OtLYoG2QmowFet5y6KNLn1pJhIaL6fB66pNOkCQCev7CIy0thfOTm0Yzb7RYjRrudaoT+2OFJmAwMt27pxmOHJjPGpxXi5JS8Gb7Zmykk1w53oMNhVn30+UAUc/4Ydg+2o91hzrFcxMZlqQVSksRxbjaAHf2Fi4HEFZLYqF4JJ2BgUIcsl0ub3QS3Nd1bZdwXRiLFdSN0bRVyZpaL/Pqa8EXy9nDJRry2tI8DpDcij02uYvdg+1WNkjuVYEW7jlpAgl5BxhaCiCRS2D2QjtDEm2U9MwRTEsef/+goRCX28+f12yaMLYZgMsj5wXuH2hFLSjg3F0AyJeEXJ2Zx145efPF9+7AcTuD+7x0CANxQJEJ3W02wmgxqfu1CIAbO5RTB7CyXcCKpXvJ7HGZV3LWRY2+b3NdjSa0SXYOgO62wmQ26m3mF+s9878XL6HZZ8bbd/Tn37exvw+kZOYf/iSNTuH27F5+6cwv80SR+emy6pHWdmpELarIzK4wGhjdv8+LZswuQJK6mNu4ebIfHYcZyKE+EXqKgj/vCiCRSajl+Pq4d7sBgh12tkPUpVaJrFT7GGAY1ueiixYPelYL2airTcknfnq+HSzbCcsmJ0DVCunf46tktgPw3FhlbtcpBB0jQK4q4TNd6d0LQS20Ypce3XhjD4fEV/N1796KvzYbnzulXHo4tBDHS5YDZaFB3+I9PruKVSz4sBuN4x95+7Bpow/23bcK4LwyryYBd/YVFgDGWkQ4m7JatvW7EklLGJqs2QmeMqZfk2jd4j9uGOX9U08el/GjG67JgsMOu22skX3HRhC+Mp8/O44MHhtU5mVp29bdhcjmCn5+YwZw/ht+5fggHNnZiS48LD78yXnRNnHOcmvHnLVC5c0cPlkJxHJtaxQnFG9810IYOuwUr2RG6IuilNkUT496KlesbDQwfPDCM315YwqXFEFbCcXSsMcNFMNzpUAVd7Dds1ovQFUG3mQ1wa2wVUTsAZGbGFGKDIuh9Wb2HtBXHewc7SjyDyiGEvFY56AAJekU5MeWHzWzAJs0L022V3zDZGSGlcmE+gH/4t3O455pevOvaAbxpazdeuLCoOw3n0mIIm7rlN9NolwNtNhOOTq7iZ8em4bAYcceOHgDAA3dtxcZuJ27Y4NEVt2y05f9iQ3Rbr/w8Ic2VRzieFnRATl3sdFoyyqD72uU0NdW/XUM087l7d+CLv6vfcj+7P4jgBy9fgYExfOjGEd2f26nYFf/tF2fRbjfjzp09YIzh928cwdGJlaI9VxYCMayEE2pjsmxu2+qFgckNwU5Mr2JjtxNtNjPaHWYdDz19NVQKp2cCYCwz9zsf798/DJOB4ZFXx+VN6TVmuAhEcRHnckXwQLtN7a2uRUToPW5bxgexqB0AStsQBeSsnK/9/vX44IHMv6XNbFSHSlzNDBeBCCbIQ28STihNmbSX3CJCz7YnSuVLT52H3WzE37x7DxhjuG2bF6uRBI5mZUCkJI7LS2FsUvxbxhj2DnXgyPgyfnFiFnfv7FVTyWxmIx7/45vx1Q9dX9IatPMaZ1ZkQd+qbHxpzysUT2ZkWfzJnVvxt+/Zk/FYvW5bRnOtcvq4CLb2uvN6/+JNpY1uo4kU/tfBCdxzTW/ejB6R6TLuC+Od+wbUytn3Xj8Em9lQdHNUDAHZlKcy0eO04LoRD545O48TU+lI3qMr6OVtip6Z9WNjl7NgDxRBT5sNd+/sxaOHJjHvj625qEgw5HEgFE9hJZzAhfmgbnQOyFGrgSEjBz29Jvm2UlIWBfft6dcdSOF1W+FxmDMGdVwtxAeTaBRWC0jQK4SkbIjuyZoVmBb0tVkuV3whXD/Sob54b93SDcaA585l+uhTyxHEk1LG1cGeoXacmQ1gOZzAO/Zm+sYep6Xk/ONMyyUKh8Wo5oAHlH4unHNE4ik4rWlRuXa4I2cYrtgcOzPjh7PMPi6l4LSaYDcbM6LbUzN+rIQTeOe+wbw/19dmU+2H37kh3c6/3W7GO/cN4Ikj0wX/hqKga5M3vyjdsd2LY5OrmFqJqFk2HY5cy0UI+ko4gViy+IbsmRI2RLV86MYR+EJxjC2G1pzhIhDCOe4Lqz179DAZDehx23JaNANpISxU5Voq+4Y7cOeO3ppMC6IIvYm4rHR3026IAoBrnZuii4F4hmXhcVqwb6gjIwUOSPdw0UaI+5TLTrfVhNu2rb1dsVeZ15hMSZj1R9DXblNboworKZ6SkJR4wTxoIB2NnZrxV60jXfbAYFGgs7OA6DHGcO1wB7b1utTfm+D33jCMSCKFX56cy/vzlxaDsJoMGCiQ0y8sLwDqB3+73YxoQsrouOgLxWFVrLClAkVSgGx5XVkKY0eRDVEtt27pVjcWO5zr89CFoL96yYdIIqWb4SL4h9/dhwfu3ppz+5YeFzZ0OdBmW38xzoPvvxb/+P7aTMAUQr6WfaFKQYJeIcSGaPY0b6tJbhO6FstFHqKcO/3ktm1eHJ3ILFQRPVw2ZkToHQCAt+zqXVckrJ3XOL0SxUC7PSd7JxyTBclR5LJfNGqaXI6s278ttF6tXXFuLgib2aBu0ubjwfdfix987Mac6O76EQ9GOh144shU3p8dWwhhY7ezYMbIrv421XJIWy7y70DYLimJYyWSwFZlj6KYjy5K4HcU2RDVYjAwfODAcMbzrxVRLSpSMsW69bh1a7euz//A3Vvxk0/euq511ANv3NSFN4x6St7crQYk6BXi5LQfFqNB9wXttppy+p6UglymzXPSoN68TR6d9cKFdLbL2GIQbpsp49iBdhv+8t/twp/elRsVlYO4lJwPxDC7GkWfZuNLfFCFE6UJeq8mTa1a6V3drsxq0XNzcj/2Yul5nU5LTm4zIEfv775uEL+9uJgxE1XLpcVQ0TcyYwzv3DeAawbaVO9a2DzCdlkOx8E5sL1XFvyigj4rrj7Ka//6/v3DGOl05FiE5dJuN6PNZlIrPYt1N9TDajIWbdfbCOwZasePPn5zSXsZ1YIEvUJcnA9ik9cJs053N7dOznYpCNsge/Nn35DcZ/vZc+n+IJcWQ9jkdeVkEPzRLRvzTlEvFfH8oltef7tNtZICaoSebp1biA6HGRbld1TNCF3bz+X8XBDb8ni7pfLuawfAOfDj13Oj9ERKwrgvXFJk9vn7duLHn7xF/V70/BC56MI/396nROhFNkbPzPjhspoyJjeVQrfLiuc+ewdu0RmeXC5DHocaeKy1LwxRGUjQK8RiKK676w4Abpt5TZuiQpSyLReT0YBbt3TjuXOLahXj2EIIm6t0qSc2rU5O+yFxuYdKdjqmGG5RLEJnjKk+ermtc0ul22WFLxRHQhkeMeuPYmsJKX2F2OR1Yd9wB/71SG6R0YQvjKTE82a4aDEaWEYWlIjURbWo8MyFNbFYJEI/PRvA9j53zXqHAGkfffMaonOispCgVwhfKJY3Bc9lNa0pDz09/SRX+N593SBm/VHc90/P4/nzC5hZjVbNuxPPf3xqBYCcU2wzG2A0MHVqUSheWoQOpAtCqmW5iA9WXyiubohuK+Dtlsp7rxvE6Rm/WsgjEBkua/n9q5ZLODNC72uXs24KReicc5yZ8Zfln1cD4aMX8s+JqwMJeoXwBeN5I85SLJcj48tYyRp0IKIzPeG755o+/P8fuxHxpIQPf+tVAPlzoNeL3WKE22pSN3772uXiELct/UElrhS0aYv5EKlrVbNcNMVF53QGLqyVd+zth9HA8ERWlC42pLN7uJSCEPRlVdDlv3mn01J0WMfMahT+aBI7yvTPK81wpxyhr8U/JypLwwn6b87O4+4Hn9WdNl4rookUQvFU3r4kLpupYNpiJJ7C+//5JXzrhcw2u4vBuDrcV4+bt3Tjl//+NnzoxhG4rCbsq2L/CtFUC0hX/bms6Q+qUImWC5BOXaxW2mK3WyvoAdjNxrI9Zj26XFa8eZsXP359KmMYxthiCB6HeU1FOnZlWLLYFBVl/x6HRXeeqxZxpbCzxhG6qPDcVuN1EA0o6MmUXGJciSkvlUItY88TcbbZzAV7uZyfDyCR4ur4KoE8oq1w8ySX1YS/fc8eHP+rt6qXvtVAiKTdbES7spHn0mTvlLopCqQtl7VUiZaCGqEHYzg/H8DW3uIZLqXy7usGMbMaxatKVgcg99BZq93FmDwMYVVjubTbzTArnSkLWS5irF2thfTN23rw0IdvwE2bumq6DqIBBT2dXdE4gu5WIvR8I87OKKlnU1ltcReDuTno+ah2ZZzwpfvb0704tJZLqZuigFxB2uO2ZrTYrcZaheWSr3pxLdy1owdWkwG/ODGr3iYyjNaKx2FRh00safqrdBexXM7MBjDksVekIGc9GA0Mb72mrybVmUQmjSfoSv7zWptdVQNxmVxoU5TzdK52NiKXeGo5U9AXgvGatuLUIqLePk0bVLfNnC4sKmNT9MZNXXj1C3erkX6lsZllz//ifBALgVhFNkQFTqsJt2/34ucnZiBJHMFYEvOB2Lo2pLUNunyasXxetxXheCqjAZqWCyUMtSBai4YT9Er0F6802o0sPUSZfL7URSHoc4GoOpACkDdFa9mKU0s6Qk970bKHLp9TOJ6C2chK6t54NfC6rfjtRbnwqhIbolru29OPOX8MRyaWcWkdG6KCDrtZHbSs7YDozdM5UjC9EqmqzUY0HvXx7isDNUKvI0Ev1tu7WMfFM7MBWE0GcJ4ewMw5ly2XGk4/0SLERTuoQLvZG46nYK9wo6310O1Kb+JW2mO+c0cPLEYDnjw+q/bQ2dhdOculSxOhA/p90UOxJFYjCfTrDMomWpfGE/R1tqOtBr5QHCYDQ5td324otOalYAyLwRhu3ixvKAkfPRhLIpaU6sdycetYLposl7Bmnmg9INbrsprUYcyVwm0z47Zt3fjFiVlcXAiBsdJ7eevRoVguksSxHM60XAD9CF0MGqlE9g7RPDScoFtNRliMa2t2VS2WgnF4nJa8m0JtBVroiuZKd+7sBZD20fNVidaKTV4njAaW0bHQZTWpU4tC8VRNe1hkI8RwS4+rKpt1b9vdj6mVCH52dBqDHfZ1NT9rd5gRS8qzWlMSzxV0nQh9WulLn6+/O9Ga1E9IVQZy1kj9ZLloL5P1cIkyeR2bSPjnd2yX29uKCF1cZlcrV7tcNnQ58fp/fIu6HwCkraRQLIlwLAlnCRuiVwtxZVPJDVEtd+/shcnAMLYYwpu2rq8fiqgzEEO+RT2Dx2GB0cB0I3QxKHyALBdCQ8NF6IDi3dZRhO4LxQpWPRby0M/OBtDplGdkdrus6ht1KZi/SrRWuLPS41zqZm8yZ/xcrRHRbaU3RAXtDrPa2Gq9PUxEg66LiqCLimN58LBFX9BXo2AMugMjiNalMQXdWrjy8mpTbDZjoalFZ2YD2N7rVieoiwh9QbFc6iXLRQ+1hW4sUXeC3qMIXbUEHQDu2yNPY1pvDx3ROlaMsdNe7eUr/59eiaDXbdPt7km0Lg35atCWnNcDxSwXp8UExnJz5yWJ49xcQJ3WPthhS1sugRgYq16/k0qgppBGkwjHk3DU0abom7Z040u/tw+3VqA9bD7etqcf917Thzu29xQ/uACq5aI0+dL+zfNVi86sRijDhcihIQXdXaQ3ytUknpQQiCYLtoI1GBhcFhP8WYI+uRxBOJ7SCLod0ysRNWXR47BktFqtN7Q1AeF4Co46Sls0GQ14z3VDVW0r22Yz4xsfvmHdszBFg64x1XLJEnTdCD2KAcpwIbKoX7UoQD1F6CJ/uLOI163XoEs0VxKCPtBhRzQhYSkUV8r+6zc6B5AxtSgUq6+0xUaiwy7/nadWIjmDs8V8VM7TbSM455heiVQ8HZNofBpT0OsoQk8XFRUWX7mFbqaHLjJchM8rcoqnVyJYDMbrJmUxH9qpRZFEfXnojYTdYlQLy7IDA6/LikSKq5WkgLxnE0tKlLJI5NCQgu62mesmy6VYYy6BPLUoK0KfC2C4065GuuISemo5UlZjrlohphYth+JIpDgJ+joQtku2ddetU1w0o1QTk+VCZFNU0Bljw4yxZxhjpxljJxljD+gcwxhj/8QYu8AYO8YYu746y5VxWU2IpyTEkvrNrq4mS0ofl2IRul5mztnZgDoMGEiP8ppaiWAxUP+CbjMbYDIwzAdkgSmlMRehj7Bdsl9Hev1cKAedyEcpEXoSwGc45zsBvBHAJxlju7KOeRuArcq/+wF8vaKrzEKbXVFrRIRerAAoe2pRLJnCpcVQxviwdrsZDosRFxeCCMVT6HbXt4fOGIPLZlJ7plCEvnbSEXqWoOtUi6YFnSJ0IpOigs45n+GcH1a+DgA4DWAw67B3Afgel3kZQAdjrL/iq1VYT4OuqZUITk37ix9YIr5QHAaWLg7JR7agX1kKIyXxjDmMjDEMdthxdEIe9VbvETog/y3m/UqETpuia0YIenaE3qtMd7q8mB5+MrMahcVkqNqAEKJxKctDZ4yNArgOwCtZdw0CmNB8P4lc0Qdj7H7G2EHG2MGFhYUyl5pGm11RLv/4b2fxJ48cXvNzZ7MUisPjKDxVCBAeenpja1yZTrShK7MoZaDDrvZ3qeeiIoHLqonQ6yhtsdEQueieLJF228zY2d+GVy4tqbdNKRkuNFCCyKZkQWeMuQA8BuDTnPPsEFfvlZUznodz/hDnfD/nfL/X6y1vpRrW03FxORRXM1PKhXOOP/vRUTx1ak69TTuQoBBuTSMrABhXZqKOZE3tGfTYkVImGzVChO62mVQ7wFHCgGhCn/Y8lgsA3LSpCwevLCOqDEiZWY1ShguhS0mCzhgzQxbzhznnj+scMglgWPP9EIBpneMqgrtAs6tihGIpBKKJjLzeUrkwH8Sjhybxw1fH1duKlf0LXFmDOcZ9YbitJngcmVaNth1qvXvogBxBig+gemrO1Wjk2xQFgJs3dyGelHBkfAWA7KGTf07oUUqWCwPwLQCnOecP5jnsJwD+UMl2eSOAVc75TAXXmUG6QrH8jovBWBIST0+pL4enTsuR+aHxZXU+6FIopnbHK0T21KJxXxjDnY6cy2atoNdz2b/ApfHNaVN07XgKROgHNnXCwICXLi4imZIw549ShguhSykR+i0APgzgTsbY68q/+xhjH2eMfVw55kkAYwAuAPifAD5RneXKuNaR5SIiZH+k/A+DX5+eBwCshBNqZ7xSI/TsjovjvnCO3QLIlgsg91C3mupfIMXfAqBN0fWwf9SDN4x6sEVnRmibzYw9g+14aWwJc4EYJE590Al9ir4DOecvQN8j1x7DAXyyUosqRrrL31osF0XQowkMoPQ3xWIwhsPjy3jXtQP48evTOHhlGZu8LqxEEgX7uAjcmo1cSeKY8IVx547cpk7iUrpeRs8Vw62N0GlTdM1s6XHjRx+/Oe/9N23uxjefH8PFeTmQoAid0KMhK0WtJgPMRrbOCL28n33mzDw4Bz526yZ0OS04eHkZy+E4OC9eVARkWi4LwRhiSQnDOhF6r9sKo4E1xIYokL7yAGhTtJrcvLkLSYnjJ0flrSny0Ak9GlLQGWNr6omeSEmIKVkm5Vouvz49j742G3YPtuH6DR4cuuIruewfyNwUzZfhAshdAoc89oZpvCSulkwGBksdd4ZsdPaPemA2Mvz8uLw11d8grw/i6tKwpqfLVn7HxZDmA8CvM2wiH9FECs+dX8B7rhsEYwxvGPXgqVNzanOt0iL0tOVyReSg6wg6ADz04f15B07XG2Jqkd1ipLzoKuKwmHDtcAdeu7wMt82UMz2KIIAGjdABeU5nuYKujejLidBfHltCOJ7C3cog5xs2dAKAmo9erHUukDm1aNwXhoHlv2ze3udumE0vEaFTymL1uWmzPKxjkOwWIg8NK+hrGRQdiqVTFcv5MPj16XnYzUbctLkLALB7sA0WkwHPnJGzXkqxXKwmIyxGAwKxJCZ8YfS322ExNeyvX0V8UFHKYvW5aZP8+iO7hchHw4ZVbqsJc0qXv1LRfgAUslxWwwn82aNHYTEa0NNmxc9PzOJNW7vVwQNWkxH7htrx2uVlAOmy7aJrVmyifCmLjYgq6LQhWnWuG+mA3WxsmtcOUXkaVtBdNhMuLpRruaQj9EJZLkcnV/DUqTn0t9sQiCYRjCXxrmszW9PsH+3Ea5eX0W43lzyo12UzIagI+p3rnENZLwjLhVrnVh+b2Ygfffwm9FGETuShYd+Fa8lyEZuiBlY4Qhdj5X7wsRux2etCIiXliPb+DR4ApW2ICtw2E+YDUSwEYuueQ1kvuMhyuarsHmyv9RKIOqZhTdy1ZLmID4Aet62goKvpiIqVoheB36AIejnl+W6rWW3d2yyXzaKvDm2KEkTtaVhBz+5eWAqiEKm/w1bQchE9ztsL9DjvcFiwZ7A9p/1tIVw2E/zKGppF0MXUIjtF6ARRcxo2rNIOueg0lRYlC8ulv91WcMiFr8Qe59//6IGS/XMgs6qyWQSdMYb+Dhv62sjXJYha07CCLgorgtFkybZHMJ6E1WRAp9OiRsp6LIfjOYMG9OgoMbtF0Kas2W01qRNqmoHH/vhm1XohCKJ2NKygq0MuyshFD0aTcFlNaLOZ4Y/IPdH1qht9objqn1cScVWh1za3kelxU3ROEPVAQ3voQHktdEOxJJxWE9rsZiQljkhCvyd6qS1xy0VYLs1itxAEUV80rKBnTwAqhWAsJQu6Yn3k2xj1hRIlWS7lIta8oUlSFgmCqC8aV9Ct5Qt6KJaEy2pUG1/ppS5yzrEcjqPTWXlPWPj+em1zCYIg1kvjCrpNiHI5EXraQwf0G3T5o0mkJF5yOX85iDRIitAJgqgGDSvo6qDoNXjo2ePgtIiiolLmhJbLzZu78Lfv2YObla55BEEQlaRhs1xEQUs5HRfVCF2JlPUsFyHo1YjQzUYDPnTjSMUflyAIAmjgCJ0xpja7KhU1y6WA5bJcxhQigiCIeqJhBR2QN0ZLHRQtSRyheAoujeWi57/7wtWL0AmCIKpJwwt6qRF6KJ5Uf8ZmNsJqMuhG6NX00AmCIKpJQwu6PLWoREFXeqE7lXTHNrtZ10NfDsVhNRlgN1OzKYIgGouGFnSXtfQWumLz1KlM1mmzmXQLi0SVaDOV5hME0Ro0tqDbzCVH6GJakfDP80bo4Tj55wRBNCQNLejuMoZciNa5YhCD22bW3RRdCsXJPycIoiFpbEG3mkrOQxeRvOqh20wI5ElbpAidIIhGpKEF3WU1IZqQkEgVn1oksmFcRTZFq9VpkSAIoto0tqArfnioBB9dTVsUHrrNDH8kCc65ekwiJcEfTVKEThBEQ9LYgm7N35MlG2G5pCN0E+IpCTHNTNJlpaiokzx0giAakIYW9EJNtrIJRpMwGhisJvmU9cr/l0Py19WYVkQQBFFtGlrQXaLjYimWSywJp8Wo5penG3Slf1ZtzFWFXugEQRDVprEFXY3Qi2e6BGMpdcAEAE0/F02EHqbGXARBNC5FBZ0x9m3G2Dxj7ESe+9sZYz9ljB1ljJ1kjP1R5Zepj0up+gzH9WeDapE7LabL+fUslyXRaZEsF4IgGpBSIvTvALi3wP2fBHCKc74PwO0A/pExdlUU0aEUCYXjuZbLBx56Cd98fkz9Pqi0zhW023M7Li6rlgsJOkEQjUdRQeecPwfAV+gQAG4mm9Mu5djSm5SvA1H1KRpvaXl9YgUvXFhUvxfDLQR6EbovFIfbZoLZ2NBOFEEQLUollOsrAHYCmAZwHMADnHPdSh/G2P2MsYOMsYMLCwvrfmK7RVgumZ8fKYkjmpBwcSGo3hbKFnSdqUXycGiKzgmCaEwqIej3AHgdwACAawF8hTHWpncg5/whzvl+zvl+r9e77ie2mAywGA0IZXnoQuAnlyOIJuT7QlmWi1X5WW3HRR+V/RME0cBUQtD/CMDjXOYCgEsAdlTgcUvCYTUinJW2KCwYzoHLSyEAQCArQmeMoc1uysiQ8YXi6KIInSCIBqUSgj4O4C4AYIz1AtgOYKzgT1QQp8WUE6GHNBbMxfkQOOc5WS6AUv6ftSlKG6IEQTQqpmIHMMYegZy90s0YmwTwlwDMAMA5/waAvwbwHcbYcQAMwOc454t5Hq7iOCzGHA89rNkkvbgQRDQhQeLpQiSB22bK3BQlD50giAamqKBzzj9Y5P5pAG+t2IrKxGE15WS5aCtHxxaCCCgtdl3ZEbqm42I4nkQ0IZGHThBEw9Lw+XlOvQhd+b7DYcbFhVDOPFGB3HFRFnR1ODRF6ARBNCgNL+gOS26ELjz1PYPtuLgQVHuh5wi63aR66KIxF3noBEE0KkUtl3rHac2N0EV/9D2D7Xj+/KKaj+7WidCXQ3E8+NQ5GJSZ0J3UmIsgiAal4QXdoZflogj63qEOAMDRyRUAuRH63bt68dz5RXzl6fOQlDkXXpetquslCIKoFg0v6E5Lbh66aNa1Z6gdAHBsclU+NkvQ3zDaiZ8/8CYEY0kcm1jBaiSBkS7HVVg1QRBE5Wl4QXdYTQgnUuCcq73OQ7EkLCYDBtptcFlNODktC7rLqn+6LqsJN2/pvmprJgiCqAZNsClqBOdANJFuHxOKy1WhjDFs9jrV+0T/dIIgiGak4QXdqTTo0laHhmMpOJTbN3td6u0Oc2YeOkEQRDPR8IKu9kTXpC4GY0m1te7mHlnQnRYjDCKVhSAIoglpeEEX/VkyIvR4Sr19U7dTOY7sFoIgmpuGF3S9qUWheLpVrojQyT8nCKLZaXhBVyN0jeUSiiVVD31DlwMGlj/DhSAIolloeEHXjdBjKTVCt5qMGOl0qJ46QRBEs9LwKqc3VzQcT2YI+J/fs0ON2AmCIJqVhhd0hzV3rmgollJvB4C37+2/6usiCIK42jS85aJG6Eq5fzwpIZ6S4CKLhSCIFqPhBd1mNoAxqP1cIoqwO2gTlCCIFqPhBZ0xljFXNKhYL07yzAmCaDEaXtCBzLmiIlKnQiKCIFqNphB0p2auqIjUnVaK0AmCaC2aQtC1EboYbuGgTVGCIFqMphB0p2auqBB0qgwlCKLVaApBd2jmioppRVRIRBBEq9EUgp6R5UKbogRBtChNIegOi1HNPxeROgk6QRCtRlMIutNqUvuhCy/dTtOJCIJoMZpC0B0WozqxKBRLwm42wkjTiQiCaDGaRtDjKQnxpIRQPEV2C0EQLUmTCLos4JF4Sm6dS0VFBEG0IE0h6Nq5ovK0IorQCYJoPZpC0LVTi0KxFFwUoRME0YI0haBr54qG4xShEwTRmjSFoDvUIRdJBGPkoRME0ZoUFXTG2LcZY/OMsRMFjrmdMfY6Y+wkY+zZyi6xOGJqUTiWQjieooHQBEG0JKVE6N8BcG++OxljHQC+BuCdnPNrAPxuRVZWBo6sTVFKWyQIohUpKuic8+cA+Aoc8iEAj3POx5Xj5yu0tpJR54rGUgjFU9SYiyCIlqQSHvo2AB7G2G8YY4cYY39YgccsCxGhL4fjSEmcInSCIFqSSiifCcANAO4CYAfwEmPsZc75uewDGWP3A7gfAEZGRirw1DIOpW/LQiAGgOaJEgTRmlQiQp8E8AvOeYhzvgjgOQD79A7knD/EOd/POd/v9Xor8NQyJqMBVpNBFXQHRegEQbQglRD0HwN4E2PMxBhzALgRwOkKPG5ZOK0mVdBpWhFBEK1IUeVjjD0C4HYA3YyxSQB/CcAMAJzzb3DOTzPGfgHgGAAJwDc553lTHKuFw2LEQjCmfk0QBNFqFBV0zvkHSzjmiwC+WJEVrRGnxYSplYj8NUXoBEG0IE1RKQrImS7q+DkqLCIIogVpGkHXijiV/hME0Yo0jaBrfXNqzkUQRCvSNIKu9c0py4UgiFakaQRdROiMATZz05wWQRBEyTSN8okI3WkxgTEaEE0QROvRNIJuV8r/aUOUIIhWpWkEXQg5pSwSBNGqNI2gi8wWB0XoBEG0KE0j6BShEwTR6jSNoIsIncr+CYJoVZpG0EVkTo25CIJoVZpG0IV3TkVFBEG0Kk0j6OkInQSdIIjWpGkEXVgtlIdOEESr0jSCLjZDKUInCKJVaRpB9zjM+MxbtuG+PX21XgpBEERNaJpwljGGT921tdbLIAiCqBlNE6ETBEG0OiToBEEQTQIJOkEQRJNAgk4QBNEkkKATBEE0CSToBEEQTQIJOkEQRJNAgk4QBNEkMM55bZ6YsQUAV4oc1g1g8Sospx6hc289WvW8ATr3cs59A+fcq3dHzQS9FBhjBznn+2u9jlpA5956596q5w3QuVfq3MlyIQiCaBJI0AmCIJqEehf0h2q9gBpC5956tOp5A3TuFaGuPXSCIAiidOo9QicIgiBKhASdIAiiSahbQWeM3csYO8sYu8AY+4tar6daMMaGGWPPMMZOM8ZOMsYeUG7vZIw9xRg7r/zvqfVaqwVjzMgYO8IY+5nyfUucO2OsgzH2KGPsjPL3v6kVzp0x9u+V1/oJxtgjjDFbs543Y+zbjLF5xtgJzW15z5Ux9nlF884yxu4p9/nqUtAZY0YAXwXwNgC7AHyQMbartquqGkkAn+Gc7wTwRgCfVM71LwD8mnO+FcCvle+blQcAnNZ83yrn/t8B/IJzvgPAPsi/g6Y+d8bYIIA/BbCfc74bgBHAB9C85/0dAPdm3aZ7rsr7/gMArlF+5muKFpZMXQo6gAMALnDOxzjncQA/BPCuGq+pKnDOZzjnh5WvA5Df1IOQz/e7ymHfBfDumiywyjDGhgC8HcA3NTc3/bkzxtoA3AbgWwDAOY9zzlfQAucOefSlnTFmAuAAMI0mPW/O+XMAfFk35zvXdwH4Iec8xjm/BOACZC0smXoV9EEAE5rvJ5XbmhrG2CiA6wC8AqCXcz4DyKIPoKeGS6smXwbwWQCS5rZWOPdNABYA/H+K3fRNxpgTTX7unPMpAP8AYBzADIBVzvm/ocnPO4t857pu3atXQWc6tzV1fiVjzAXgMQCf5pz7a72eqwFj7B0A5jnnh2q9lhpgAnA9gK9zzq8DEELz2Ax5UfzidwHYCGAAgJMx9ge1XVXdsG7dq1dBnwQwrPl+CPJlWVPCGDNDFvOHOeePKzfPMcb6lfv7AczXan1V5BYA72SMXYZsq93JGPsBWuPcJwFMcs5fUb5/FLLAN/u53w3gEud8gXOeAPA4gJvR/OetJd+5rlv36lXQXwOwlTG2kTFmgbxR8JMar6kqMMYYZB/1NOf8Qc1dPwHwEeXrjwD48dVeW7XhnH+ecz7EOR+F/Dd+mnP+B2iNc58FMMEY267cdBeAU2j+cx8H8EbGmEN57d8Fed+o2c9bS75z/QmADzDGrIyxjQC2Ani1rEfmnNflPwD3ATgH4CKAL9R6PVU8z1shX1YdA/C68u8+AF2Qd8DPK/931nqtVf493A7gZ8rXLXHuAK4FcFD52z8BwNMK5w7gPwE4A+AEgO8DsDbreQN4BPJeQQJyBP7RQucK4AuK5p0F8LZyn49K/wmCIJqEerVcCIIgiDIhQScIgmgSSNAJgiCaBBJ0giCIJoEEnSAIokkgQScIgmgSSNAJgiCahP8DuXTh7V5pOSMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, nrgz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfec3d41",
   "metadata": {},
   "outputs": [],
   "source": [
    "kde"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}