pydhamed/cg-rna/cg_RNA_duplex_formation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    },
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1 toc-item\"><a href=\"#Free-energy-profile-for-RNA-duplex-formation\" data-toc-modified-id=\"Free-energy-profile-for-RNA-duplex-formation-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Free energy profile for RNA duplex formation</a></div><div class=\"lev1 toc-item\"><a href=\"#Prepare-calculation\" data-toc-modified-id=\"Prepare-calculation-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Prepare calculation</a></div><div class=\"lev2 toc-item\"><a href=\"#Transition-counts\" data-toc-modified-id=\"Transition-counts-21\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Transition counts</a></div><div class=\"lev2 toc-item\"><a href=\"#Load-bias-potentials\" data-toc-modified-id=\"Load-bias-potentials-22\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Load bias potentials</a></div><div class=\"lev1 toc-item\"><a href=\"#Solve-DHAMed-equations\" data-toc-modified-id=\"Solve-DHAMed-equations-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Solve DHAMed equations</a></div><div class=\"lev2 toc-item\"><a href=\"#Use-all-data-points\" data-toc-modified-id=\"Use-all-data-points-31\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Use all data points</a></div><div class=\"lev2 toc-item\"><a href=\"#Block-average\" data-toc-modified-id=\"Block-average-32\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Block average</a></div><div class=\"lev2 toc-item\"><a href=\"#Vary-observation-interval\" data-toc-modified-id=\"Vary-observation-interval-33\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Vary observation interval</a></div><div class=\"lev1 toc-item\"><a href=\"#References\" data-toc-modified-id=\"References-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>References</a></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from pydhamed import count_matrix\n",
    "from pydhamed import run_dhamed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Free energy profile for RNA duplex formation\n",
    "\n",
    "We have run coarse-grained simulations of RNA duplex formation with the oxRNA model (https://dna.physics.ox.ac.uk/index.php/RNA_model_introduction). \n",
    "\n",
    "# Prepare calculation\n",
    "\n",
    "## Transition counts\n",
    "\n",
    "Load trajectory data. Each time point needs to be assigned to a structural/conformational state. I use the number of native base-pairs to assign structural states for the coarse grained RNA. Second but lust column contains the state assignment. The state assigned trajectory is used to determine the transition counts, which serve as input for DHAMed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "traj = np.load(\"us_step-e4.npy\")[:,-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 8.,  7.,  8., ...,  0.,  0.,  0.])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "traj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "c_l = [count_matrix(traj,n_states=9)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "200001"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "traj.__len__()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "## Load bias potentials\n",
    "\n",
    "Bias potentials were used to sample unfavorable structures. For each structural state and for each window (or simulation run) we need to specify the bias that was acting on the simulations system. The bias is saved as $w_i^\\alpha = \\exp(-\\beta u_i^\\alpha)$. To get the correct form of the bias we take the log to obtain $u_i^\\alpha = -\\ln(w_i^\\alpha) \\, k_BT$. The bias is in units of $k_bT$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "v_ar = np.genfromtxt(\"wfile.txt\")[:,1].reshape((9,1))\n",
    "v_ar = -np.log(v_ar) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9, 1)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_ar.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "\n",
    "\n",
    "For the simulation of RNA duplex formation we have nine structural states and the simulation was run all at once (in a single simulation run or window). The bias array thus has the shape 9x1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Solve DHAMed equations\n",
    "\n",
    "Solve DHAMed equations for Umbrella Sampling simulations of coarse grained RNA duplex formation.\n",
    "\n",
    "## Use all data points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36\n",
      "loglike-start 307329.010043\n",
      "307329.010043\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 288165.228180\n",
      "         Iterations: 22\n",
      "         Function evaluations: 27\n",
      "         Gradient evaluations: 27\n",
      "time elapsed 0.690969944 s\n"
     ]
    }
   ],
   "source": [
    "og = run_dhamed(c_l, v_ar, numerical_gradients=False, jit_gradient=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36\n",
      "loglike-start 307329.010043\n",
      "307329.010043\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 288165.228180\n",
      "         Iterations: 23\n",
      "         Function evaluations: 300\n",
      "         Gradient evaluations: 30\n",
      "time elapsed 0.00719094276428 s\n"
     ]
    }
   ],
   "source": [
    "og = run_dhamed(c_l, v_ar, numerical_gradients=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "## Block average\n",
    "\n",
    "Estimate uncertainty in the DHAMed free energy profile by splitting the data into five blocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32\n",
      "loglike-start 59535.0856402\n",
      "59535.0856402\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 55154.290427\n",
      "         Iterations: 21\n",
      "         Function evaluations: 24\n",
      "         Gradient evaluations: 24\n",
      "time elapsed 0.00709700584412 s\n",
      "32\n",
      "loglike-start 54945.0688078\n",
      "54945.0688078\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 50890.763573\n",
      "         Iterations: 21\n",
      "         Function evaluations: 24\n",
      "         Gradient evaluations: 24\n",
      "time elapsed 0.00628900527954 s\n",
      "33\n",
      "loglike-start 30329.4839929\n",
      "30329.4839929\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 27634.820755\n",
      "         Iterations: 21\n",
      "         Function evaluations: 25\n",
      "         Gradient evaluations: 25\n",
      "time elapsed 0.00654888153076 s\n",
      "34\n",
      "loglike-start 52453.6162341\n",
      "52453.6162341\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 48500.807605\n",
      "         Iterations: 21\n",
      "         Function evaluations: 25\n",
      "         Gradient evaluations: 25\n",
      "time elapsed 0.00642395019531 s\n",
      "33\n",
      "loglike-start 53979.5557127\n",
      "53979.5557127\n",
      "Warning: Desired error not necessarily achieved due to precision loss.\n",
      "         Current function value: 49892.483234\n",
      "         Iterations: 20\n",
      "         Function evaluations: 85\n",
      "         Gradient evaluations: 73\n",
      "time elapsed 0.0166161060333 s\n"
     ]
    }
   ],
   "source": [
    "n_blocks = 5\n",
    "bl_l = []\n",
    "for bl in np.split(traj[1:],n_blocks):\n",
    "    c_l = [count_matrix(bl,n_states=9)]\n",
    "    bl_g = run_dhamed(c_l, v_ar)\n",
    "    bl_l.append(np.exp(-bl_g) / np.sum(np.exp(-bl_g)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "bl_ar = np.column_stack((bl_l))\n",
    "bl_ln_ar = np.log(bl_ar)\n",
    "\n",
    "# minmum and maximum values of the PMF (from the five reconstructions) to indicate the uncertainty.\n",
    "min_pmf = np.min( bl_ln_ar - bl_ln_ar[-1], axis=1)\n",
    "max_pmf = np.max( bl_ln_ar - bl_ln_ar[-1], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9, 5)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bl_ar.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9,)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dh_log_err = np.std(-np.log(bl_ar), axis=1) / np.sqrt(n_blocks-1)\n",
    "dh_log_err.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAADRCAYAAABreLGdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8W9WZ8PGfLFnWYnlf4i3OfrI7IXsChATowEuHMi2l\n0NJSKBQ6MHTaTuedd1oKLR06bad0p+uktCxtgSm0pZNSskM2spA468ni2Ikd75tkybIsS+8fUoyd\n1U5kX9l6vp9PPrHv1b338Y3z6Oicc59jCofDCCGEGF5JRgcghBCJSJKvEEIYQJKvEEIYQJKvEEIY\nQJKvEEIYQJKvEEIYwGJ0AIPV2Oi5rLlxmZkOWlt9sQ5nxJL70Z/cj3PJPenvcu5Hbq7LdKF9CdPy\ntVjMRocQV+R+9Cf341xyT/qL9f1ImOQrhBDxRJKvEEIYQJKvEEIYQJKvEEPs3Xd38eUv/99+21at\n+jl/+MPLAOzevZOHH36ARx75NPff/wleeunFc86xatXPufPOD/Loow/x8MMP8J3vfBOPx3PB8z/1\n1FfZuvXt3u/37dvLNdcs4PjxY73bVq9+nZtuWkEwGOzd5vF4WLFiCatXvz6gn62i4jj/9E8PDui1\noj9JvkIMA5PpgoPefPvbT/Hkk//Jj370c5555pds3Lielpbmc153xx138YMf/JQf//gXTJ8+g8cf\n//cBnR9gzZo3KC0dx5o1b/TbnpGR0S9Jb9iwlvz8MQP9sQZ0bXF+knyFMJjH48Hr9QJgtVr58Y9/\nQVZW9kWPufnm99PTE+Tgwf0AXKw6YSgUYuPG9Xzxi19izZq/9du3ePEy3nzzvYS8fv0a5s9f1Pv9\nq6++wmc+8ykeeeTT/P73LwDQ2NjAQw/dxxe/+Flef/21wf2woteIm+crxHDwB4LUNHkpynFisw7t\nf5P773+QBx74BHPnzmfBgkXceONNuFyuSx43adJkTpyooLCwiD17dvPoow8BkUR88mQVK1ZcD8CO\nHdsZN248ZWVzyMhI58CB/cyYMRMApaby298+R2dnJ35/J93d3WRnRxJ/dXU169ev4Sc/+W8APvOZ\n+1ix4gZeeeX33HDD+7j99jt54YVfc+zY0aG4LaOeJF8hzuIPBHny1zupbfZRkO3gsXvmD0kCPvNx\n/bbbbmf58pVs376VTZs28JvfrGLVqucv2foNBoMkJUU+vM6dO48nn/zP3n1PPfXV3q/ffPOv3HDD\n3wFwww1/x5o1b/QmX4AlS65m06b1dHR0cO21K+joiPQll5eXU119ikcffYhwOExnp5/a2loqK0+w\ncuWNvdfdvn1rDO5G4pFuByHOUtPkpbY58iRTbbOPmibvFZ0vIyMTj8fdb1tbWyvZ2TkAdHV1kZmZ\nxU033cJTT32bhQsX88472y553sOHDzFlytSLvqarq4vNm9/i5Zd/y333fYzXX/8jGzas7feaFStu\nYP36NWzcuI7rrlvZ24VhtVpZuvRqfvCDn/LDH/6MZ599kbKyOUCYpKTIG0coJIsxXC5JvkKcpSjH\nSUG2A4CCbAdFOc4rOl9JyViamhqpqakGoLW1ld27dzFrVhnV1af41Kc+TmdnJxDpn21qaqKgoOic\n8/Tt1/3jH/9AenoGEydOuui1N29+i3nzFvDrX/+OVate4LnnXmLs2HHs3r2z9zVTp06jrq6Onp4e\ncnPzerfPmDGD3bt30tXlJxwO8/3vf4dAIEBJSSmHDh0E4N13d55zTTEw0u0wjPyBIPuO1DCpJJvM\n9FSjwxEXYLNaeOye+THr87VYLHzlK0/yrW/9B+FwmHA4zOc//69kZmaSmZnJ3Xffw2c/+xnsdjuB\nQIBrrlkebWH298orv2PjxnUEAgHGji3liSe+fslrr137Bu9//wf6bbvllr9nzZo3mDWrrHfbwoWL\nyczM6ve6goIC7rjjLh5++NOYzWauuWY5VquVD3/4Lr7ylX9j06YNl0z+4sJMI20Nt8strJOb66Kx\n0RPrcAasbz+iy27m7uV5zJ05GYvFmPc/o+9HvJH7cS65J/1dzv2QwjpxoG8/oqezhxMtFtZs3Y8+\nduKi04SEEKOTJN9hUpTjJC8jBYB0p5XMNBs2Vy6nWmHdlr2crK4xOEIhxHCSPt9hYrNa+PCyHI6c\nbCMn3Qo9AbDYsCRbITmXQ9U+TlTvZdrEYvJyLz7FSAgx8knLdxh1BXqYWprBjVfPZXJBCoGOBnp6\nIs/Vp9gchG057D7SwJad+3uf2xdCjE6SfIdRc1sH2RmRWQ7jxhZzw9Iy8hwBOtsbe/t9bc40ApZM\n3t5Tyc7yQwQCASNDFkIMEUm+w6SlpYWuQIDCMfm925KSkpg1bTIrFk3FSTv+jvbefXZXFp4eF+u2\nHWTfoaOEQiEjwhZCDBFJvsOkpr6Z7Iy0804ts9lsLJw7g4UzCkjqaqKrM/JElclkwpaWS6MvhTc3\nl3O0onKYoxZCDBVJvsPE4wvgclov+pqszEyuXVTG9LGpBDsaCHZHuhySzGZSXLlUNoVYt2UPNafr\nhiNkIcQQkuQ7TNo7ushItQ3otSVFhVy/bA7F6WE63e/1BydbU0hy5LK/0s1b75TT3NIylCELMSp1\ndfl57rlfGR2GJN/h0NbWhtfXSXHhwItUm0wmpk6ZwA1LZpBm9tDpae3dl+Jw0mPN5p2DtWzbtZ8O\n75UVfhHG2759K6+99j+GxlBXV8v73recRx99iEce+TR33nknb721gbq6Wq65ZkFv7eAzPvWpj/dW\nT+t77D/904M8+uhD/PCHT1/wWrt377zo/ssRDofPOeezz/6St96KVIo7Q+vD/VbvgEjZzQcfvJdV\nq37eb/umTRv44AdvYePGdUCkktzPf/4ML774HKtXv84772zjhRd+fVnxyjzfYVBT10hGWio228Ba\nvn0lJyczb/Y0PB4P5Ycr8XRbsTkiMybsqRn4gbd2HycvPZnZ0yaRnJwc4+hHl3A4fE6FMQCrNYTb\nfXnT+1yutCtezWHRoiVXdHyslJaW8oMf/BSI3JNbb/0ATz/9QwoLi1iz5g2mT4+Uoqypqcbr7bjg\nsZdy5Mhh1q59k/vu+zRO53t1TqqqKtm0aT0f//i9g4rb7XazevWf2bPn3d5tO3e+A8A111zH0aNH\n2Lv3XXJycikpGcuuXTv6Hb9gwSKamhr5y1/+xH33fbp3e1eXn+LiEpYvXwnA17/+ODfeeBPLll0D\nRJZnOnmyalCxniHJdxi4vd24HBfv770Ul8vFsgWzOF1Xz+GKWkKWdCzWyDntrmzagyHWbT1ISb6T\nqZMn9NZ5Ff15PG7e3H4Mu6N/pbJUZwsd3q5Bn6/T5+XGRZNIS0u/6OtWr36dd9/dRXt7G5WVJ3jg\ngc+wZs0bVFZW8vjjT1JZeYKKiuNMmDCR8vI9tLW1curUSe666+PccsutgzpXaek4nnjiy3R1+fH7\n/Xzuc19k6tTpPPvsL7HZbNx55908++wvSU5O5mMfu+eCMaenp5OTk4vVmsKMGbPYsWM74XAYk8nE\nmjVvsHDhYvx+/6Dv2cmTVUybNpMVK27gj3/8Ax/96Cd69+3evZPJk9Wgz5mWlsZHPvIxNm9+q3fb\nvn17e0tuTpmi2LVrB3PnzsNqTaG29jQdHR2kpr6X+M1mMyUlpVRVVVJaOo6Ojg5crrTe/eXl5VRU\nHOtNvACzZpVx4MC+QccL0u0wLNy+AGlXmHzPKByTz4olZYzNNuH3NPZOQUtKSiIlLYfajmTWbNnL\nicpTMbneaGR3OHE4Xf3+OFPTztk2kD9nJ/GLqamp5pvf/C533/1Jnn/+13zjG9/h7rvv4c03/wrA\nmcZzRcVxvvGN7/DUU9/mlVd+N6hzrVnzBi0tLdx66218//s/4cEHH+b5558F4O67P8n69WupqDjO\ntm1buPPOu885b98yI9XV1bS3t9HT04PFYmH69Jm9pSjffnsjS5Ysu+CxF7Nnz27KyuZw++0f4dVX\nX+n9HX7nnW38+c+v0dBQf9417Aaib52U1tYW7HY7AHa7g5aWFubOnYfT6aStrfW80zeXL1/B+vVr\nADhx4jjjxo3v3bdz505KS8efc8ytt/7DZcUqLd8h1u5ux+P1U1Q4IWbnNJlMTJk4jgmlQfYdPEpd\nexC7K1IO0Gy2YE7N41hDF5Wn9zBtYiFj8vMucUYxHKZOnQZAdnYOEydOwmQykZWVzb59e/u9bubM\nWQDk5ub3ru02mHNlZmbyq1+t5cUXn6O7O4DdHqlNbLFYeOCBz/CP/3g/3/729zCbzeec9+TJKh59\n9CFCoRAul5Mvf/lrva9bseIG3nzzr2RlZZObm4/NZj/vsWdaxwsWLDqn+8DtdpORkQlAUVExU6Yo\n1q79GzfeeBMLFy7mf/7n9/2S2YkTFezYsf283To33/z+fi3Xs4VC4d5PgKFQD2Zz5Otx48bzrW99\n7zxHmJg/fyHPP/8sn/zk/QMueOUYxBtwX5J8h1hNbROpThtpfT6+xIrFYmHu7Gl0eL2UHzqBu8uC\nzRlZ+yvZmgLWXPYcb8V1sp4ZU0rISM+IeQxi4Pomu75fn/2fvP++yN+vvfYKa9e+SWZmFkuWLLvo\nuV566bfk5eXz2GNf4/DhQzzzzPd797e0NJOWlkZDQ/15Y+zbb3umhGJdXS0A8+cv5Omnv0l2dg7X\nXXf9RY+9kM2bN/G+993c+/0dd3yUH/3ou9x44020traQlZXT7/Xjx09g/PiBN1z6JumsrKzeIvVe\nr7c36V+MxWIhP38MR4/qfn3RAAsXLuSVV84dFN21awfz5i0YcIy91xr0EWJQ3N4uXPahHQRLdTpZ\nOn8m9Q2NHDxeQ9CcFkm+gM2RSjewtbyG7NQaZk2d0PtRTMSX8ze0Ihtvu+12brvtdiDS53sx7e3t\nTJoUKXK+adN6uru7Aejo6ODll3/Hz372LF/4wiMsWbLsnFbbxRp7FouFOXOu4i9/+RMvvvgKWh8e\n8LFA7wyDvm8WZWVzCYXC7Nmzm46ODqZPn8HhwwcZP34CKSm23pbv2UwmEzfddMs5C432fSObPXsO\nhw8fZMmSZRw8eID58xdePMDovb722hWsWvVzvvGN71Bbe7r3nDNnzmTChEls3vxWb79vQ0M9qamX\nXuz0fIY9+SqlZgKvAU9rrZ9RShUDzxHpf64FPq617h7uuIaK2xdgTJZjWK6Vn5dLfl4uFZVVHD3V\nSLIji6ToL7rdlYEP2LDzKAWZVlZcM3dYYhIDd/4JE4OfRXHTTbfw9a9/hXXr1vKhD93B2rV/Y/Xq\n1zl06AB33vkxMjMzuf32O/nZz37M5z73rwOI4T0rVlxPW1vbeT9qX+rYdevW8MwzP+AXv/hJv+1e\nbwcvv/xbPvGJT3HkyGGKiopJSYnMDBpoy7ezs5M///lVTp6s5KWXXuTWWz/IvHkL2LZtC+vXr8Fk\nMrFw4eILHr99+1Zeeum3WK1WlixZxtGjGp/Px2uv/Q81NdVs3LiO22//AF/60hOsWvVzqqpOkJOT\nS0ZG5kXPezHDupKFUsoBvA4cAcqjyXcV8LrW+g9KqW8CFVrrn13oHCNpJQu3x82arZqr544nLzfn\n0gfEUE9PDwf0cWqa/dhd/UtUdgWC+HztzB6fzcTSc9cKG83c7vYLzHZIGdLZDiORrGTRX6xXshju\nlq8fuBn4tz7brgMejH79R+BzwAWT70hSc7oRm9VMbs7w1+c1m83Mnj6FyZ2d7D1UQavPhD01ne5g\niL9uP0W7N8DG8jbuurqNhXOnX/E81ZHC5UrjxkXnrjuWk+Oiqeny5/kKMVjDmny11iGgS6l+8/ic\nfboZ6oCC4YxpKLV7u3A5rIYmNrvdzuKrZtDc0sL+I6doclto90ZqRrh9QSrbLPi27mVB2WRSnVe2\nSu9IYDKZzttKTU93EQjIzEsxfOLht61vN4LprO9HNI+vG5cjPp44y87KYvniMuZPSiPNHvlnT3da\nyUpzYHLk8vbuY5ysPm1wlEIkjniY7dChlErRWncBRUQG3S4oM9OBxXLu/MSByM29vFHJy+HxeDBZ\nkpkyqWBYr3spubnTKZs9iXVbD+Lv6sGVERkMdKUWU93aQXdPJQuvmpGQT8jF079TvJB70l8s70c8\nJN81wIeAF6N///ViL25t9V3WRYZ78ODw0Qo6O7uwWhxxOWjx9yvncOx4Nbv3n6AnOT2ylhwWKpuC\nVP7pbRbOnozLdeEJ7KONDC6dS+5Jf5c54HbBfcPavFFKXaWUWg/cA3xWKbUO+CrwSaXURiATuLwS\nQXGmvaMLl9183uLp8SI9LZ3rlpSRnxqksyNSNc1stmB25rF5TwWVJ2VFZSGGynAPuO0GVpxn1/uG\nM47h4PEGyEhLMTqMSzKZTMyaNpnitjZ2HzgBtizMZgs2Vza6xktjywHmzZ6WkN0QQgwl+R81BLxe\nL75uyHSNnCfJMjMyWLl0Drk2P35vZC25FLsTT4+L9Vv34j5PGUYhxOWT5DsEqk/Xk2QKD6p4ejww\nmUzMnqFYMG0MQW8jPT1BksxmzM48tuytoqLypNEhCjFqSPIdAu0dXdiTTaSkxH+3w/lkZ2Wxcsls\nspJ9va1gW2oWx2r9bN+9n56eHoMjFGLkk+Q7BNy+wBUXTzdaUlISc2dP46opeQQ6GgiFQljtTryk\ns25rOW3tbUaHKMSIJsk3xnw+H/6g6ZIrFY8Uebk5XL9kNumWDvxeN0lJSSSn5rFt3ymOnbi85VOE\nEJJ8Y+7U6XrAREnh6ClgbjabmTd7GmWTsgl4Iq1gW2oWx+sDbN2575zFCIUQlybJN8baO/wkJ/Vc\ndo3PeFaQn8fKJbNINbnp8nWQYnPQmZTB+m37aGltvfQJhBC9JPnGmMcbiJt6DkPBYrGwcO4MZo5L\nJ9DRiMlkIjk1j3f213DkeKXR4QkxYkjyjSG/309n0DTiB9sGoqhwDCsWzcAeaqOr04vNlUlVUw9b\ndpZLN4QQAyDJN4ZOna7DbLFSkJdldCjDIjk5mcXzZjK12Inf3YjFmoI/KZP12/bT3NJidHhCxDVJ\nvjHU5vETDvrJzkqM5HtGaUkRKxdPJ6WnlWDAT3JqLjsO1XL46AmjQxMibknyjSG3t4s0hyVhVoXo\ny2q1snT+LKYU2PG7m7A5MzjVGmbzjvLeBRyFEO+R5BsjXV1dkfm9Q7xScbwbV1rM8gVTsASaIBwm\nYMli3bb9NDZJN4QQfUnyjZHq07WYk23kZI2+KWaDZbfbuXphGePzLHR1NJPiymO3ruOQPm50aELE\nDUm+MdLi9tMT8DEmP9/oUOLGpPGlXDNvMiZ/I0nJNmrcSby1fQ+BQMDo0IQwnCTfGPF4u3CmJGE2\nX94SR6OV0+Fg+eI5lGYnEezqIGjNYcM7B2lobDI6NCEMJck3BgKBAL5uEmJ+7+WaMnEcy8rGY/I3\nkWRN590jDezef4Tjp9vxB2ResEg8A1rJQilVDHwBuAkojW6uJLLe2ne11qeGJLoRovp0LVa7i8w0\nafVejMuVynVL5nDoSAXHaoP8Zl0zbl81+Zl2Hr93ATZr/C65JESsXbLlq5S6D3iTSLL9EJAL5AC3\nR7e9oZS6d+hCjH8t7k66OzsoKhhZxdONMm3KBEoLc3D7IlPQ6ls7efeQFGoXiWUgTY1ZwCyt9dmf\nDQ8CB5VSPwX+M+aRjSBubwBbcnjEFk83wpTSPAqyTlLb4sNlN9Pi6WJn+SGumqlkvTiREAaSfLcD\ndymljmitt5+9U2sdAD4f88hGiEAgQGcgTG669PcOhs1q4bFPzqemyYvdHODA0WraApms37qXeTPH\nk5GeYXSIQgypgTQxPgw8f77Eq5RK+CZKzek6UhxppDul1TtYNquFiYXpFObnsnLJbNLMXoJJTrbv\nq0Yfk0eTxeg2kOT5N611+AL7HoxlMCNRi9tHl99HUUGu0aGMaGazmYVzZzC12Ek41MPJ5pDMCRaj\n2kC6HW5RSqVfYN/NwE9iGM+I4/YGSDaFRmXxdCOUlhSRk53Jzr1H8JvTWL/9ILOnFFKQP3pWBhEC\nBtbydQD2C/xJ6EIG3d3d+AIhXHbp740lp8PB8iVzKEoPETYls/dYM++WHyIUChkdmhAxM5CW78ta\n65+db4dSKqEfU6qprSPFkU6qQx4SGArT1STyc1vYfbCK5q4MNmwrZ/6sCaS50owOTYgrNpCW74qL\n7PvvWAUyErW0eQl2d1OYn210KKNWdlZWZDDO4qPb5GDznkqOVlQaHZYQV2wgyfcNpdQ9SqmFZ+/Q\nWvuHIKYRw+0LEA76yMrMNDqUUa3vYJwpHOZEQ5AtUidYjHCX7HbQWv9qOAIZaYLBIF5/iHRnYhZP\nN0LfwThvOI112w5QNqWQMTIYJ0agQT1Mr5R6H/BB4KvAU8AY4Hta6zeGILa4VlNbR4ozHZddpkIN\npzODcQf1MU42JrPnWBNjGloom6nkTVCMKIN9SEIBTwIPAZ8F/g+wINZBjQTNbR2EQiFys2XwxwjT\n1SQWTC/AHO6isTOFDVv34vF0GB2WEAM22ORbpbWuAf6itXZHH74oH4K44p7H102330N+nnzkNcqZ\nwbj0ZB8Bk4PNeyqoqJQCPWJkGGwNvy8ppe4GspVSzYCJSInJP8U8sjjW09OD199Dqs0sxdMNdmYw\nrupUDYdOdHK0rou6pnIWzpmOxSIlKkX8Guxv59e11n/uu0Ep9f4YxjMi1NTWkWxPx2XrMjoUEVVa\nUkRuThY79mi8oTTWbd3HnKkl5OXmGB2aEOc1qG6HsxNvdNvrsQtnZGhu68CUZCIzzWF0KKIPh93e\n+2RcyGTl3SMNlB88Qjh8odIkQhhnUMlXKTXhrO8v9gDGqOX2BvB72ykqlOLp8ejMYJwp5Keuw8KG\nrXvo8HqNDkuIfgY74PbUmS+UUnOAn8Y2nPh3pr/XngxWq9R0iFfZWVlcv7SMjOROAiYnb+8+RmVV\nQq92JeLMYPt8W5VStxNZ3WIm8P3YhxTf6urrsdjTSLNKf2+8O3swTtd0Ute0j/ll02QwThhusH2+\nnwGOAA6t9YeAhJvX09jiwWJJJi3VZnQoYoBKS4q4dsEUbCYf7qCDdVv30djUYnRYIsFd8u1fKfWP\nwCeBMzPYTYBLKbUamA0UDVl0ccjj66bT76ZoSumlXyzixpnBuMiTcVZ26TqKGpuZNW2y0aGJBDWQ\nz17ZwDKt9TlVTJRS/xD9e9H5lhkabUKhEJ7OIMlmE6mpqUaHIy5D3zKVte4UmrfuYUHZFHKRYvhi\neA2k22Ex8IRS6nx1Ezcppf4DeCy2YcWnuvp6LDYXLntC15Af8c48GZdp9eMP21m3/Qjrth3GH5C6\nzGL4DKTl+/fAF4gsE18JnALCRJ5sKwH+C7h1iOKLK43NbqwpqbhssqLCSGc2m1kwdwZHK07xoz/X\n4dneQobzGF/71GJSHbIYqhh6AykpGQK+rZR6mkgRnRIi/b4ngR1a656hDTF+uDu76Qx4KCrNNzoU\nESNJtjQ8nZFf4TZvD3/eUM7KeePIz5cFUcXQGvB8m2iS3Rb9k3BCoRAdnUEId5MpxdNHjaIcJwXZ\nDmqbfWS5kkl1prD7aCNj6psom6lIShrsVHghBkYmOw5QQ2MjZmsqDnNCL94x6tisFh67Zz6+YBiH\nxUSgy8eufRU0+lJZt2UvV80YJyuViCEhb+sDVN/UhtVmJ80hT7WNNjarBVWahc1qIc2VxnVLyijO\nhGA4me37a9h/+KjUhxAxd8nkq5Safdb3CblapMfXTZe/k9zsdKNDEUPMZDIxXU1i8exirEnd1LTC\nRinWLmJsIC3f7531/ctDEUg8C4fDuH1Bero6yM+TgZhEkZGewcqlcyjMAH8ombf3VMjKySJmBtLn\ne/bCWDFdKEspdRXwR+Bo9NzlWuvPxvIaV6qhsYkkqwN7UpcMwCQYk8nErGmTKWxpYc+hKo7XBaht\n2MOCMoXdbjc6PDGCDST5nt3ZFevOr1TgZa3152N83pipb2zFZnfiMsv83kSVnZXFiiUZlB84Qr3b\nyqYdR5g6PofSkoR6ul7EUDw041zEuDUda25fgGB3gKx0p9GhCAMlJSUxZ9ZUrlJjMCeFOFjlZuvO\nfQSD8mScGLyBtHyXKqX6Vi/Li35vAsJa67FXGEMqcLVS6i+AE3hCa73hCs8ZM+FwGE9nkFConaKC\nGUaHI+JAbk4WK5dmsGe/pt6dxLpt+5k1uZCCfFlMVQzcQJKvGuIY9gBf1Vq/rpSaDKxRSk3UWsdF\nc6KxqZmkZAcp4U6Sk6Wmg4hISkriqtnTqG9opPxIDXuONlJT38xV8mCGGCDTpeYvKqWuvdh+rfWm\nWAaklNoO3KG1rjrf/mCwJ2yxDN+KwbvLD9PgteGyeFi2YNawXVeMHD09PWzffYC6thDJpi4Wz5lM\nbk6W0WGJ+HDBLtWBtHw3AIeBd4DQWScLA1eUfJVS9wKpWusfKqXGAHlAzYVe39rqu6zr5Oa6aGz0\nDPq4quo2AuY0UtPCl3V8vLrc+zFaXen9mDxuPKn1Dew76mH1Jk1xTgozp07CZIrr4YyLkt+R/i7n\nfuTmXrhU6UCS73IixdSvBv4CPK+13j2oCC7uVeCF6PJEVuCheOlyCIfDeHwBQqZ2SqZOuPQBIqEV\n5OeRm53FrvLDnGwM0NS6l3mzJpDmSjM6NBGHBlLV7C3gLaWUHfgQ8K1oC/VF4IULdQ8MlNa6Dbjl\nSs4xVJpbWjAlO7GGOnE4ZJl4cWkWi4VFV82k+nQtB47Xs3lPJeMLXEydNN7o0EScGfDIgNa6U2v9\nPHAT8APg88CuoQosHtQ2tGBzOHE5ZKBNDE5xYQErF88gw5HEsRovm7btwdfZaXRYIo4MOPkqpaYp\npf4LqABuBh4ECocqsHjg9gYIh8OSfMVlSU5OZsm8mcyakImvGzbtOEJlVbXRYYk4MZAFND8N3Etk\ncO05oExr3TrUgcUDty8A5g6KJowxOhQxgpWWFDEmL4dd5Zp9J5o53djCgrJpMnUxwQ1kwO2nROou\nnAbuAD6s1HtTf7XWK4cmNGO1tLSAxQFBHxnpGUaHI0a4lJQUli6YzYnKk+iTrazffoCZEwsoLJBV\nURLVQJKy+q9jAAAR+ElEQVTvmZGCMHH+GHAs1dQ3Y3ekYumOi4kXYpQYP24sY/Jz2bXvCLsO11JT\n38RVs6ZiNg/f3HURHwaSfFuBLwHTgLeA78bLVLCh5PEFwOwgTfp7RYzZ7XauXljGsRNVHDnZwrqt\n5cyZWioPZiSYgQy4PUOkxftzIgn48SGNKE64fd10dfrIz5EuBzE0Jo0vZfkChdViYlt5JeUHNJ1d\n3Rw/3S7L2CeAgbR8x2mt7wZQSq0G1g5tSMZra2sjbEohFPCSmzvR6HDEKOZ0OFi+eA5HjldyqLKF\nVWs24+4MUZDt4LF75mOzyjKLo9VAWr7dZ76IrmA86hezqqlrxJ6ahtNmliIpYlhMmTiOCWMLcXdG\nakbXNvuobpBli0azgWSWoS6mHnfc3sj7jfT3iuE0uTSHguzIk5QuWxJHKk5wuq7e4KjEUImHer5x\nx+0LYLYFyMpKNToUkUDOLGNf0+Qly2nm4JEKdh48TW5NI/NmTcFqlZWzR5N4qOcbV9rd7fRgJehr\np1CKp4thZrNamFgYWSF72YLZVJ2q4XBlA2u3HmBSSSaTJ4wzNkARMwMprHNFhXNGmuraSH9v2N8i\nTyAJw5WWFFFUkM+e/Uc4XNXG6Ya9zJ0xXiqljQIymnQWjzeAyWTC5ZCPeCI+WCwW5s+ZzqKZxXQH\nQ2zaVUH5wSOEQrKg60gmyfcsbl+AUE8PGS6b0aEI0U9ebg4rl5ZRmu/kVIOPdVv2Ut/QaHRY4jJJ\n8u3D7XETDCfT2dFOsTxzL+JQUlISs6ZN5pp5E0i2JLHjYA3bd++nu7v70geLuCLJt4+a043YU9NJ\nMYew2+1GhyPEBaW50li+uAw1NoOWjiBrtx7gROXJSx8o4oYk3z7avV2YTCZSZX6vGCEmTxjHioVT\nSXOY2V/Zylvv7MXjkYczRgJJvn14fN2Ew2HSnClGhyLEgNlsNpbOn8XcKXn4ukK8tfsY+w8f5VIr\nkwtjSfKN6ujw0B224Pd6KCrIMTocIQatuLCAG5bOJj8zhar6TtZt2UNDY7PRYYkLkOQbVX26AXtq\nOqZQF+lp6UaHI8RlMZvNzJs9jUUziwHYcaiGnXsPEgxKlbR4I8k3qr13fq/094qRLyc7i5VL5zA+\n30FDezfrtu2n6lSN0WGJPiT5Rnm8AQBJvmLUMJlMTFeTuGbuRFIssK+ihc07yunweo0OTSDJFwCv\n14u/x4y/00t+TqbR4QgRUy5XKssXz2FGaToef4i3dx3joD4uA3IGk+QLVNc24HClEwr4yM3JNjoc\nIYbE+HFjWbloOulOMyfqO1m/dS/NLS1Gh5WwJPkC7R1+kpKScNktUjxdjGpWq5Ul82Yyd0ouPSHY\nvr+GXeWH6OnpMTq0hCOZBnBH+3vl4QqRKArH5LNyySwKMizUtwdZs6WcUzWnjQ4roSR88vX5fHQF\nk+gOBMjOkOLpInGYzWbmzp7G4pnFWM0m9lW0snlHOZ2dnfgDQXRViyzkOYQSfnW+U6frsbsy6HQ3\nUzhmptHhCDHssjIzuW5JBvroCSrrPKzZrtl4wEOLp1sW8hxCCd/ydUf7ex0pSVgs8gsmEpPJZGLq\nlAlcPW8S3YEgLZ5IlbTaZh81TTI1bShI8vVF5/c6pXi6EKlOJ//nujnkpEX+P7jsFo4fP0HlyWqD\nIxt9Erqp5/f76ew24ejpISNNiqcLAZF15L52/2K83WFqq09R2+jm0KkOTlTvZXxxNuPGFhsd4qiQ\n0Mn31Ok6HKkZ+DraKJ4xxehwhIgbNquFkiIX2c5kZkwJc+xEFZW1fg6e9EgSjpGE7nZo8/hJMpux\nWcLYbNLyFeJ8TCYTkyeM44alZUwpshMKhzlQ5Wb9lr3SHXEFEjr5ur1dALJYphADYDKZmDR+HCuX\nljGtxEkoDAclCV+2hE2+XV1d+LtNkeLpknyFGDCTycSEcaVcv6yM6WNdABysamf9lr1SOW0QEjb5\nnqqpxZaaTmeHm0Ipni7EZRlXWsKKpWVML43UwN5/ok2S8AAl7IBbq8eP2ewiiYAUTxfiCo0bW8y4\nscWcrK7h+Kkm9p1opeJUExNKsiktkYG580nYlq/nTH+vXeo5CBErY4uLWLGkjDkTszElwb6KVtZv\n2cvJamkJny0hW76BQABfN6QixdOFGArFRQUUFxVQU1vHsaoG9h5r5vjJJiaOzWFscZHR4cWFhEy+\n1adrsadm4O/0MqYoy+hwhBi1igrGUFQwhtq6eo5W1UsS7iMhk2+LuxOz2UVXoIXcnMlGhyPEqFcw\nJp+CMfnU1zdwpKqevUcbEz4JJ2TydXsDJNkj/b0mk8nocIRIGPn5eeTn51Hf0MjRqrpIEj7VxMSS\nxEvCCZd8A4EAnYEwTrv09wphlPy8XPLzcmlobOZI5WnePdKQcEk44ZJvzek6UpzpBAJd5OS6jA5H\niISWl5tNXm42jU3NHKms5d0j9f2SsD8QpKbJS1GOc9TVFB5dP80AtLh9WCwuvG2NjMkfa3Q4Qggg\nNyeb3JxsmltaOFxxmj26nsOVjawp76DJHRiVRd0Tbp7vmfXanDYpni5EvMnOymLZ/JksmzsOX1eY\nJnfk/2tts4+t7x6ls7PT4AhjJ6GyT3d3N96uEKl2SJV6DkLEreysLG66Jo1tR7fT0NZFusOM2+Pj\nfzfuJd1lJ81hJd1lo7ggD6fTaXS4lyWhkm9NbR02Zzo9PUEy0+1GhyOEuAib1cIT9y3q7fNNSTbT\n2tpKXWMrHl+AU/Vu9hyqIs3lJCs9lYzUFIoK8khNHRkL4cZF8lVKPQ0sBkLAP2utdw7FdVravFiS\nXXjdzRRJ8XQh4p7NamFi4Xu1V7KyssjKeu/BqO7ubmrrG2ht66Cxzcf+oztIMiUxJi+HDJeNksI8\nUlPjc2Dd8OSrlLoWmKS1XqqUmgb8ikgijjm3LwA2sJmleLoQo0FycjJji4vou6hGR4eHmrom3N4u\nNmzdQ5e/k+zsbPKzMygqzCXNlWZcwH0YnnyB64HXALTWh5RSGUqpVK11RywvEgwG8fpDOG1SPF2I\n0Sw11YWadKa1O5VQKERjYxMNLe1s370fn9eLLSWFgoJ8SovHGFbVMB6S7xigbzdDQ3TbsVhe5GT1\naVKc6YTDYVmpWIgEkpSU1Ptk3axpkXICfr+f03X17N2v8Xk7sFjMZGVlM2FcMRnpGeecwx8Ioqta\ncFhMMZvuFg/J9+zne01AONYXaWrpwJJsxedpo2SyzO8VIpHZbDYmjCtlwrhSAMLhMG1tbZysruWw\nPobFYsFms1JUUIDd6eLJX++kttkX0/nG8ZB8a4i0dM8oBOou9OLMTAcWi3lQF+jsClLV4CPFmYkV\nM+PHF15epKNMrjzh14/cj3Ml0j3Jy0tjypT3GmbBYJDaugYOHKuittkHROYb+4JhSoqu/L7EQ/L9\nG/AE8Aul1FygRmvtvdCLW1t9gzq5PxDsfddKdzZy28J0Ghs9VxTwaJCb65L70Ifcj3PJPQFbiosZ\nk+wUZDf2tnwdFtOA78vF3rwMT75a661KqV1Kqc1AD/BwLM9f0+Ttfddq9wboCkoVMyHEwNmsFh67\nZz6+YHjU9fmitf73oTp3UY6TgmwHtc0+XHYLMybmDdWlhBCjlM1qoaQotp8E4iL5DqUz71rb9h3H\n3e6jMD/X6JCEEGL0J1+IJOCibBuWULcUTxdCxIWEqmomKxULIeJFwiTfrq4ucrISZ9qMECK+JUzy\nDQa8jMnPNzoMIYQAEij5Tp08AbN5cA9nCCHEUEmY5FtSUnzpFwkhxDBJmOQrhBDxxBQOx7yGjRBC\niEuQlq8QQhhAkq8QQhhAkq8QQhhAkq8QQhhAkq8QQhhAkq8QQhggIaqaKaWeJrIcfQj4Z631zksc\nMqoppb4FXA2Ygf/UWr9qcEiGU0rZgAPAV7XWvzE6HiMppT4GfBHoBh7TWv/V4JAMo5RyAr8BsoBk\n4Gta67/F4tyjvuWrlLoWmKS1Xgo8APzI4JAMpZS6DpgevR83A98zNqK48RjQZHQQRlNKZQFfAZYC\n7wduMzYiw30SOKy1XgF8GPh+rE486pMvcD3wGoDW+hCQoZRKNTYkQ20k8ksE0Ao4lFIJXeRYKaWA\nqcBfjI4lDtwAvKm19mmt67XWDxkdkMGagOzo11lAY6xOnAjJdwz9b1gD/VdLTiha67DWujP67QPA\n/2qtE/0xx+8AnwcS+k0oahyQpJT6nVJqo1JqpdEBGUlr/XugVCl1FNgA/Euszp0Iyffs/1AmINGT\nDUqpDwD3Ao8YHYuRlFIfB7ZorauimxI9AZuAIuCjRH4/fmVsOMaK9n9Xaa0nE/kUHbNuy0RIvjX0\nb+kWAnUGxRIXlFJ/B/w/4CatdWKvDQ63AB9QSm0F7ge+nOCtvXoib0YhrXUF4FFK5RgdlIGWAW8A\naK3LgSKlVEzyZiIk378BtwMopeYCNVprr7EhGUcplQZ8C3i/1rrd6HiMprW+U2u9SGu9BPgl8KTW\nep3RcRnob8BKgGjSdWqtE3kg8hiRmVIopUoBj9Y6FIsTj/qpZlrrrUqpXUqpzUAP8LDRMRnsI0QG\nEF6KDrSFgU9orauNDUvEA631aaXUK0qp9YCdBO+WAn4GrFJKbSAyNfPBWJ1YSkoKIYQBEqHbQQgh\n4o4kXyGEMIAkXyGEMIAkXyGEMIAkXyGEMIAkXyGEMIAkXzGklFKlSqmQUuqus7afGMCxH4v+na+U\n+n0MYgnF6umkWFNKvaiUKjA6DjF84vIXUYw6R4DHo7VRz7joBHOllJlIaUOi1bU+EoM44nZSu9b6\no1rrWqPjEMNn1D/hJuJCLfBXIsn0//bdoZRyEClWnQm4gJe11t8G/ptINam/Enmq6G3g74A/aK2n\nRo8tBrYBJUTKZD4CdAFtwKe11q1nxWEC/k0pdQ2QS+TJvoNKqduAfwU6ifyf+LjW+qRS6rPAxwAv\n4APu1lq3KqUeiV6vG6gEHtZad/X5mUqBtcD/AnOIJP07tda1SqmHgE9E4/QDH9Fau6OfBK4HriFS\nRzcDeDp6T/4F8BBpLN2rta4cxL0XcUpavmI4hImUbbxFKTX5rH15wKta6+uJrK7xpWi95ceBBq31\nTWfOobU+CPiUUjOj2+4AXiRShevfgeu11jcCbwFfukAsh7XWNwPPAF+NbssA7ojGsJr3Hqn9KnBL\ntJD294BCpdQC4B+01su11jcAbiIFec42Hliltb6WSA3lL0S324Abo+esAu7uc4/OKANu1lqvjv5c\nD0dj+9fozypGAWn5imGhtQ4qpb4I/IDIChpnSjc2ANcqpf4RCAApRIpWX8hviRRK2k+kTsX9wBKg\nAHgjWq/CClyoT3lN9O8tvJcQG4DfRPuD84Gt0e2/jJ7zFSIt8qPRn2GiUmpd9GdwROM+W7PWek/0\n683AZ6NftwCrlVIhoBQ4Hd3et5Tlbq11MPr1r4BfR2N4VWu94wI/lxhhJPmKYaO1Xq2Ueij6Mf9M\nS++fAavWehmAUupSKwX8FvhfpdSzQIrWep9SajywXWt96yWODRNZxw+idZ2VUhbgd8AcrXWFUuph\nYF403n9RSpUQKTv5qlLqX4h0F/xJa/3oJa7V91PlmWsVAf8FTNNaNyulvn2BY3uTudb6+0qpF4Gb\ngJ8qpX6ptf7FJa4tRgDpdhDDoW+r7nPAN4i0cCHS0qwAUErdSqSSVgqRJGk7+xzR6mvNRBZ4fD66\nbwewUCmVHz3P7Uqpv79AHNdHv74a2EekTzUJqI4uovkBIEUplaGUehyo1lr/lEg3xQIirdibzwwe\nKqU+o5RadJ5rZSmlZve5VjmRLhZPNPFmAe/rcx/OoZRKUkp9A3BrrZ8j0g2y+EKvFyOLJF8xHHr7\nM6MFul8hknQBVgF3KaXWEPkY/kL0Tw1wWim1A3DSv0/0BeBTRPp7ic4S+CzwerT0331EBuLOFgRm\nRAfxPg08ER2U+zWwnUir+ltE6tleD6QCO5RSbxJp/f5Ca70L+DGwQSm1CVgO7D3PtWqAe5VSa4ks\nRvldYA+wTym1DfghkQHIe5VSyzjPTIxo3dgmYEv0/nyOSMtZjAJSUlKIGIvOdnhba11idCwifknL\nV4ihIa0acVHS8hVCCANIy1cIIQwgyVcIIQwgyVcIIQwgyVcIIQwgyVcIIQwgyVcIIQzw/wHgcq3E\nISXOpgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f090fbcf910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5,3))\n",
    "\n",
    "# error bar is barely visible\n",
    "#plt.errorbar(rc, rna_pmf, yerr=dh_log_err, fmt=\".\", capthick=1,capsize=4\n",
    "rc = np.arange(0,9)\n",
    "rna_pmf = og*-1\n",
    "rna_pmf -= rna_pmf[-1]\n",
    "plt.plot(rc, rna_pmf, \".\", label=\"US DHAMed\")\n",
    "\n",
    "plt.fill_between(rc, min_pmf, max_pmf, alpha=0.4, label=r\"min-max PMF $\\Delta t=10^4 \\, \\mathrm{MC}$ \")\n",
    "\n",
    "ax.set_xlabel(\"Native base pairs\")\n",
    "ax.set_ylabel(\"PMF ($\\mathrm{k_BT}$)\")\n",
    "ax.set_ylim(-1,13)\n",
    "ax.set_xlim(-0.4,8.4)\n",
    "ax.set_yticks([0,5,10])\n",
    "ax.set_xticks([0,2,4,6,8])\n",
    "ax.tick_params(direction='out')\n",
    "\n",
    "plt.legend(borderaxespad=0.1)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"rna_pmf.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "The figure reproduces the free energy profile presented in Figure 7B of the manuscript [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "## Vary observation interval\n",
    "\n",
    "Here I vary the observation interval. Using a sliding window to determine transition counts I can analyze trajectories even at long observation intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "c_l_lag_e7 = [count_matrix(traj,n_states=9, lag=5000)] # 10^4 * 5000  -> 5x10^7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36\n",
      "loglike-start 1134671.76583\n",
      "1134671.76583\n",
      "Warning: Desired error not necessarily achieved due to precision loss.\n",
      "         Current function value: 1072622.665740\n",
      "         Iterations: 20\n",
      "         Function evaluations: 82\n",
      "         Gradient evaluations: 70\n",
      "time elapsed 0.0179169178009 s\n"
     ]
    }
   ],
   "source": [
    "og_e7 = run_dhamed(c_l_lag_e7, v_ar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAADRCAYAAABreLGdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4W+WV+PGvLFmWZcv7EttxnP3NvpCFLISQsAwMDGVa\nylKgLVAKHSgt7TAzv2kptHTotB3oTjcmhbK1wBRaoGHJzpKEJMZx1jeL9323ZMm2vOj3h2RjJ3Zi\nJ7KvbJ3P8/iJreXeoxv76NV733uOyefzIYQQYnRFGB2AEEKEI0m+QghhAEm+QghhAEm+QghhAEm+\nQghhAEm+QghhAIvRAQxXba3rnNbGJSbaaWz0BDucMUuOR39yPE4nx6S/czkeqakO02D3hc3I12Ix\nGx1CSJHj0Z8cj9PJMekv2McjbJKvEEKEEkm+QghhAEm+QghhAEm+Qoywjz/ex7e//e/9btuw4Xf8\n5S8vA5Cbu5d7772L++77Ml/60ud56aUXTtvGhg2/46abPs3999/DvffexeOP/xCXyzXo9h977Lvs\n3Pl+788HDuxnzZplnDx5ove2jRvf4Mor19HZ2dl7m8vlYt26lWzc+MaQXltBwUm++tW7h/RY0Z8k\nXyFGgck06Elvfvzjx3j00f/ml7/8HU8++RTbt2+loaH+tMfdcMPN/Pznv+FXv/o9c+bM5eGH/3NI\n2wfYtOltcnIms2nT2/1uT0hI6Jekt23bTHr6hKG+rCHtWwxMkq8QBnO5XLjdbgCsViu/+tXvSUpK\nPuNzrrrqGrq6Ojl8+CAAZ6pO2N3dzfbtW3nwwW+xadM7/e5bsWI17777SULeunUTS5de2Pvzq6++\nwle+cif33fdl/vzn5wGora3hnnvu4MEHv8Ybb7w2vBcreo25db5CjIY2byfldW6yUmKwWUf2z+RL\nX7qbu+76PIsXL2XZsgu5/PIrcTgcZ33e9OkzKCwsIDMzi7y8XO6//x7An4hLSopZt+5SAPbs2c3k\nyVNYuHARCQnxHDp0kLlz5wGg1CxefPFZWltbaWtrpaOjg+Rkf+IvKytj69ZN/PrX/wvAV75yB+vW\nXcYrr/yZyy67guuvv4nnn3+GEyeOj8RhGfck+QpxijZvJ48+s5fKeg8ZyXYe+sLSEUnAPR/Xr7vu\netauXc/u3TvZsWMbf/zjBjZseO6so9/Ozk4iIvwfXhcvXsKjj/53732PPfbd3u/fffctLrvsHwC4\n7LJ/YNOmt3uTL8DKlRexY8dWWlpauPjidbS0+OeS8/PzKSsr5f7778Hn89Ha2kZlZSVFRYWsX395\n7353794ZhKMRfmTaQYhTlNe5qaz3X8lUWe+hvM59XttLSEjE5XL2u62pqZHk5BQA2tvbSUxM4sor\nr+axx37M8uUr+OijXWfd7tGjR5g5c9YZH9Pe3s4HH7zHyy+/yB133MIbb/yVbds293vMunWXsXXr\nJrZv38Ill6zvncKwWq2sWnURP//5b/jFL37L00+/wMKFiwAfERH+N47ubmnGcK4k+Y6CJUvmkZYW\nN+DXkiXzzr4BMaqyUmLISLYDkJFsJysl5ry2l509ibq6WsrLywBobGwkN3cf8+cvpKyslDvvvI3W\n1lbAPz9bV1dHRkbWadvpO6/717/+hfj4BKZNm37GfX/wwXssWbKMZ575Exs2PM+zz77EpEmTyc3d\n2/uYWbNmU1VVRVdXF6mpab23z507l9zcvbS3t+Hz+fjZzx7H6/WSnZ3DkSOHAfj4472n7VMMjUw7\njILS0hJqapy0eTs5cKyc6dnJJMbHApCWFmdwdOJUNquFh76wNGhzvhaLhe9851F+9KP/wufz4fP5\n+MY3/o3ExEQSExO59dYv8LWvfYXo6Gi8Xi9r1qwNjDD7e+WVP7F9+xa8Xi+TJuXwyCPfP+u+N29+\nm2uu+VS/266++p/YtOlt5s9f2Hvb8uUrSExM6ve4jIwMbrjhZu6998uYzWbWrFmL1Wrls5+9me98\n5z/YsWPbWZO/GJxprPVwO9fCOqmpDmprXcEOZ0jS0uIoKWvonUd0RJu5dW0ai+fNIDMziZoa59k3\nEmRGHo9QJMfjdHJM+juX43Gmwjoy8h0lfecRXa1dFDZYqN/5yTIhWSspRHiROd9RkpUSQ1pCFADx\nMVYS42zYHKkAbPlwPyVl5UaGJ4QYZTLyHSU2q4XPrk7hWEkTKfFW6PKCxQaAOSaVI2UeCsv2M3va\nRNJSz7zESAgx9snIdxS1e7uYlZPA5RctZkZGFN6Wmt77omx2fLYUco/V8OHeg73X7QshxicZ+Y6C\n7OxJg65qSJ+QQWtzLba4FEwmE7aYOLzA+3lFpMZbWDBrGlardXQDFkKMOEm+o2DfvoM0NDSwZfdh\nrrt8BRZL/8Pe1tZG/pGTNLSYsMXGAxDtSMLV5WPLrsNkpcYwV03rvZpJCDH2yV/zKCmvric5Ie60\nxAtgs9lYvnguy+dmENFeR3ur/4oqk8mELS6VWk8U736Qz/GColGOWggxUiT5jhKXx4sj5szTB0mJ\niVx84ULmTIqls6WGzg4vABFmM1GOVIrqutnyYR7lFVWjEbIQYgTJtMMoaW5pZ0Z20tkfCGRnZTIx\nMwN9vJDCqlpsDv98cKQ1CqypHCxyUlBWw5zpE0lOGto2hQglb7zxGoWFBXR0dHDdddczdeo0o0Ma\ndeZHHnnE6BiGxePxPnIuz4uJicLj8QY5mqFpampCF1awdMHMAacdBmIymUhJTmRyZjJN9VU0Oz1E\nRkUDYIm04jPbKS6vo762hoT4mGGflDPyeISikT4eVVWVfPrT17Bnz27efPNv/O1vrxIfH09UVBRX\nXbWeFStW9aurcOedt3Ho0AHWrLmk33P//vfXeeutNzlx4hgXXrjytP3k5u7llVf+NOB9w9X3mJw8\neYKkpCTKy8uw2WyYzQN38vX5fPzylz/pt/+nn34Kl8vJ9u1bWLhwMQDJySlccsmlmEwRvZc579mz\nm+997yFqa2tYvHhJ7/N37NjGAw/cy4QJGUyePAXwV3R76qnfoPVRyspKaWxsYNu2zSxYcPpl2cFy\nLr8jMTFR3x3sPhn5joLyqloS4mKx2WzDfm5kZCRLFszG5XKRf7QIV4cVm91fFyI6NoE24L3ck6TF\nR7Jg9nQiIyODHL0IlpycHH7+898A4HQ6ueOOW3jiiV+QmZnFpk1vM2eOv8hSeXkZbnfLoM89k2PH\njrJ587vccceXiYmJ7XdfcXERO3Zs5bbbbh927F/96t1YrVZuvvlWbrzxlgEf43Q62bjxdfLyPu69\nbe/ejwBYs+YSjh8/xv79eSxcuIjk5BRyc/cyffqM3scuW3YhdXW1vPnm37jjji/33t7e3sbEidms\nXbu+97bvf/9hLr/8SlavXgP42ySVlBQP+3UZSeZ8R4HT3YHDfn7LxRwOB6uXzWfB1AS6PbV0ej95\nB452JNPcGcuWnYc5rE/Q3d19viGHrdGqQBcXF0dKSipWaxRz585nz57dvVXLNm16m+XLVwx7myUl\nxcyePY916y7jr3/9y2n35+buZcYMdU7xfv3rD/LaaxsHTbzgf0033ngLMTGfVIE7cGB/7z5nzlTk\n5u7pvW/nzg9ISUnttw2z2Ux2dg7FxUUAtLS04HD0X6Z59OhhCgpO9CZegPnzFzJlytRzem1GkeQ7\nCpweL3HnmXx7ZE5IZ93KhUxKNtHmqu1NtBEREUTFpVDZEsmmD/dTWFQalP2Fm54KdAN9lZaWnNe2\n+9awqqysoLm5ia6uLiwWC3PmzOst8/j++9tZuXL1oM8dTF5eLgsXLuL662/k1Vdf6fcm/NFHu3j9\n9deoqakesD/c2Wh9mJ073+fFF58DoKuriy1bNvHMM//Lxo1v8PjjP6SiojwQ6yfBNjY2EB3tny6L\njrZTX//Jvtva2gbc19q169i6dRMAhYUne6caenz8cS45OVNOe9611/7zsF+XkWTaYYQ1O5txudvI\nygzeu7LJZGLmtMlMzenkwOHjVDV3Eu3wn3gzmy2YY9M4UdNOUUUes6dlMiE97SxbFKOhpKSY+++/\nh+7ubqKibHz729/rnTtdt+4y3n33LZKSkklNTcdmix7wuT1FmJYtu7Df9IHT6SQhIRGArKyJzJyp\n2LLl3d4OFsuXr+D//u/P/RJUYWEBe/bsHrCo01VXXUNq6ietjO677wFMJhMVFRV89NEu4uMTWLNm\nLVu3bqazs4N16y7tLQ7fV3e3r3d9end3F2bzJ+O9b37z3097PJhYunQ5zz33NF/84pfO2JvuVHb7\n+dVdHm2SfEdYeWUdsTE24hzBr9trsVhYvGA2LW43+UcKcbZbsMX4/2B6VkbknWzEUVLN3JnZJMQn\nBD0GMXQDzdtWVVUCsHTpcp544oe9J6KG8ty+PvhgB1dccVXvzzfc8Dl++cuf9CbfxsYGkpL6J8cp\nU6YO6aP63//+Ot3dXVxzzXVERUVx8uQJbr75VgAOHcrnxhtvITPzk+LvfZN5UlJSb6F4t9vd+wZx\nJhaLhfT0CRw/rk+btwZ/66K33jq9tf2+fXtYsmTZWbcfKmTaYYQ53e04okf2JFhsTAyrls5j8Yxk\nfK21dHjbe++z2WPpiExiZ345H318qPcPQYy+Mw3iLBYLixZdwJtv/o2LLlpz2v1nem5nZydAvxUI\nCxcuprvbx/79/pNfhw4dZM6cuRw9epj2dv/H/cLCAl566cXTvl5++U/9aovExyewevXFgP/NYsaM\nmRw5cojm5qbeBp49+/HH+kmwCxYs4uRJf4PNw4cPMXfu/MFfiP/ZAFx88To2bPgd06ZN7y1A32PW\nrNlMnTqdDz54r/e2mppqYmPP3nQ0lIz6yFcpNQ94DXhCa/2kUmoi8Cz+N4JK4DatdcdoxzVSnB4v\nE5Lso7Kv9LRU0tNSKSgq5nhpLZH2JCICf5DRjgRuv+1SaqsrBnxudvYk9u07OCpxhquzlWxet+5S\nmpqaBvz4fKbnbtmyiSef/Dm///2v+93udrfw0ksvsHDhYlJSUjl27ChZWROJivKvuhnqyHfVqot4\n5ZU/Y7fbSU1NY+nS5Tz99FMkJSUzb94CduzYRnx8Aq2trbz++quUlBTx0ksvcO21n2bJkmXs2vUh\nW7duwmQynfFE4u7dO3nppRexWq2sXLma48c1Ho+H1177P8rLy9i+fUvviodvfesRNmz4HcXFhaSk\npJKQkHhOJymNNKqdLJRSduAN4BiQH0i+G4A3tNZ/UUr9ECjQWv92sG2MpU4WTpeTTTs1Fy2eQlrq\n6fNhI6mrq4tD+iTl9W1EO/wlKm+4YjYvvXOEdm8nHk8zC6YkMy3H/3ExLS3OkI4aoaLn9+NMxyHc\njpF0suhvrHeyaAOuAv6jz22XAHcHvv8r8AAwaPIdS8orarFZzaSmjH59XrPZzII5M5nR2sr+IwU0\nevy/Ax2d3by1u5Rmt5ft+U3cfFETyxfPGfX4QtWZKtBlZ08a5WjEeDaqyVdr3Q20K9VvrWFMn2mG\nKiBjNGMaSc3udhx2q6EtgqKjo1lxwVzqGxoAqGtw0uz2rxF2ejoparLg2bnfsPhCjUy9iNESCifc\n+k4jmE75eUxzeTpw2EPjirOeGhBLp8cRF+3/b4+PsZIUZ8dk9y90LykbeD5YCBF8obDUrEUpFaW1\nbgey8J90G1Rioh2LZeDrys+m77rFkeZyuTBZIpk5PWNU93s2qy6cw8IF09my8zBt7V04Ej45GVjW\n2EFHVxHLL5gblrWDQ+n/KVTIMekvmMcjFJLvJuAzwAuBf98604MbGz3ntJPRPnlw9HgBra3tWC32\nkDpp0RPLP61fxImTZeQeLKQr0l/A3dtloaiuk6K/vc/yBTNwOE5fYzleycml08kx6e8cT7gNet+o\nJl+l1AXA40AO0KGUuh64BXhGKXU3UAw8M5oxjZTmlnYc0eYhVzEbDWc6mZSSOgHwXyFHTBof5BUw\nKyeZyZOyBny8EOL8jPYJt1xg3QB3XTGacYwGl9tLQlyU0WH00/dk0qnv4o1NTeQeKgRbEmazBZsj\nGV3uprbhEEsWzA7LaQghRpL8RY0At9uNpwMSHdFnf3CISExIYP2qRaTa2mhzNwMQFR2Dq8vB1p37\ncbrCZ32rEKNBku8IKKuoJsLkY2LmBKNDGRaTycSCuYplsyfQ6a6lq6uTCLMZc0waH+4vpqDo/Kp6\nCREK2tvbePbZPxgdhiTfkdDc0k50pImoqNCadhiq5KQk1q9cQFKkp3cUbItN4kRlG7tzD9LV1WVw\nhOPP7t07ee21/zM0hqqqSq64Yi33338P9933ZW666Sbee28bVVWVrFmzjMOH+6+BvvPO23jsse+e\n9tyvfvVu7r//Hn7xiycG3Vdu7t4z3n8ufD7fadt8+umneO+9bfzxjxt6b9P6aG89jB579uzm7rtv\nZ8OG3/W7fceObXz601ezffsWwF9H43e/e5IXXniWjRvf4KOPdvH88+d2mip0zgaNI8Gs32uUiIgI\nFi+YTU1tHXm6FIs9BWt0DO7uaLbszGfZ/Cljskqaz+fDNcAUitXajdN5bmf2HY64876QJhhtf4Kh\nb/U0q7Wba6/9VNC7bcDgHTfOtdvG0LpofExKSmqgjsmefs8faheNYHbQkOQbZB6Ph7ZOE1ln6VQ8\nVqSlpnBpUiJ5h45R6wJbTBwRsWnsOlDK9Kxmpk/JMTrEYXG5nLy7+wTRpxSviY1poMXdPsizBtfq\ncXP5hdOJi4s/4+M2bnyDjz/eR3NzE0VFhdx111fYtOltioqKePjhRykqKqSg4CRTp04jPz+PpqZG\nSktLuPnm27j66muHta2cnMk88si3aW9vo62tjQceeJBZs+bw9NNPYbPZuOmmW3n66aeIjIzkllu+\nMGjM8fHxA3bbMJlMvd02BiuIfiandtz43Oc+33vfuXbb6Omi0bfS2YED+5k5cxbg76Kxb98eFi9e\ngtUaRWVlBS0tLcTGfpL4+3bRyMmZfFoXjfz8/AE7aBw6dGDY8YJMOwRdaUU1YCI7c/wUMDebzSxZ\nMJuF05Pxumro7u7GFpvEyWovO/ceOO0jXKiLtsdgj3H0+4qJjTvttqF8nZrEz6S8vIwf/vAn3Hrr\nF3nuuWf4wQ8e59Zbv8C77/qXtvcMngsKTvKDHzzOY4/9mFde+dOwtrVp09s0NDRw7bXX8bOf/Zq7\n776X5557GoBbb/0iW7dupqDgJLt2fchNN9162nb71tkqKysLercNGLzjxvl22/DHMHgXjYaGBhYv\nXkJMTAxNTY0Dtts6UxeNvXv3BrWDhox8g6y5pY3IiK4xV1t0KDLS00hNTiL3gKbRYyHKHktrt42t\nuw6wZO5kkhLPXig7nM2aNRvwd+6dNm06JpOJpKRkDhzoX1tj3jx/zdvU1HTcbvewt5WYmMgf/rCZ\nF154lo4OL9HR/qsYLRYLd931Ff7lX77Ej3/80wE7EPfttuFwxAS12wYM3HFj8+Z3uPzyK8+p20bf\nkeupBuuiMXnyFH70o58O8Ixz66Jxrh00JPkGmcvtDZl6DiPBYrGwfPFcyiuqOFRQRWRMCpGxaXx0\nsJypmc3MnDbZ6BBDVt9k1/f7U//I+9/n//e1115h8+Z3SUxMYuXK1Wfc1ksvvUhaWjoPPfQ9jh49\nwpNP/qz3/oaGeuLi4qipqR4wxr7ztj1rwYPVbQMG77hx+eVXnle3jR4j2UVj+fLlvPLK6SdFz7WD\nhiTfIGpra6O100TaGD/ZNhRZmRNIS01mX76m2WvF5kikuK6NusZ8li+aE1JX9o0VAw+0/Dded931\nXHfd9YB/zvdMmpubmT59OgA7dmylo8NfNLClpYWXX/4Tv/3t03zzm/excuXq00ZtQ+228cILr6D1\n0SE/F87ccSMvL5eWlpbebhtTpkwlKsrWO/I9lclk4sorr8bh6P8J89QuGkePHmblytUcPnyIpUuX\nnznAU7po/OAHj1NZWdG7zXnz5vV20OiZ9z2fDhryFxJEpRVVmC1WMtKSjA5lVERGRrJiyTyKS8s5\nUlhLlCOFNp+VrbsOcsGcSb2V1MTQDLxgYvirKK688mq+//3vsGXLZj7zmRvYvPkdNm58gyNHDnHT\nTbeQmJjI9dffxG9/+yseeODfhhDDJ8612wacuePGyy+/yOc/f+c5d9sYjS4a11//qaB20BjVThbB\nEMqdLPbsP0ptUxtXXbzQ0Bq+QxHs4+H1etmbr2nptGG12WlzNzE5zc6sGaefoDCS09k8yGqHqBFd\n7TAWSWGd/sZ6J4txzeluJ85uCfnEOxKsViurls6nqLiMo8V12OJSKG1sp35PPssXzSYyMjTmwR2O\nOC6/cPppt6ekOKirO/d1vkIMlyTfIGlvb6et00RKXGgkGaNMzplIeloy+w4co8MXizcyiS27DnLB\n7BxSU4yfhjCZTAOOUuPjHXi9svJSjB75bQuSsopKzJE2UpLG3xKz4YqOjuai5QuZkmahvaWeKEca\nubqKI/qk0aEJETIk+QZJg7ONLq+HCenpRocSMqZPyWHNkhmY2mqJiLRR7ozgvd15eL1eo0MTwnAy\n7RAkLnc7MVERAy5cD2cxdjtrVyzi2MkiCipdmGNSWLR4HnW1VQM+3n/dvTSxFOOfJN8g8Hq9eDog\nI3H8r+89VzOnTSYjrYV9B09QV1vFs6/uIjXeQnxSOlkpMdis/l/FwTptCDHeDCn5KqUmAt8ErsTf\nAgigCH+/tZ9orUtHJLoxoqyiEmu0g8Q4GfWeicMRyyUrFwHg7ezkj1vqcXrKSE+M5uHbl/UmYCHC\nwVnnfJVSdwDv4k+2nwFSgRTg+sBtbyulhlf/bZxpcLbS0dpCVsbYKp5upJzMFJwe/5VX1Y2tfHxE\nCrWL8DKUocZ8YL7W+tTSVYeBw0qp3wD/HfTIxhCn24st0jdmi6cbYWZOGhlJJVQ2eHBEm2lwtbM3\n/4jRYQkxaoaSfHcDNyuljmmtT7vIWmvtBb4R9MjGCK/XS6vXR2q8zPcOh81q4aEvLqW8zk202cuh\n42U0ef2FT5qam8ZkoXYhhmMoS80+Czw3UOJVSoX9UrXyiiqi7HHEx8iod7hsVgvTMuPJTE9l/coF\nxJn95RN3HyhDnyg0ODohRtZQRr7vaK0Hq6dwN/DrQe4LCw1OD+1tXWRlTDI6lDEjO3vSoKsaMjIy\n8XV3UVLfTU19HhcunoPVKp8qxPgzlOR7tVJqsKohVxHmydfp9hJp6h6XxdNHytnW8bo9HvbuP0ab\nOY6tuw+zYGYmGenjpzOIEDC0aQc7ED3IV1gXMujo6MDj7cYRLSOzYIqx21m7chFZ8d34TJHsP1HP\nx/lHBmz7IsRYNZSR78ta698OdIdSqi7I8Ywp5ZVVRNnjibWPrR5mY8UcNZ301AZyDxdT357Atl35\nLJ0/lTipIibGgaGMfNed4b7/DVYgY1FDk5vOjg4y05ONDmXcSk5K8p+Ms3joMNn5IK+I4wVFRocl\nxHkbSvJ9Wyn1BaXUaT04tNbD7xs9jjg9XnydHmkcOcLMZjPLF89l1sQYTD4fhTWdfLgnv7c9jhBj\n0VmnHbTWfxiNQMaazs5O3G3dxMeEZ/F0I+RkZ5GSnMje/cdw++LYsusQC2dmMkFOxokxaFgX0yul\nrgA+DXwXeAyYAPxUa/32CMQW0sorq4iKiccRLeURR1PPybjD+gQltZHknahjQk0DC+cpeRMUY8pw\nL5JQwKPAPcDXgH8Eht8zeRyob2qhu7ub1GQ5+WOEOWo6y+ZkYPa1U9saxbad+3G5WowOS4ghG27y\nLdZalwNvaq2dgYsv8kcgrpDn8nTQ0eYiPU0+8hql52RcfKQHr8nOB3kFFBRJgR4xNgy3ht+3lFK3\nAslKqXr8fa1zgL8FPbIQ1tXVhbuti1ibWYqnG6znZJy/fX0rx6vaqarLZ/miOVgsUqJShK7h/nZ+\nX2v9et8blFLXBDGeMaG8sorI6HgctuG3GhcjIyc7i9SUJPbkadzdcWzZeYBFs7JJS00xOjQhBjSs\naYdTE2/gtjeCF87YUN/UginCRGKc3ehQRB/26OjeK+O6TVY+PlZD/uFj+HyDlSYRwjjDSr5Kqamn\n/HymCzDGLafbS5u7maxMKZ4einpOxpm626hqsbBtZx4tbrfRYQnRz3BPuD3W841SahHwm+CGE/p6\n5nujI5FqWyEsOSmJS1ctJCGyFa8phvdzT1BUHNbdrkSIGe6cb6NS6nr83S3mAT8Lfkihraq6Gkt0\nHHFWme8NdaeejNPlrVTVHWDpwtlyMk4Ybrhzvl8BjgF2rfVngLBb11Pb4MJiiSQu1mZ0KGKIcrKz\nuHjZTGwmD85OO1t2HqC2rsHosESYO+vbv1LqX4AvAj0r2E2AQym1EVgAZI1YdCHI5emgtc1J1syc\nsz9YhIyek3H+K+Os7NNVPLh+JdVVlQM+Pjt70lnrDgtxPoby2SsZWK21Pq2KiVLqnwP/XjhQm6Hx\npru7G1drJ5FmE7GxsUaHI85B3zKV1VWVvPTXHSxbOJPJOenU1rp6HzdYpw0hgmUo0w4rgEeUUgPV\nTdyhlPov4KHghhWaqqqrsdgcOKLDuob8mNdzZRxAmy+aLbuPsWXXUdq8UpdZjJ6hjHz/Cfgm/jbx\nRUAp4MN/ZVs28D/AtSMUX0iprXdijYrFYZOOCmNdz5WJ0ybY+eXrVbh2N5AQc4Lv3bmCWLs0QxUj\nbyglJbuBHyulnsBfRCcb/7xvCbBHa901siGGDmdrB61eF1k56UaHIoIkwhaHq9X/K9zk7uL1bfms\nXzLZ2KBEWBjyeptAkt0V+Ao73d3dtLR2gq+DRCmePm5kpcSQkWynst5DkiOS2Jgoco/XAv7/84iI\n4S6FF2JoZLHjENXU1mK2xmI3h3XzjnHHZrXw0BeW4un0YbeY8LZ72HegAIAtH+7ngrmTpVOJGBGS\nfIeouq4Jqy2WuEiZ7x0vsrMnDbqqIX1CBp2+SHYfLCc7tY65aroUaxdBddbPVEqpBaf8HJbdIl2e\nDtrbWklNjjc6FBEk+/YdpKbGSU2NE5/P1/t9TY2TA/maFQsmYo3ooLwRtkuxdhFkQ5nQ+ukpP788\nEoGEMp/1FrSEAAAV60lEQVTPh9PTSVd7C+lpqUaHI0ZJQnwC61ctIjMB2rojeT+vQDoni6AZyrTD\nqZ+1gvrZSyl1AfBX4Hhg2/la668Fcx/nq6a2jgirneiIdjkBE2ZMJhPzZ88gs6GBvCPFnKzyUlmT\nx7KFiujoaKPDE2PYUJLvqcVQg10cNRZ4WWv9jSBvN2iqaxuxRcfgMMt8b7hKTkpi3coE8g8do9pp\nZceeY8yakkJOdlhdXS+CKBSGcQ6CPJoONqfHS2eHl6T4GKNDEQaKiIhg0fxZXKAmYI7o5nCxk517\nD9DZKVfGieEbysh3lVKqb/WytMDPJsCntZ50njHEAhcppd4EYoBHtNbbznObQePz+XC1dtLd3UxW\nxlyjwxEhIDUlifWrEsg7qKl2RrBl10Hmz8gkI12aqYqhG0ryVSMcQx7wXa31G0qpGcAmpdQ0rXVI\nDCdq6+qJiLQT5WslMlJqOgi/iIgILlgwm+qaWvKPlZN3vJby6noumKfkvIAYEtPZ+lsppS4+0/1a\n6x3BDEgptRu4QWtdPND9nZ1dPotl9DoG5+YfpcZtw2FxsXrZ/FHbrxg7urq62J17iKqmbiJN7axY\nNIPUlCSjwxKhYdAp1aGMfLcBR4GPgO5TNuYDziv5KqVuB2K11r9QSk0A0oDywR7f2Og5p/2kpjr6\nlQwcquKyJrzmOGLjfOf0/FB1rsdjvDrf4zFj8hRiq2s4cNzFxh2aiSlRzJs1ti/MkN+R/s7leKSm\nOga9byjJdy3+YuoXAW8Cz2mtc4cVwZm9CjwfaE9kBe4JlSkHn8+Hy+Ol29RM9qypZ3+CCGsZ6Wmk\nJiexL/8oJbVe6hr3s2T+VOIcUhtYnG4oVc3eA95TSkUDnwF+FBihvgA8P9j0wFBprZuAq89nGyOl\nvqEBU2QM1u5W7HZpEy/OzmKxcOEF8yirqOTQyWo+yCtiSoaDWdOnGB2aCDFDPjOgtW7VWj8HXAn8\nHPgGsG+kAgsFlTUN2OwxOOxyok0Mz8TMDNavmEuCPYIT5W527MrD09pqdFgihAw5+SqlZiul/gco\nAK4C7gYyRyqwUOB0e/H5fJJ8xTmJjIxk5ZJ5zJ+aiKcDduw5RlFxmdFhiRAxlAaaXwZux39y7Vlg\noda6caQDCwVOjxfMLWRNnWB0KGIMy8nOYkJaCvvyNQcK66mobWDZwtmydDHMDeWE22/w112oAG4A\nPqvUJ0t/tdbrRyY0YzU0NIDFDp0eEuITjA5HjHFRUVGsWraAwqISdEkjW3cfYt60DDIzpCtKuBpK\n8u05U+AjxC8DDqby6nqi7bFYOkJi4YUYJ6ZMnsSE9FT2HTjGvqOVlFfXccH8Wb095UT4GErybQS+\nBcwG3gN+EipLwUaSy+MFs504me8VQRYdHc1FyxdyorCYYyUNbNmZz4P3fY6KioGXt2dnT2LfvoOj\nHKUYaUM54fYk/hHv7/An4IdHNKIQ4fR00N7qIT1FphzEyJg+JYe1yxRWi4mKinI2bd1DcWk9O/NK\nKSlr6C3sXlpacvaNiTFnKCPfyVrrWwGUUhuBzSMbkvGamprwmaLo9rpJTZ1mdDhiHIux21m7YhEA\nhVVuNmz6AGdrNxnJdh76wlJsVun0NV4NZeTb0fNNoINxsOv5hpzyqlqiY+OIsZmlSIoYNVMnZeJs\n9deMrqz3UFYjbYvGs6FklpEuph5ynG7/+43M94rRNCMnhYxk/5WUDlsExwoKqaiqNjgqMVJCoZ5v\nyHF6vJhtXpKSYo0ORYSRnjb25XVukmLMHD5WwN7DFQB4vV6sVqvBEYpgCoV6viGl2dlMF1Y6Pc1k\nSvF0McpsVgvTMv0dslcvW0BxqX8FxOadh5ienciMqZMNjE4E01AK65xX4ZyxpqzSP9/ra2uQK5DE\nqMnOnkRa2sDVzyZOzCY1PoqjxU1U1Oxn8dwpUiltHJBTqadwub2YIqKJtctHPDF6hrKOt6a2jnxd\nyo59BUxKszNv1nQ5ITyGyf/cKZweL91dXSQ4bEaHIkQ/aakprF+1kJz0GEprPGz5cD/VNbVGhyXO\nkSTfPpwuJ52+SFpbmpko19yLEBQREcH82TNYs2QqkZYI9hwuZ3fuQTo6Os7+ZBFSJPn2UV5RS3Rs\nPFHmbqKjo40OR4hBxTniWLtiIWpSAg0tnWzeeYjCIrkSbiyR5NtHs7sdk8lErKzvFWPEjKmTWbd8\nFnF2MweLGnnvo/24XHJxxlggybcPl6cDn89HXEyU0aEIMWQ2m41VS+ezeGYanvZu3ss9wcGjxzlb\nZ3JhLEm+AS0tLjp8FtrcLrIyUowOR4hhm5iZwWWrFpCeGEVxdStbPsyjprbe6LDEICT5BpRV1BAd\nG4+pu534uHijwxHinJjNZpYsmM2F8yYCsOdIOXv3H6azc9xXgR1zJPkGNLu9mEwm6dcmxoWU5CTW\nr1rElHQ7Nc0dbNl1sPdqOREaJPkGuNxeAEm+YtwwmUzMUdNZs3gaURY4UNDAB3vyaXG7jQ5NIMkX\nALfbTVuXmbZWN+kpiUaHI0RQORyxrF2xiLk58bjaunl/3wkO65NyQs5gknyBssoa7I54ur0eUlOS\njQ5HiBExZfIk1l84h/gYM4XVrWzduZ/6hgajwwpbknyB5pY2IiIicERb5Fp5Ma5ZrVZWLpnH4pmp\ndHXD7oPl7Ms/QldXl9GhhR0prAM43V4sMcjFFSJsZE5IJz01hfxDx6hs7mTTh/nMmZrOdddeMWjP\nOGnkGVxhn3w9Hg/tnRH4vF6SU6R4uggfZrOZxQtmk9PYyP4jxRwoaKS0tITi4mpM5kg8nT7sFlNv\nH7nBSl6KcxP2ybe0oppoRwKtznoyJ8wzOhwhRl1SYiKXrExAHy8EYNNuzfZDLhpcHdLIcwSF/QSn\nMzDfa4+KwGKRXzARnkwmE7NmTgWgw9tJg8tfJa2y3kN5nSxNGwmSfD2B9b0xUjxdCIB/vGQRKXH+\nvwdHtIWTJwspKikzOKrxJ6yHem1tbbR2mLB3dZEQJ8XThQB/H7nvfWkF7g4flWWlVNY6OVLqr5RW\nVFLG5EkTDY5wfAjr5FtaUYU9NgFPSxMT5840OhwhQobNaiE7y0FyTCRzZ/o4Uehv5Xi4xEVh2X6m\nTEyWJHyewjr5NrnaiDA7sFl82Gwy8hXiTI08s7MnMTMrmuKKdg4VOyUJn6ewTr5OdztmuwOHNMsU\nAujfyDM11UFtreu0x0yb7KOwuITC8iYOSxI+Z2F7wq29vZ22DpO/eLokXyGGzGQyMXVyDpeuXsic\nSQ4ADhc3s/XD/VI5bRjCNvmWlldii42ntcVJphRPF+KcTM7JZt2qhczJ8dfAPljYJEl4iMJ22qHR\n1YbZ7CACrxRPF+I8TZ40kcmTJlJSVs7J0joOFDZSUFrH1OxkcrJlOmIgYTvydbnbAXBESz0HIYJl\n0sQs1q1cyKJpyZgi4EBBI1s/3E9JmYyETxWWI1+v14unA2KR4ulCjISJWRlMzMqgvLKKE8U17D9R\nz8mSOqZNSmHSxCyjwwsJYZl8yyoqiY5NoK3VzYSsJKPDEWLcysqYQFbGBCqrqjleXC1JuI+wTL4N\nzlbMZgft3gZSU2YYHY4Q417GhHQyJqRTXV3DseJq9h+vDfskHJbJ1+n2EhHtn+81mUxGhyNE2EhP\nTyM9PY3qmlqOF1f5k3BpHdOywy8Jh13y9Xq9tHp9xETLfK8QRklPSyU9LZWa2nqOFVXw8bGafkl4\nyZJ5476oe9gl3/KKKqJi4vF620lJdRgdjhBhLS01mbTUZGrr6jlWVMnHx6o5WVpHaWkJNTVO2ryd\nlNe5yUqJGXdF3cMu+TY4PVgsDtxNtUxIn2R0OEIIIDUlmdSUZOobGjhaUAHAO+/nsSm/hTqnd1wW\ndQ+7db5Ot79+b4xNiqcLEWqSk5JYvdTfUcbT7qPO6f97raz3sPPj47S2thoZXlCFVfLt6OjA3d4N\nQKzUcxAipF25Zj5pCVEAxNvNOF0e/r59PwAffXwEfaIQt3vsdtkIq6FfeWUVtph4uro6SYyPNjoc\nIcQZ2KwWHrnjwt4536hIM42NjQD4gNJqJ3lHiolzxJAUH0tCbBRZGWnExo6NRrghkXyVUk8AK4Bu\n4Ota670jsZ+GJjeWSAduZz1ZUjxdiJBns1qYlvlJ7ZWkJP9FURcung34P81WVtfQ2NRCbZOHg8f3\nEGGKYEJaCgkOG9mZacTGhuaJdcOTr1LqYmC61nqVUmo28Af8iTjonB4v2MBmluLpQoSysxV17xEZ\nGcmkiVn0LSXc0uKivKoOp7udbTvzaG9rJTk5mfTkBLIyU4lzhMZqCcOTL3Ap8BqA1vqIUipBKRWr\ntW4J5k46Oztxt3UTY0OKpwsR4s5nHW9srAM1vWe0O4vu7m5qa+uoaWhmd+5BPG43tqgoMjLSyZk4\n4YxVDUdyvXEoJN8JQN9phprAbSeCuZOSsgqiYuLx+XzSqViIMBIREdF7Zd382f5yAm1tbVRUVbP/\noMbjbsFiMZOUlMzUyRNJiE/ofW7f9caeTh92iylo641DIfmeen2vCf98+nkLh6tkhBDDZ7PZmDo5\nh6mTcwDw+Xw0NTVRUlbJUX0Ci8WCzeYfpLV5O3n0mb1U1nuCut44FJJvOf6Rbo9MoGqwBycm2rFY\nzEPacGlpCT6fj9b2Tt7YmkdUTCK+tgY+dfkyTCYTqWF+hVu4v/5TyfE4XTgdk7S0OGbO/GQ+ubOz\nE4BDJ4qprPcA/vXGnk4f2Vn+43I+xycUku87wCPA75VSi4FyrfWgi/caGz3D2nhpeWPvu1Z8TC3X\nLY/vbQo4UHPAcDFYc8RwJcfjdHJM/OZOzyEjubZ35Gu3mIacQ86UnA1PvlrrnUqpfUqpD4Au4N5g\nbr+8zt37rtXs9tLeKVXMhBBDZ7NaeOgLS0+b8z1fhidfAK31f47UtrNSYshItlNZ78ERbWHutLSR\n2pUQYpyyWS1kZwX3k0BIJN+R1POutevASZzNHjLTU40OSQgxRgx1vfG5GPfJF/wJOCvZhqW7Q4qn\nCyGGrO+KqGDPgY/r5DuS71pCCHE+xnXy7fuudVAfI9ERS1ZmpoERCSGEX9iUlOz0upmQnm50GEII\nAYRR8p01Yypm89AuzhBCiJEWNsk3O3vi2R8khBCjJGySrxBChBKTzxeUGjZCCCGGQUa+QghhAEm+\nQghhAEm+QghhAEm+QghhAEm+QghhAEm+QghhgHFd26GHUuoJ/O3ou4Gva633nuUp45pS6kfARYAZ\n+G+t9asGh2Q4pZQNOAR8V2v9R6PjMZJS6hbgQaADeEhr/ZbBIRlGKRUD/BFIAiKB72mt3wnGtsf9\nyFcpdTEwXWu9CrgL+KXBIRlKKXUJMCdwPK4CfmpsRCHjIaDO6CCMppRKAr4DrAKuAa4zNiLDfRE4\nqrVeB3wW+FmwNjzuky9wKfAagNb6CJCglIo1NiRDbcf/SwTQCNiVUmFd5FgppYBZwJtGxxICLgPe\n1Vp7tNbVWut7jA7IYHVAcuD7JKA2WBsOh+Q7gf4HrIb+3ZLDitbap7VuDfx4F/B3rXW4X+b4OPAN\nIKzfhAImAxFKqT8ppbYrpdYbHZCRtNZ/BnKUUseBbcC/Bmvb4ZB8T/2DMgHhnmxQSn0KuB24z+hY\njKSUug34UGtdHLgp3BOwCcgCPof/9+MPxoZjrMD8d7HWegb+T9FBm7YMh+RbTv+RbiZQZVAsIUEp\n9Q/A/wOu1FqHe2/wq4FPKaV2Al8Cvh3mo71q/G9G3VrrAsCllEoxOigDrQbeBtBa5wNZSqmg5M1w\nSL7vANcDKKUWA+Vaa7exIRlHKRUH/Ai4RmvdbHQ8RtNa36S1vlBrvRJ4CnhUa73F6LgM9A6wHiCQ\ndGO01uF8IvIE/pVSKKVyAJfWujsYGx73S8201juVUvuUUh8AXcC9RsdksBvxn0B4KXCizQd8Xmtd\nZmxYIhRorSuUUq8opbYC0YT5tBTwW2CDUmob/qWZdwdrw1JSUgghDBAO0w5CCBFyJPkKIYQBJPkK\nIYQBJPkKIYQBJPkKIYQBJPkKIYQBJPmKEaWUylFKdSulbj7l9sIhPPeWwL/pSqk/ByGW7mBdnRRs\nSqkXlFIZRschRk9I/iKKcecY8HCgNmqPMy4wV0qZ8Zc2JFBd68YgxBGyi9q11p/TWlcaHYcYPeP+\nCjcREiqBt/An03/ve4dSyo6/WHUi4ABe1lr/GPhf/NWk3sJ/VdH7wD8Af9Fazwo8dyKwC8jGXybz\nPqAdaAK+rLVuPCUOE/AfSqk1QCr+K/sOK6WuA/4NaMX/N3Gb1rpEKfU14BbADXiAW7XWjUqp+wL7\n6wCKgHu11u19XlMOsBn4O7AIf9K/SWtdqZS6B/h8IM424EattTPwSeBSYA3+OroJwBOBY/KvgAv/\nYOl2rXXRMI69CFEy8hWjwYe/bOPVSqkZp9yXBryqtb4Uf3eNbwXqLT8M1Gitr+zZhtb6MOBRSs0L\n3HYD8AL+Klz/CVyqtb4ceA/41iCxHNVaXwU8CXw3cFsCcEMgho18ckntd4GrA4W0fwpkKqWWAf+s\ntV6rtb4McOIvyHOqKcAGrfXF+GsofzNwuw24PLDNYuDWPseox0LgKq31xsDrujcQ278FXqsYB2Tk\nK0aF1rpTKfUg8HP8HTR6SjfWABcrpf4F8AJR+ItWD+ZF/IWSDuKvU/ElYCWQAbwdqFdhBQabU94U\n+PdDPkmINcAfA/PB6cDOwO1PBbb5Cv4R+fHAa5imlNoSeA32QNynqtda5wW+/wD4WuD7BmCjUqob\nyAEqArf3LWWZq7XuDHz/B+CZQAyvaq33DPK6xBgjyVeMGq31RqXUPYGP+T0jva8DVq31agCl1Nk6\nBbwI/F0p9TQQpbU+oJSaAuzWWl97luf68Pfxg0BdZ6WUBfgTsEhrXaCUuhdYEoj3X5VS2fjLTr6q\nlPpX/NMFf9Na33+WffX9VNmzryzgf4DZWut6pdSPB3lubzLXWv9MKfUCcCXwG6XUU1rr359l32IM\nkGkHMRr6juoeAH6Af4QL/pFmAYBS6lr8lbSi8CdJ26nbCFRfq8ff4PG5wH17gOVKqfTAdq5XSv3T\nIHFcGvj+IuAA/jnVCKAs0ETzU0CUUipBKfUwUKa1/g3+aYpl+EexV/WcPFRKfUUpdeEA+0pSSi3o\ns698/FMsrkDiTQKu6HMcTqOUilBK/QBwaq2fxT8NsmKwx4uxRZKvGA2985mBAt2v4E+6ABuAm5VS\nm/B/DH8+8FUOVCil9gAx9J8TfR64E/98L4FVAl8D3giU/rsD/4m4U3UCcwMn8b4MPBI4KfcMsBv/\nqPpH+OvZXgrEAnuUUu/iH/3+Xmu9D/gVsE0ptQNYC+wfYF/lwO1Kqc34m1H+BMgDDiildgG/wH8C\n8nal1GoGWIkRqBtbB3wYOD4P4B85i3FASkoKEWSB1Q7va62zjY5FhC4Z+QoxMmRUI85IRr5CCGEA\nGfkKIYQBJPkKIYQBJPkKIYQBJPkKIYQBJPkKIYQBJPkKIYQB/j8gMGSszLk13wAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f090fe6f690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5,3))\n",
    "\n",
    "# error bar is barely visible\n",
    "#plt.errorbar(rc, rna_pmf, yerr=dh_log_err, fmt=\".\", capthick=1,capsize=4\n",
    "rc = np.arange(0,9)\n",
    "rna_pmf = og*-1\n",
    "rna_pmf -= rna_pmf[-1]\n",
    "plt.plot(rc, rna_pmf, \".\", label=\"US DHAMed\")\n",
    "rna_pmf_obs_e7 = og_e7*-1\n",
    "rna_pmf_obs_e7 -= rna_pmf_obs_e7[-1]\n",
    "plt.plot(rc, rna_pmf_obs_e7, \"s\", mew=1, mfc=\"None\",\n",
    "        label=\"PMF $\\Delta t=5x10^7 \\, \\mathrm{MC}$ \")\n",
    "\n",
    "plt.fill_between(rc, min_pmf, max_pmf, alpha=0.4, label=r\"min-max PMF $\\Delta t=10^4 \\, \\mathrm{MC}$ \")\n",
    "\n",
    "ax.set_xlabel(\"Native base pairs\")\n",
    "ax.set_ylabel(\"PMF ($\\mathrm{k_BT}$)\")\n",
    "ax.set_ylim(-1,13)\n",
    "ax.set_xlim(-0.4,8.4)\n",
    "ax.set_yticks([0,5,10])\n",
    "ax.set_xticks([0,2,4,6,8])\n",
    "ax.tick_params(direction='out')\n",
    "\n",
    "plt.legend(borderaxespad=0.1)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# References\n",
    "\n",
    "[1] Dynamic Histogram Analysis To Determine Free Energies and Rates from biased \n",
    "Simulations, L. S. Stelzl, A. Kells, E. Rosta, G. Hummer, J. Chem. Theory Comput.,\n",
    "2017, http://pubs.acs.org/doi/abs/10.1021/acs.jctc.7b00373"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": true,
   "user_envs_cfg": false
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "12px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": true,
   "toc_position": {
    "height": "830px",
    "left": "0px",
    "right": "1932px",
    "top": "107px",
    "width": "212px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}