{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Utils: Compute and plot free energy surface\n", "\n", "Authors: Ioannis Galdadas & Luigi Bonati" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/luigibonati/mlcolvar/blob/main/docs/notebooks/tutorials/utils_fes.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following example we show how to use the function `compute_fes` from `mlcolvar.utils.fes` to calculate and visualize the free energy surface." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from mlcolvar.utils.plot import paletteFessa\n", "from mlcolvar.utils.io import load_dataframe\n", "from mlcolvar.utils.fes import compute_fes\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# Load COLVAR file containing collective variables (and bias information)\n", "colvar = load_dataframe('data/muller-brown/biased/opes-y/COLVAR')\n", "\n", "# Simulations parameters\n", "temperature = 300 \n", "kb = 0.0083144621 # Boltzmann constant in kJ/(mol·K)\n", "kbt = kb * temperature \n", "\n", "# note, for a toy model you should use directly kbt = temp\n", "kbt = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate statistical weights in case of a biased simulation (uncomment one of the options)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# (1) unbiased simulation\n", "#bias = None\n", "\n", "# (2) reweight for a single bias\n", "bias = colvar['opes.bias'].values\n", "\n", "# (3) reweight multiple bias potentials (e.g. OPES and harmonic walls)\n", "#bias = colvar[['opes.bias','lwall.bias','uwall.bias']].sum(axis=1).values\n", "\n", "# (4) reweight all field *.bias in the COLVAR\n", "#bias = colvar.filter(regex='.bias').sum(axis=1).values\n", "\n", "# calculate the weights\n", "weights = np.exp( bias / kbt )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate and plot 1d free energy surface. For the full list of options see the documentation. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG3CAYAAABYEDo3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAABcSAAAXEgFnn9JSAACEFElEQVR4nOzdd3hc5Zn4/e+Z3iSNei+23HsBGwMugKnBgSSkLRtCSQLZFBKyu9lsILDZZDe/FBIWkjdlEwikU5ZggwHTDLhi3LslW9Xq0kjT63n/GPtYwnKTRhqN5v5cl+E8zzlz5p7xWLrnqYqqqipCCCGEEGlIl+wAhBBCCCGSRRIhIYQQQqQtSYSEEEIIkbYkERJCCCFE2pJESAghhBBpSxIhIYQQQqQtSYSEEEIIkbYkERJCCCFE2pJESAghhBBpSxIhIYQQQqQtSYSEEEIIkbYkERJCCCFE2pJESAghhBBpy5DsAMaToqIivF4vFRUVyQ5FCCGESCkNDQ3Y7XZaW1tH9XmlRSiBvF4v4XA42WEIIYQQKSccDuP1ekf9eaVFKIFOtgTt27cvyZEIIYQQqWXmzJlJeV5pERJCCCFE2pJESAghhBBpSxIhIYQQQqQtSYSEEEIIkbYkERJCCCFE2pJESAghhBBpSxIhIYQQQqQtSYSEEEIIkbYkERJCCCFE2pJESAghhBBpSxIhIYQQQqQtSYSEEEIIkbYkEUoBsXCIWCiY7DCEEEKIcWdcJEJdXV0UFBSgKAqTJk0667VPPPEEixYtwuFwkJOTww033MDGjRtHKdKhcb2+mu61T0syJIQQQiTYuEiEvvGNb9DZ2XnO6772ta9xxx13sHfvXlauXMmiRYtYt24dy5Yt4/nnnx/5QIco6ukl0ttNqKUh2aEIIYQQ40rKJ0Kvv/46v//97/n85z9/1utee+01HnnkEXJzc9m1axfPP/88L7/8Mm+//TZ6vZ477rgDl8s1OkEPQbSvl2BTXbLDEEIIIcaVlE6E/H4/d999NzNmzOCf//mfz3rtww8/DMD999/P5MmTtfolS5Zwzz334HK5+O1vfzui8Q5H1O0i1NpMLBxKdihCCCHEuJHSidB//Md/cPToUX75y19iNBrPeJ3f7+eNN94A4JZbbjnt/Mm61atXj0ygCaBGIkQ9vYRam5IdihBCCDFupGwitHv3bn7yk59wxx13sHTp0rNee+jQIYLBIPn5+ZSVlZ12fsGCBdo9xyJVVVHDYaJ9LkJNx5IdjhBCCDFuGJIdwFDEYjE+97nP4XQ6+eEPf3jO6xsa4oOMB0uCAOx2O06nk56eHtxuNxkZGQmNd6hCLY10PP2/uN/fAKqKYjETbGkiFgqiM5mTHZ4QQgiR8lIyEXr00Ud57733ePzxx8nNzT3n9R6PBwCbzXbGa+x2Oy6X67wSoZkzZw5aX1tbS3V19TnjOV+K0UTfxte1cizgJ9LVRu/6tWQtvQad5cyvRwghhBDnlnKJUENDA/fffz/Lly/n9ttvT3Y4I8qYV4ixqIzwiXFBOqOZUHsLajSCGg6RufRaDBlZSY5SCCGESF0plwh96UtfIhQK8ctf/vK8H+NwOADw+XxnvMbr9QKcV7fYvn37Bq0/U0vRcNhnLsB1IhGK9vVgn3cJwYZa/OEwsUiYrKXXYszOS/jzCiGEEOkg5RKhNWvW4HQ6ueeeewbUBwIBAJqbm1mxYgUAf/nLXygqKqKiogKApqbBZ1x5vV5cLhfZ2dljZnzQSfaZC3G9/gIA4Y5WFKMJS9VkAg21BI4egmgU5xU3YnDmJDlSIYQQIvWkXCIE4HK5WL9+/aDnAoGAdu5kcjR16lTMZjMdHR00NzdTWlo64DHbt28HYM6cOSMY9dDYZszXjmMBHzFPH/qMLCyVkwk21hJqbcRfe4CMhZclMUohhBAiNaXc9HlVVQf9c+xYfFp5dXW1VldVVQWA1WrlyiuvBODpp58+7Z7PPPMMAKtWrRqdF3EBDM4c9JlOrRzqaAFA0esx5hYScfcRbmtGVdUkRSiEEEKkrpRLhIbqvvvuA+B73/seR44c0eo3bdrEr371K5xOJ3fddVeywjsrY0GJdhzuaNWOdXYHasBPpM9F1NOXjNCEEEKIlJY2idDKlSu599576erqYt68edx8883ccMMNLFu2jEgkwuOPP47T6Ux2mIMy5hdrx+GOFq31R9Hp0dlsRL192swyIYQQQpy/tEmEAH72s5/x+OOPM336dNatW8emTZtYuXIlb7/9NjfffHOywzuj/i1CajBAtLdbK+vtmUQ9fYTampMRmhBCCJHSUnKw9GCqqqrOa5zM7bffnnLrD5nyCtHZM4h53QAE6mtwOOMLSeodmYQbOuItRdEIin7c/JUKIYQQIy6tWoRSlSEnD1NRuVYONtSiRiMAKGYLAFGPm3BnW1LiE0IIIVKVJEIpwJCTj7lyIujif11qKEjoeHz/NEVR0NsziHr7CLUdT2aYQgghRMqRRCgFGLPz0Tuc6DOztbrAscPasd5mJ+b3EXV1JSM8IYQQImVJIpQC9JlO9FYbhpx8rS7c0aJNmVcs1viGrP0GUQshhBDi3CQRSgGKTochOxdDbgE6m0Or9x+J73mmM1tQwyGiPi+xgD9ZYQohhBApRxKhFGHIKcBgc2AqKtPqAnVHiAX8KDo9itFELBiQViEhhBDiAkgilCKMOXnorDZ0GVkoZmu8MhbFX7MfONEqFPQT6e1JYpRCCCFEapFEKEUYcvLRWWwQDmGdNF2rD9QeIBYOaeOEotIiJIQQQpw3SYRShM5qR293oJjMmIorUAxGANRImGDdEXRm64muMWkREkIIIc6XJEIpQlEUjLmF6Kx21HAIy8Sp2rlg0zF0J1uE+npQY7EkRiqEEEKkDkmEUoghrxC93UHU58FcPlGrj3R3EAuHQY0RDQZkJ3ohhBDiPEkilEKMeUXorHZifi+6zGx09gztXPh4fXzAtIwTEkIIIc6bJEIpxODMQW9zgKKDUBBzWZV2LthchyLjhIQQQogLIolQClF0Ogy5+eht9nj3WGmVdi7S1Q4gK0wLIYQQF0ASoRRjzCtCZ3MQ83vRO3MHdI9FejqJBf1EpUVICCGEOC+SCKUYY14hepuDmNcDcFqrUCwYIOLpIxYKJilCIYQQInVIIpRijDn56Kw2VDWGGg4N2HIj3NWOotOhBgPSKiSEEEKcB0mEUoxiMMZXmbbaifk8GHILtMUViUVRQyFiQT8RV1dyAxVCCCFSgCRCKchUWIo+I4uouxdFp8OYX6ydi3r7TgyYlhYhIYQQ4lwkEUpBpuIK9I5Mol43aiyGsbBEOxdxdccTIWkREkIIIc5JEqEUZMjOxZDpRDGZifk8mApLtXMxr5uIp1e22hBCCCHOgyRCKUhRFEzF5Vr3mN6ROWAafazPRSwYIOrpTWKUQgghxNgniVCKinePZRH19KGq6sBWIZ9XxgkJIYQQ50ESoRRlKihB78hEVVXUoB9jv0Qo4nYR9XuJyjghIYQQ4qwkEUpRisGAqbAEfUYWkT4XpoIS0OvjJ6NRwu0t0iIkhBBCnIMkQinMUlGNwZkbnyGm18WToRPCHa1EejpRVTWJEQohhBBjmyRCKcxUWokxOxfFaCLq7sVUUqGdi/Z2E/G4iXndSYxQCCGEGNskEUphik6PZcJUjNl5RLo7MRWVa+fUSJhwezPhE7vSCyGEEOJ0kgilOMvEqRiycogF/aAoGHILtHPB5npJhIQQQoizkEQoxeltDkyllRicuYS72jAVn+oei3S0EO5sS2J0QgghxNgmidA4YJs6B2NeIVF3L4acfK0+5vcRrD+CGgknMTohhBBi7JJEaBww5hVirZ6OsaCEqNs1IBny1x4k3N2RxOiEEEKIsUsSoXHCPvvi+PR5RcGQV6jVh5rrCLcfT2JkQgghxNglidA4obNYsc++GHNJJYreALr4X60aDuHevinJ0QkhhBBjkyRC44hl4lRs0+dhmTAVfUaWVu/dtYVYwJ/EyIQQQoixSRKhcURRFBwLLsU+fS7WybO0+tDxeno3vJbEyIQQQoixSRKhcUZRFBwLLydr2bXoLLZ4parS+/ZagscbkhucEEIIMcZIIjQOKTodGRcvJ2PRcq0ucPQg7m1vo0YiSYxMCCGEGFskERqnFEUh/9P3DBg07du3k8CxQ0mOTAghhBg7JBEax4zZudimz9fKgdoD+A7vQY1FkxiVEEIIMXZIIjTOZV//ce046uklcOwIwfraJEYkhBBCjB2SCI1zjgWXYsjO08r+w7vxHdqNGoslMSohhBBibJBEaJzT6XRkXHKFVg53thFqbSLc0ZLEqIQQQoixQRKhNJC98iYUgzFeiMXw1xwg1NqU3KCEEEKIMUASoTRgKqnEMmmGVg41HSN4vDGJEQkhhBBjgyRCaUBRFLKWX6+VYwEfgaMHiHr6khiVEEIIkXySCKUJx9xLMOQVaeXA0UOEWqRVSAghRHqTRChNGLJzsU2bo5VDrU0Em+qSF5AQQggxBkgilCYURSHzsqtBp49XRCN4dm9FjYSTG5gQQgiRRJIIpRHrpBmYSiq0cuDoAekeE0IIkdYkEUojhowsHLMv1srh9hY8u99DVdUkRiWEEEIkjyRCaSbrig+hmC3xgqri2blZFlcUQgiRtiQRSjPm8olY+60pFKw7jP/g7iRGJIQQQiSPJEJpRlEUnFfdpJWj7l58h/cQ7u5IYlRCCCFEckgilIYyFi9H78zRyr6Du3FveYtYMJDEqIQQQojRJ4lQGtLpDWQsuEwrh9uP42+opW/TG6ixaBIjE0IIIUaXJEJpynnlKu1YDQUJ1h0hUHcE797tSYxKCCGEGF2SCKUpy8RpGIvKtLIaDBA8XofvwE4ZLySEECJtSCKUphSdjoyFp7rHQm3N6EwWQm1NeLa9K11kQggh0oIkQmksc+l1KEZTvBCNoIaCRN29BI834D+0J7nBCSGEEKMgZROhhx9+mI9+9KNMnjyZrKwszGYzlZWV3HbbbezZc+Zf4k888QSLFi3C4XCQk5PDDTfcwMaNG0cx8rHDUlaFecIUrRyor8FYUEKotRHfoT3EwqEkRieEEEKMvJRNhP7rv/6LtWvXkpOTw1VXXcWHPvQhLBYLTz31FAsXLmTNmjWnPeZrX/sad9xxB3v37mXlypUsWrSIdevWsWzZMp5//vnRfxFJphiMZFy8TCtHe7tRo1FQVcLdHQSOHkpidEIIIcTIMyQ7gKH6+9//zsKFC7FYLAPqf/GLX/ClL32Jz33uczQ1NWEwxF/ia6+9xiOPPEJubi6bNm1i8uTJAGzatIkVK1Zwxx13sGLFCpxO52i/lKSyT5tLtzOXqKsLiK80bZ04jUhXG/4je7FOmo6iT9mPiRBCCHFWKdsidNlll52WBAH80z/9E9XV1bS1tbF//36t/uGHHwbg/vvv15IggCVLlnDPPffgcrn47W9/O/KBjzHGghIsE6dp5WDjMRSbg1goRLirnWDD0SRGJ4QQQoyslE2EzsZoNAJgMsUHAvv9ft544w0AbrnlltOuP1m3evXqUYpw7NBnOrFNmQl6fbwiGiHUXIcxJ59wVzu+I3tld3ohhBADqKqK/8g+wt2dyQ5l2MZdIvTUU09x6NAhJk+erLX8HDp0iGAwSH5+PmVlZac9ZsGCBQDs3p1+m48qioKpqAxjQYlWFzh6CL0zh5jXTbiznUhP6n/QhRBCJE64oxX3jk14dm5KdijDlvKDP370ox+xb98+vF4vBw4cYN++fZSUlPDnP/8Z/YlWjoaGBoBBkyAAu92O0+mkp6cHt9tNRkbGqMU/FhjzS7BUTSHc0gjEB01H3b3oM7KI9nYTOHYYY05+kqMUQggxVkTdvajhEEQiyQ5l2FI+EXrllVd4/fXXtXJlZSVPPvkkCxcu1Oo8Hg8ANpvtjPex2+24XK7zSoRmzpw5aH1tbS3V1dUXEv6YYCwoxlRYgs5qJ+b3AvFWIdu0OYRaGgk2HsUxdzGKIeU/LkIIIRIg5nPHE6FxIOW7xl577TVUVaWnp4e3336byZMns3z5cr7//e8nO7SUoXdkos/IwlBQrNUFG4+C0QSqSqS3m2BzXfICFEIIMaZEvZ5xkwiNm6/4TqeTpUuX8tJLL7FkyRIeeOABrrnmGi6++GIcDgcAPp/vjI/3euMtIefTLbZv375B68/UUjTWKYqCMb8Yc0kloea6eFNnNEKo6Rh6Zw4RVzeBuiNYKiclO1QhhBBjQNTrJhYaH4lQyrcIfZDRaOSTn/wkqqpqs8AqKioAaGpqGvQxXq8Xl8tFdnZ22o0POslUWIoh04khO0+rCxw7jCErh2hfD+G2ZqI+bxIjFEIIMVbEfOOnRWjcJUIAeXnxX+YdHfFd1KdOnYrZbKajo4Pm5ubTrt++fTsAc+bMGb0gxxhTYSl6ewaGzGytLurqIuZ1ozNbibhdhJqPJTFCIYQQY4EaCRMN+CURGsvWr18PoA1ctlqtXHnllQA8/fTTp13/zDPPALBq1apRinDs0VmsGPIKMeQWoM86lQz5jx5Cn5VNtNdFoFESISGESHfRE61B42WNuZRMhDZs2MDLL79MLBYbUB8Oh3n00Ud56qmnsFqtfPKTn9TO3XfffQB873vf48iRI1r9pk2b+NWvfoXT6eSuu+4anRcwRpmKytA7MjHmFmh1wcaj6Kw2ot4+wh0tRL3uJEYohBAi2aKe8TNjDFJ0sPSRI0e44447yMvLY+HCheTm5tLZ2cmePXtoaWnBYrHwxBNPUF5erj1m5cqV3HvvvTzyyCPMmzePq6++mlAoxLp161BVlccffzzt9hn7IFNROXpHJorZimIyo4aCJwZN18WTIXcvwaZj2KambxeiEEKku/E0PghStEVo+fLl/Pu//ztTp05l9+7dPP3002zYsIGcnBy+8pWvsGfPHj7xiU+c9rif/exnPP7440yfPp1169axadMmVq5cydtvv83NN988+i9kjDFk52LIyEJntmAqrdTq/bUH0Gc4ifT2xKfVCyGESFtRr2fczBiDFG0RmjBhwpDXCbr99tu5/fbbExvQOHFyu41AfS2K3kjw2BFAJeZ1E/P7iPk88S033L0YMrKSHa4QQogkiHndqOFgssNImJRsERIjx1Rcjj4jCzUcxFh8qmsxcOwQOqudqLs3vtaQEEKItBT1uVHHUYuQJEJiAFNROYYTs8bM5RO0+nBbM4pOR6Svh2BTXZKiE0IIkWzjaVVpkERIfIBiMGAurUSflYOiN6DPytHOBY/Xx7vHutqIevqSGKUQQohkiIWCxIIBYpFwskNJGEmExGnMFdUYsrKJul1Yp8zS6kNN9aDTxWePNdcnMUIhhBDJoM0YGydrCMEwBks/+eSTCQvitttuS9i9xPAZ84sxOHNQjCb0WTno7BnEvG5AJdLVgd7qINhch23q7GSHKoQQYhSNt24xGEYidPvtt6MoSkKCkERobFF0OszlEwk21RFzu7BNmY1nx0YAQq2N6BwZhDtaifo86G2OJEcrhBBitEQ9fcRCQRR9Sk46H9SwXsncuXO56aabhvz4559/nt27dw8nBDFCLJWT8O3fib+zFXPFJHQHbMQCPojFiHR1EPX0Emyqw9av60wIIcT4FnW7UENBdGbzuOkdG1YiNG/ePB588MEhP76urk4SoTHK4MzFVFpBuLOVqKsL2/R5WqtQpKuNUFszIUmEhBAirUT7eokFAygmC2owkOxwEmLIg6UzMzOx2WzDenKr1UpmZuaw7iFGjm3aPAy5BUR6ezCWVqFznPq7Ch47QqijhajPk8QIhRBCjBZVVYm4e4gFA+jMlmSHkzBDToRcLhePPfbYsJ78F7/4BT09PcO6hxg5xvwizEXl6DOdRF2d2Gcu0M5F3S6CzcdkcUUhhEgTsYCfWCAAkTA6kznZ4SSMTJ8XZ6QoCtZpczDmFRJxdWHILcDgzNXO+w/sJtBwLIkRCiGEGC1RtyveLWY0gTJ+0ofx80rEiDAVl2OpnISxoIRQcz3WGfO1c1G3C+++94n6vEmMUAghxGiI9rlQQ+OrWwwkERLnoCgKGRcvw1xahc5qh0gYY36xdt63+z0CdYeTGKEQQojREHWfHCg9frrFYBiJkF6vH/Ifg2H8rD+QDnRmC5mXXIG5tBJVVdE7T227EfX04npjDWoslsQIhRBCjLRInws1GBx3LUJDzkjKy8sTtqCiGPuMuQVkLb0W3bZ3CLY0Em47TrQvPtDds30DwcZjWCqrkxylEEKIkXJyjJA+Ow//kb3xWWR9PaiqmtL5wJATobq6ugSGIVKBqaiM7Ks/gmfnZhQUXK+/AKjEfB66X/wzJf90f7JDFEIIMQJioSBRn5dYOIgaChA4ehCAYP0R8m7+DBiMSY5w6GSMkLggOrOFzMUryPvIbVinzdHq+7a8RaijNYmRCSGEGCnxgdJBdHoDUXevVm9w5qKkcBIEkgiJITLmFVJyz7fgxH4zajBAx19/k+SohBBCjARtoLTZQqS3W6s35OQnMarESHgitHv3bu6++25mzJhBVlYWWVlZzJgxg3vuuUe20xhnTEVlZC29Rit7TowfEkIIMb5E+k6tKB11nUqEjDkFSYwqMRKaCD3yyCNcdNFF/O///i8HDx7E7Xbjdrs5ePAgv/71r7nooot45JFHEvmUIsnyPv550MU/RrGAj87nn5QZZEIIMc5EejqJBXyntwjlSouQZt26dXz961/HZDLx9a9/nR07dtDT04PL5WLnzp184xvfwGw2c9999/H6668n6mlFkply88m4aKlW9u7YhL9mfxIjEkIIkUiqqhJxdRML+OLlUFA7Z5SusVMefvhhDAYDr776Kj/+8Y+ZO3cuWVlZZGZmMmfOHH70ox/x6quvotPp+MlPfpKopxVjQO5Hb9eOo30uetevRY2EkxeQEEKIhIm6e4n5fRCNxv9/gmIyozNbkxhZYiQsEdq6dSvLly/n0ksvPeM1S5YsYcWKFWzZsiVRTyvGAGvVZMxVU7SyZ/dWaRUSQohxIuLqIhbwoTNbifbrFtNZ7UmMKnESlgj5fD7y88/dRJafn4/P5zvndSK15N74Ke040tGCe/tGYv2aT4UQQqSmSE8nUb8PncU6YHyQJEIfUF5ezqZNm4hEIme8JhKJsGnTJsrLyxP1tGKMyLx0JfoMp1b27n4P/5F9yQtICCFEQkR6TrQIWW1EXJIIndFNN91EfX09d955Jy6X67TzfX19fP7zn6ehoYGbb745UU8rxghFpyPz0qu0crDpKP7De1CjZ06MhRBCjG3xgdLxREjRG4l53do5nW18JEIJ2/30W9/6Fs899xx//OMf+fvf/851111HVVUVAPX19bz88sv09fUxceJEvvWtbyXqacUYknPjp+h57XmIRiESwXtgJ44Fl2KpnJzs0IQQQgxB1NOnDZSOBvoNazEYxs0u9AlLhHJycnjnnXe4++67efHFF3n66adPu+ZDH/oQv/rVr8jOzk7U04oxxJRfjHXKbPwHdgIQOHoIf80BSYSEECJFnVw/SGe2EnV1afUGR1ZKb7TaX8ISIYCSkhJWr17NsWPHePfddzl+/LhWf/nllzNhwoREPp0Yg5xXfVhLhGKePvw1+3HMXzIu1poQQoh0oyVCFiuB+matXu/MTWJUiZXQROikCRMmSNKTpjIXLafdmat9c/Af3kug9oAkQkIAaiyGd+/7BBtqUAxG9PYMbDMXyL8PMWZFujuI+n3oHZmEu9q0ekN2XhKjSizZdFUklM5kxjFvsVYOtTbhr69BPctsQiHSQSwUpO/dV/Hs2Ih37/t4dm3FvWMTrrdeJNTafO4bCDHK1Fg03iLk94Kqogb82jmDtAid2fr161m/fj0tLS0Eg4OvI6MoCr/97W8T/dRijHBeuYret1+BWBSiEQI1Bwi1NGAun5js0IQYNjUWjSf2qorObDmvx8QCflxvryVYX0OorRlTYRmK0Ui4q53A0UMQjZJ56VWYS6tGNnghLkDE1U3U7wNVJdLXo9XrM7LQjZOB0pDARKi7u5uPfOQjvPvuu6iqetZrJREa36zV0zGXTyRYfwSAQN1hgo3HJBESKU2NRvDs2Eyg7vCJjYVVDBlOTMXlWKqmYHDmDPq4WChI7zuvEDh2mHBnG5bKSegsNiA+/TjU3EDg2GEAspZei6mobGRfRywKgKLTj+jziNQX6Won5veis9oJd7Ro9YbcwiRGlXgJS4S+/vWv88477zBz5ky+8IUvMHHiRBwOR6JuL1KIojeQsXiFlghFXV34jx3CcfFSdEZTkqMT4sJF/T76Nr5GoL6G0PF61GgUFAW91Uag/gi+g7sxl1VhnTQTY2GJNpsm3NGKe/sGAvU1p5Kgfq1IiqLDVFpJqLmeQH0N6HQ4l12PMb8oofHHwiECtQcItTQS6elEVVXsMxZgnTobRScjJMTgwt0dxPw+9FY7/s52rd6YJ4nQoFavXk1ZWRmbNm2SBEiQeelKul/4o7Zbsb9mP6HjDVgqJyU5MiEunHvzG/iP7CPc2Ya5tAqd3QExlajXTbSvB3/HPsIdLQTqjmDIzsPgzEENhwi1NBFqbybqcWOpqB60K01RFEylFQQbjxGor8H19loyl1yFuaRiyPGqqooaCRPp6STcdhz/0YOEu9rjC+P5vCh6PVFPH8HmOjIvuxq91Tact0eMU5GudqI+LwZnDlFPr1ZvlBahwUWjUZYsWSJJkADAVFiCpXoavn3bAQg2HiPYeFQSIZFyIn09hNqaCXe0YJkw9VQyowdDphNDppNYKEikq51A3WF0xxtQjCZQVaI+N/oMJ9aJ01AMZ/5xqyg6zGUTCDYdI3D0EGokQsbFS7FOnHb+cbq6CTTUEmquI+rzxscyBfxEfR4iri6IxTDmFqArKCHm9xFsPBZfLVinJ2vF9dJVJgaIBXxE3L2oQT/R/qtJW6zo7A5iXk8So0ushCVCCxcupLW1NVG3EylOURQyl1ylJUIxTy/+mgNkLFo+rgbZifEvUFdDpLcHvT3jjIOjdSYzpuJyjIWlRD19EI2AosOYX4TOYj2v51F0OszlEwi1NBI4dhg1GiHq7sM++6Kzdl+p0QieXVvxH95DpLeHaG8PsVAQNRpBZzCiszsw5haiz3RqXXY6swWdzU6g7jABsxX9zi1kLLj0wt8cMW6FuzriibLJTKS7Q6s35BaOm4UUT0pYInT//fdz3XXX8fLLL3Pdddcl6rYihTkuXor+r78m6o43qfpr9xNqbcJSUZ3kyIQ4P2osRrChhkhv93l1Byg6HYZM55CfT1F0mIorCHe2Ejh2CDUSJup2nfELRNTrpm/zmwQbagkdr0dntWPILUBnsaIYDGdt5dGZzJiLKwg116GYzBiyc7FOmDrk2MX4Eu43UDpQe1CrN+YXJzGqkZGwROjKK6/kT3/6E7fddhs33HADV199NaWlpejO8E1m2bJliXpqMUYZs3KwTp2NZ9u7AISa6gg210kiJFJGuL2FiKsbNRRE78gcledUFAVTfjE6k4VgQy2xYIBon4uMS67EmHNqEbvg8QbcW9cTPF5PpLsTc2nlBceoz8jCEMgn2HgUjzG+wKOpoCTRL0mkoEh3OzGfD8VoIOp2afWmwvH3+UjoOkIejwej0chTTz3FU089ddZro9FoIp9ajFFZy67TEqGY34vvwC4yFy1H0Y/IouZCJFSg/giR3m4MmdmE21vw1+wjFvSjoKBYbBjzizAVlmDIGnzq/HAYsrLRmcwEm47h83mIuF2YCssw5uYT7uog1NJAsLkeYlEsE6YMucvZkFdILBgg0HgMZdMbOK+4cVitWiL1qdHIiRWlPSjBU62KOpsDnT0jiZGNjIT9NnriiSe46667UFWV+fPny/R5AYBjziIMuQVEuuJTL/1H9hFqOz6sGTFCjAY1FiN0vIFwZxuRzrYB66jEdRFubcS3J77KrmXiNMzlE1AMxsHvp6rEvG4ifS50ZjOG7PxzTl3XWW1YJk4l1NKE//BeQi1N6KxWYj4vUa8bQ1YuxsJiFOXcU+BVVSUW8BHzedHbHOhOzBRTFAVTSQXB+hqCDbX0vr2WrGXXSzKUxiI9XScWUoxPoT/JWFAy7sYHQQIToR/+8IeYzWZeeuklVqxYkajbihSns9iwz1xI79trAQg11xFsOiaJkBjzop5eIn4vwboj8S0GziLi6sKzfQPe3VsxV1RjKipDn5WNGokQdbsINdcTam1CDYe0xygmM6aSSuwz52sLLA5G0Rswl1URC4eIurqJBfzo7BkYC0pOawVSo1GiXjdqKEAsGEQNBYh6PfGp8z2d8UHcJ+iduZjLJmCdNANFr8dcPoFAfU18cUdVJWv59Rgys4f47olUFu5qI+bzoLPaCB/Zp9WPx24xSGAiVFdXx/LlyyUJEqfJuuJD9L7zcnyvmlAQ764tZFy0VBZyE2NapKeLUEPtgCTIkFcYn9KuKERc3YTbj8cTjBPUSJjA0YMEjh4c7JYDqKEgwbrDhJrrsM9djLmi+qzftnVGE7ozLLQY7u4gcOwwoaZjqJHweb2+qKsLn6uLYP0RHAsvw5hbiKVyEoH6WvwnV7pe8SEMGVnndT8xfoQ724j6vKDGUIOn9hczjtPxYwlLhEpLS7HZZFEucTrb5FmYisoItTQC4D8cX3zOVFia5MiEOLNwTwf+ft+GjUVlZF66UktWzGUTgIVEersJHD1EsKH2vJIQnc1BzO8DNQaAGg7h2fYOwaZjOOZfit5mP+8Y1UgY7+6t2hYd50MxGAZsghx199L71ktYJk7DNmvhiWSohkDdEXj7ZZwrbkA/DseFiMGpqhrfWsPnIdZ/k9XsvHG79EnCEqHbbruNn/70p3R3d5OTk/iBgyJ1KQYDjvlL6D6RCAVb6gkcOyKJkBjTPO9vJNZvITn7zAWDttgYsnJwzF+CffZFBBuPEmppItzTEd+pW9Ghs1ox5hZiKq2Md2cZTcTCIYL1Nfj2va8lJeHWJlzr/g9L9XQsE6eit519jGW4qw33e+8MiBEARUFntqKYzehMFhSzBYMzF2NuAfqsbBSDkai7l0DtgQEtV4GjBwm1NGCffTHmiokETyRDve+8gvPKVcP+JaiqKlFPH9HebhSTBVPB+JuGPR5E3b1EvR7UcIhIz8DxQQOdfU/RVJKwROjf//3f2blzJ1dccQWPPPIIy5cvH5eDqsTQZF1xI90vPxvfkT4Swf3+O2QsWiqr2YoxSVVV3O+t18qm0koMztyzPkYxGLFMmIplwtT4xtPRCOgNg/4c1BlNWCfNwFRSgWf7RsJtzfHnjYTxH9qN/9AejAXFmEorMRWVD2glioWC+A/vwX9oL/1/GemzcrBOmYW5pOKMA7ZPMmQ6ccxfgrmiGs/7G7Tp0TG/D/fW9RhyC7DNmE+0t5tg41HcW94i87Krh9ydHentoW/ja0Rc3cSCfnRmK5mXXomlcvKQ7idGTrirLT5bzGQh3Nmm1ZuKBn5xjfo88VWmHanfWpiwRGjy5PgHur6+nquuugqj0UhRUdGg6wgpikJtbW2inlqkAHNpVfxbZl18I1Z/zX5Crc0yaFqMScGGGqKubq1smz7vgh6vKAqcIxkB0NscZF52NcGGGry7tvYbTK0Sbj9OuP04Xjahz8hCZ3OAqhLuaNW61U48G9Zpc7BNn3vBXyyMuQU4r/ow/sN78B3Ypd030tVO3zuvYCqtJBYMoJgs6DOdOOYuvqD7A4Ram+jb9AaBxlqiru749iM6Hej1KEaz/AwYYyKdbcR8XmKhAMTinwfFaMKQU6Bdo6oq0d4eTKWVmMsnJivUhEnoYOn+QqEQDQ0Nibq9SHGKTkfGxcu1RCjc1oy/Zr/8EBRjkuf9jdqxzp4xIusEnaQoCpbKyZiKKwjUHSFQe4CYb+A+TlF3r7ZCe386RyYZFy3FmFtw2rnzfn69Htv0eZjLJuDds41Qy6mf26HmehSzhVgkHP9l6My9oP0Cg03H6N34OsGGWlBVrFNmgU5P6HgDwfoa+nQ6cq67ZdQWqxTnFu5oJerzEu1zaXXGgpIBrYEnJxDoM52Yi1P/Z3jCEqFYLHbui0Racy6/nq7nn4x/643FcG97h6xLV551M0ohksG7f4d2fHKlZTUWJerpiw901ulQ9IZ414DFlpAZkDqTGduUWVgnzyTS1U6wuY7Q8YbTkiIAxWzBUjUZ27S55+wGO1/6jCwyL72KUPtxvLu3Eu3tAUANBggeOwSxKDqTCX2mE2N23jnudioJCtQdQWcyYSqp0NY7MpWUE6yvJdzVjr/2II65ixLyGsTwRH1eIu5eYn6vtvYbgKmobOB1vT3os7Ixl1aNi5/fQ34FgUAAi2XwDQiFGIwxrxBL9XT8B3cB8cGZweP1suWGGFNUVR0wiNhYUEy4o5VwV1s88bFlQCxGLOQj0t2JGg6hszswZGWjz8ga9rg3RVEw5hVizCuEuYuJevoId3fEv0BEo+izsjHmF4/Y8hOmghKMV32YQM0BvHveA1WFWIzgscMoRhN6mwPnVR8+4wa0AMHm+nhLUN2R+Ia0JRUDxkopig5DbgHh9uME6g5jn7VAVpsfA8IdLcR8HhSVActG9B8fpKoxIn0uLJWTxs3P7iF/8vLy8li5ciWrVq3ixhtvpLDw3BsSCpF56VVaIhTpbMO7f/u4+cckxodgc/2AVhjFaCLi6sJSPR2jMxdjURmcWKU50t1B1OMm4nYRcXUTam3C4MxN6FRjvSPzvLuOVDVGzOtFjYRQoxFQdCh6fTw5O5E4qdEoxKKokQhqLL7VkaLo4pu0GowoRlN84PekGeizsnFveQs1FAQgcHgfOosVvSOTzMuvGTQZC7U00rfxNYL1R1C0lqDTB4zrHZmEW5uIuLoINh2TgdNjQLizlajPQzTg0+r0ztwBC37GPG4UgwGDM2fcrCs05ERoyZIlrF27lhdeeAGdTsfChQu56aabuPHGG5kzZ04iYxTjSOalV9Px519p3za82zfhXHb9OacKCzFavLs2a8d6RxYRVzeWykk45l2Cbfq8Ab/UVVUl2uci2HSMYEO8qye+uOEh9I5MDDkF6K0jv76aGgkT7mwn0tuNYjCiM5tBb4gvYhqNQDSKemIgtKLTo+j18Rlt+njrlRqNEAsG4glUOBxfD0nRxbvrZi7Et3fbiYHcKr4929CZrRhy8rHPWjggjlBLI70b1hGor0ExGDGVVJ5x9rCiKBicuUR6ugjUHpJEaAwId7QS83qI9vVodR+cLRbu6cSQnYe5fOK4WRR3yInQunXrcLvdrF27ltWrV7N27Vruv/9+HnjgASoqKli1ahWrVq1ixYoVGI2J6cMWqc/gyMA6bQ7eHZsA4oND645gnzE/yZEJEefdtVU7VswWTCUVWKfMOi0JghO/zLOyMWRlY5sxn1BLI/6affE9yro7CTbWojsxyFifma0lHomiqiqRnk7CHS3oM51YKidhcOZgcOaimMwQixILhVAj4XhCpKrojGYUkwnFZEZnNMXvE43G9yHz+4j5vUSDAdRwCDUUJNLXg6l8AsH6GohGQVXxbHsXncmMcmJcE0CgoRb3ljcJ1Nei6HSYSs+cBJ2kd+YSrt1PqL2ZiKsbg1PWoEuWqN9HpLeHqNdDpKdLq+8/PigWDBDze7GUT8Q6aUYywhwRw+qUzcjI4BOf+ASf+MQniMVibNiwgb///e+sWbOGxx57jJ///Oc4HA6uu+46Vq1axQ033CCLLQqcS6/TEqFoXw/ePe8N+ktGiNGmxmIDVpPWZ2Vjyi8i46LLz/n5VBQFc0kF5pIKwj2d+I/sI9hQS6S3O75dR1sz+gwnhpw89NbzXz36TGIBH6GWRlQVzBWTMJdUYJsxH1NR6bDHKamRCFGvO75qdt1hQscb0Vts+A7uik+pVmP0bXqdWCRMtPcaYqEgwcZjBBpq+o0JOndrgc5oRG/PINLnIthcJ4lQEsXHB52YNq/2nzafr10T6e7A4MzFVFY1rmb6JWx0mk6nY+nSpSxdupQf//jHHD58mBdeeIEXXniB5557jqeffhq9Xs+SJUu0LrSpU6cm6ulFCnFcdDn6zGyt+dW3fwfhjlZZaVYkXailccAgUXNZFZaqyRc8kNeYnYdx0XJicxYRqD9C4NgRwt0d8fEwjUfRGc0Y84vQ2TMu+AtALBgg3NVO1N2LMb8QU0Ep9tkXYZk4LWFdFYrBoLV0WSqqCXe0orc7UIwmvLu3agOoPVveIuLqwpCdT8zbhyHDibGoDEVRUKMRQscbCLUfJ9LThRoKYCwowTZ1Dvp++5fpM7KI9LkItTZhn7kgIfGLCxefNu8ZMD7OWFiqJbRqNEKkrwfLxGlYJ81MVpgjYsSG6U+ZMoV//ud/5p//+Z/p6elhzZo1rF69mldffZV3332Xf/3Xf+UHP/gB//Iv/zJSIYgxSmcyY5+1kL6NrwEQbDhKoP6IJEIi6QLHDmnHitmCMacAc9WUId9PZ7FimzoH65TZRLo78NceiLcS9XQSamlEMRox5haic2SeMSFST2xWHPW6ibp7ifl9GLJzsVRPwzphCva5l4z4OCRjfhHOlTfHEzeDEc+Ojdpie4FDe9A5MrBOiG9GG6yvIdzZSuh4Q78FIuOC9TUE62uxTpmJbdZFKIqCzp5JrLUpvr9VwI/OYh3R1yIGF+5oJertI9JvIdH+3WKRni70tgxM+cUYz7D5b6oalfmK2dnZfOYzn+Ezn/kM4XCYt956i9WrV5OdnT0aTy/GoKwrV2mJUMzvxbv7PTLmL0nYmihCDEX/zUv1GVmYikoTsvu6oigYcwsw5hYQnX1xfIuMmv2Eu9oJtR2Htmb0jsz47BydEp+eHwwSC/rj6xYBersDgzMHfdkEzGVV2KbOGdVfSHqrDefy6+Nx2hy4N72ubTIb87jjU+3Pi4r/8F70mdlYKiehMxpRTGaiXg+htuYLWrBRJEZ8fFA3UVfPgN3mT+4HeXIsmqmkEuukGeNuGMOoL9xgNBq5+uqrufrqq0f7qcUYYp8xD2N+UXy7AMB/eA/B5jqZOSKSyl+7Xzs25hVhmTD01qAz0VttOOYuxjp5Fv4j+wgcPRj/JeTzEHG7QFVRdDoUkxlDVg66ojJ0FhvG3HxMxeWYSioTkpwNhaI3kHHxMgxZ2eisdtxb3xywFclp9AbMZVUY84pAVfEd2q1tEuvduRljfhF6mwO9I5Oop5dQa5MkQkkQbj9OzOclGjjVLWzIztNa56KePtDpMGTnjostNT5oyInQULfP0Ol02O12aQ1Kc4pOj2Ph5fS8/AwQX8rfL1NoRRKpqkqwrkYrm4rKMZdOGLHn09vsOOYuwjZjHqHj9US6O4n0xpMKxWBEZ7VrM8AMzpwxs+CgoijxcT6OLAzObHz7dhBqazrRihVA0esx5hZgyCs6bQNYY2EJrnXPx2exRcJ4tr1L5tJr0dszCbU0Em5rRo3Fxs207FQR7mgh6nUT87i1OuOAbrFODDn5WKqmjIuVpD9oyK+oqqpqWM1jdrudFStW8O1vf5vFiy98Iz+R+pzXfISeV//vxOJuYby7tpC5eDl6e+rvZixST7j9OLF+C8lZJozOD32d0RT/ApBiXwLMpZU4l10H0SiK0YjBmXvOPc/0Ngf2eZfg2fYOEP8FHGppxFRcjhoNE3H3EunpHNbeaeLCqKpKuL0l3irpdmn1J7vFYsEAMd+JKfPV05IU5cgactpdUVEx5D8FBQV4vV7WrFnDsmXLePfddy/ouX0+H88//zx33XUXU6dOxWKxYLfbmTt3Lt/97nfxeE7fm+ekJ554gkWLFuFwOMjJyeGGG25g48aNZ7xejBxzcQWWqlM//OMDSY8mMSKRzvqPD9JZbZhLK5MYTWow5hXhmL8Ec9kEwp1tRL3ucz7GXFGNsfDUIn2Bmn0oioLenkHU6ybU2jSSIYsPiHndRPpc8b3FVBUAxWLVps1HerowOHMwlVSMqynz/Q35684Hd5u/UF6vl8cee4xvfetbfPe73+XVV18978f+6U9/4vOf/zwA06dP58Mf/jB9fX1s3LiRBx98kD//+c+sX7+egoKB3yq+9rWv8cgjj2C1WrnmmmsIBAKsW7eOV199lWeeeYabb755WK9JXBhFUci8/BptX6dIZyu+w3uxTZ+b5MhEOhowUPo8WjdEnKV6OuHujvhaQk11mCuqzzqLTVEUbNPm0tvWDMRnK0VcXejsGcS8bsKdbaMVugBC7ceJ+txEfaeSWHNxfFsUNRYj0tuNpWoylonjszUIhtEiNFx2u51vfvObLF68mPfeO9/ZBnFGo5EvfOEL7N+/n/379/O3v/2Nl19+mUOHDjF//nwOHjzI1772tQGPee2113jkkUfIzc1l165dPP/887z88su8/fbb6PV67rjjDlwuV+JeoDgvzhXXo+u3uJx3z3tE3L1JjEikq/5T5405BRgkEToviqKQsfAyrNXTMRWVEmysJdpvLabBGHILMPTbwd5/ZB96myM+YLy7XdsDTYy8cHsL0b5eon0urc5UUgHEB0krJlO8RegDO9CPJwlLhP785z+f13WqqnLXXXdp5cmTJ9Pbe2G/+D772c/yq1/9iunTpw+oLy4u5uc//zkAzz33HKHQqTUsHn74YQDuv/9+Jk8+1R2zZMkS7rnnHlwuF7/97W8vKA4xfHpbBrZ+22sE62sINR1LYkQiHamqOiARMpVWJmzT1HSg6A1kXnoVlqopmApKCTYcJdhcT9TvRT3R3TLgekXBOvnUonzBxqPxzWBV4ls8dHeMZvhpS1XV+DitE4PdIT5Q35gfX9Mt4urC4MzFUjl5XA9gT9gru/3223nllVfOed0dd9zBE088oZWffPJJYif+AhJh7tx4t0owGKSrK75fit/v54033gDglltuOe0xJ+tWr16dsDjE+XNe8SHtOOp24d23PYnRiHQU6e4g2q8lsv8vaXF+dEYTWUuvxT7n4vhaMwYDoaY6/Id2E2yuIxbwD7jeVFp1qjVYVQk21KC3x1uFwh3SPTYaor09RDx9RHo6tTpjURmKXk8sHCLm88RXFx+BZSTGkoQlQnl5eXzsYx9j06ZNg55XVZXbbruNJ598kksvvTRRT3uao0fjg22NRqO2r9mhQ4cIBoPk5+dTVnZ6896CBfFl3Xfv3j1icYkzcyy4FH3WqT2GvPu2E+nXTCvESOs/PkixWLFUjL+1UkaDzmwhc/EKsq/5CBkLLsM2cz6W6ukoJjOB+iMEm+u0FiJFpxvwCzZYX4NitRPzeQl3tibrJaSVUGsj0T4XUXefVmc+2S3m6kaf4cRUWDZuB0mflLBE6JVXXsFkMrFq1Sr27t074JyqqnzmM5/hD3/4A5dffjlr165N1NOe5pFHHgHguuuuw2yON22fXPNosCQI4uOVnE4nPT09uN3nnvUgEkvR6XHMPbWEQrChlkCjzB4To6d/t5ghOw9jbmESo0l9xpx8Mi+9itwP/yM5134snhRNmY0aDhNubdKSIXPFqcUTo54+CAWI+uIDptUE9hSIwYVamwm1NEA0Eq/Q6TAWlcVXkj7ZLTYhtZZ1GIqEJUKzZs3ixRdfJBgMcu2113LsWHycRzQa5R/+4R/405/+xLJly1i7di0OhyNRTzvASy+9xG9/+1uMRiP/+Z//qdWfnE5vs515JoPdHm+iPZ9EaObMmYP+qa2tHeYrSF/Oq2/WjmNeN7697ycvGJF2Akf7DZTOzkMna1klxMmtRTIvuYLMy67GXFEd39W+ux2IbxvSf5uQYEt8rErU5yHi6kpW2GkhFg4R7mwl3HZcqzMWlKIzmuKrfysKhuycEV1UdKxI6OinJUuW8Mwzz9DZ2ck111xDY2Mjn/70p/nrX//KihUreOmll7SEI9EOHjzIP/7jP6KqKj/60Y+0sUIiNdgmz9QG6AF4971P1Hf2mSdCJIr/xBIOQHz13HG2l9JYYC6tJGPBEswV1YQ727TFK/u3CoWajqGz2OLdYx3SPTaSwu3H4wtY9huGYC6rAk4Mks7OxVIxaVyuJP1BCX+F1157LU8++SS33nor06ZNw+/3c8UVV7BmzRqs1pHZVbi5uZnrrruOnp4e7rvvPu69994B50+2QPl8vsEeDsTXNQLIyDj3N8F9+/YNWj9zpgywHA7HxUvpeelvAIQajhI8Xo9t0owkRyXGu4irm2jvqf2yLNXTz3K1GA7LpJnxhRf7XITajmOpnISptAp2boZoBPXEAN2o0Ui4owWmzk52yONWqLWJUHO/bjFFF1/hOxIm6nFjKiof94OkTxqR+XCf/OQn+fnPf47f72flypW89NJLI5YEdXd3c80111BfX88dd9zBj3/849OuqaiID/5qahp8xVKv14vL5SI7O/u8EiExMnKu/Zh2HAv48G6XFb/FyAvU9RsobTLLitIjSFEU7LMvxphfRCwUJOrpQ2c0DnjPw13tRL0ewh2tsp7QCFFV9dT4oBOMhSXoTGYirm70jgyMBcUYnLlJjHL0DLlFaOLEc8+qMBgM7N+//7T1fhRFSch4Go/Hw/XXX8/+/fv56Ec/ym9+85tBm7SnTp2K2Wymo6OD5uZmSktLB5zfvj0+XXvOnDnDjkkMnamwFFPZBG0dIc/ureR+5DZ0RlOSIxPjWf8ZYwZnLvos2RB6JOkdmVgnzSDS00Wo/TgWewbm8okEG+K/E8IdLRhy8uLjhHq6ZIXvERDtcxHp7R4wDstcWqUNkjYVlWOZODWJEY6uEd9i4/jx4+e+aAiCwSA33XQTW7du5dprr+XPf/4zer1+0GutVitXXnkla9eu5emnnz5t1elnnonvgL5q1aoRiVWcv6wlV9Hx9P8C8UXWQscbsFROOsejhBg6f81+7diQk48hw5m8YNKEbfo8AnVH4q0/fS6MBSUoZgtqMACxGGogQNTrJtzeIonQCAg2HSPYeHRgt1hJRXyQNGBw5mApT58lJIbcNRaLxYb1Zzii0Sif/vSneeONN1i6dCnPPfccJtPZWw3uu+8+AL73ve9x5MgRrX7Tpk386le/wul0DljxWiSHc+VNoMQ/lmowQN/W9UmOSIx3/VuETMUVaTE4NNl0ZgvWKbMw5hfF1wxSFMylVdr5sKuLmNdNqH1kvkinM1VVCTYeI9hUp9UZC4rRmcyEuzsw5ORjmTAFxWBMXpCjLCX/xT/22GP83//9HxBfyPGf/umfBr3uxz/+MXl58f1sVq5cyb333ssjjzzCvHnzuPrqqwmFQqxbtw5VVXn88cdxOp2j9RLEGRiysrFUTyNw4lu6d8cm1I/djqIbvLVPiOGIevq0qdwAlippfRwt1kkz8B/aQ7ijhai7F3P5RG0D5mhPJ5He7vh6QtEIij4lf1WNSdHeHsLdbUR6Tm1jYi6bQCwYIObzYi6bgDXNJqmk5Kerp6dHOz6ZEA3moYce0hIhgJ/97GfMmzePxx57jHXr1mEymVi5ciUPPPDAiK52LS5M1vLrtUQo2HiUUGuzttqpEIkUqDvVOqwYjJjKxv+aKWOFzmTGOnkm4a52wp2tmKumoLM5iPni675FenuI+TyEu9oxFZQkOdrxI9h4lGDDMYic7BZTMJVUEOnuwODMwVxaOe5Xkv6gIXeNdXd3n3U6+vnw+Xx0d3ef+8IPeOihh1BV9Zx/qqqqTnvs7bffzrZt2/B6vfT09LB27VpJgsaYrMuvhRPfANVwiL5NbyQ5IjFe+Q7u0o71WdkYnTlnuVokmnXyTAy5+RCLEfO6MZefSkSjfa74KtPtLUmMcHxRVZVg0zFCzXVanbGgBEVvINLbjSE7L+1ag2AYiVB+fj5f+cpXhvXkX/rSlygokIFwYiC91YZt2qkZfJ733x10B2shhst3YKd2bMwtRJ8pM8ZGk85swTZpBsb8YsLtLZjKTg3QjXndhDvbCbU0JjHC8SXS00Wou4Nw/26x0ioivd3orHaMeYUYC0vPcofxaciJ0MlWl+GSX3BiMFlX3KgdBxtqCfcbxyFEIqiqOmCPMWNhCXrZWmPUWafOOTEzTEVRFPSZTu1cuLWJcFcbkd6eMz5enD9/zT5CDTUDusWMxeVEujsw5uRjnTQjLVdVH9YYoXfffZc777xzWI8XYjCZFy2l1WBEjYRRI2H63nmFvJtvS3ZYYhwJtx8n1m8bF0vlZBTdiKwxK85CZzJjnTqbcE8n4bZmzOUT8e2Lr+0WdfcSdnUTbKjFMPuiJEea2qJeN8H6GgINp9bwMxaUQDgEKPENVqvG/wargxlWIlRTU0NNTc2wAkjH7FOcm85ixT53EZ73NwDg3rpeEiGRUP27xXRWG6bisuQFk+ask2cSqNlPuLMVve3UptyxgI9QaxOBhlpssxbK74th8B/eS6ijhWivS6uzVE4i3NOJITcfy8SpaTVlvr8hJ0JvvvlmIuMQ4jRZy2/QEqFAfS1hV7cMZhUJ492zTTvWO3MxZOed5WoxknRGE7bpcwl3dxJsPoYhJ59Id3wcS6SjlYiri3BHK6aC4nPcSQwmFvDhP3qQ4LHDcGLbEsVgxJBbQLDxKJbyiVjTeI+9ISdCy5cvT2QcQpzGMe+SU6vNRiP0vvUieTd/JtlhiXHCX3tAOzbmFWFwSiKUTJbq6fhrDxLp6cDgzNUSoWhvN+GeDoINtZIIDYGqqnh2v0eku2PAlhqmsglE+1wYsnIwpeGU+f6kQ1yMWTqTmYyFl2vl3g3rkhiNGE+iAT/htlOrFpuKyjBkOZMXkEDR6bHPWYSxsBTFaIYTi6iqkTDBuhoCdYeJevqSHGXqCdQewH9kH4GmOqLuXq3eXD4hPmU+Jz+tW4NAEiExxmVedrV2HDpeT8R14etOCfFBvn3vg3piqx+dDnP5RFm9eAwwFZdjKZ+IsaB4wB5jke52Qu3H8ezemsToUk+o/TieHRsJNh5FDZxa909nz0AxGNFZbBhz8tNyynx/kgiJMc0xdxE6qz1eiEbpef3vyQ1IjAue7Zu0Y32GE2NufhKjEScpioJ97iJMBcXo+40HjPm8BBuPEairIdTalMQIU4OqqvgO78X11ksEGo6iGIwD1mOyVE0m0tMV31esenraD0KXREiMaYrBSMaiU+PR+t6V7jExfP1njBly82Wg9BhiyMzGOmUWlqrJ6PqNW4n2uQi3NePZuRk1Ek5ihGObGovifu9t3FvXE6g9iKLTo4ZDqKFg/AK9AWNxOWosGt9lPk2nzPcniZAY8zIvvUo7DrU1Ee7pTGI0ItVF+lwDvh2biitkoPQYY5sxH1NByYA9xiJdbYQ6Wwk21+HeJqvND0aNhOl7dx2+fTsI1B3GkJOHsaRiwMQAS9VkYp6++LpBFdXoTOYkRjw2SCIkxjz7zIWnvhnGYnS/9LfkBiRSWt/mN/qND9JjKiqXgdJjjM5owj7rIqzT5qBYrFp9pLON4PEG/DX7CdTsS2KEY48ai9H77jp8R/YSbK7DXFqFMSefcEsjsX6DzC0TphJ192LIzsUycVoSIx47JBESY55iMJC5eIVWdr/3dvKCESnPs32jdqzPcmLMzpWB0mOQuXISpqIyLFVTtbpobzfEYgQbj+HesVn2IesnULufQEMtoZZGzBUT0TsyUaMRvP0GmJtKKlEjIfSOTEz5RdIlfMKIJ0I1NTVs3ryZxkb5wIqhy1xyqnss3NZMqN/UZyHOVywcItB//aBc+WUwVik6HfaZC7FOnjFg/7HA0YOgUwg21NK76XXC3R1nvkmaiHr68OzZRuh4PcaCEvQnJpj4Du4m5vPEL1IUrNPnxQdJZ+fJIOl+hpwItbW18be//Y0NGzYMen7Dhg1MmzaNqVOnctlll1FVVcUll1zC/v37hxysSF+2GfNO7QyuqnS/9NfkBiRSknfP+wPWUjGVVsgifWOYqaQi3ipUPR1O/NJWgwECtQdRY1GC9bX0vvMKwaZjaTtmSFVV3NveJdTShGIwYHDmAvF92vyH92jXWSfNRGeIt3wasnKwlE9MSrxj0ZAToT/84Q98+tOfHjSxqamp4brrruPIkSOoqkpOTnwa5NatW7nqqqvo7pa1YMSFUXR6Mi+5Qit7dmw6y9VCDM695Q3tWDGZMeYWxDeeFGOSoijYZy3EXF6NqaRCq4/29RBsqCXqc+Ov2Y9r/Vp616+lb/Ob9G15C3/tgbRZfDHc2kSwuY5ITwem4goURSEWCsbHwsXiY+F0Vhu2GfOI9HRiyMnDUjU5bfcVG8yQE6H169djNBr5xCc+cdq5hx56CK/XS1lZGbt27aKjo4OOjg5uvPFG2tvbeeyxx4YVtEhPA7rHOloINNclLxiRcqI+D74Du7SywZmLMTsfndmSxKjEuZgKS7GUT8AycTqm4nKtPtLZhv/QHiK93fhrD+DZtYW+LW/Rt+UtXG+9RNeav+B+7+1xPdVeVVW8B3YS7mzD4MxFZzLHp89vfoNon0u7zj53MWosRtTrwZCVG29hE5ohJ0IHDhxg/vz5ZGVlDagPh8M8//zzKIrCD3/4Q2bPng1ATk4Ov//977HZbLz00kvDi1qkJevU2afGc6gqPWufSW5AIqX49m0n3N5vW43iCozSLZYS7HMWYcovQp9bgKn4VMuQGgoSOLyXwJF9BBtqUSNhdAYjke4O/DX78Ox+j57XXiDS25PE6EdOuK2ZUEsjUbcLY24BaiyGe+vbhDtatWusU+dgLq0i0t2JPisbc2klhn5jrsQwEqH29nYmTjy9j3H79u34fD4sFgs33XTTgHPZ2dksWrSIw4cPD/VpRRpTdDoyL7lSK3t3bU5iNCKVqLEorvVrtbJiNGEsqRiwTo0YuwxZ2Viqp2EqKMFQWIJ93iXQb6afGg4Raq7D897bePduw5Cdi7lsAqHWJnyH99Dz+t8J1Nck8RUknqqqePfviLcGZeWCTod7y5uE+rWUm8snYpu5ADUWJdLbhfHEStJioCEnQsFgEK/Xe1r99u3bAZg3bx4Wy+lNzkVFRYM+TojzkXnpqUQo3NFKoPFYEqMRqSLYVIf/yF6trM/KwWB3YMgrTGJU4kLYZy6Ib71hsqAYTTivWoW5ctKAhAgg0t2B6401hFoasUyYghoK4q/ZT++GdXh2bEI9MW4m1YXbjmutQYbcfNxb1xM63qCdNxaW4lh4OYqiEHF1x/cVyy0Y0L0o4oacCBUVFbFnz57T6t966y0UReGSSy4Z9HEej0cbPC3EhbJMmomh32aMPa8+l8RoRCpQVRX3lvVE+3WPGIvKMOYUoDOakhiZuBA6i42sy6/BXFkNqkqkpwvrlNlkX/cxHPMvHbgMQiyKd9cWvLvfw1Q+Ab09g8CxQ3h2b6VvwzpiJ7ebSFGqquI7sINwVxuGrBz8B3YRaq7XzpuKy8lcchWKXo+qxoh0tWPMK8Q6aYZMmR/EkBOhpUuXUldXx69//Wut7vDhw7zwwgsAXH/99YM+bvfu3ZSWpvdOt2LoFEUh89KVWtmzffDlG4Q4KfyBXct1NjumgmKMhdItlmqMeUVkXXY1lspJ6Kw2Il3tBI4dJhYOYamehmPRMhTzqZWog3WH8W57F2NeEaaicoKNR/Ed2k3v+rUpnQyFO1oItjQS7XURcfcNaO00FpaScckVKHo9AFFXN4rRiDG3AMuEqWe6ZVobciJ03333odPp+OIXv8iyZcv42Mc+xqJFiwgGg8yYMYOVK1ee9pg9e/ZQX1/PRRddNKygRXrLXHKqeyzS1U6g4WgSoxFjnW//ToJ1p8Yl6jOcGDKdmMsmJDEqMVSmojKcV32YrCVXYZ9zMbZpc7FMmII+w0nM5yVj8YoBrcbBxqO4t7yF3pGBuWISobZmAnVH6Nv4Omo0ksRXMnS+/fGZYorBiG/ve1q9PiuHjMVXoOjiSZCqxgh3tmHMK8Y2bQ6KQVZQH8yQ35X58+fz61//mi996Uu8++67Wn1BQQF//OMfB33MyWnz11577VCfVggsE6ZizC/SZka43lxN0WfvTXJUYiyKuLpwb3sbNRyKVyg6jMXlGHMKMGRlJzc4MWTGnHyMOflAfLXwqKcP355t+C02gk1Hsc1cgP/ALsIdLQCEjtfTt/F1MpdciaWimkB9LYrBgM5qI2PR8pTqLgq1NhNsaSDi6ibUfExbK0ixWMm8dCU646n1gaKubhSDEWNeAZYJsq/YmQxri4077riDmpoafv3rX/P973+fJ598ksOHDzNnzpxBr583bx4//elPufrqq4fztCLNKYpCxuJ+iytueyeJ0YixzHtwN75Dp8YyGnLyMOUVYq6QVXXHC53RhDE7j8zLr8ExbzGWyslEujuxzb4YY1GZdl24rZm+DetQDEbMpVUEm+vwH96L/+Cus9x9bFEjETw7NhJubSbqdg0Y95ax8HL0Nnu/a8OE2lswFpZgnSqtQWcz7HempKSEz33uc+d17Re/+MXhPp0QAGQsWk73mj8DEO5sO9H8KzOAxClRrxvvzk0Ddt42ZOWgz8jCLNsLjDuKTod9ziLUSBg1GiHYeBT7vMX49hi0KeXhjlZ6332VzEtXYiooIdh4FMVoRJ+Vg7nfytXn6+Q4I53JnMiXcka+g7sItTYR7m4n1HRqxqxlwhRM/ZI+gFBrE4asbMwllVgnyZT5s5Hd50VKsk6aMWDvsZ431yQ3IDHm+A/vxbd/h1bWZ+VgzC/ClF+C3p6RxMjESFEUBcf8S7FWx1ehDjU3YJ93CeaKau2aSFc7fe+8gs6egd6RSbCpDveWtwh3tp7lzqeosSj+owdxrV9L19//QNcLf8S97R0i/VZyHgnhrnZ8B3YRPN5AuL0FTuytprM5sM1eNODaSF8PMb8PU1EpGRcv1cYMicENORG68847+d3vfjfouRdeeIGdO3cOeu7BBx9k4cKFQ31aIYD4t7+Mi5dqZffmN5MYjRhr1EgYz/7t8V8YJxiy8+Lf/KVbbFyLL7x6BZaKSRhz8gg112GfvwTLhCnaNRFXF33vvoo+Jx9QCRw7jGv9WkItjWe8rxqLEWw6Rs+r/0fvO6/i2bER38Hd+A7txv3eO/S88iyh1qYReU3B4w241r9EsOkoUXcv0d5T+3U65i8ZOC7I5yHU0oiptBLb9PnaJqzizIbcNfbEE08A8YTog26++WZuv/32QROlhoaGMyZJQlyIjEUrcL0eX64h1NJIxN2LISPrHI8S6SDYXI9v73atrLPa0GVkYXDmDGgdEOOTYjCSednVxEIBYjUHCDXXYZt3Cej0BGoPABDp6cS98TUyLl1JuK2ZwLHDEItiLq/GMnEqepsjfl2fi0h3B4H6I0R6ugi3txD1eeIb9uYXo0YihDtaiHrdoDeQfeWNw0o+1EiYiKuLWChEzOsm1NZM8HgDwYZaYpHwgJWjTaVVA7rEon4vwcZjmEoqsU6chm3GvCHHkU5k9JRIWfaZ89FZbMQCPohF6X37ZXI/9MlkhyXGAH/tQYINtVrZmF+MKScfS8WkURvPIZJLb7OTddnVxIIBAkcPEW5pxDYn3oWkJUPdHfFk6LKVRDra8B3eS6i9BX/NftDrUVDiyVTAT7S3h1goiCE7F2tJhbZODyYzuoqJBOtrCTYepffdV3FeuUpLpM5XuKsdf81+gs31xLxu1GgENRIm6vUQ9fSis2UQbjqmzYBUDEbsc091iUX6XPGWoOIyrBOnxRdUlC6x8yKJkEhZisGIfcES3BtfB6Bv42uSCAliAR/urevh5BoxegM6mwODM1cGjaYZgzOXrEtXokYiBI4dItLZin3uYlBVAkcPAieSoXfXkbF4BcacfMI9nfF9yWIxUGMoJjM6ixVDTh76DCeK7vQRJYqiw1w+gUDdEQINtfS+/TLOFTegs9jOGaOqqvgP7caz+z3CHS1EertRVFCMRtAb0NscGCuqCR1vJNyv68026yL0VjuqqhLpaifc1Y65fAKWqilkXnqVzBK7APJOiZSWecmVWiIUbKglFvCjs1jP8SgxngUaagestGsqLMXgzNX+L9KLqaiMzIuXxWeSHTuMojfEN21VVQLHDgHxbjLXa3/HPmcR5opqlCHMIFP0BswV1QTrjhCoO0Lv26+QtezasyZDsVAQ99b1+I8dItRUh85iw1w2AZ3FNmBto0hP54DV0Y2FpVgmTkVVVcJtx4l6+rBMmIJ9xnzscxdJS9AFkkRIpDTHnEUoBmN8ymw4RO+GdWRf9eFkhyWSyLN9IzGvWyvrHZkYcvKwSGtQ2rJMmEIs4IdYjEDdYRS9Hvv8JaAoWsuQGgnj2b4B795tmMsnYiwoxphbiM58+ubhZ6IzmjBXVBOoPwLEBy7bZ1+EZcKU05KTcHcn7s1vEGg6Rrj9OKaCEvTO3NMWdwx3tdO3YZ3WwqmYzGQsvByIj42MBXxYqiaTcdHlWCfPHPJ7lM4kERIpTWexYpt9Ed4dmwDok0QorUU9fXj3bNPKemcOOrsDY3Ye5tLKJEYmks06bQ6xUAA1FiVYX4sJsM+7BENuAd4dm1AjYQDUUJBA7QFtHJHOZkef4YwPuD85vkw98R9FhyEnH1NRmTZmSGe2YKmcTLC5jkhfD5E+F96972MurUSfkYWi0xNsaSDU0kiotZmY1x3fO+0DLUexYAB/zX78R/ad6uZVFDIuWorOaiPU3nIiCZpC5uIVWKomj8K7OD5JIiRSXuYlV2qJkP/oQWLhkOwqnqaCzfUEG0/tPWfIysGQnYelajKKXn7cpTNFUbQFFxUUAg01GGMxLBXVGHML8R3YQbDpGESjAx4X83mJ+bxnv7fRhLlqMrbp89EZjfFkaMJUIj2dBBtqCekNBJuOoRhNKDpdfAB0nwu9IwPLxKnaZzPq8xBsqCXcfpxwVwfE+sWi05Gx+ApMxeVEeruJuDqxTJhK5qLlkgQNk6KqJ1ZlukA6nW5Y+7NEP/BhGw9mzow3S+7bty/JkaSXqKePw1+4Udtzp+QrD5F12emb/orxr+0PP9dWHEfRYZ0yG/vM+eTc8An0jszkBifGBFVV8e7agnfv+wQbatHZ7JiKylH0emKhIMGmY4TbWwh3tqIGAxd0b53NjmPBZZgKSwc8X8znIerpQ41EIBZFMVsxOHO0FqZITye+AzsJtTRxorlpAMVoImPRckxFZUT9PoINNZgrqnHMW4JjzsXDej/GkmT9Dh3WV6Qh5lAptcGdGPv0jkysU2Zrewb1rn9JEqE0pEYiePsNKDVk52HIzcdUXC5JkNAoioJ97mJ0ZiuKyUzoeAOBowcx5hWiz8rBOnEa1onT4gmM30vE1U3M6yYW9KOGTm7eG/8dFgsG4osonui6ivm89L37KubKSdjnLEJnMqMoCnp7xqCrmauRMN592wnUHGDQBMhkxjp5JpaJ09CZzKiRCKGmY5gK41Pk7bNkceJEGHIiFDvx7VuIsSDr8mu0RMh3aDexYOCCBjmK1BdoOkaw32JzOnsGxuw8LBNl120xkKIo2KbPxZhfhHvrekJtzYQ7Wwl1tGJw5mBw5qIzmdHbHOdcD0iNhPGd3Lz1RONAsL6GUGsT1imzsE6chmIwnva44PEGvDs3E/MP7HbT2TMwV1RjKijBkJOvTddXVZVgcx06RybmsioyFi0fdCq/uHDSaS7GhcxLrqD18Z9CNIIaDND7zitkr7wp2WGJUeTe8tap8R16A4acPAzOHEzF5UmNS4xdxrxCsq/5CP6jh/Af2Uu4q4OIq5PA0UPo7Q4MuQXnTIQUgxH7jPmYSyvxvL+BSE8nAGowgG/PNvwHd2MqLseYX4xiNBHzewk21RHpahtwH50jE8ecRRiLyk7rNVFVNb79hxrDXFZF5qUrZWHQBJJESIwL+hM/RDw7NgLgeuslSYTSiKqqePe8p5UNWTnxmWLlE+VbszgrxWDENmUW1upp8W6yY4cJtjTGBzo316EzWTAVlp5zfTJDVg5ZKz6Ev2ZffLPfE0m5Gg4RbKgdsNL5wAB0WKfOxjZtzqAD+lVVJdzaNGCGmCEre9ivW5wy5J8QV155JT/84Q8HPdfQ0EB3d/eg54QYKVlXfEg7Dhw7RMTdm8RoxGiKuLoHzBbTWa3os7IxV0xKYlQilSh6A+byiWQtu46c6z9OxsVLsU2Zjd6eQaD+CMHjDdoU+zPeQ6fDNmU2Odd9HOuUWXCOmYqGvEKcK2/CPnPB4EnQib3Foj4PlsrJZF6yAnNp1XBephjEkFuE3nrrLaqqqgY9N2HCBG6//XZ++9vfDvX2Qlwwx7xLTu09Fo3gev3v5N18W7LDEqPAs2PjqRk+ioIxvxhjdj6GbFlJWlw4Q6aTjIWXY5s6B+/ebQSOHSHUcRx/7QGMuQUYsvNP7TU2CJ3Fin32xVinzSPc3kyopZGouxc1GgVFwVRQjKm0CkN23qCTh9RYjEhvN+GOVvQZWVirp+G4aCmWSpkmPxJGpGtMVdUhzygTYqh0JjMZi5bT+/ZaAHrffkUSoTTheX+DdqyzZ2DMLcRSWS0zVMWw6B2ZZF5yJdZJM/Ds3EyopYlQRwvhrn3xWYknBlWfic5oxFxadd6tOGokQqSnk3B3BzqLFXNpFaaScjIuXoYxOy9Br0p8kIwREuNK1oobtEQodLyBYEsT5uKyJEclRpIajcR3Cz9Bb89An5mFuaI6iVGJ8cSYV4TzqpsINtTgO7CLUEcLke5OAscOoRhN6G2O+MrTZguKyXLB49LUSJhwVzsRVxf6k7PG8gpPTJ2fKouBjjB5d8W4Yps2F31WDtHebkCl55VnKLr9a8kOS4wg3+G9J/6+40wFpRjziwddt0WIoVIUBUvlZMzl1QSb6wjWHYl3eXk9RP0eIr09qMEAaiSMzmxFZ3Ogz8xCZ7WfsWVSjYQJd3cQ6elEn5GFZcJUTAUlWKfOxlw2QQb6jxJJhMS4ouh0ZC27ju7VfwLiU6oLP3uvdJGMY+5Nb2jHitmKsagEc9mEJEYkxjNFp8NSPhFL+URiAR+hjlYinW1EXF1E3b1E/T5iPi9Rn5tQcz0A+ows9I7M+HpCikIs4CPq6SPq7osnQFVTMBWVYp+5cNDp82JkSSIkxh1nv0Qo0tOJ7+Bu7NPnJjkqMVK8+7Zrx3pHBnpHFuayquQFJNKGzmLDUj4RyidqdVGfl3BXG6GWRoLNdUT7XEQ9fYQ6WiASQSW+Mave7sBYUIqpqBTb1NmYSiolAUqSpOw1pigKkUhkSI8dy2SvsbGj9hu3at/GMi65krKvfTfJEYmREO7ppOafPgpqfKV767R5ZC27luwrVyU5MiFOTH9vaSTUdpxwRwtqJIwai6G3Z2AsKMFcUo4ht1ASoBPSaq8xmVEmRlrWig/R8cdfAODdvYVYKITOJDvSjzd9G17TkiD0BkwlFbLOihgzFIMRc/lEzP1ajMTYM+SRWLFYbFh/hBhJzqXXghL/eMd8XnrXv5TkiMRIcL//rnasd2RiyHJKt5gQ4oLIkHQxLhmcudjnXKyVe15/IYnRiJEQi4QJHD2olY35RRjzimS2mBDigkgiJMat7Ktv1o6DDTUEWxqTF4xIOO+OTQNWk44vXFeZ3KCEEClHEiExbjnmLUHnyIwXYjG6X/pbcgMSCdW3+dS0eZ3NgSE7D5OMDxJCXCBJhMS4pRgMZC27Xiu7t7xFLDr+ZiumK+++HdqxIScfY04+hoysJEYkhEhFkgiJcS175U3acbSvB/eG15IYjUiUYFMdUVeXVjaXTcAk3WJCiCGQREiMa+aSCqxTZmnl7leeTWI0IlFc77ysHStmC6aCYpk2L4QYEkmExLiXfcMntePA0UMEmuqSF4xICO+OTdqxITsXgzMXfVZ2EiMSQqQqSYTEuJd50VL0J8eOqDG6V/8xuQGJYYl4+gg2HtPKpuIKTKVVsjqvEGJIJBES455iMOC86tRYIfd7bxMLBJIYkRgO96Y3Bq4mXVqJuaQiuUEJIVKWJEIiLWRf8xHQnVppuuulvyY5IjFUfVvf0o4NWdkYnbkYcguSF5AQIqWlbCL0/vvv84Mf/ICPfvSjlJWVoSjKeTWNP/HEEyxatAiHw0FOTg433HADGzduHIWIRTIZc/LJWLRcK/e8+pxs9ZKC1EgE/6E9WtlYVIqpuEK6xYQQQzasTVeT6T//8z/5+9//fkGP+drXvsYjjzyC1WrlmmuuIRAIsG7dOl599VWeeeYZbr755pEJVowJeR+9HffmNwGIurroe/dVnMuuS3JU4kJ4925DDQXjBUXBUl4ts8WEEMOSsi1CS5Ys4YEHHuCFF16gpaUFs9l81utfe+01HnnkEXJzc9m1axfPP/88L7/8Mm+//TZ6vZ477rgDl8s1OsGLpLBUVGObuUArdz3/ZBKjEUPRu2Gddqy3Z2LIycNYUJTEiIQQqS5lE6FvfvObfPe732XVqlUUFZ37B+HDDz8MwP3338/kyZO1+iVLlnDPPffgcrn47W9/O2LxirEh76O3a8eh4w24d24688VizPHu2aYdGwuK491iOn0SIxJCpLqUTYQuhN/v54034vsS3XLLLaedP1m3evXqUY1LjD7bjPmYJ0zRyu1/+HkSoxEXItB4dMBq0qaKiTJbTAgxbGmRCB06dIhgMEh+fj5lZWWnnV+wIN5dsnv37tEOTYwyRVEo+NTdWjnUVIf7/XeTGJE4X33vvKIdKxYrpvwSTEWn/3sWQogLkbKDpS9EQ0MDwKBJEIDdbsfpdNLT04Pb7SYjI+Os95s5c+ag9bW1tVRXVw8vWDHi7HMWYZkwlcCxQwC0//EXZCy8PMlRiXNxbzuVsBrzCjEVlaIYjEmMSAgxHqRFi5DH4wHAZrOd8Rq73Q6A2+0elZhE8iiKQsE/3KOVQ8cb6N30RhIjEucS6mwjdLxeK5tltpgQIkHSokUo0fbt2zdo/ZlaisTYY5t1EZbq6QRqDwDQ8edfkrloGYpe/kmMRb3r12rHitGEqaQcU3F5EiMSQowXadEi5HA4APD5fGe8xuv1ApyzW0yMD4qiUPDpU61C4fbjdL/yXBIjEmdzcv0nAGNeEaaCUnRmSxIjEkKMF2mRCFVUxGeWNDU1DXre6/XicrnIzs6WRCiN2GYuwFI9TSt3v/RXIu7eJEYkBhPu6STYWKuVTRUTMZdWJjEiIcR4khaJ0NSpUzGbzXR0dNDc3Hza+e3btwMwZ86c0Q5NJJGiKOT/wz9p5UhnG90v/iWJEYnB9L37qnasGE2YSyoxyfggIUSCpEUiZLVaufLKKwF4+umnTzv/zDPPALBq1apRjUskn2PmAiwTT7UK9b79MoHm+rM8Qoy2vo2va8eGvEJMhaXorWee+CCEEBciLRIhgPvuuw+A733vexw5ckSr37RpE7/61a9wOp3cddddyQpPJFFB/1ah7g56XnkWNRpJYkTipEhvD4G6w1rZXD4Bc1lV8gISQow7KZsIvfjii1xyySXan1AoBDCg7sUXX9SuX7lyJffeey9dXV3MmzePm2++mRtuuIFly5YRiUR4/PHHcTqdSXo1IpnssxZgrjq12rR763p8NQeSGJE4qW/zG6CqACgGI+ZSSYSEEImVsnOFOzo62LJly2n1/es6OjoGnPvZz37GvHnzeOyxx1i3bh0mk4mVK1fywAMPcOmll454zGLsKvj03TT+9zeA+M70vetfxFY9TRbsS7LefuODTnaL6SzSLSaESBxFVU983RLDdnIdoTOtMyTGtqP/+lmCDfHZSYbsfEq/9h/YpsoA+mSJuHs58oVVoMYAcFy8jLybP4O1enqSIxNCjIRk/Q5N2a4xIRIt9yO3aceRng76NqxDjYSTGFF6c29dryVBGAyYKybIatJCiISTREiIEzIvuRJjQYlWdu/YjL9mfxIjSm/9p80b84owF1eis1iTGJEQYjySREiIExRFIfvaj2rlSEcL7m3voMaiSYwqPUU9ffgO7dHK5opqzJWTkhiREGK8kkRIiH6yr/ko+qwcrezZvolg49EkRpSe3O9vgJMJqN6AuXIS5uKK5AYlhBiXJBESoh+d0UTW8uu1cuh4Pe5tG5A5BaOrb0O/brH8IizlE1EMKTvJVQgxhkkiJMQH5N74aXQ2h1Z2v7eecPvxJEaUXqI+L959O7SypXISlgrpFhNCjAxJhIT4AEOmk4xFy7VysKEWz67T16wSI8OzYyOcXNlbb8BSPQ1jQXFygxJCjFuSCAkxiNyb/hHFbIkXYjH63l0nO9OPkr6Nr2nHxoJiLJVTUHTyo0oIMTLkp4sQgzAXl+OYt0Qr+48ewLdvexIjSg+xgA/vrq1a2VI5GYvMFhNCjCBJhIQ4g7yP3nZqi41IhJ43XiAWCiY3qHHOs3PzqUUs9XpsU+dgcOac/UFCCDEMkggJcQbmiknYZi3Uyv6Du/HXyPYpI6lv0xvasbGgBMvEqUmMRgiRDiQREuIMFEUh7yOfBZ0eADUUpOeV/0ONxZIc2fgUCwbw7NiklS1VUzBXVCcxIiFEOpBESIizsE6ZhW36XK3s3b2FYENNEiMavzy7tqCe7HrU67HPvgi9VXaaF0KMLEmEhDgLRVHI/9TdoMT/qcT8Pjpf+KMssDgC3Jvf1I5NBSVYJkxJYjRCiHQhiZAQ52CbPBP7nIu1smfbu4SO1yUvoHEoFgrGt9U4wVw5WXaaF0KMCkmEhDgPhZ/9KuhPjRXqePp3SY5ofPHueQ816I8XdHocC5agM5mTG5QQIi1IIiTEeTCXVJK5+Aqt7Nn2LsHjDUmMaHzp2/yWdmwsLME6QWaLCSFGhyRCQpynws/ei3KilUKNhGn/4y+SHNH4oEbCeLa9rZUtVZMxyU7zQohRIomQEOfJkJWNc+XNWtmzYxOBeplBNlzePduI+X3xgk6HY+HlstO8EGLUSCIkxAXI//hd6Kz2eCEWpf2PP09uQONA35a3tGNjgXSLCSFGlyRCQlwAvdVG7k3/qJW9e7bhO7Q7iRGlNjUSwf3eeq1sqZyEqag0iREJIdKNJEJCXKDcGz+FPuvE/leqStvv/ye5AaUw7573iHk98YKiw3HxMhS9dIsJIUaPJEJCXCDFYKTg1i9q5cDRg/RtfSt5AaWw3nde1o6NhSXYqqcnMRohRDqSREiIIci6/FpM/Rb8a3/q58RkD7ILEvX7cL/3jla2TpqBsbAkiREJIdKRJEJCDIGi01F4x9e1crijBderzyUxotTj3roeNRyKFwxGMi5ehnJig1shhBgtkggJMUSOWQuxzVyglTue+R2xYCCJEaWW3rfXasfm0kqsE6clMRohRLqSREiIYSi64+ugKADEPH20//XXSY4oNYS7O/Dt36GVrdPmYsgrTGJEQoh0JYmQEMNgLptA5mVXa2XXa38n3NORxIhSQ89rfwdVBUCxWMlYeBnKiYRSCCFGkyRCQgxTwa3/hGI0AfENWVt/99MkRzS2qZEwrnX/p5UtE6ZiqZiUxIiEEOlMEiEhhsmYnUf2tR/Typ73N+Dd934SIxrbeje/SdTdGy8oChnzlmBw5iQ3KCFE2pJESIgEyL/lDnQ2R7wQi9L+x/+PWCiY3KDGqO41f9aOjcXl2GbMS14wQoi0J4mQEAmgs9jI+8htWjlw9CA9r7+QxIjGJl/NfoJ1R7SyfcYCzJXSLSaESB5JhIRIkJwbPoEhJ18r96x9mnBnaxIjGnv6b1Krz3TimLcYncmcxIiEEOlOEiEhEkTRG8j/5Oe1crj9OJ0v/Ak1Fk1iVGOHd/8O/Ad2aWXrtLlYZO0gIUSSSSIkRAJlLbseU0mFVnZvegPfQdmdXlVV2n7/iFbWOTKxz74IQ3ZeEqMSQghJhIRIKEVRKLj1n7Ry1O2i++VniPp9SYwq+fo2vEawvkYr22YtxDZltqwdJIRIOkmEhEiwjIWXY5k0Qyt7d27Gs31DEiNKrmjAR9uTp1qD9JlO7NPmYi6fmMSohBAiThIhIUZA4We+DLr4Py81FKTn1f8j3J2eK063/u+Pifa5tLJtziJs0+ai6OTHjxAi+eQnkRAjwDZ1DvZ5l2hl/5G99L77KuqJbSXShWfPe/RtWKeVjUVl2CZOw1whrUFCiLFBEiEhRkjhrV9CMVvjhVgM15trCB2vT25Qoyjq99Hyy/8+taeY0YRj7iXYZsxH0emTHJ0QQsRJIiTECDGXVpK1/HqtHG5ppPvlZ9JmOn3Lr39ApKtdK1smz8RcWS0LKAohxhRJhIQYQfmf+DyG3AKt3LfhNXyH9iYxotHRu+kN3Jvf1Mr67DxsU+fgmLdEZooJIcYUSYSEGEEGRwZ5H71DK8d8Hjqf/R2xYCCJUY2sSJ+Ltt89fKpLzGTGMWcR1knTMebIukFCiLFFEiEhRphzxfXYZszXyr592+l5c00SIxpZx3/xPaJul1a2TJyGuawK+6yLkheUEEKcgSRCQowwRW+g6K5/QWe1xytUla7/+z2hjpbkBjYCXG+9iHfnZq1syCvEOmU2joWXozNbkhiZEEIMThIhIUaBubSC3I9+VitHe3to+91PUWOxJEaVWJHebtqe/B+trLPYsM++GOuk6Zj7bTsihBBjiSRCQoyS3A99ckAXmWfHRrpffjqJESVW8yMPEvN54wVFwTJ5BuayKhz91lMSQoixRhIhIUaJotNT8pWH0DsytbqOv/4G/7HDSYwqMXpeex7f/h1a2VhYinXSTDIuXobOZE5iZEIIcXaSCAkxiozZuZR86f5T228EAzT/9H4iPk+SIxu6UFcHbU89ppV1Ngf2uYuxz5iHqaAkiZEJIcS5SSIkxChzzL+U3FW3auVw+3Ga/t+/EIum3kKLqqrS/LP7UU8uB6DTYZuxAEv5RGwzFyQ3OCGEOA+SCAmRBPmf+gK2ftPJ/Yf20PLz/0xiREPTvfZvBI7s08rm8olYq6eReckVso2GECIlSCIkRBIoikL5v/4/TKWVWl3fxtdo+e2PkxjVhQk2HqXjT7/UyrqMLGxzFpNx0eUDxkEJIcRYJomQEEmiM5mp/M6j6LNytDrXuudpffynSYzq/ET9Php/+E3USDheodeTcdEyHDPmYS6XneWFEKlDEiEhksiQlUPlg4+iz8jS6npeeZbW3z+SxKjOTlVVWn7534T7LQhpm7EA29TZ2OcuSmJkQghx4SQREiLJzCWVVDzwP+j6dSf1rH2a1t//LHlBnUXHM7/DveXUhqqmE2sFybggIUQqkkRIiDHAUlFN5Xce/UAy9MyYGzPU89rzdD37uFbWZzrJvOwaspZeg85iTWJkQggxNJIICTFGWCqqqbj/kQHJkGvd8xz/xfeJjYGtOFxvvUjrbx/WyorFRtYVN+JcfgOGzOwkRiaEEEMniZAQY4i1ajKVD/zPgDFDvW+v5fj/PJjUdYa6Vv+Jll/+N6jxhEwxmcm++mZyrr4ZY05e0uISQojhkkRIiDHGUjmJiu88ir5fK4t785s0/ejfiEUjoxpLLODn+P/3X7T/8RdanWI0kX31zeSu+jTGvKJRjUcIIRJNEiEhxiBL+UQqHnwMvTNXq/Pu3ETj979ONOAflRj8Rw9y7N8/R+/6l7Q6ndVO7k2fIe+Wu6Q7TAgxLqRdIuT3+/nOd77DlClTsFgslJSUcOedd9Lc3Jzs0IQYwFJaSeV3HsWQk6/V+fbvoOE/7yXS1zNizxv1umn53U+o+/bnCR2v1+oN2XkU/OM/kfeRz6C32kbs+YUQYjQpqqqqyQ5itAQCAa644go2b95McXExS5cupa6ujq1bt5Kfn8/mzZuZOHHoi8HNnDkTgH379p3jSiHOX+B4PY3fv49IV5tWZyoup+RLD2CdNCNhzxMLBuha82e6V/+J2AdanSyTZ1Lw6XuwTZ+HoigJe04hhDgpWb9D06pF6Hvf+x6bN29myZIlHD58mL/+9a9s2bKFn/zkJ3R0dHDnnXcmO0QhTmMpqaTiWz/GkH9qPE6opZGG73+drjV/Rh3muKGIq4vWJx/lyN2r6Hz6twOSIJ3VRubyGyi5+9+wz5gvSZAQYtxJmxahUChEQUEBvb29bN++nfnz5w84P3fuXHbv3s22bdtYuHDhkJ5DWoTESAoeb6DpJ98i1Fw/oN5cOYmCf/gi9jmLLihRCdQdoWv1n+jb/AZ8cEaaTodl4nSyLr+azMuuxtBvFpsQQoyEZP0ONYzqsyXRhg0b6O3tpbq6+rQkCOCWW25h9+7drF69esiJkBAjyVxSQfm//5TW3/wQ784tQPw7TLC+hsb//gbGwlIyLrmCzMUrsFRNQdENbPBVY1ECx47g3b2V3g3rCDUdO/1JdDpMxRU4FlyKY84ibDPmoejT5seEECINpc1PuF27dgGwYMGCQc+frN+9e/eoxSTEhTLlFlD61YfoeuFPuN5cQ7S3WzsXbmum++9/oPvvfwC9AWN2HvpMJ4rBQNTrIdx+HDUcGvS+isGIuXwC1unzsVZPwz5nkbQCCSHSQtokQg0NDQCUlZUNev5kfX19/aDn+zvZfPdBtbW1VFdXDzFCIc6P3uYg/xN3YZ9zMT2v/h++fe8TdfcOvCgaIdzZSriz9ez3cmRhnjAF65TZWCqrsU2ZjcGZM4LRCyHE2JI2iZDH4wHAZht82q/dbgfA7XaPWkxCDJWi02OfMR/rpBkEag/i3r4B34GdhNuOE3W7zvJABUN2HsaiMmxTZmOZOBVTYSnm0koUg3HU4hdCiLEibRKhRDrTQK4ztRQJMVJ0JjO26XOxTZ9LpK+HcEcrobZmQs31hHs6UUNBFIMRncWKpXISlkkzMGbnobPYUAzyz18IIdLmJ6HD4QDA5/MNet7r9QKQkZExajEJkUiGzGwMmdlYq6cnOxQhhEgZabOOUEVFBQBNTU2Dnj9ZX1lZOWoxCSGEECK50iYRmjt3LgDbt28f9PzJ+jlz5oxaTEIIIYRIrrRJhC677DKysrKora1l586dp51/5plnAFi1atUoRyaEEEKIZEmbRMhkMvHlL38ZgC996UvamCCAhx9+mN27d7N8+XJZTFEIIYRII2kzWBrg/vvv57XXXmPjxo1MnjyZpUuXUl9fz5YtW8jPz+d3v/tdskMUQgghxChKmxYhAIvFwptvvskDDzyAzWbj+eefp76+nttvv53t27cPa+d5IYQQQqSetNl0dTTIpqtCCCHE0CTrd2hatQgJIYQQQvQniZAQQggh0pYkQkIIIYRIW5IICSGEECJtSSIkhBBCiLQls8YSKCMjg3A4THV1dbJDEUIIIVJKbW0tRqMRt9s9qs8rLUIJZLfbMRqNyQ5jXKitraW2tjbZYaQlee+TR9775JH3Pjn6v+9GoxG73T7qMUiLkBiTZE2m5JH3PnnkvU8eee+TYyy879IiJIQQQoi0JYmQEEIIIdKWJEJCCCGESFuSCAkhhBAibUkiJIQQQoi0JbPGhBBCCJG2pEVICCGEEGlLEiEhhBBCpC1JhIQQQgiRtiQREkIIIUTakkRICCGEEGlLEiEhhBBCpC1JhIQQQgiRtiQREknl9Xp56qmn+MpXvsLixYsxm80oisJDDz00rPuuXr2a5cuXk5mZSWZmJitWrODFF19MTNDjyIYNG7jhhhvIycnB4XCwaNEinnzyyQu+zxNPPIGiKGf886lPfWoEoh/b/H4/3/nOd5gyZQoWi4WSkhLuvPNOmpubL/hePT093HvvvVRWVmI2m6msrORrX/saLpcr8YGPA4l676uqqs76uT548OAIvYLU9P777/ODH/yAj370o5SVlWnv01CN1ufekNC7CXGBjhw5wm233ZbQe/7sZz/j61//OgaDgZUrV2I2m3n11Ve58cYbefTRR/nyl7+c0OdLVc8++yyf/OQnicViLFu2jLy8PF5//XU++9nPsnv3bn784x9f8D3nzp3LvHnzTqtfvHhxAiJOHYFAgCuvvJLNmzdTXFzMTTfdRF1dHY8//jhr1qxh8+bNTJw48bzu1dnZyZIlS6ipqWHixIncfPPN7Nu3j0ceeYS1a9eyadMmcnJyRvgVpY5Evvcnffaznx20PisrKxEhjxv/+Z//yd///veE3GtUP/eqEElUU1Oj3nXXXeovf/lL9f3331e/+93vqoD64IMPDul+Bw8eVPV6vWo2m9WNGzdq9YcOHVJzc3NVg8GgHjlyJEHRp66uri41MzNTBdRnn31Wq29tbVUnTZqkAuqbb7553vd7/PHHh/X3Nt58+9vfVgF1yZIlqtvt1up/8pOfqIC6fPny877XrbfeqgLqRz/6UTUcDmv1X/nKV1RA/exnP5vAyFNfIt/7yspKVX5Nnr8f/OAH6gMPPKC+8MILaktLi2o2m4f8/o3m517+hsWY8t///d/D+oX6xS9+UQXUe++997RzDz/8sAqoX/7yl4cX5Djw//7f/1MB9aabbjrt3HPPPacC6o033nje95NE6JRgMKhmZWWpgLp9+/bTzs+ZM0cF1G3btp3zXsePH1d1Op1qMpnU1tbWAecCgYCan5+v6vV6ta2tLWHxp7JEvveqKonQcA01ERrtz72MERLjyslxQLfccstp507WrV69elRjGovO9j596EMfwmKx8NprrxEIBEY7tJS3YcMGent7qa6uZv78+aedv5DP4csvv0wsFmPp0qUUFhYOOGc2m1m1ahXRaJSXXnopMcGnuES+9yJ5RvtzL4mQGDdcLhcNDQ0Ag/4QLC8vJy8vj/r6evr6+kY7vDFl165dACxYsOC0cyaTiVmzZhEIBDh8+PAF3ff999/nX/7lX7j77rt58MEHWb9+fULiTSVne2/71+/evXtU75UORur9+tGPfsQ999zDvffey69//Ws6OjqGF6g4q9H+3MtgaTFunEyCsrOzsdvtg15TVlZGZ2cn9fX1zJ49ezTDGzP6+vro7e0F4u/HYMrKyti2bRv19fXMmTPnvO+9Zs0a1qxZo5W/+93vsnz5cv7617+e9s1uvDr5OTzbewtQX18/qvdKByP1fv3rv/7rgPLXv/51Hn30Ue68884hRCnOZbQ/99IiJMYNj8cDgM1mO+M1JxMkt9s9KjGNRSffJzjze3Wh71NxcTEPPfQQO3bsoLe3l9bWVl544QWmTZvG+vXrufHGG4lGo8MPPgWc63N4Ie9tIu+VDhL9fn34wx/mueeeo76+Hp/Px969e7nvvvsIBoN87nOfS9gMKTHQaH/upUVIDMtHPvIRDhw4cEGPefLJJ1m0aNEIRZQextr7fu2113Lttddq5czMTFatWsUVV1zBwoUL2bZtG3/729/49Kc/PSLPL8RI+J//+Z8B5ZkzZ/KTn/yEadOm8YUvfIFvfvOb3HTTTUmKTiSKJEJiWI4dO8ahQ4cu6DE+n29EYnE4HOe8v9frBSAjI2NEYhgtw3nfT75PJ+syMzNPuzZR75PD4eCrX/0qX/7yl3nllVfSIhE61+fwQt7bRN4rHYzW+3XXXXdx//33c+jQIerq6qiqqhrW/cRAo/25l0RIDMvOnTuTHYKmoqICiK9G6vV6Bx0n1NTUBEBlZeWoxpZow3nfMzMzycrKore3l6amJmbMmHHaNYl8nyZPngxAS0vLsO+VCk5+Dk++hx90Ie9tIu+VDkbr/dLpdFRXV9Pe3k5LS4skQgk22p97GSMkxg2n06n9A9qxY8dp5xsbG+ns7KSysnLQVpB0MnfuXAC2b99+2rlwOMzevXuxWCxMmTJl2M/V09MDcMYB7OPN2d7b/vXnMwg9kfdKB6P5fqXb53o0jfbnXhIhMa586EMfAuCZZ5457dzJulWrVo1qTGPR2d6nNWvWEAgEWLlyJRaLZdjP9eyzzwJnngo73lx22WVkZWVRW1s7aMvdhXwOr7vuOnQ6He+88w7t7e0DzgWDQVavXo1er+eGG25ISOypLpHv/dns27ePQ4cOYbPZmDZt2rDuJU436p/7hCzLKESCnO/K0lOnTlWnTp2qNjU1Dajvv8XGpk2btPrDhw/LFhv9nGmLjba2trNusXGm9/2//uu/1I6OjgF1oVBIfeihh1RAtVqtpz1mPDu5zcOll16qejwerf5M2zw8+uij6tSpU9V/+7d/O+1eJ7ca+NjHPjZgq4GvfvWrssXGIBL13r/44ovq66+/ftr9d+3apU6fPl0F1K9+9asj8hrGi3OtLD1WPveSCImku/nmm9XFixerixcvVsvLy1VALS0t1epuvvnm0x4DqIB67Nix086d3ErDYDCo119/vXrTTTepVqtVBdT/+Z//GYVXlBqeeeYZVafTqYqiqFdccYV6yy23qE6nUwXU++67b9DHnOl9B1Sz2axedtll6qc+9Sn1hhtuUEtKSlRAtVgsA5KtdOD3+9XFixergFpcXKx+4hOf0Mr5+flqbW3tgOsffPDBM/5w7+joUKurq1VAra6uVj/5yU+qs2bNUgF18uTJaldX1yi9qtSQqPf+ZH1lZaX64Q9/WP3Upz6lLlq0SDUYDCqgrlixQvX5fKP4ysa+NWvWaD+3Fy9erCqKogID6tasWaNdP1Y+95IIiaQ7uZ/Pmf5UVlae9pizJUKqqqovvPCCunTpUtXhcKgOh0NdunSpunr16pF9ISno3XffVa+77jrV6XSqNptNveiii9QnnnjijNef6X3/zne+o1599dVqRUWFarVaVYvFok6aNEm9++671YMHD47wqxibfD6f+sADD6jV1dWqyWRSi4qK1Ntvv11tbGw87dqz/UJQ1XgL3le+8hW1vLxcNZlManl5ufrVr35V7enpGdkXkaIS8d5v3LhRvfPOO9XZs2drrck5OTnqihUr1N/85jdqJBIZpVeTOk7uOXi2P48//rh2/Vj53CuqqqqJ6WQTQgghhEgtMlhaCCGEEGlLEiEhhBBCpC1JhIQQQgiRtiQREkIIIUTakkRICCGEEGlLEiEhhBBCpC1JhIQQQgiRtiQREkIIIUTakkRICCGEEGlLEiEhhBBCpC1JhIQQQgiRtiQREkKMGq/Xy8MPP8wVV1xBYWEhJpOJ7OxslixZwne+8x0aGhoAOHz4MIqikJGRgc/nO+d9b7jhBhRF4dFHHx3plyCEGGdk01UhxKjYuHEjH/vYx2htbcVms3HJJZdQWFhIb28v7733Hh0dHZjNZtasWcPKlStZvHgxW7du5U9/+hOf/vSnz3jf9vZ2SktLATh+/Dj5+fmj9ZKEEOOAtAgJIUbczp07ueqqq2htbeWb3/wm7e3tvP766/zpT3/ixRdfpLW1lWeffZaysjKampoA+MxnPgPAH/7wh7Pe+y9/+QuRSITrrrtOkiAhxAWTFiEhxIhSVZU5c+awd+9eHnroIR588MEzXtvb20tjYyOzZs2is7OTkpISVFU9a0vPokWLeO+99/jLX/7CJz/5yZF6GUKIcUpahIQQI+rll19m7969lJWV8e1vf/us12ZlZTFr1iwA8vLyuPbaa4lEIvz1r38d9PojR47w3nvvkZmZyYc//OEhx/jQQw+hKApPPPEEW7Zs4dprr8XpdJKZmcnVV1/N5s2bz/teO3fuxGw2k5OTo7Vu9ff5z38eRVG4++67hxyvECJxJBESQoyoF198EYCPf/zjGAyGC3rsubrHTtbfcsstWK1W4FRSc/vtt19wrBs3bmTZsmU0NTVx/fXXM3XqVF577TWWL1/Oq6++el73mDdvHt///vfp6enhs5/9LP0b3Z9//nn+93//lylTpvDwww9fcHxCiMSTREgIMaJ27twJwIIFCy74sR/+8IfJyspiy5Yt1NTUnHb+j3/8I3AqYRqu3/zmN/zLv/wLe/fu5c9//jPvvfcev/jFLwiFQtx+++34/f7zus83vvENrrzySt544w1+8pOfANDS0sLnP/95jEYjf/zjH7Hb7QmJWQgxPJIICSFGVFdXF8CQBjJbLBZuueUW4PRWoU2bNlFbW0t5eTnLly/X6vPy8pg6dSrFxcUX/HyVlZVai9JJX/ziF1m8eDEtLS08++yz53UfRVH4/e9/T3Z2Nt/+9rfZsWMHt99+O52dnTz44INcdNFFFxybEGJkSCIkhBjT/vEf/xE41fpz0snyrbfeOiBx+fKXv8zBgwf57//+7wt+ro997GODdt+dnL7/zjvvnPe9ysrK+NWvfkUoFGLFihW8+uqrXH755fzbv/3bBcclhBg5kggJIUZUbm4uAB0dHUN6/PLly6moqKCmpoYtW7YADBhAnahuMYi3CA2mqqoKiK9TdCE+/vGP85GPfIS+vj5sNhtPPfUUer1+uGEKIRJIEiEhxIiaN28eANu3bx/S4xVF4dZbbwVOdY+9/PLLdHZ2smDBAmbMmJGQOEdCS0uL1ork8/nYv39/kiMSQnyQJEJCiBH1oQ99CICnn36aSCQypHucbPX561//SiQS0RKik91miVJfX3/W+pKSkvO+l6qq3HHHHXR2dvLpT38avV7PnXfeOeSWMSHEyJBESAgxoq677jpmzpxJU1MT3//+9896bV9fH/v27Tutfvr06SxYsICOjg6effZZXnjhBfR6/Vm33hiK5557jmg0elr9X/7yFwAuv/zy877Xo48+yiuvvMJll13GU089xbe//W3a2tq46667EhavEGL4JBESQowoRVH4wx/+gMVi4aGHHuJb3/oWXq93wDWqqvLCCy9w0UUX8d577w16n5OtQl/+8pfx+/1cffXVFBUVnXbdY489xrRp0/jWt751wbHW1dXxH//xHwPqfv3rX7Np0yYKCwv52Mc+NuDcVVddxbRp09i6deuA+v379/PNb36TjIwMbVzQAw88wKJFi1i9ejW//OUvLzg2IcTIkERICDHi5s2bx2uvvUZhYSE/+MEPKCgoYOXKldx6663ceOONFBcXc9NNN9HY2Eh5efmg9zjZvdTZ2QmceZB0Z2cnhw4doqWl5YLj/PznP88PfvADZs2axT/8wz+waNEi7r77boxGI0888QQ2m23A9bW1tRw6dAifz6fVhUIhbr31VgKBAI899hgTJkwAwGAw8Ic//AG73c43vvENDh06dMHxCSESTxIhIcSouOyyy6ipqeHHP/4xF198Mbt37+Zvf/sbGzZsoKqqigcffJAjR45w1VVXDfr4wsJCrrnmGgAcDgc333xzwmO89NJLWb9+PUVFRaxZs4YDBw5w1VVX8dZbb3Hddded1z3uv/9+du7cycc//nFuu+22AecmT57Mww8/jM/n49ZbbyUcDif8NQghLoxsuiqESHsPPfQQ//Ef/8Hjjz8+pK05hBCpS1qEhBBCCJG2JBESQgghRNqSREgIIYQQaUvGCAkhhBAibUmLkBBCCCHSliRCQgghhEhbkggJIYQQIm1JIiSEEEKItCWJkBBCiP+/3ToQAAAAABDkbz3IRRFsiRAAsCVCAMCWCAEAWyIEAGyJEACwJUIAwJYIAQBbIgQAbIkQALAVXtPAuqt6hPUAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cv_name = 'p.x'\n", "cv1 = colvar[cv_name].values\n", "\n", "fes_params = {\n", " 'blocks': 3, # Number of blocks for error analysis\n", " 'bandwidth': 0.01, # Kernel bandwidth (sigma) for density estimation. if 'range', the bandwidth is expressed as ratio of the range of values (e.g. here it is 1% of it)\n", " 'scale_by': 'range', # Method to scale the bandwidth\n", " 'temp' : temperature, # temperature for proper energy scaling (alternative to kbt)\n", " 'fes_units': 'kJ/mol', # units of the free energy\n", " 'weights': weights, # Statistical weights from the bias \n", "}\n", "\n", "fig, ax = plt.subplots(figsize=(4,3),dpi=150)\n", "fes1D, grid1D, bounds1D, error1D = compute_fes(\n", " cv1,\n", " plot=True,\n", " ax = ax,\n", " **fes_params\n", ")\n", "ax.set_xlabel(f'CV: {cv_name}')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2D free energy surface" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adjusting regularization (eps) to 1.0e-13 to avoid NaNs.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/lbonati@iit.local/work/code/mlcolvar120/mlcolvar/utils/fes.py:248: RuntimeWarning: invalid value encountered in log\n", " * np.log(kde.evaluate(cartesian(pos)) + e)\n", "/home/lbonati@iit.local/work/code/mlcolvar120/mlcolvar/utils/fes.py:248: RuntimeWarning: invalid value encountered in log\n", " * np.log(kde.evaluate(cartesian(pos)) + e)\n", "/home/lbonati@iit.local/work/code/mlcolvar120/mlcolvar/utils/fes.py:248: RuntimeWarning: invalid value encountered in log\n", " * np.log(kde.evaluate(cartesian(pos)) + e)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAGiCAYAAAChyG+jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABKU0lEQVR4nO3de3hU1b0+8DcJZJKQTBAiGSIhImAEJIAoMdgCSjREVKj86u2cAh6KtQe0EqsWHyUIbaFFBW0RPJaLPZSiWNEeUWjFAhWDFiQFb6mhCKGSULDkJgkE1u8POsNMMpd9Wfs67+d55nnIzJ49azaT7He+a+21EoQQAkREREQGSLS6AUREROReDBpERERkGAYNIiIiMgyDBhERERmGQYOIiIgMw6BBREREhmHQICIiIsMwaBAREZFhGDSIiIjIMAwaREREZBgGDSIiIodYtmwZCgoK4PV64fV6UVRUhLfeeivw+JgxY5CQkBByu/feey1sMZDAtU6IiIic4f/+7/+QlJSE/v37QwiBF198EYsWLcKePXswaNAgjBkzBpdeeinmzZsXeE5aWhq8Xq9lbe5k2SsTERGRKjfffHPIzz/5yU+wbNky7Ny5E4MGDQJwLlj4fD4rmheWI4LG2bNn8eWXXyIjIwMJCQlWN4eIiGxMCIHGxkbk5OQgMdGYEQItLS04deqUlH0JITqc2zweDzweT9TnnTlzBuvXr0dzczOKiooC9//mN7/BmjVr4PP5cPPNN+Pxxx9HWlqalLZqIhygpqZGAOCNN9544403xbeamhpDzkknT54UPbokS2tnenp6h/vKy8sjvv7evXtFly5dRFJSksjMzBQbN24MPPb888+LTZs2ib1794o1a9aIiy66SHzrW98y5Dgo5YgxGvX19ejatSs+nDkK6R5HFGHIJNkjU6I+3tS/UdX+Nrd019OcuFaSclzxtumfZ0R8rO69FhnNoTjW1NqGK365HSdOnEBmZqb0/Tc0NCAzMxO7Z45Chs5zUmNrG4b/cjtqampCxlFEq2icOnUKhw4dQn19PV555RX86le/wrZt2zBw4MAO277zzjsYO3Ysqqur0bdvX11t1coRZ21/SSnd00n3fyq5i7dL54iPNeU3wKvwI77xZBYAIC1dSrMM88bBoQCAm/IqLW1HOH+GD+NTjynb+IqTSK8KPzjNe31n1G4/KbFlFK+M7mrPkHhO8l9FokRycjL69esHABg+fDj+8pe/4JlnnsHzzz/fYdvCwkIAYNAg0sI3KjXiY035DYr34w8ZVvOHCFnbWhFGNp7MUh42iEiKs2fPorW1NexjlZWVAICePXua2KJQDBrkSE4PGWpChazXMCt4KA0bTfkNEasavlGprGoQhTF79myUlpaid+/eaGxsxNq1a7F161Zs3rwZ+/fvx9q1a3HjjTeie/fu2Lt3L2bNmoVRo0ahoKDAsjYzaFBcMjtgmBEs1LTB6NDBygaRMY4ePYrJkyfjyJEjyMzMREFBATZv3ozrr78eNTU1ePvtt7FkyRI0NzcjNzcXkyZNwmOPPWZpmxk0yHH0VjPMChl2CBeRmBE6lISNaFUNIupoxYoVER/Lzc3Ftm3bTGyNMgwa5BpWhww7B4to7DrAlN0nRO7AoEGOEqmaYVXIcGq4CMeIwMGqBhExaJBjROsyiUVmyHBTuAhHduDQM16DVQ0i52PQIMeLVc2QETLcHi7CMbNLhVUNIvfiMvHkCFq7TPSGjDcODo3LkBFMxvu3y1wlRGQ+Bg2iMBgwQplxLNTMf0JEzsGgQbZnZjWDASMyvcdFa1VDz9gcIrIegwbZmlkhgwFDGaOPEasaRO7DwaBkK0q+vRoRMki5Nw4O1TxAlDOGEsUfVjTIFnyjUmOGjKb8BoYMm9Bz3LR0obD7hMi5GDTIcjKqGFowZOhjVHcTu0+I3IVdJ2QpJVUMpdR8U7Y6ZFQfzlD9nH69Gg1oiX5aulLYhUIUPxg0yDJWLfVuZsjQEiiU7stOwUPPuI1wOIEXkXswaJAl3BoyZAYLNa9lh9ChNmywqkEUHxg0yHQyQobWOTKMYGa4iNUGqwOH7MoGETkfgwaZSm/I0Drpk+yQYYdwEU714QzLw4Ya0aoa7D4hcgdedUKmiTb5llNCRvXhDNuGDD+r22f1QFsishdWNMgUWmf4BOwTMJzEaZUNInIvVjTIcE4OGU6oYERiZbtZ1SAiP1Y0yFBmL+8uM2C4ASsbRGQ1Bg0yjNmrrsrgloBBRGQX7DohU9l51VW3hgy7vy+t1SsicgZWNMgQWhbBsmoKcbufiN2Ok3YRuRsrGiSdli4TK0KGkwd6qhUv75OI7IcVDTKFjJDBLhIiIudh0CCpwlUz7BIyGDDMw2nIiciPQYOkMLK7RG/IYMAgIrIOx2iQYRgy7MWsY8FqBhEFY0WDdDOiu4QBw5nUhoxoV5xwQTUid2DQIF3sFjIYMKzDSgYRhcOuE5KKISM+aQkZnD+DKD6wokGaqZmUy8iQYaeA8dWOf3W4r9s1F1jQEnMYVcVgtwmRezBokDSRqhluDRnhQoWS7ZwePGSEC1YziOIHgwZp0r6aYXbIsCJgKA0WSvfjlMAhu2oRK2SwmkHkLgwapIqawZ9uCBmywkW0fds5cJgdMiKp3X5SajuIyDwMGhRTtLEYZoYMswKGkeEi0uvZMWxYETJYzSByHwYNisqoGT/tFjLMDhfhXt+OYUMWhgyi+MXLWymieAgZX+34l+Uhw88u7ZBN78BPdpsQORsrGhSWHUKG0QHDjtxW2VAaMljNIHIvBg1SzA0hw64BI5gbwoaaKka0kMFqBpHzMWhQB1ZPK25EyHBCwLATrQNB1XaTMGQQuR+DBsXk5JDBgGE8rWMwGDKI4gMHg1KAb1Sq4om4ALkho/pwBkNGECe0fXzqMUNCBhFFtmzZMhQUFMDr9cLr9aKoqAhvvfVW4PGWlhbMmDED3bt3R3p6OiZNmoS6ujoLW8yKRtxz46BPJ5yknc7IgMFqBlFkvXr1wsKFC9G/f38IIfDiiy9iwoQJ2LNnDwYNGoRZs2Zh48aNWL9+PTIzMzFz5kzceuut2LFjh2VtZtCIU1oCBmDvkMGAEV2/Xo1S9sOQQSRfQ0Po316PxwOPx9Nhu5tvvjnk55/85CdYtmwZdu7ciV69emHFihVYu3YtrrvuOgDAqlWrMGDAAOzcuRNXX321cW8gCgaNOKFlds9gSgIGwJBBHTFgkFtlj0yBt0tnXftIaz4NPAXk5uaG3F9eXo65c+dGfe6ZM2ewfv16NDc3o6ioCLt378bp06dRXFwc2Oayyy5D7969UVFRwaBBxnFzyHBzwHDyZa5qxmAwZBABNTU18HrP/96Eq2b47du3D0VFRWhpaUF6ejo2bNiAgQMHorKyEsnJyejatWvI9tnZ2aitrTWq6TExaLicnpBhRMAA5IQMNwcMq0W7tDVat4mWAZ4MGUTn+Ad3KpGfn4/KykrU19fjlVdewZQpU7Bt2zaDW6gdg0acYsggGbRePcKAQaRdcnIy+vXrBwAYPnw4/vKXv+CZZ57B7bffjlOnTuHEiRMhVY26ujr4fD6LWsugEXdkBQzA/JDBgOF8DBhE8p09exatra0YPnw4OnfujC1btmDSpEkAgKqqKhw6dAhFRUWWtY9Bw8XUzIkBsIrhdtWHM6RdeRINwwSRcWbPno3S0lL07t0bjY2NWLt2LbZu3YrNmzcjMzMT06ZNQ1lZGbp16wav14v77rsPRUVFlg0EBRg06N8YMsjvjYNDFU1B7u82YbAgMs/Ro0cxefJkHDlyBJmZmSgoKMDmzZtx/fXXAwAWL16MxMRETJo0Ca2trSgpKcFzzz1naZsZNOKE2Uu7A+wqiQcMGUTmWrFiRdTHU1JSsHTpUixdutSkFsXGKcjjAEMG+Rm1Ki4RUSSsaLiUf3yGnpBhRcAAGDL87DiHxsaTWZpnBiWi+MSKRpxiyCAiIjOwouFCsaoZskMGZ/mMP1x9lYiUYtBwgXCzfzJkOJsdu02IiLRg0HAgrdOKywwZDBjEK06ISAmO0XAYhgz3YzWDiNyEFQ0HiRQy9M74yZBhHwwZROQ2DBoOYVTIUINXlRAvbSUitdh14mAyQobSagZDhvHMqGaYsdYJEVEwBg2HYsggoyhZ54SISCl2nTiAmlVYZS+OxpBhDidVMziHBhGpwaDhMGaGDL0YMOxDdpcJL20lIqXYdeISRoQMPdUMhgxnitZtwoGgRKSF6qCxfft23HzzzcjJyUFCQgJee+21qNtv3boVCQkJHW61tbVa2xy3tM72CZwLGAwZ9mR0twkHgBKRlVQHjebmZgwZMkT1WvdVVVU4cuRI4NajRw+1Lx3X9IYMNRgy3ENNyGA1g4iMoHqMRmlpKUpLS1W/UI8ePdC1a1fVzyN9zBqPATBk2AmrGERkF6aN0Rg6dCh69uyJ66+/Hjt27Ii6bWtrKxoaGkJu8SralON2XeqdlJPdbdKvV6OmkMFqBhEZxfCg0bNnTyxfvhy/+93v8Lvf/Q65ubkYM2YMPvzww4jPWbBgATIzMwO33Nxco5vpOEaswsouE+fSGjAAzptBRMYy/PLW/Px85OfnB34eOXIk9u/fj8WLF+N///d/wz5n9uzZKCsrC/zc0NAQl2HDX82INTlXe2Yv9c6QoY2MaobeLpJYIYPVDCLSy5J5NEaMGIF333034uMejwcej8fEFtmP1i4Thoz4YXTIICKSwZKgUVlZiZ49e1rx0o4na5E0jsWwlt5qhhkhg9UMIpJBddBoampCdXV14OcDBw6gsrIS3bp1Q+/evTF79mz84x//wK9//WsAwJIlS9CnTx8MGjQILS0t+NWvfoV33nkHf/jDH+S9C5eJVs2IxsxpxQFWM7Syeil4hgwiMpPqoLFr1y5ce+21gZ/9YymmTJmC1atX48iRIzh06FDg8VOnTuHBBx/EP/7xD6SlpaGgoABvv/12yD7oPDXrmgRjyHAGK8dlKO0qiRYyuM4JEamlOmiMGTMGQoiIj69evTrk54cffhgPP/yw6obFm3BVjPYhQ2+3CUOGtawKGWrGYigJGVznhIjU4KJqFovUTaLmShNOymV/VnWXyA4ZRERqMWhYJNo4jHAhww7VDIYMbWSFDFnTiYfDkEFERmHQsIDakKEXQ4Z1rKhkyKpiAAwZRKQfl4m3ES0Lp5nRbcKQoY3MkKG0msGQQUR2w4qGyZQM+pSJE3I5n+yQoeTSVYYMIpKFQcNEWkKGldUMhgztzB6XIXNujGghg1ecEJFaDBom0HJliYxF07hImjWsGPwZC0MGEVmFQcNgRoQMsi8r5sqQsTBarK4Shgwi0opBwyBaryxREjJYzbAnJ4YMBgwiMhqDhgGsDhl6MGRYw4gZP8OFDAYLIjIbg4ZkdggZWqsZDBnaaa1mGLVuSaSQwSBBRGbjPBomsXvIIO2cEDKIiKzCioZEsi9fBdR1lXBchvnsFjKIiOyGQUMCrQujyQwZZD4zQ4beybg4ARcRWYVBQwc9a5bIDhmsZpiLIYOISBkGDY0YMuKXWSFD5pTiAK8oISJrMGhoYJeQwYGf5jNjNVYujEZEbsKrTlTSc/mq3UIGqxnq6AkZXH2ViGRYsGABrrrqKmRkZKBHjx6YOHEiqqqqQrYZM2YMEhISQm733nuvRS1mRUMVO8yRATBkWMFOlQyuvkoUv7Zt24YZM2bgqquuQltbGx599FHccMMN+OSTT9ClS5fAdtOnT8e8efMCP6elpVnRXAAMGorZJWSQ+fSGDCXVDFmrrzJgEDlTQ0PoecTj8cDj8XTYbtOmTSE/r169Gj169MDu3bsxatSowP1paWnw+XzGNFYlBg2dzA4ZrGaYy4xKhhJ6QwYHghLJ19S/EYkZ+k6jTY1tAIDc3NyQ+8vLyzF37tyYz6+vrwcAdOvWLeT+3/zmN1izZg18Ph9uvvlmPP7445ZVNRg0dGDIcDezFkmLVs3g8u5E8aGmpgZe7/nf43DVjPbOnj2LBx54ANdccw0uv/zywP133XUX8vLykJOTg7179+KRRx5BVVUVXn31VUPaHguDhkZODBmknB0qGVrnxWCwIHIer9cbEjSUmDFjBj766CO8++67Ifffc889gX8PHjwYPXv2xNixY7F//3707dtXSnvVYNBQoP34DKeGDFYzlJEVMmSNzQA6hguGCaL4NnPmTLzxxhvYvn07evXqFXXbwsJCAEB1dTWDhhMwZLiXHaoYfsHVDIYMIvITQuC+++7Dhg0bsHXrVvTp0yfmcyorKwEAPXv2NLh14TFo2BhDhnlkhwyZ1YxgDBlE8W3GjBlYu3YtXn/9dWRkZKC2thYAkJmZidTUVOzfvx9r167FjTfeiO7du2Pv3r2YNWsWRo0ahYKCAkvazKChgtPmyiBlrAgZSgVXMxgyiGjZsmUAzk3KFWzVqlWYOnUqkpOT8fbbb2PJkiVobm5Gbm4uJk2ahMcee8yC1p7DoCGBnUMGqxmRGdFVImMGUKVXmhBR/BFCRH08NzcX27ZtM6k1ynAKcoUiVTM4IZczWRkytGA1g4icihUNg1l5GSurGR0ZNeBTTchQW81gyCAiJ2NFQwc7LffeHkNGR3YPGe1xOnEicgNWNGKItsaJTBz8aRwjL1uV2V3CsRlE5EasaGgku5ohE6sZ59kpZGgZAMpuEyJyOgYNG2A1wxhOCRlERG7GrpMo/N0m0ebPCIfzZVjL6Bk+ZYcMdpkQkZuxohFBtLEZSi5ptUq8d5u4JWRwICgRuQUrGv8ma9An58ywjp26SgBWMoiIgDgIGnoDhNpuEzXYbSKPnUKGkvEY0UIGqxlE5CauDBpawoWRgcIs8dht4rSuEkB5yOAVJ0TkBq4KGmoChtZgYefxGfHGTlUMgJUMIqJwXBM0lIYMu1Qu2G2ijx1m+QymN2QQEbmVK4KGkpBhRsDgQFDj2S1gAKxkEBFF4/igEStkyAwY7Daxjh0DBsCQQUQUi2ODhlkBg+HCevEYMjgQlIjcwrFBIxoO9HQHuw329FM6nThDBhGRQ4NGtGqG2pDBcGFP8RgyGDCIyI0cFzRkhAwjwoUdBoJ2u+YCx8+lYbd5MfxkBAwgfMhgwCAiN3NU0MgemRLxMStDBulndMAAjJsbw09tyGDAIKJ44KigEYmSkGFkwLBDNcOpzAgYgP1CBhFRvHB80HBiyDBysi6ndJ/ES8AA2F1CRPHN0UHDiSEj3pkVMADrQwYHfRIROTxoxGJUyLB7wLBjVcPuAQOQN+ATYFcJEZGfY4NGrGqG3pBhVJgwa40TO4QNM8OFn91DBqsZRBRvHBs0otEaMoyuVMTDQmpWhAvA+rkxgNhVDIYMIopHjgwa0aoZWkKG3btCtPKf9I2sbFgVLIJZHTIYMIiIInNk0JDFzIBhZTVDZjeKHYJFMDuHDAYMIiIHBg1Z1Yx4CRl+WqsbdgsWfmYsiAYoDxkMFURE4TkuaMjg1q4SJewaHNQwo4oBMGQQEcmQaHUDZLHr1OJ2qGa4id1CBhERReeoikZT/0Z4dTY5nqsZTmf0pasA1yshIpLNNRUNJcwOGaxmyGPHkEFERLE5qqIRiV27TUg/MwIGwJBBRGQUVwQNO2I1Qz87VDEArllCRKRH3AQNjs1wDjsM+PRjJYOISJ+4CRpmYjVDO6sn4ArGybiIiPRj0CBbsMMEXH6cUpyISB7HBw0OBHU2vQED4MqrRER25vigYTfsNlFGRsAAuPIqEZHdqZ5HY/v27bj55puRk5ODhIQEvPbaazGfs3XrVlxxxRXweDzo168fVq9eraGp5Ab9ejU6KmTUbj/JkEFEpIPqoNHc3IwhQ4Zg6dKlirY/cOAAxo8fj2uvvRaVlZV44IEH8N3vfhebN29W3dj22G3iHDIDBmBsyPCHCwYMIrKbBQsW4KqrrkJGRgZ69OiBiRMnoqqqKmSblpYWzJgxA927d0d6ejomTZqEuro6i1qsoeuktLQUpaWlirdfvnw5+vTpg6eeegoAMGDAALz77rtYvHgxSkpKwj6ntbUVra2tgZ8bGjqu2MqQ4Qwyw4WfUSGDwYKI7G7btm2YMWMGrrrqKrS1teHRRx/FDTfcgE8++QRdunQBAMyaNQsbN27E+vXrkZmZiZkzZ+LWW2/Fjh07LGmz4WM0KioqUFxcHHJfSUkJHnjggYjPWbBgAZ544gmDWyYfx2ecZ0TAABgyiCi+bdq0KeTn1atXo0ePHti9ezdGjRqF+vp6rFixAmvXrsV1110HAFi1ahUGDBiAnTt34uqrrza9zYavdVJbW4vs7OyQ+7Kzs9HQ0ICTJ8P/cZ89ezbq6+sDt5qampDH1VYzOFmXeWR3kfjdlFdp2sBPIiKzNTQ0hNyCq/rR1NfXAwC6desGANi9ezdOnz4d8gX/sssuQ+/evVFRUSG/4QrY8qoTj8cDj8cT9jF2mdiTURUMgAujOZVvVGrEx1hBIjfY3NIdaZ0669rH1y2nAQC5ubkh95eXl2Pu3LlRn3v27Fk88MADuOaaa3D55ZcDOPflPjk5GV27dg3ZNjs7G7W1tbraqpXhQcPn83UYhFJXVwev14vU1Mh/iMLZ3NIdaekyW0cyWN1N4qc1ZPCkp0+0QKH2Ofy/oHhVU1MDr/f836hIX7aDzZgxAx999BHeffddI5umm+FBo6ioCG+++WbIfX/84x9RVFRk9EsDYLeJkexSxQCUTytO+mkJFlr2zdBB8cTr9YYEjVhmzpyJN954A9u3b0evXr0C9/t8Ppw6dQonTpwIqWrU1dXB5/PJbLJiqsdoNDU1obKyEpWVlQDOXb5aWVmJQ4cOATg3vmLy5MmB7e+99178/e9/x8MPP4zPPvsMzz33HF5++WXMmjVLzjuIgiHDGEaNwwDUjcXw4+Jo5vCNSjU0ZFj9ekROIITAzJkzsWHDBrzzzjvo06dPyOPDhw9H586dsWXLlsB9VVVVOHTokGlf8NtTXdHYtWsXrr322sDPZWVlAIApU6Zg9erVOHLkSCB0AECfPn2wceNGzJo1C8888wx69eqFX/3qVxEvbSX7MrKCAaivYpB5rDzh+1+bFQ6ic90la9euxeuvv46MjIzAuIvMzEykpqYiMzMT06ZNQ1lZGbp16wav14v77rsPRUVFllxxAmgIGmPGjIEQIuLj4Wb9HDNmDPbs2aP2pXSxoprRr1ejay9xNTpkWMk3KpUnsQjsVFHg/xMRsGzZMgDnzqvBVq1ahalTpwIAFi9ejMTEREyaNAmtra0oKSnBc889Z3JLz7PlVSd6sctEHrMChtXVDJ7EOtISMpryO06uF43aLi1WNyjeRfui75eSkoKlS5cqnsHbaK4MGqSPG6sXTfkNMU9qDBvnqA0YasNFpOeqCR38vyJyDsMn7DKb1dUMJ5+kjRzkGYsZ/296TojxQk3IaMpvkHpM/ftTuk87desQUWSsaJCjw5FasSobWr8paznp2eUbuZkVDDWvoaTCwcoGkf25KmhYXc3wc8qgULcEjI0ns1TNo6EkbPi1P4nJ/BYdbl9mnTTNGH+hl//12OVF5GyuChp2YtewYedw8cbBoZoHhcoOG35ml+eNCh9634eV3U4cX0PkbK4Zo2GXakYwu5zU/WMv7NIeo6hdB8cpYza0Tlzlf56TQ4aaNnDMBpE9saJhMKsqG24PFZH4w4bS6obSyoYdKJ2eW9YJ1w4BI5jSrhQishdXBA07VjOCmRE23BAsZM6loaYrxUlhw8/Ib+96A4bSypIR69OwC4XIflwRNJxAdthwQ7AwmllhQ++J2eqQI6tyobbrSm31yc+oK4eIyBgMGibyhwMtgcPtwcKomUHVhg0rtH9dM4KHzPeqNmBEer4Rg3mJyHqODxp27zYJJ1xoCA4fbg8V7Vk9/bjdaJ0tU81+ZdEbMsLtS0YVilOVE9mH44OGW7glXEQLDe1DoVkBQ+2lr3aidwCkUVUamQEj3L5l/X+xG4VIuWPHjqG5uRl5eXmB+z7++GM8+eSTaG5uxsSJE3HXXXep3i+DBkmhJDRYWblwctgA1AUOo7uAjAwZwa+h5P9L6RwbfgwdRJHdd999yMnJwVNPPQUAOHr0KL75zW8iJycHffv2xdSpU3HmzBl85zvfUbVfRwcNJ3abuBG7PswTKXCYNb7EjJAR/FqywoaflbOxEtndzp07sXr16sDPv/71r9GtWzdUVlaiU6dOePLJJ7F06dL4ChpkPSeFDKdXNYJZMXDVzJBhpliXCjOIULyora3FxRdfHPj5nXfewa233opOnc5FhVtuuQULFixQvV/XzAxKRMbYeDLLspCh9HVlryQbLHiGVc4+Sm7m9Xpx4sSJwM8ffPABCgsLAz8nJCSgtbVV9X5Z0SDNnFTN8DOqqqHlRGy36oobKhZmzB7KMR/kVldffTWeffZZvPDCC3j11VfR2NiI6667LvD43/72N+Tm5qreL4MGxR01YcPoqyuCmR08nBIstIRDs6Yr52W05Cbz58/H2LFjsWbNGrS1teHRRx/FBRdcEHh83bp1GD16tOr9MmiQJk6sZgRrf/Kyw0lX60yZWl4jXpg1GRoDB7lBQUEBPv30U+zYsQM+ny+k2wQA7rjjDgwcOFD1fhk0KG7Z9aQrO3DY9X2qIavLK9Y4Dr1BhIGDnC4rKwsTJkwI+9j48eM17ZODQYlsSs8gTP9z3RAyzOQfVBp804IDR8mJbrzxRtTX1wd+XrhwYcjg0OPHj2uqaDBokGpO7zZxGjWBwc3hwqr3pSd4MGyQk2zevDnkqpKf/vSn+OqrrwI/t7W1oaqqSvV+2XVC5BBuDRBq2GEuFLUDTdmdQk4hhIj6s1aODRqcFdQarGaQ1YIDl5WhQ0vgYNigeOTYoEFEZPUlwoC6wMHqBtlZQkICEhISOtynF4MGKcZqBtmdldUOtYGDYYPsRgiBqVOnwuPxAABaWlpw7733okuXLgCgaVZQwMFB46a8SnafEFFEVoUOpYu8MWyQ3UyePDmkgvGf//mfYbdRy7FBg8zFagY5mZKBtDLDiNLqBsMG2cmcOXNw8cUXIzFR7gWpvLyVYmLIoHgQPPeIrHlIlFwOy0tgyS769++PY8fOB+7bb78ddXV1uvfLoEFEFIXewMGwQU7R/nLWN998E83Nzbr3y64TiorVDKJz9Iz5UNKVwm4UcitWNCgihgyi8LRWObROaU5kBl7eGgavPDEOQ4a7Rfu94f+9crJnKmVVg6wU6/JWv1dffVXVfh0dNEg+nmT0URJ8zTzGWoJ4uOfwcxGZ2tV2Y13+yrBBVpkyZUrIz+Eub9XC8UGDVQ15eDJRR+vnzsgwYtTvQvB++TkJT011Q+lcG0RmWrVqlSH7dXzQIP144ojN7DBr5/Dcvm38/JwnK2ywqkFu4oqgwaqGNjxBnMPPjj5qjl88fObssMIskZ24ImgADBtqxcMf/HD4GbGWjOPvhM+u0rDBLhSKB64JGkThMFi4j1O6bvRWNth9Qm7hqnk07PoHh8z1xsGhgRu5n53/r5XMtcG5NcjtXFfRYBdK/OH/NwHnPwf8wkFkL64LGgDDRjzg/y9FYrfAwcGhFO9cGTSA839keEIK742DQy35Q8z/DzKL3QJHNJEGhXKcBrW3fft2LFq0CLt378aRI0ewYcMGTJw4MfD41KlT8eKLL4Y8p6SkBJs2bTK5pee5aoxGOE74I2MVzg1B8YBjdshNmpubMWTIECxdujTiNuPGjcORI0cCt9/+9rcmtrAj11Y0grErJTKzvvXx+JMdWDXnB7tPSJbS0lKUlpZG3cbj8cDn85nUothcX9EgZRgEiEKxCkJmamhoCLm1trZq3tfWrVvRo0cP5Ofn4/vf/z6OHz8usaXqxUVFA2BVQwmjqhs87uRkZo1n4jgN59l0aDCSu6To2sep5hYAryE3Nzfk/vLycsydO1f1/saNG4dbb70Vffr0wf79+/Hoo4+itLQUFRUVSEpK0tVWreImaJByThpER2QGqwZPU/yoqamB13s+aPqXalfrjjvuCPx78ODBKCgoQN++fbF161aMHTtWdzu1iKuuE/6hUIelYyJ5lEzeFYlvVCp8o1Iltobsxuv1hty0Bo32LrnkEmRlZaG6ulrK/rSIq6ABMGxowcBBZE4XYKxZQhk4SK3Dhw/j+PHj6Nmzp2VtiLugQdoxcFC80/v5lzUluT9wMHTEn6amJlRWVqKyshIAcODAAVRWVuLQoUNoamrCQw89hJ07d+KLL77Ali1bMGHCBPTr1w8lJSWWtTkugwarGvowbBAZS836Jwwd8WXXrl0YNmwYhg0bBgAoKyvDsGHDMGfOHCQlJWHv3r245ZZbcOmll2LatGkYPnw4/vznP0vritEibgeD8ioUfTg4juKV3s++miXkAahaRt4fNniVinuNGTMGQoiIj2/evNnE1igTlxUNP54oiUgLM7pQ/JryG1Sv8MoKB9lJXAcNwBlh46a8Ske0MxxWjcitzAwbwPnAobZbhchqcR80AHuHjeC2OTlwELmR3gHSWi95VRM4GDbIagwa/2bHE3ikNtmxrUTxTG/YMDpwsCuFrMSgEcROJ/BYbbG6rVa/PpHdWFXdANQFDiKzMWi046QTqN27Ujg+g+KRVdUNQPkcHERmYtAIw+oTuNrXNrutdg438aj6cEaHG1lLxkBRI7tTGDbITHE7j4YSaubaiHbyVfNHR+tJ3P88o6sIDBn2ECtMRHu8X69Gaa8jg5r2OImMuWaCw4aSuTeCRVoNlshsjgoaJSnH4U3tpKu0qFb7PxRaVjYN3jZSEJB1AlfyWjL2TdaQceK3W8WjfXvcFDxkTmzn/7unJnBECxtcfp7M4qig4Tc+9ZipYSOY3j8aZp6slbxWrDCitb0cnyGf3QKCUfzv0y2BQ/YsumoDB8MGWc2RQQMI/SWzKnS4AasUzhAvISOYmwKHEVP2qwkcDBtkJVcMBh2fekx1/yWRU8RjyAjmlvdvVJVP6RetaANEOc8GGckVQcOPgcM+WCmRwy0nWb3ccjWNncMGwJVgyRiagsbSpUtx8cUXIyUlBYWFhfjggw8ibrt69WokJCSE3FJSUjQ3WAl/4GDoIHIXho3IZIUNPwYOkkV10HjppZdQVlaG8vJyfPjhhxgyZAhKSkpw9OjRiM/xer04cuRI4Hbw4EFdjVaDgYOcyg0nVSO4obphh7DBwEFmUR00nn76aUyfPh133303Bg4ciOXLlyMtLQ0rV66M+JyEhAT4fL7ALTs7W1ejtWDgMB+7T8hITg8ceqcsj0TtEvRKMWyQVqqCxqlTp7B7924UFxef30FiIoqLi1FRURHxeU1NTcjLy0Nubi4mTJiAjz/+OOrrtLa2oqGhIeQmCwOHuRg2tHHyCdRsTj9W/sAhM3SoDRusbpCRVAWNY8eO4cyZMx0qEtnZ2aitrQ37nPz8fKxcuRKvv/461qxZg7Nnz2LkyJE4fPhwxNdZsGABMjMzA7fc3FwAQPrnGdJmumPgIHIPp4cNP5mhQ+1l/1x6noxi+DwaRUVFKCoqCvw8cuRIDBgwAM8//zzmz58f9jmzZ89GWVlZ4OeGhoZA2AAQEjbUlP7CceJ8HNECkh3fg5qp3Im0qj6c4Yo5N/y0zELc3saTWZqmLgcQ80udP2xwDg6KRVXQyMrKQlJSEurq6kLur6urg8/nU7SPzp07Y9iwYaiuro64jcfjgcfjUbQ//y+D3sABnD+B2/FkDSibmCfcNnZ4PwwbRNroDRxawgagPnC0xwBCfqqCRnJyMoYPH44tW7Zg4sSJAICzZ89iy5YtmDlzpqJ9nDlzBvv27cONN96ourHRGFXlCKZmJj6ZJ3e9XTx2CR8MG2Q0t1U1gumZXVTLOil+Whdni9a9whASX1R3nZSVlWHKlCm48sorMWLECCxZsgTNzc24++67AQCTJ0/GRRddhAULFgAA5s2bh6uvvhr9+vXDiRMnsGjRIhw8eBDf/e535b6TIDKrHMHU/JK231btid3o8SN620dE5tM7lbnWwCF7JdjgEMLQ4X6qg8btt9+Of/7zn5gzZw5qa2sxdOhQbNq0KTBA9NChQ0hMPD/G9F//+hemT5+O2tpaXHDBBRg+fDjee+89DBw4UN67iEBmlUOvaFUFOwxKVdsGrcGEVQ0ympurGrJoXQkWiN2VohbHerhfghBCWN2IWBoaGpCZmYm/PXgdMjznspHWUc9WBw430Ro2GDRic8tVFFZxe9CQfdm4rC87MkKIjMDR2NqGS596B/X19fB65QYj4Pw56bZXHkdyF30zXZ9qbsHL/2++YW21A8eu3hr8YVQTOuxU5XC68anHVIcNhgwi/axeej6SSH9T1QQQVjjcxxWLqtVuP6npQ5le5ZVeBrSC/3245f0QUWxWzyqqhn+ODs5EGp9cETT89AYOJ52ko7WZwYPiXbx0PTkpbPhxYrD449iuk2i0dqsA+voYje6K0dq29s+zqsuI3SbK9evVGDcnS9JHdjcKYM5AdTXzdLAbxdlcGTSC6QkdaoX7hZFxUpddlTCqnURkDRmziIajdbIvNZQEDoYNZ3N90Ajm/6CaWY5TW02wqqvD6PDBagaR8YJ/z2SFDjPCBhB7rg6GDeeKq6DhF+7Dalb4MDJIyK7epFd5o4YNTvJFZF8yqxzBv+uRJvuTEUYYNtwpLoNGOFo+vHYYqBSp3WZ2GZEx/HNBcKwG6SG7WyXSF4xI98uehZRhw3kYNHRo/2E384Su5hdNybYMI/bFwEEyGDWOIxats5AybLgHg4ZERlcRjPzFYgXE/mLNdskgQkoYcZWKEmoDB8OGezBoGIS/AGS29kGEwYMisSpsAOoCh5KwAfDvrd25asIuIjqvX6/GwI2oPauvBFM6mFzJ1W++UamBG9kPgwZRHIi30BEv71MvN4UNP9+oVGSP1LfQGcnFoEEUZ9weOtz6vozipLDBiQWdiWM0iOJYuJOyU8d2MGBoZ+WYDUDdpGBKpy4n+2DQIM3ULBN/U16l5d+cSBmtJ2yzAwqDhVxOChtAaHcKQ4e9MWhQVLFmByXyi3bibx9CGBLsyWlhw49VDntj0CDTsKoRvxgsnMOI9VLU0LO2SqzLYckaDBqki5ruE4Bhg8hJtPyuylxbRWt1o6mxTXcb7Gr79u1YtGgRdu/ejSNHjmDDhg2YOHFi4HEhBMrLy/HCCy/gxIkTuOaaa7Bs2TL079/fsjbzqhMynZWlWSIy1hsHh0r7MrHxZBYXb2ynubkZQ4YMwdKlS8M+/vOf/xzPPvssli9fjvfffx9dunRBSUkJWlpaTG7peQwaFFOsUqSWbx0MG0TuJrNyycBxXmlpKX784x/jW9/6VofHhBBYsmQJHnvsMUyYMAEFBQX49a9/jS+//BKvvfaa+Y39NwYN6sCs6Xxvyqtk4CByMdndpP7A4cbQ0dDQEHJrbW1VvY8DBw6gtrYWxcXFgfsyMzNRWFiIiooKmc1VhWM0SJFYV5+oHasRjOM2iNzLqCtZ9AwaleXvX6ajU6q+ac/bTp47Defm5obcX15ejrlz56raV21tLQAgOzs75P7s7OzAY1Zg0CBp9IYNQPs3IKV/yBhoiMxnZNgAtHXf2k1NTQ283vPd1B6Px8LWyMWuE1LMjMvG1HanGL09EclhZMh3Q1eK1+sNuWkJGj6fDwBQV1cXcn9dXV3gMSuwokFS6alqBDM6DOitoBCRev7fN7d2pVitT58+8Pl82LJlC4YOHQrg3NiP999/H9///vctaxcrGhRWpAGhSqoaTvplZ3WDyHwyL4EN5obKRixNTU2orKxEZWUlgHMDQCsrK3Ho0CEkJCTggQcewI9//GP8/ve/x759+zB58mTk5OSEzLVhNlY0SDUl05LLqmyYgYNRiaxhxCykG09m4euW01L2ZUe7du3CtddeG/i5rKwMADBlyhSsXr0aDz/8MJqbm3HPPffgxIkT+MY3voFNmzYhJSXFqiYjQQghLHt1hRoaGpCZmYm/PXgdMjzMRmbyjYo8olrJGihOCRt+DBxE1tMbOr5uOo3/GvEa6uvrQwZYyuI/J1357M8kXHVyErvuf8SwttoBz9oUVe32kxHDhtLKBuCcwBH8B46hg8gaVq+3QnJxjAbpovRKFCeN2/DzX6HCP3RE1mHgdz5WNCimaFUNQPlS8k4at9Fe+7DBP35E5jHyahUyHoMGKSIzbADO6UqJhJfHul/14QxV2/fr1WhQS8jPqIm/yFgMGqSYrLABOLu6EYxjOpxHbYAwe78MLNExbDgPgwapIjtsAM6vbvjxMll7MypgyKannfESUhg2nIWDQUm1WKu7qp2qfHzqMemDRY3YpxIcPGpPTgkZelUfzuhwcyuGeudgRYM0kVnZ8NNT4YgUKqyqmrC6QXYRLWw4vQLCyoYzMGiQZkrCBqBsYq9gRlQirAgcHDBqD27+Vq+XG7ppGDbsj10npEusbhTAnFVflbKiS4V/BMmN7NRFY9TaKSQHKxqkW6zKBqC9umGU4LBhRpWD1Q1rWH0CjCftj7UVFQ9WN+yJQYOk8Fc2lAQOu4QNPzNDBwMHxYvg4GFm6Hjj4FCcam4B8Jppr0nRMWiQVE6sbgQL161iRPgI962L4YPcyh867DKug8zFoEHSKQkbgL0DRzCzKh6RSr5KAki0cjEDDNkFA0d8YtAgQyjtSgE6DhbVEzyMDi9WXr2i9/kMHGQXDBzxhUGDDKW0uhFMxlUqMsNLOE6c1ZSBg+ym+nAGw0Yc4OWtZLja7ScVXQZrpPQqryGX2Vox+6heHJVPdmKHy2PJWAwaZBo7BQ6ZocOq6c714FTpZDcMG+7FrhMyXaywobarRSvZ3StOXJGW3SlkJxy74U4MGmQ7eqseWoNKuCqHmeu1WMmNy93zG7JzMXC4C4MGuY6aK15i0btei9MCB8AqB9mHVZN+kVwMGuRadgocgPNChxurHORcrHI4F4MGuV77rhg9wUPPPB1mzTpqBFY5yC5Y5XAeBg2KO5HGgKgJILLWbGkfPuwePFjlIDvhPBzOwKBB9G9qKx9GzELqpOBh5yoHB4LGD3ap2B+DBlEEalakBYyZ9lzW/BxmrtFix+BB7sfqhn0xaBDF4IZF4sy8CkbJRGAMI2QEhg17YtAgUkDrInF2Cx12uezWzt0u5GzVhzPQ+4IWq5tBQTgFOZEKaicTC57y3Ii1VrTyT5tu9dTpnAadjPD3L9OtbgIFYUWDSCU983MYvaqsFlbP9XFTXiUrG0QuxqBBpJGMCcHsFjysCh0MG0TuxaBBpFO47hRZ661YGTysrnQQkTswaBAZQMakYED4hd7UkjmxGAMHEanFoEFkIpnrrygl8yoYI6sc7D4hcicGDSILBFc8rAgdbp8+vV+vRs4OSmQTDBpEFlN7yWwwO40FiXaprJ1CCBGZi0GDyMGMGAti1FTqZocNVjWI7IETdhG5UO32kyE3NYyaYEzJ5GCcwIvIfVjRIIoDWseE2Hn9FiJyBk0VjaVLl+Liiy9GSkoKCgsL8cEHH0Tdfv369bjsssuQkpKCwYMH480339TUWCLST0u1w8xp1GVWNbjAFpH1VAeNl156CWVlZSgvL8eHH36IIUOGoKSkBEePHg27/XvvvYc777wT06ZNw549ezBx4kRMnDgRH330ke7GE5F+artX9KzfYvXaKkRkPtVB4+mnn8b06dNx9913Y+DAgVi+fDnS0tKwcuXKsNs/88wzGDduHB566CEMGDAA8+fPxxVXXIFf/vKXuhtPRPJoGc8BdAweMqoeHKtB5B6qgsapU6ewe/duFBcXn99BYiKKi4tRUVER9jkVFRUh2wNASUlJxO0BoLW1FQ0NDSE3IjKHnoGkftHChtlVDXafkJvMnTsXCQkJIbfLLrvM6mZFpWow6LFjx3DmzBlkZ2eH3J+dnY3PPvss7HNqa2vDbl9bWxvxdRYsWIAnnnhCTdOIyCBqwobMycc4UyhReIMGDcLbb78d+LlTJ3tf12HLy1tnz56N+vr6wK2mpsbqJhGRSjKqGrK6UFjVILtrX8VvbW2NuG2nTp3g8/kCt6wse0+IpyoGZWVlISkpCXV1dSH319XVwefzhX2Oz+dTtT0AeDweeDweNU0jIhuo3X5ScVVD6SResiobnMCLZPvXzhNISm7RtY8zp849Pzc3N+T+8vJyzJ07N+xzPv/8c+Tk5CAlJQVFRUVYsGABevfurasdRlJV0UhOTsbw4cOxZcuWwH1nz57Fli1bUFRUFPY5RUVFIdsDwB//+MeI2xORe5hxOawarGyQXdXU1IRU8mfPnh12u8LCQqxevRqbNm3CsmXLcODAAXzzm99EY6N9P9uqO3bKysowZcoUXHnllRgxYgSWLFmC5uZm3H333QCAyZMn46KLLsKCBQsAAD/4wQ8wevRoPPXUUxg/fjzWrVuHXbt24X/+53/kvhMisqX0Km/ECb/MrmoA58MGqxtkJ16vF15v7GBeWloa+HdBQQEKCwuRl5eHl19+GdOmTTOyiZqpDhq33347/vnPf2LOnDmora3F0KFDsWnTpsCAz0OHDiEx8XyhZOTIkVi7di0ee+wxPProo+jfvz9ee+01XH755fLeBRE5lhVhA2DgcKt+vRpxqrkFu6xuiEm6du2KSy+9FNXV1VY3JSJNQ1VnzpyJmTNnhn1s69atHe779re/jW9/+9taXoqIHCbcOI1oVQ3AurABhHanMHQ4Wzx2jTU1NWH//v34zne+Y3VTIrLlVSdE5D6yxmsYOZlXv16NITciu/nhD3+Ibdu24YsvvsB7772Hb33rW0hKSsKdd95pddMiYtAgIukizb1hp4m8lGDYcI54+b86fPgw7rzzTuTn5+O2225D9+7dsXPnTlx44YVWNy0ie8/yQURxxcoulEg4lsP+4iVkAMC6deusboJqrGgQkSG0VDUAdZN5mbkmSjydzJyE/y/2x6BBRKaTOb8Gw0b84v+HMzBoEJHtqB2vwbARf/j/4BwMGkRkS1rChlmBg1elWIvH3lkYNIjIVVjdcDcec+fhVSdE5DrBYcPoq1N4VYp5GDKciUGDiFzNHzoYOJyLAcPZGDSIyHTRpiP3UzKfhhrtu1SMCh7tT4oMHtoxYLgDgwYRGaL9eid2Y3alIxjDR2QMF+7DoEFEcc3MWUb9uJBbKIYLd2PQICJTWdFtEosVYcMv3rpaGCriD4MGERGsDRvBop2InRZCGCoIYNAgIgqwS9iIROmJ2+xAwkBB0TBoEJGtmN1t0p7dw4YSsU784YIIwwIZhUGDiKgdN4SNaBgqyEycgpyIKAwzpzIncjMGDSKiCBg2nOemvEqM673P6mZQEAYNIrIVtau2Go1hwxnMXL2X1OEYDSIHUDPLZu32kwa2JD65fcyGkzFc2B+DBpEN6Zm+u/1zGTzkYNiwFwYM52DQILIRI9YHCd6nHUJHepVX0eygdmTW+igUHsOFMzFoEFnMzMXH/K9lh8ARzfjUY5bPpxENA4d5GC6cj0GDyEJWrXBqRuCo3X7S9iu46sXAIR+DhfswaBBZwC4nYN+oVEuqG0q6T+xe1QgWfHJk6IiOQSL+MGgQmUxGyGh/kk6v8mrel1Vhw62inUjjJYQwTFAwBg0iE2kJGUoGTobbRk34MKorRW/3iZOqGkrEOgE7JYgwSJAaDBpEJlFzwpVxVYZ/H2oDh1nVDaVXn7gtbETDEzi5EWcGJTKB0pDRlN8g/dJPtfuUPX5ERnCx22yhZF/jU4+hJOW41c2gIAwaRAZTEzKMpCZwmDVYVc/YEqJg41OPMZDaFIMGkYGUnLCNqGLEej2zRatqKA0bPIlQOAwY9segQWQQpSHDCkpe1y6X4AbjCYWA8+GCnwdnYNAgsojV03CbHTZkVDUAho14FBws+P/vPLzqhMgAsU7QVocMv6b8BtuMk1CzBko8XYkSLxgg3ItBgyjOxQobMi95lTktuf/ExMBhfwwR8Y1Bg0gyo6sZkU6sev6Y26WyoWVlVxnVjfbHjuFFP4YL8mPQIDKR1pCh5MTn38aIP/BmVjW0hg1AXUCIdpyCH7N76ND7/y3z/dkhXKRXeXG2+bTVzaAgDBpEEsm+UkPLSUBr4DCzqmFE2ACMOdHZpdph1Ek83H7VvkerA4YdqnEUGYMGkUnUnjj1ntA2nsyy/ASgh9awYbRIxzT4/8vJxx1wRvsZLpyDQYPIhmR9a1YbNqJVNWSvg6JkYKhdw0Y4Tjg5OxVDhbMxaBBJEu2kaeXJ0s6VDbeFDdKGQcLdGDSIbMbugw9lY9iILwwV8YdBgygO2LmqASgPG4B9JjsjdYwOGMHdeo2tbYa+FqnDoEFkI3aoZpg5TiOY0sm8WN2ILtYJ3exjJztgGPX5c5KlS5di0aJFqK2txZAhQ/CLX/wCI0aMsLpZETFoEJFtqAkbfvEcOrScxM06djICBkNFRy+99BLKysqwfPlyFBYWYsmSJSgpKUFVVRV69OhhdfPCckTQEEIAAJpYDiMbS4sySVBTo7LP7tctxk001NCm/Pcn2oRHRpelG//YiOyRKcqf8OH5YNLUv1HVa6V/nhH2frX7MVqkdjZA5+dFx7GLxN9WtW2re69FyusD588V/nOHUc6ebpW2j4aG0NDn8Xjg8Xg6bP/0009j+vTpuPvuuwEAy5cvx8aNG7Fy5Ur86Ec/0t0eQwgHqKmpEQB444033njjTfGtpqbGkHPSyZMnhc/nk9bO9PT0DveVl5d3eN3W1laRlJQkNmzYEHL/5MmTxS233GLIe5XBERWNnJwc1NTUICMjAwkJCYqe09DQgNzcXNTU1MDrtf8oZ6e1F3Bem53WXsB5bXZaewHntdlp7QXMb7MQAo2NjcjJyTFk/ykpKThw4ABOnTolZX9CiA7ntnDVjGPHjuHMmTPIzs4OuT87OxufffaZlLYYwRFBIzExEb169dL0XK/X65hfRsB57QWc12antRdwXpud1l7AeW12WnsBc9ucmZlp6P5TUlKQkqKiiy+OJVrdACIiIootKysLSUlJqKurC7m/rq4OPp/PolbFxqBBRETkAMnJyRg+fDi2bNkSuO/s2bPYsmULioqKLGxZdI7oOtHC4/GgvLw8bD+XHTmtvYDz2uy09gLOa7PT2gs4r81Oay/gzDbbVVlZGaZMmYIrr7wSI0aMwJIlS9Dc3By4CsWOEoQw+PofIiIikuaXv/xlYMKuoUOH4tlnn0VhYaHVzYqIQYOIiIgMwzEaREREZBgGDSIiIjIMgwYREREZhkGDiIiIDOPYoPGTn/wEI0eORFpaGrp27Rpz+9OnT+ORRx7B4MGD0aVLF+Tk5GDy5Mn48ssvQ7a7+OKLkZCQEHJbuHCh6e0Fzk1LO2fOHPTs2ROpqakoLi7G559/HrLNV199hf/4j/+A1+tF165dMW3aNDQ1Nelur5Z9f/HFFx2Onf+2fv36wHbhHl+3bp0lbQaAMWPGdGjPvffeG7LNoUOHMH78eKSlpaFHjx546KGH0KZikTJZ7f3qq69w3333IT8/H6mpqejduzfuv/9+1NfXh2wn8xgvXboUF198MVJSUlBYWIgPPvgg6vbr16/HZZddhpSUFAwePBhvvvlmyONKPtd6qGnvCy+8gG9+85u44IILcMEFF6C4uLjD9lOnTu1wLMeNGyetvWrbvHr16g7taT9DpdHHWG2bw/2OJSQkYPz48YFtzDjOZBHLVlnRac6cOeLpp58WZWVlIjMzM+b2J06cEMXFxeKll14Sn332maioqBAjRowQw4cPD9kuLy9PzJs3Txw5ciRwa2pqMr29QgixcOFCkZmZKV577TXx17/+Vdxyyy2iT58+4uTJk4Ftxo0bJ4YMGSJ27twp/vznP4t+/fqJO++8U3d7tey7ra0t5LgdOXJEPPHEEyI9PV00NjYGtgMgVq1aFbJd8Hsys81CCDF69Ggxffr0kPbU19eHvK/LL79cFBcXiz179og333xTZGVlidmzZ5ve3n379olbb71V/P73vxfV1dViy5Yton///mLSpEkh28k6xuvWrRPJycli5cqV4uOPPxbTp08XXbt2FXV1dWG337Fjh0hKShI///nPxSeffCIee+wx0blzZ7Fv377ANko+11qpbe9dd90lli5dKvbs2SM+/fRTMXXqVJGZmSkOHz4c2GbKlCli3LhxIcfyq6++0t1WrW1etWqV8Hq9Ie2pra0N2cbIY6ylzcePHw9p70cffSSSkpLEqlWrAtsYfZzJOo4NGn6rVq1SfOJu74MPPhAAxMGDBwP35eXlicWLF8tpXBhK23v27Fnh8/nEokWLAvedOHFCeDwe8dvf/lYIIcQnn3wiAIi//OUvgW3eeustkZCQIP7xj3/oaqesfQ8dOlT813/9V8h9ADqsPiiD1jaPHj1a/OAHP4j4+JtvvikSExND/pgvW7ZMeL1e0draanp723v55ZdFcnKyOH36dOA+Wcd4xIgRYsaMGYGfz5w5I3JycsSCBQvCbn/bbbeJ8ePHh9xXWFgovve97wkhlH2uzWxve21tbSIjI0O8+OKLgfumTJkiJkyYoLttkahtc6y/IUYfYy1tbm/x4sUiIyMj5Euc0ceZrOPYrhMZ6uvrkZCQ0KErY+HChejevTuGDRuGRYsWSSmRq3XgwAHU1taiuLg4cF9mZiYKCwtRUVEBAKioqEDXrl1x5ZVXBrYpLi5GYmIi3n//fV2vL2Pfu3fvRmVlJaZNm9bhsRkzZiArKwsjRozAypUrISRM56Knzb/5zW+QlZWFyy+/HLNnz8bXX38dst/BgweHrJhYUlKChoYGfPzxx5a0N1h9fT28Xi86dQqd6FfvMT516hR2794d8hlMTExEcXFx4DMY7j0Fbw+cO1b+7ZV8rrXS0t72vv76a5w+fRrdunULuX/r1q3o0aMH8vPz8f3vfx/Hjx/X1Va9bW5qakJeXh5yc3MxYcKEkM+hkcdYT5uDrVixAnfccQe6dOkScr9Rx5ms5dopyGNpaWnBI488gjvvvDNkNcH7778fV1xxBbp164b33nsPs2fPxpEjR/D000+b2r7a2loACLscsP+x2tpa9OjRI+TxTp06oVu3boFt9Ly+3n2vWLECAwYMwMiRI0PunzdvHq677jqkpaXhD3/4A/77v/8bTU1NuP/++y1p81133YW8vDzk5ORg7969eOSRR1BVVYVXX301sN9w/w/+x8xub7Bjx45h/vz5uOeee0Lul3GMtSxJHelYBX9m/fdF2kYrGUtoP/LII8jJyQk5iY4bNw633nor+vTpg/379+PRRx9FaWkpKioqkJSUZHqb8/PzsXLlShQUFKC+vh5PPvkkRo4ciY8//hi9evUy9BhrbXOwDz74AB999BFWrFgRcr+Rx5msZaug8aMf/Qg/+9nPom7z6aef4rLLLtP1OqdPn8Ztt90GIQSWLVsW8lhZWVng3wUFBUhOTsb3vvc9LFiwoMM8/Wa1Vyalbdbr5MmTWLt2LR5//PEOjwXfN2zYMDQ3N2PRokURT4JGtzn4JD148GD07NkTY8eOxf79+9G3b1/V+zPrGDc0NGD8+PEYOHAg5s6dG/KY2mNM5yqZ69atw9atW0MGV95xxx2Bfw8ePBgFBQXo27cvtm7dirFjx5rezqKiopAFtEaOHIkBAwbg+eefx/z5801vj1orVqzA4MGDMWLEiJD77XacSR5bBY0HH3wQU6dOjbrNJZdcous1/CHj4MGDeOedd0KqGeEUFhaira0NX3zxBfLz801rr3/J37q6OvTs2TNwf11dHYYOHRrY5ujRoyHPa2trw1dffRVxyWClbday72CvvPIKvv76a0yePDnmtoWFhZg/fz5aW1vDLrpkVpuD2wMA1dXV6Nu3L3w+X4cR9f5lmsPt14z2NjY2Yty4ccjIyMCGDRvQuXPnmO8p2jEOR8uS1D6fL+r2Sj7XWulZQvvJJ5/EwoUL8fbbb6OgoCDqtpdccgmysrJQXV2t+wQoY9nvzp07Y9iwYaiurgZg7DHW2+bm5masW7cO8+bNi/k6Mo8zWcziMSK6qRkMeurUKTFx4kQxaNAgcfToUUXPWbNmjUhMTJQ2+lntYNAnn3wycF99fX3YwaC7du0KbLN582apg0G17nv06NEdroSI5Mc//rG44IILNLfVT9bxePfddwUA8de//lUIcX4waPCI+ueff154vV7R0tJienvr6+vF1VdfLUaPHi2am5sVvZbWYzxixAgxc+bMwM9nzpwRF110UdTBoDfddFPIfUVFRR0Gg0b7XOuhtr1CCPGzn/1MeL1eUVFRoeg1ampqREJCgnj99dd1t1cIbW0O1tbWJvLz88WsWbOEEMYfYz1tXrVqlfB4POLYsWMxX0P2cSbrODZoHDx4UOzZsydw+eSePXvEnj17Qi6jzM/PF6+++qoQ4lzIuOWWW0SvXr1EZWVlyCVU/isH3nvvPbF48WJRWVkp9u/fL9asWSMuvPBCMXnyZNPbK8S5S9S6du0qXn/9dbF3714xYcKEsJe3Dhs2TLz//vvi3XffFf3795d6eWu0fR8+fFjk5+eL999/P+R5n3/+uUhISBBvvfVWh33+/ve/Fy+88ILYt2+f+Pzzz8Vzzz0n0tLSxJw5cyxpc3V1tZg3b57YtWuXOHDggHj99dfFJZdcIkaNGhV4jv/y1htuuEFUVlaKTZs2iQsvvFDa5a1q2ltfXy8KCwvF4MGDRXV1dcjnuK2tTQgh9xivW7dOeDwesXr1avHJJ5+Ie+65R3Tt2jVwBc53vvMd8aMf/Siw/Y4dO0SnTp3Ek08+KT799FNRXl4e9vLWWJ9rrdS2d+HChSI5OVm88sorIcfS/3vZ2NgofvjDH4qKigpx4MAB8fbbb4srrrhC9O/fX1fI1NPmJ554QmzevFns379f7N69W9xxxx0iJSVFfPzxxyHvy6hjrKXNft/4xjfE7bff3uF+M44zWcexQWPKlCkCQIfbn/70p8A2+PdcAkIIceDAgbDbBz9n9+7dorCwUGRmZoqUlBQxYMAA8dOf/lTKB11te4U4983k8ccfF9nZ2cLj8YixY8eKqqqqkP0eP35c3HnnnSI9PV14vV5x9913h4QXPWLt239Mg9+DEELMnj1b5ObmijNnznTY51tvvSWGDh0q0tPTRZcuXcSQIUPE8uXLw25rRpsPHTokRo0aJbp16yY8Ho/o16+feOihh0Lm0RBCiC+++EKUlpaK1NRUkZWVJR588MGQy0nNau+f/vSniJ/jAwcOCCHkH+Nf/OIXonfv3iI5OVmMGDFC7Ny5M/DY6NGjxZQpU0K2f/nll8Wll14qkpOTxaBBg8TGjRtDHlfyudZDTXvz8vLCHsvy8nIhhBBff/21uOGGG8SFF14oOnfuLPLy8sT06dM7zFthZpsfeOCBwLbZ2dnixhtvFB9++GHI/ow+xmrbLIQQn332mQAg/vCHP3TYl1nHmazBZeKJiIjIMHE9jwYREREZi0GDiIiIDMOgQURERIZh0CAiIiLDMGgQERGRYRg0iIiIyDAMGkRERGQYBg0iIiIyDIMGERERGYZBg4iIiAzDoEFERESG+f+kujjDc0PV7wAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cv2d = colvar[['p.x', 'p.y']].values\n", "\n", "# or equivalently\n", "# cv1 = colvar['p.x'].values.squeeze() \n", "# cv2 = colvar['p.y'].values.squeeze()\n", "# cv2d = np.stack((cv1, cv2)).transpose()\n", "\n", "fes_params = {\n", " 'blocks': 1, # Number of blocks for error analysis\n", " 'bandwidth': 0.015, # Kernel bandwidth (sigma) for density estimation\n", " 'scale_by': 'range', # Method to scale the bandwidth. \n", " 'kbt': kbt, # Thermal energy for proper energy scaling\n", " 'weights': weights, # Statistical weights from the bias\n", "}\n", "\n", "fes2, grid2, bounds2, error2 = compute_fes(\n", " cv2d,\n", " plot=True,\n", " plot_max_fes=40,\n", " **fes_params\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "mlcolvar111", "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.10.14" } }, "nbformat": 4, "nbformat_minor": 2 }