pmda-rmsf-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PMDA: RMSF example "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import MDAnalysis as mda\n",
    "from MDAnalysis.analysis import align\n",
    "from MDAnalysis.tests.datafiles import TPR, XTC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create RMS-superimposed trajectory\n",
    "\n",
    "We first need to fit the protein to a reference structure. We choose the mean structure as the reference (as opposed to the first frame of the trajectory)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = mda.Universe(TPR, XTC, in_memory=True)\n",
    "protein = u.select_atoms(\"protein\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need to center and make whole: this trajectory\n",
    "# contains the protein being split across periodic boundaries.\n",
    "# Fit to the initial frame to get a better average structure\n",
    "# (the trajectory is changed in memory)\n",
    "prealigner = align.AlignTraj(u, u, select=\"protein and name CA\",\n",
    "                             in_memory=True).run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ref = average structure\n",
    "ref_coordinates = u.trajectory.timeseries(asel=protein).mean(axis=1)\n",
    "# Make a reference structure (need to reshape into a\n",
    "# 1-frame \"trajectory\").\n",
    "ref = mda.Merge(protein).load_new(ref_coordinates[:, None, :],\n",
    "                                  order=\"afc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "aligner = align.AlignTraj(u, ref, select=\"protein and name CA\",\n",
    "                          in_memory=True).run()\n",
    "# need to write the trajectory to disk for PMDA 0.3.0 (see issue #15)\n",
    "with mda.Writer(\"rmsfit.xtc\", n_atoms=u.atoms.n_atoms) as W:\n",
    "    for ts in u.trajectory:\n",
    "        W.write(u.atoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the RMSF "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pmda.rmsf import RMSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/oliver/anaconda3/envs/pmda/lib/python3.6/site-packages/MDAnalysis/coordinates/XDR.py:195: UserWarning: Reload offsets from trajectory\n",
      " ctime or size or n_atoms did not match\n",
      "  warnings.warn(\"Reload offsets from trajectory\\n \"\n"
     ]
    }
   ],
   "source": [
    "u = mda.Universe(TPR, \"rmsfit.xtc\")\n",
    "protein = u.select_atoms(\"protein\")\n",
    "calphas = protein.select_atoms(\"name CA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmsfer = RMSF(calphas).run(n_blocks=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3.4120324 ,  2.80786034,  2.54288758,  2.13817685,  2.27161787,\n",
       "        2.05617269,  2.23263611,  2.26324985,  2.6685174 ,  3.34765102,\n",
       "        3.34649963,  3.67809519,  3.12530468,  3.04820756,  4.06848446,\n",
       "        4.28387404,  3.50432274,  3.65109967,  3.88320942,  3.57783562,\n",
       "        3.5559858 ,  4.1139274 ,  4.19462169,  3.99780905,  3.86728712,\n",
       "        3.57739603,  3.11901273,  2.71796965,  2.1907306 ,  2.01977274,\n",
       "        1.55588573,  1.9707523 ,  2.35237367,  2.05900165,  2.04981109,\n",
       "        2.69986978,  2.78669618,  2.51341295,  2.94705014,  3.39237424,\n",
       "        3.39615728,  3.47870721,  2.10949452,  1.63478931,  1.46610329,\n",
       "        2.08843289,  2.54029465,  2.16283894,  1.98612016,  3.00872536,\n",
       "        3.26406271,  2.78350623,  3.15157422,  4.09545959,  4.30529783,\n",
       "        4.05972764,  3.42932333,  2.75299481,  1.90037993,  1.56063035,\n",
       "        1.18161541,  1.20062227,  1.24029507,  1.20542892,  1.12901317,\n",
       "        1.24455601,  1.44645618,  1.80180757,  1.9404226 ,  1.97164966,\n",
       "        2.56329004,  3.04080912,  2.98180739,  3.0766044 ,  3.48753088,\n",
       "        3.64188448,  3.1457597 ,  3.58785793,  4.03127864,  3.31054317,\n",
       "        2.71834203,  2.45252201,  1.9749031 ,  1.97444495,  1.5650119 ,\n",
       "        1.18783697,  1.04361209,  1.00633681,  1.02301773,  1.28392244,\n",
       "        1.41744775,  1.17559251,  1.4332404 ,  1.79447553,  1.69970472,\n",
       "        1.62779012,  2.10749398,  2.31673317,  2.08195313,  2.48872005,\n",
       "        2.26479541,  2.66914862,  2.60074942,  3.02414758,  2.84180965,\n",
       "        2.37831783,  2.62804072,  2.44469618,  2.85970076,  2.81443738,\n",
       "        3.12594609,  3.10757295,  2.98239164,  3.31337701,  3.60417109,\n",
       "        3.38519409,  3.49163139,  3.95066449,  3.99368091,  3.87733486,\n",
       "        4.24801468,  4.91079151,  5.14602984,  5.87738503, 29.93689205,\n",
       "       35.44181935, 35.640867  , 32.98859581, 35.23634362, 22.65857828,\n",
       "        6.77082372,  6.24465234,  6.276922  ,  5.69802383,  5.75754529,\n",
       "        5.6902914 ,  5.9221854 ,  6.55007787,  6.70941175,  7.37004744,\n",
       "       34.16201472, 31.34996441, 31.74424566, 20.91568832, 20.66101542,\n",
       "       31.69069741, 34.00941606, 33.78754933, 20.3886708 ,  8.93876787,\n",
       "        8.39040492,  8.16684692, 21.38551723, 21.66163617, 36.28303378,\n",
       "       36.94809815, 30.82556538,  4.43303086,  4.45646918,  4.29140925,\n",
       "        3.82407131,  3.63712906,  3.54261973,  3.34451417,  2.89724049,\n",
       "        2.73667817,  2.61035716,  2.26938328,  1.90595577,  1.79565567,\n",
       "        1.66453859,  1.50176738,  1.21944205,  1.20492936,  1.3635034 ,\n",
       "        1.29675512,  0.98971101,  1.13806773,  1.6069187 ,  1.61052174,\n",
       "        1.73121472,  2.05321741,  2.34673194,  2.56467655,  2.71605102,\n",
       "        2.99706657,  3.26533731,  3.48917845,  3.65729177,  3.27499154,\n",
       "        3.0463187 ,  3.18937707,  2.67626899,  2.97969328,  2.9177627 ,\n",
       "        3.21926125,  3.38790406,  3.58477594,  3.84035115,  4.20389765,\n",
       "        4.67678676,  4.94923976,  5.23180414,  4.88679123,  4.41375615,\n",
       "        4.80224295,  5.03888067,  4.44825968,  4.17962361,  4.80101035,\n",
       "        4.69129912,  3.86170048,  3.95110073,  4.83380157])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rmsfer.rmsf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11848bef0>]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXmYHGW1/z+nunuW7NtkD2SDQBYSIIQdQtgDCCgIiAgCN6CAevXqD1QUuXr1iiwKV7wgS9jlQhDEAAmIQiAJmSxk38i+T/ZMJjM93fX+/qiunp6e3nu6Z2pyPs8zz/RUV1e9FJ1vnfq+55xXjDEoiqIo3sdq6QEoiqIozYMKuqIoShtBBV1RFKWNoIKuKIrSRlBBVxRFaSOooCuKorQRVNAVRVHaCCroiqIobQQVdEVRlDaCv5gn69Gjhxk4cGAxT6koiuJ55s6du9MYU5Fuv6IK+sCBA6msrCzmKRVFUTyPiKzPZD+1XBRFUdoIKuiKoihtBBV0RVGUNoIKuqIoShtBBV1RFKWNoIKuKIrSRlBBVxRFaSOooCuK0gRjDFPmbWJvTbClh6JkgQq6oihNmPnFLr7/6uf8fdHWlh6KkgUq6IqiNOEvlRsBqK4NtfBIlGxQQVcUpRF7a4K8s3gbADXBcAuPRskGFXRFURrxt4VbCYZsAA7Vq6B7iaI251IUpfWzaXcNpX6LdiU+aoJquXgJjdAVRWlEyDYEfBbtSvwcCtotPRwlC1TQFaXI7DtUz7qdB1t6GEkJhW18llBe4uNQvUboXiKtoItImYh8JiKfi8gSEflFZPt9IrJZRBZEfiYWfriK4m3mrt/DhQ9/xJcem4ExpqWHkxAnQpeI5aIeupfIxEOvAyYYY6pFJADMEJF3Iu89bIz5XeGGpyhti9uer2TXwSDGwN6aejqXBxABEWnpoUUJhY0ToQdU0L1G2gjdOFRH/gxEflpnaKEorZj6sM3O6iBH9ewAwI4Dddz/9lJufnZOC4+sMSHb4LecSdFDKuieIiMPXUR8IrIA2AFMN8bMjrx1p4gsFJGnRaRrwUapKG0AN9od2L09ADsO1LJg41427jnUksNqQsi2CfgcD12zXLxFRoJujAkbY8YA/YFxIjISeBwYAowBtgIPJvqsiEwSkUoRqayqqmqmYSuK93DFcVCPiKDvr2PL3kPUh1tXJknIdi0Xv0boHiOrLBdjzF7gn8BFxpjtEaG3gSeBcUk+84QxZqwxZmxFRdpFqxWlzXKwLhKhRwR9895D7DhQRyjcuhzMUNiOpC36qNHCIk+RSZZLhYh0ibwuB84DlotIn5jdrgQWF2aIitI2cCP0ig6ltC/xsXDTPoDWF6FHJkXVQ/cemWS59AEmi4gP5wbwqjHmbRF5XkTG4EyQrgNuK9wwFcX7uBF6u1IfPTuVsWDjXiB/Qa8P29QEw3QuD+Q9RohMivosykt81IVswhELRmn9pBV0Y8xC4PgE228oyIgUpY3iRujtS/z07FjK2khxUb6Wy7OfrOPPM9Yw655zk6Y/Lty0l35dyuneoTTt8UK2jT8SoYPTz6VDqXYJ8QJaKaooReJgxL5oH4nQXert/CL0quo6tu+v40Bd4owU2zZc+8QsLn10Bku37E97vFDY4LeE8hJHxDXTxTuooCtKkaiJCG67SITukm+EHradz+88UJfw/V0Hg9QEw2zbX8stk9PnvDuWi9AuEInQ1Uf3DCroilIkqiOC3r40TtBtk1cbADvy2aokgr51n5PnfnTPjmzfX5v2eG5hUXnEctFqUe+ggq4oRcIVxnYlPnp2auxl1+cRpduRCL2qOpmgOyI+oFs7bNOwfzKctEVRQfcgKuiKUiQOBkOU+C0CPoueHcsavRfKw0d39TmZ5bJ1rxOhH9GtHZDes3ezWlzLpVZz0T2DCrqiFImaujDtI1Gva7n0jkyO1oeawXJJEaGX+Cx6RZ4K0nn29WEbf6QfOmiE7iVU0BWlSBwMhqIiObBHe64ZO4CLR/UG8st0Se+h19K7cxkBn/PPPZ2gOx56rOWiWS5eQQVdUYpETV2Y9qWOSAZ8Fv991XEc3asjkF+mi3svSDUp2rtzGX6fk6Oe7ubhpC1aDXnoGqF7BhV0RSkSsRG6iz9SgZlPtWg4EqHvrA4mfH/L3lr6di7Db2UaoTcuLFLLxTuooCtKkagJNkToLiV+559gPoKeynIJ24bt+2vp06W8IUJPc65wJA+9PKZSVPEGKuiKUiQO1iWK0CNRc5pUwlS4aYg7q+uapCTuqq4jZBv6di4jEBH0dOeqj1SKlvgsLFEP3UuooCtKkagJNu2JkmnUnApXn0O2Yd+h+kbvbYnkoPfuXB5juWQSoVuICO1K/Gq5eAgVdEUpEjXBUNSXdglEBT33CD0cU2X6k78u4r63lkT/dnPQ+8RE6OnOVR+2o95+eYlP89A9hLZQU5QiUV0Xon1chN6QSph7hB7bNmDqom0Mrmgf/XtPjROxd+9QwrZItJ6uiMnt5QJOVatG6N5BI3RFKQJh21BbbzeJ0F0bJK8I3TbR6Ns5lt3kdYnPirF3kp/LGONYLpFxlQdU0L2ERuiKUgRie6HHEmgmD31g9/aM6NuJtTsPsn1/Q7aLe9xApOUApH4acCdMXctFVy3yFhqhK0oRiDbmiktb9Lsim0+lqG0o8Vs8cu3xjOrfudHNIRgboVvps1zcVrzuuJxJUc1y8Qoq6IpSBA7WpYvQ8+vlYkVWKgr4rKiIQ0OPmIDPiop0qqcB9z1X/MvUcvEUKuiKUgRiW+fGkml/lVSEDVgRAS7xWU08dJ8l+CxpyENPca6GCN3Zt9Rv5ZUjrxSXtIIuImUi8pmIfC4iS0TkF5Ht3URkuoisivzuWvjhKoo3ORizuEUszVH6b4zBXcM54LMIhhoLuivkDUVMqSL0xh66SPr+6UrrIZMIvQ6YYIwZDYwBLhKRU4C7gQ+MMUcBH0T+VhQlAeki9Lx6udgGX8RyKfFb2KYh0g6G7eg5MrF34j10nyXR1gJK6yetoBuH6sifgciPAS4HJke2TwauKMgIFaUNcDAysRhfKRq1XPIp/Y/z0KHhBlEftimJbMtkAtb9nC8SoVsijQqXlNZNRh66iPhEZAGwA5hujJkN9DLGbAWI/O6Z5LOTRKRSRCqrqqqaa9yK4ilq6pwIvTw+D7050hZtiLgp0SjcnRitD5moyDfYO8kF2r2xuMexRMgjAUcpMhkJujEmbIwZA/QHxonIyExPYIx5whgz1hgztqKiItdxKoqniU0fjCXQDIVFsRF6tHtjqCFCD/gbR++pJ0XdCN3Z1xLUcvEQWWW5GGP2Av8ELgK2i0gfgMjvHc0+OkVpI7ii6GajuLhim0/pf9iYqEXSYLk09dD90W6L6SdFA5HjqYfuLTLJcqkQkS6R1+XAecBy4C3gxshuNwJvFmqQiuJ13EwRd/LSpaH0P79KUcnAQ8/kaSB+UlREyGNoSpHJpPS/DzBZRHw4N4BXjTFvi8hM4FURuQXYAFxdwHEqiqdxNdSKE/RmKSyyDb5o2qLzoi5quZimEXoWhUU+Sy0XL5FW0I0xC4HjE2zfBZxbiEEpSlvDjdCtuGdiEcFvSX6l/zEeeqm/aYQezUPPYIGL+MIiS9Ry8RJaKaooRcBN/fPFeejgiGdelaK2iXrz8ZZLMBSTh56BveM+KTRKW9TCIs+ggq4oRSA6KSpNBT1gNe6/ki3G0KhSFOI89EjUblmCJZmV/rvHsUTQAN07qKArShGIWi4JBD3vCD1Blksw0pQr1kN3zmVRnyrLxW5cWOSz0AjdQ6igK0oRcAPwRJZLwGfl7aFLNA+9caFSrIcOTjpiqptHKJq22BChq4fuHVTQFaUIhKOWS9P3nIZauYumMQ3pkE089Jg8dHAi9FRZLuG4CN3SPHRPoYKuKEXAtp2OiJLMcskjQg/bpomHHgw1zUN33peUWS7RwqJIVO8TQR0X76CCrihFIDa1MJ6Az8rLQ7dNQ5aLOwGaqJcLOIVMmfVDbyj9Vw/dO6igK0oRCMeIbjx+S/JszhXTyyWu9D+2lws4VkrKSdG4wiJ3zEZtF0+ggq4oRcCO6VkeTyBulaGsj52Fhx5Ik1GTqLAodrvSulFBV5QiELYTZ7iA66Hnl7YY3z43trCoJH5SNGXaYuPCIve36rk3UEFXlCJgxywTF0++EbqJXeDC33RStLGHLqmbc0XG4aYtug8VmuniDVTQFaUI2Ck89IAvtcimI5zEQw/bBtsQZ7mkTlt0nxR8MVku7viV1o8KuqIUgXAKD93JPMnTQ29SKWpHo/7YSdF09k59gsIid/xK60cFXVGKQOoI3cq7fa57r/BZ4mSyhO2EqyQFrNT2TqLCImf8OQ9PKSIq6IpSBGy76eIWLoE8C4ts0zj6dywcO7oMXeNK0dRZLu6NJZq26HroquieQAVdUYpAOMWkqD/PCD0+xz3gc7o3usd0i42i50rTD92Shsi8IctFBd0LqKArShGw7RSWS76FRaYhGwUci6U+HOOhN7JcJPWKRbYdrRKFhlYFYRV0T6CCrihFILbFbTx5l/7b8ZaLRX3IRD302G6L6SyXcNhEF4iGBptI9dwbZLJI9AAR+VBElonIEhH5bmT7fSKyWUQWRH4mFn64iuJNYqs54/H78o3QG/eJCfglYrk0nRRN1w89ZDe+8bgvNcvFG2SySHQI+IExZp6IdATmisj0yHsPG2N+V7jhKUrbIDYTJZ58CouMcXLNE3roocarD0EG/dDtxoVIlnroniKTRaK3Alsjrw+IyDKgX6EHpihtibCdynLJvfTf1Vkr3kMPNaQtBuImRVMWFoXjI/SIoOf+AKEUkaw8dBEZCBwPzI5sulNEForI0yLStZnHpihthnCK9rn+PCL06OLTMccu8cdPisalNKYpLIqN0N2XOinqDTIWdBHpALwOfM8Ysx94HBgCjMGJ4B9M8rlJIlIpIpVVVVXNMGRF8R4m1aRopL9KLi1qo4tPx1ku9WGT2ENPU5Uatu3EEboKuifISNBFJIAj5i8aY6YAGGO2G2PCxhgbeBIYl+izxpgnjDFjjTFjKyoqmmvciuIpUlkubppgLhOPrhVixRUWBZOkLaYtLLJNtHVu7HG1sMgbZJLlIsBTwDJjzEMx2/vE7HYlsLj5h6cobYOwSbz8HDQIbi4+up1grVJ3kjWYaFI0TZZLOGyiVaLOcbX030tkkuVyOnADsEhEFkS2/Ri4TkTGAAZYB9xWkBEqShvAyRVP/J7rcQfDNmUBX1bHjXroMSJc4rMaNecqiW3OlUGWi99K4KGronuCTLJcZgCJvopTm384itI2SWm5RLbnUlxkIsG2xBcWJbVcLEK249cnemII2abRJKp66N5CK0UVpQikXCQ6klaYSwvdhiyXhm1OlotJWvoPye2dpGmLKuieQAVdUYpASkGPWBzBHAQ9WZZLMGQTDDf10N0J2GRPA6G4Xi66BJ23UEFXlCKQOssld8ulYVI0Ng+9cfvcRv3Q3TVHk0yMhuImRd3DqofuDVTQFaUIhOPK82OJRs05lGMmTluM89DjJkUh+c3DSVtsGqHnkiOvFB8VdEUpAsYkz3IpcaPmPCL0GA1uUliU2HJJfPMI23bCtEWN0L2BCrqiFIHYhZzjcdMEcyn/d4U2Pssl1kP3N/LXXcsl+aSo5qF7FxV0RSkC4RQLXPjziNBdJ6RRL5dIpWgwZFPisxqJvXvzSBahh5pUijq/NcvFG6igK0oRiF/3MxbXEskpQo9muTRsc5ecq60PN8oph/Q3j1A4vrBI0xa9hAq6ohQB25C2sCinXi4JslzcG8TBulCj1rnOuVJPwIbs+CwX9dC9hAq6ohSBVAtcRNMWc2rOlVzQa4LhRhOijc6VNEJvbLk0ZLlkPTSlBVBBV5QikGpNUTdqDueStuh66LETn5Go/GAw1CgHHRomRZNWisYVFukSdN5CBV1RikDYTu6h+/Lo5RKORugN29w0yJq6BB56JpOiWvrvWVTQFaUImJSFRYXx0KvrQkktl+SToqbRpKgKurdQQVeUIpAqQncj4lRLwyUjcem/66E3FfRAiqrUJz9aQ3VdiP5dy6PbtJeLt1BBV5QiEDamUWphLL7m9tDdLJdgOCruLslK/2et2cWvpi7jklF9uPG0gdHt6qF7CxV0RSkCdspK0fw9dGnkoUci9LpEk6KJc96/qKoG4N5Lhzdun6t56J5CBV1RioCdKsslDw/dJFixKDZCj23MFXuu+CyXULTVbuP91UP3FiroilIEUvVyiWa55CDo4QR56D07lUZfN5kUTdI3xv3bH7e/L1pYlPXQlBZABV1RioBtSNucK7csF+d37LGP7tWR+y4bDpA8Dz3O3nFvJvERuntYjdC9Qdo1RUVkAPAc0BuwgSeMMb8XkW7AX4CBOItEf9UYs6dwQ1UU7+IscJH4PTdCz6WXS0OWS+PtN50+iH5d21HRsbTR9mS91928dH/czG00y0UnRT1BJhF6CPiBMeZY4BTgDhEZDtwNfGCMOQr4IPK3oigJcLJckjXnyj8PPZE/f/7wXowZ0CXhuYJxEXp9Wg8966EpLUBaQTfGbDXGzIu8PgAsA/oBlwOTI7tNBq4o1CAVxeuYFN0Wm8NDl2SNYuIIJKkUDdk2PkuaHMcN2MNquXiCrDx0ERkIHA/MBnoZY7aCI/pAzySfmSQilSJSWVVVld9oFcWjZLLARW5ZLs7vZBk08bh9XuLtnfiFLVzcm5AuQecNMhZ0EekAvA58zxizP9PPGWOeMMaMNcaMraioyGWMiuJpjDHOpGgS0XU355flktn+gSSl//Vh0yQjxjmuts/1EhkJuogEcMT8RWPMlMjm7SLSJ/J+H2BHYYaoKN4mWs2ZJEIXEfyWJG2YlfrYTdMWUxFIkrbodFlsegxLS/89RVpBF8dUewpYZox5KOatt4AbI69vBN5s/uEpivdJtJBzPH6fNFtzrlRYlnPzaJqHbppkuDjHjZxHFd0TpE1bBE4HbgAWiciCyLYfA78BXhWRW4ANwNWFGaKieJuoLZLCF/FbVm4LXGTpoYNz84i3XEJhu0mGS+xxNQ/dG6QVdGPMDCDZt+Xc5h2OorQ9MomifVZuEXq2Hjo41aPBULzlYhJbLq6HroLuCbRSVFEKjCu6yTx0cBp0JVvnMxXRm0UWil7is5qcqz5sR/31WKxolkvWQ1NaABV0RSkw0fL8FKLrsySnbouu0GbqoYMTodeHmjbnShyhO781y8UbqKArSoGxoxF68n2cCD13yyVV9B9PwN90UjRk2wknRdVD9xYq6IpSYMIZ2CJ+n5VXlksWek7AsggmyHIJ+JvKgVs5qlku3kAFXVEKjG2nnxTNNUJP1cslGQGflTBCDyQ5hs8SzUP3CCroilJgMkktdLJccl+CLisP3d/Ur69P4qE7x9YsF6+ggq4oBcYVw1Q+d66Tog057pl/JuBrark4eeiJD2KJqIfuEVTQFaXA2HZ6n9vvy81yMVlWikIyyyVxcy732OqhewMVdEUpMNFMlAJUiiZagi4dgQSVoo7lklgO1EP3DiroilJgMpm49OfpoWeVtpgoQk9S+g/Ok4XmoXsDFXRFKTCZlv7n4qFH0xaz9NCb9HKxEzfncsemHro3UEFXlALjBsMpI/Q8uy1mE6GXJIjQ68OJ2+eCTop6CRV0RSkwmTTQ8lkW9Tl56ESOna2HnlkvF/fYObRqV1oAFXRFKTCZWC6BnD307NMW/T6L+vhui2ny0HUJOm+ggq4oBSaTSdGcPfScslwsgk2yXJLnoefa2lcpPiroilJgMlrgImcP3fmdnYfetFVv2jx01XNPoIKuKAUmsyyX3JpzhXNpzpXUcknioVvabdErqKArSoGJZrmkac5Vn4OHbozBkoauiJkQ8DdNW6y3k+eh+zTLxTNkskj00yKyQ0QWx2y7T0Q2i8iCyM/Ewg5TUbxLJhOXfksI59jLJRv/HJwJ2GDYjk50hm2DMSTNQ3eyXFTQvUAmEfqzwEUJtj9sjBkT+ZnavMNSlLZDwwIXqT30XBeJzmb5OSA6+emez01hTJrlYokuQecR0gq6MeYjYHcRxqIobZJMFrjINZPEjlgu2eAuZOFm1bjCnsxysbT03zPk46HfKSILI5ZM12YbkaK0MTJpoJVrcy7bNllluEBDhO620A25EXoKy0U9dG+Qq6A/DgwBxgBbgQeT7Sgik0SkUkQqq6qqcjydongXVwvTL3CRW5ZLth56SSQSd60Wd4I0eYSugu4VchJ0Y8x2Y0zYGGMDTwLjUuz7hDFmrDFmbEVFRa7jVBTPkslCzn6raTl+JpgcPHQ3PdE9n5uTru1zvU9Ogi4ifWL+vBJYnGxfRTncCWeS5ZJjYZGT5ZLdZ1zLpT4U8dAjEXrywiL10L2CP90OIvIyMB7oISKbgJ8D40VkDGCAdcBtBRyjoniaTMrzfREP3RiTVU65bUxWC0RDg7Xi5r27kXrSJei0fa5nSCvoxpjrEmx+qgBjUZQ2SSaLRLvRsW0giZWd9NjZ3ADAaZ8LsZZLJEJXD93zaKWoohSYcIYLXABNeqykI58sF9dyqU+T5eITIYciVqUFUEFXlAITLSzKIELPtuNiLnnobiTekLaYOstFpOGmpLRuVNAVpcBkssCFP656M+NjG5N1lktTyyV9lov2Q/cGKuiKUmAy6bboRujZZpMYk10vdGhaKRrNQ0/RPlezXLyBCrqiFJhMF7iA7D30sJ1LlktchO6mLabIcsmhb5jSAqigK0qByWiR6BwjdNuYrHqhQ4NX7nro9Xaa5ly6BJ1nUEFXlAKTySIUvjwmRXPOcomfFE2R5aKWizdQQVeUApNJ+9z4lraZHzsHD72JoKeO0EWXoPMMKuiKUmCy8dDD2XroOWS5RCtF3UnRNO1zfVbDTUlp3aigK0qByWiR6OikaLZZLtnnoTdJW0xXWKSl/55BBV1RCkxmi0Tn5qHnleUSis9ySW65aGGRN1BBV5QCk9Ei0b5cs1yy7+Xib2K5OAMsSVZYJLoEnVdQQVeUApPZItHupGiWvVyMyaqZF8RE6G63xVDqSlFtn+sdVNAVpcBktEh0Xr1c8mvOlbbbonronkEFXVEKTDirLJfsPfRss1x8luCLWSGpPk0euiWiWS4eQQVdUQqMK4apvG43Os46D92kbvqVjIBPMs5D92keumdQQVeUAmOb1NE5OCsWQQ6TojlkuYATjTeU/qdZgs7S9rleQQVdUQpMOIPyfFdMs10oOhcPHZyOi65fHwrb+C1J+gRhibbP9Qoq6IpSYGzbpMxwgdzTFsM5tM+FOMvFNkntFtD2uV4iraCLyNMiskNEFsds6yYi00VkVeR318IOU1G8S9hOH0UXs1IUnEyXqOUStpNOiIJbKZr9OZTik0mE/ixwUdy2u4EPjDFHAR9E/lYUJQHOws+F8dBzqRQFp4ioPmq5pI7QRbSXi1dIK+jGmI+A3XGbLwcmR15PBq5o5nEpSpvBzqCBVq4Rei6VouBYPNHSf9tOWlQEbpaLCroXyNVD72WM2QoQ+d0z2Y4iMklEKkWksqqqKsfTKYp3ySSKbujlkuWkqJ19P3RwLBe3KrU+bJIuPwfuikUq6F6g4JOixpgnjDFjjTFjKyoqCn06RWl1hDPwuXPPQ08/4ZoIx0OPyXJJEaFbmofuGXIV9O0i0gcg8ntH8w1JUdoWJoPUQn+ueeg5pi2W+Kyo5VKfNstFPXSvkKugvwXcGHl9I/Bm8wxHUdoeWVkuOVWK5pKH3rhSNH2Wiwq6F8gkbfFlYCYwTEQ2icgtwG+A80VkFXB+5G9FURIQzmCZuIZFonMpLMp+TH7LarSmaOosF8dy0eKi1o8/3Q7GmOuSvHVuM49FUdoktilchJ5Lcy5wPPTYJejSZbmAm36Z9amUIqKVoopSYJzCotT7RBeJzrJ9rsnRcunbpYw1O6vZVV0XsVxSe+iA2i4eQAVdUQpMJnno7tu5ROi5pC1+49SB1IVsnvlkXVrLxYraQSrorR0VdEUpMHYGzblEBL8luXnoOfwrHtqzAxeN6M3kmevYXROMPiEkwn0C0AC99aOCrigFJtPyfJ8lueWh5xChA9w5YSi19WFW76hO2jrXGZfzW4uLWj8q6IpSYMJ2ZuX5fksIZ70EXW4eOsCIvp154ZaT6VweoGu7kqT7WdFJURX01k7aLBdFUfLDGEMKRyOK32fl5qHnkrcY4eTB3fnXD8dT4k9vuWhxUetHBV1RCkwmC1yAE6GHcvDQcwzQo3RJEZ1DQ0ql6nnrRy0XRSkwmeaK+6zsF5LItTlXNrhD1yyX1o8KuqIUmEwnLv2WZJWHbtuG+hwLi7LBPb5WirZ+VNAVpcDYdvoFLgB8vuwi9I9X7yQYshnVr3M+w0uLezPSLJfWjwq6ohSYcIa54n4ru0nRVz7bQNd2AS4Y0SuP0aUntvRfad2ooCtKgbEzzETJZlK06kAd05du5ysn9KfU78t3iClxHy40y6X1o4KuKAUmnKGHfkS3dny2djeHguGU+322djdX/+lTbGO4dtwRzTXMpDRkuaigt3ZU0BWlwGRa/HP7+CHsrA7y4uz1SfeprQ9z+wtzCRvDczefzNCeHZpzqAmJeugaobd6VNAVpcBkarmcNLAbpw7uzsPTV3L2Ax/yf5Ubm+zzzuKt7D4Y5NdXHscZR/UoxHCbYGkeumdQQVeUAuO0z80stfAnlxzL6AFdCNuGh6avjC5C4fL8zPUM7tGe04Z0L8RQE6Ltc72DCrqiFBg7w9J/gJH9OvPSv53CL740gq37apm6aGvU6nh74RbmbdjL9accWfDc81h82svFM2jpv6IUkH2H6llTdZDTh2Znj5wzrCeDK9pzz5RFfP/Vz+nftZwNu2s44YguXHvSgAKNNjGiHrpnyEvQRWQdcAAIAyFjzNjmGJSitBXeXbyVYNjmS6P7ZvU5yxJ+esmxTP50PcN6d+SLHdWcNqQHP7t0OOUlhU1TjMcXrRQt6mmVHGiOCP0cY8zOZjiOorQ5/jp/C4N6tOe4/tlXc044phcTjils0VAmaC8XOFgXojzgK6rVlQvqoStKgVi9o5pZa3fLd2dHAAAZA0lEQVRx+Zi+GfVDb61El6A7TEP0xZv3ceqvP+DOl+c16WcTCmfXHbPQ5CvoBpgmInNFZFJzDEhR2gJLt+zn2idm0rk8wNVji+t5NzcNS9AdfoL+RVU1Nz79GSHbMHXRNp79dF30vWc/WcvI+95j5he7Uh7DGENNMFTgkTrkK+inG2NOAC4G7hCRs+J3EJFJIlIpIpVVVVV5nk5RvMHdUxbis4TXbj+Nfl3KW3o4eeGLToq28ECKzOa9h7jhz7MRgb/ddQbnHduTX09dzsbdNTzx0Rfc97el1NbbPP6vLwBHuO/962J+8bclrN91MHqcldurGXP/dP6xfHvBx5yXoBtjtkR+7wDeAMYl2OcJY8xYY8zYioqKfE6nKJ5gyZZ9LNy0j9vPHlKUSs5CczjmoYdtw6TnKjlQF+K5m09mSEUHfnnFKCwL/v0vC/jtuyuYOKo33zvvKD5aWcWq7QeYt2Evz89azzOfrOOyR2dwsM6Jyv+xfAfBkM2xfToVfNw5C7qItBeRju5r4AJgcXMNTMmfw3kSqyX5y5yNlPgtrjy+X0sPpVlwPfTvvDyfu16ez/b9tS08osLzxvzNLNmyn19dOYrhfR0h7t25jH87czCV6/fQtX0J/3XlKL5x6kBK/RYPv7+SyZ+uo2Opn/+94UT214b4eJXjSHy4fAfD+3SiT+fCP6nlE6H3AmaIyOfAZ8DfjTHvNs+wCs/Hq6q48OGPmL9hT8r91u48yJR5m9h9MFikkeVPXSjM3a8vZMz906JfKqU4HKit5435m7l4ZO+0S7t5hWP7dOLikb054YiuTFuyjQse/ojNew+19LAKxqFgmAenrWB0/85cdlyfRu/ddvYQLhjei0euGUOXdiV0a1/CXROGMnXRNt76fAtfObE/5x7Tky7tAry3ZDt7a4JUrt/NhGN6FmXsOactGmPWAKObcSxF493FW7njpfmEbcOLszdw/BFdE+732D9W8btpKwHo3amMh64ZzWlDitM/Ix/ufGk+05dup3enMm55tpJnbz6pxcYdCts8++k6rji+Hz06lLbIGArN+l0HWb2jmmP6dOKhaSs5WBfiljMGtfSwmo3O5QEe//qJAKzafoBLHp3Bg++t4KFrxrTwyArDr99ZxtZ9tfz+2uObZCd1KPXzxDcal9vccc5QDtSFeH7mem449Uj8Potzj+nF9KXbOGNoD2wDE44tjqAflmmLz89azxHd2jFxVG+mLdlGMGRTFwrzq78v5bW5mwCnwu9P/1rD+GEVPHfzOEoDFl97cjY3PDWbXdV1LfxfkJy9NUHeX7ad284ezHvfO4uKjqX8/v1VLTaeF2at55d/X8afP17bYmMoJC/OXs/ZD/yTWyZXctZvP+T1eZu4c8JRHNe/S0sPrSAc1asjt5wxiCnzN7Ng496WHk6zM33pdp6buZ5bzhjEuEHdMvqMiHDPxccy797zGVLhzJlcMKIX+2tD/OSvi+jRoYTRRfo+HBaC/n+VG5m7fjfg2BGV6/YwflgFV53Yn/21If66YDPffGYOT368lh++9jnTlmzjhVnrqa4L8cMLh3HW0RW8890z+fHEY/hs7W5ufnZOdMKjtTFj9U6MgQuG96ZzuwDXn3IEs9fuZk1VdVHOv2N/LR+u2IExhl3VdTw03XnCeXPB5jaxQEJNMMS+mnrAuXn+9t0VnDyoG3+ZdAq3nDGIr5zQn7smDG3hURaWb48fQs+OpVz/5Cxemr2B2vrU/duLTdg2TJm3ibU7D6bfOYaPVlZx50vzGNmvEz+6aFjW5y0LNFTwnnVUBacM7sbFI/vwzE3jMuq22Rx4QtDrQmGWbtmf02ffmL+JH762kO++soBQ2GbBhr3UhWxOHdydM4ZW0LHMz49eW8icdbv5rytHMapfZyY9P5cH3lvB2UdXMKKvU+HXrsTPpLOG8NjXTmDR5n3c+9fmm/+1bcOB2npWbT/As5+sZdOempyP9fHKnXQs8zM6Upl41Yn98VvCK3OatmJtbv4yZwMTHvwX33xmDjNW7+R301ZSEwzznXOPYuu+Wmav3V3wMRSSvTVBvvTYJ5z1wIe8u3gb9765hAO19dx/+UhOHtydH088lge/OppApp24PErHsgBv3HE6I/p25sdvLOKkX77PczPXtYo89fW7DvLlxz/l+69+zi3Pzkm7WAg46YbPz1rPrZMrGVzRgeduPjnvVaDKS3y8MulUHr5mDKNyqBLOFU8057rn9UX8a2UVn9w9odFdMBV//ngNby7YwsrtB+jXpZxNew7x9sKtrNt1EBE4eVB3SvwWd54zlCVb9vPv5x/NoB7tuWhkb6bM28Sug0GuOrF/k+OeP7wX3xo/hP/58AtuOPXIpP57ptSHba7535nM29Dw+LpuVw33fWlERp83xrBg415G9euMzxI+WlXFGUN74I+ISs+OZZw/vBevzd3EDy44umDLlS3evI97pizi5EHd2bC7hp+8sZiNe2q45fRBfOvsITz18RremL+JU4vY9rU5qQuF+bfnKtmwq4a+Xcq4/YW5gBOtDuvdsYVHV3z6dSnnlUmnMGvNLv700Rp+9uYSlm09wK+/PKrFxjR3/W5unVyJAe48ZyiPfbia/353ecp/S9v31/LTvy5m+tLtnH10Bb+/doynJ7M9IehfObE/U+Zv5u2FWxOKbDw1wRCPvL+Kio6lTBzVh3suPoavPzWbP/xjFSU+ixF9O9G5XQBwZq1j6da+hFvPHJzy+N8eP5RXKzdx/9tLef320/Lq7zD503XM27CXSWcNZkhFeyZ/uj7jpxFjDP81dRlPfryWuyYMZeKoPmzdV8t3zm2c73/duCN4Z/E2pi3ZzmVZNolKhdvf4lB9mHumLKJ7h1L+dMOJTFuyjR++tpAeHUr4znlHUV7i45Lj+vC3z7fyk4nDo9e+NbNpTw0vzNpAWcDi7KMreLVyI3PW7eHR647n7GEVTJm7iTOProh6pocjliWcNrQHpwzuzq+mLuOpGWu59Lg+WXeWbA7mrt/DDU99Rq9OZTxz00kM7NGe6roQz366DtsYtuytpTRg8dh1DROdG3bV8OXHP+FAbYgfTzyGW88Y3Op7taTDE4J+2pDuHNWzA898spavnNAvbV+MdxZto7ouxNM3nRSd2Ph/Fx3Dt16cRzBkc9vZqQU7He1L/fzowmH88LWFvDJnI187Obd1HXfsr+WR91dxzrAK7rn4GESExZv388Z8x29O9+Wa/Ok6nvx4LT06lPLnj9fyr5VVdCz1c27cjPoZQ3vQv2s5r8zZ0GyCvnL7AS59dAblAR/1YZuaYJg/Xn8CncsDXHl8Pz5bu5uJo/rQqcwR75tOG8SrlZt4ec4Gbo+7iRaTN+Zv4ld/X8aPLjqGq0/s3+S7tO9QPRj4xtOfsW7nQWwDj0Qmlb81fkj0+t10etvJYskXyxJ+eOEwpi3dxi/+toS37zqTEn/xbKclW/Zx0zOf0bNjKX+57RR6diwD4KeXHEvItnlu5noCPqE+bLj6xP7MWrObjbtrWLZ1P/Vhw9/uOoOje7WNpyxPCLqIcONpA/npXxfz/rIdnD88dQe6Vys3MrB7O04a2GCHnHtsL+bfez7zN+xl9ID8Pa2rTuzPG/M38+upyzjnmIqMigZq68O8tWALM1bvZNJZg3n6k7UEQzY/u2xEVFhG9O3E87PWs3FPDUd2b5/yeH9ftJWR/Trx2HUncN5D/2Lhpn08cs2Y6BfaxbKEa08awO+mrWTdzoMM7JH6uJnwm3eWU+q3uGhEbwJ+4crj+3Pikc719vssHri6cUbr8L6dOHVwdyZ/uo5bzhhUdJ+5Pmyz/1A99721lNr6MD96bSE79tdy54Sjovu8v3Q7tz5Xic8SBHjp305xFpyYvZ4d++v4jwuynyg7XCgL+Lj3kuFMen4ulz06g59fNpxTh3QveFOyldsP8I2nPqNjqZ8Xbj250Xff77P4z8tHMnFkH4b27MCVf/yUu16ez4HaEN3bl1BdF+K5m8e1GTEHj0yKAnzlhP4c26cTd740j0+/SN6td/Hmfcxeu5urxw5o8mVqX+rnjKN60LEs/0d+EeE3Xz6OetvmgfdWZPSZ+95awo9eX8jURVv52pOzmDJvM7eeOYhBMQLrVqUtSWC7rKmqZmckZbImGGLBxr2cMbSCgT3a89NLjuWuCUO5Ikl1onM9nFVv8mXWml38Y/kOvj1+KP991XH88opRUTFPxaSzBrN1Xy0/fWNxUTJe3ErZf67YwTH3vsu5D/2Lg3Uh3rzzdC4a0ZvHPlzNzC92cf2fZ/HsJ2u5983FDKlozzdOPZInbxzLKYO706HUmQz/6aXDi5ap4FUuGNGbP39jLNV1Ib7259lc+8Qs54mnANSHbR6ctoJLH52BCLxw68n079quyX4iji3Us1MZd00YyoHaEF8+vh+VPz2PRfddyMmDvTmnkwzPCHp5iY8XbhnHkd3bcedL8xNWbobCNv/v9YX06FDK108+suBjOqJ7O64/+UjeXLClUTMel/W7DkbHua/GqSD86tj+vP/9sykN+OjTuYw7zmmc4nZ0r474LGnio1cdqOPyxz7hq3+ayaGgk3pZHzbRScabTh/ED1JEkL06lTFmQBemL9uR7382j/1jNb06lfLN0wdm9blzjunJdyYM5S+VG/ndtMxugtlSEwxhjGHFtgOMuX8aP39zMT95YzH9upQzun8XfjzxWI7p3YmfXHIstg3XPTmLWWt2c9/flrJtfy0PXD2an182gnOGFacQpK1x3vBefPCDs7nvsuHMXb+He6YsbPbsl32H6vnmM3N49B+ruXhkb96+60wGZzCX8dWxA3jx1pP5zVeOQ0SKagsVC09YLi7dO5Ty2NdO4JI/fMx3X5lPKGwY1rsj3x4/hP/5cDVzN+xhyZb9jpdbpIm3284azPOz1vP791fx4FdHR58Ktuw9xCV/mEGXdgGmfPs0pi7cSl3I5hunDmRgj/a8972zCIVt2pc2/l9QFvAxtKIDS7bsa7T9wWkrqKkPs2bnQX79zjLal/rxW8LYDCJjl/OO7cUD761g+/5aenUqS/+BBCzftp8Zq3fyo4uGZZxxFMu/n380m/fW8uTHa7j2pCM4onvTqCpX3lm0lf/4v88ZO7Ab+w7VUxeymTxzPSLw2u2nNXqKGNCtHd8517m5PHXjSXy2djcicEKeWUuK8x2+6fRB1IZsfvPOcl6cvYGvn9I8AVZdKMytk+ewYONeHrjquKxaE1uWtMiEbTHxlKCDE8F+a/xQ/vDBKrq3L2Hmml08P2s9VuQf4w8vHMbFI3sXbTw9O5Vx02kDeeKjNew6GOT3146hc3mAn7yxiLBt2FUd5KrHZ3KoPsyofp0Z2c/x77u1T54aNaJvJ/65soq6UJhSv4/l2/bzl8qN3Hz6IGxjeOaTdZQFLEYP6NLkhpAKV9A/WLYjp4ncXdV1PP7PLygLWFx3Um4TwSLCjy4axt8XbeH+t5dy29mDOX5Al2iaZa68vXALd740n6N6duDjVVXOZOY1Y7CNIRiyE1pCd044ijvOGYqItCkftbUw6czBfPrFLu5/eyknHtk1726DxhjueX0Rc9bt4bGvHc+lxzVfxlZbwXOCDvCdCUMZN7AbJw3qyruLt/Hmgi388MJhRWlPmYi7LzqGfl3K+c+3l/Kbd5Zz8uBufLiiinsvHc7RvTrw0PSVHIoU2GTClSf0Y8r8zTw/cz23njmYZ2aso9RvcdeEoXQo9dOlvIRH/7Eq64Y/R/fqwIBu5bz1+WauG9d0jiEZy7ft57+mLuejlU6jrxtOOZKuKW5I6ejVqYxJZw3hDx+s4v1l27lkVB/+5/oTsjrGnoNBDtSGOKJ7O2rrw/x66nJG9uvE6986jVlrdrN4876MVgry8kpCrR3LEh766mgu/v3H3Dq5kv+48GguPa5vzhPif/znF0yZv5nvn3+0inkSpJjVXWPHjjWVlZVFO1+xuf9vS3n207W0L/VzVM8O/N/tp+U8kXbDU7NZuGkfb95xOhf9/iOuGNOP33zluOj7e2uCdCj1Zx3ZPj1jLfe/vZSfXnJs2nx7gI27a5j4h4/xWcI3TxvEiL6dOHtYRd5ZKsYY5m/cy7uLt/HER2v44/UnMHFUn/QfxOmGd9ljM1i9o5rR/TvTq1MZ05Zu58VbT27zj9ReZN6GPdz9+kJWbq+mX5dyvjV+CNeNOyKjfxurdxxgx4E6pi7ayguzNnD5mL48cs2Yw+5GLCJzjTFj0+6ngt587K0JcvYD/+RQMMzU757B0J65P8Yv3bKfLz02A8sSgiGbv3/njGgbgnwwxnDb83P5x/IdvPat0xgzIHnToGDI5ponZrJ6ezVTv3smA7o1n9/tEgrbXPnHT9my9xDTv392Uitqx4FaurUrwe+z+PEbi3j5sw3ccvogPlu3m0Wb9zH+6Aqe+WaT9VWUVoJtGz5csYPH//kFlev3MHpAF/736yfSu3PyuZx3F2/l9hfmAVDis7hqbH9+dunwnOZuvI4KegtRuW43NcEwZx2d/+pMn2/cyz1TFtGzUynPNqNY7TtUz8WPfERZiY+/33UmABt219C+1EefzuX4LGHL3kPc9fJ85q4vvF+5fNt+Lnt0BheP7MMfrju+0Xs7q+v49dTlvDF/EycN7MapQ7rzyPuruO2swdwz8Vjnv6emnrISq2BtDZTmwxjDW59v4SdvLKZP5zIm3zyOzXsP8cB7Kwj4hItH9uH6k4+gPmw476F/UR7w8fPLhjO0V4cm9RWHEyrobQhjTLM/Yn6yeifX/3k2ZQGL2vqGxSI7lvoZ2a8zc9fvwe8TfnvVcUXxK//wwSoemr6Sm08fxO3jB1PRoZSd1UG+9uQs1u+uYeLI3ry9cCsh23DZ6L48/NXReU+kKi3HzC92ceMznxEMOd+9Pp3L6FjmZ+X2au6++Bj2HAzyvx+tYfLN4zi7GYIjr6OCrqTlr5Ge1t3bl3BE93YcCob5fNNe5qzbw2lDunPrGYObNa0wFfVhm5+9uZhX5mzEGCgP+KgLhSnxWzxz0zhOHdKdGat2Mn/DHr41foiKeRtg8eZ9fLZ2NwGf8OUT+lMe8HHHS/N4Z/E2AK4Y05dHrj0+zVEOD1TQFU+yYtsBPl5VxdZ9tbQr8XHhiN7RVE+l7XMoGOZ301Zw6uDunHtsz8Nu8jMZmQq6J9MWlbbLsN4dD8t2tIpDeYmPey8d3tLD8Cx5PbeKyEUiskJEVovI3c01KEVRFCV7chZ0EfEB/wNcDAwHrhMRvbUqiqK0EPlE6OOA1caYNcaYIPAKcHnzDEtRFEXJlnwEvR8Qu1Dlpsi2RojIJBGpFJHKqqqqPE6nKIqipCIfQU80/dwkZcYY84QxZqwxZmxFheaTKoqiFIp8BH0TENu7sj+Q/+oJiqIoSk7kI+hzgKNEZJCIlADXAm81z7AURVGUbMk5D90YExKRO4H3AB/wtDFmSbONTFEURcmKolaKikgVsD7Hj/cAki8mquj1SY5em9To9UlOa7k2Rxpj0k5CFlXQ80FEKjMpfT1c0euTHL02qdHrkxyvXRvtcKQoitJGUEFXFEVpI3hJ0J9o6QG0cvT6JEevTWr0+iTHU9fGMx66oiiKkhovReiKoihKCjwh6NqmtzEisk5EFonIAhGpjGzrJiLTRWRV5HfXlh5nsRCRp0Vkh4gsjtmW9HqIyD2R79IKEbmwZUZdHJJcm/tEZHPk+7NARCbGvHc4XZsBIvKhiCwTkSUi8t3Idu9+d4wxrfoHp2jpC2AwUAJ8Dgxv6XG18DVZB/SI2/Zb4O7I67uB/27pcRbxepwFnAAsTnc9cFo9fw6UAoMi3y1fS/83FPna3Af8R4J9D7dr0wc4IfK6I7Aycg08+93xQoSubXoz43JgcuT1ZOCKFhxLUTHGfATsjtuc7HpcDrxijKkzxqwFVuN8x9okSa5NMg63a7PVGDMv8voAsAynY6xnvzteEPSM2vQeZhhgmojMFZFJkW29jDFbwfmiAj1bbHStg2TXQ79PDneKyMKIJeNaCofttRGRgcDxwGw8/N3xgqBn1Kb3MON0Y8wJOKtF3SEiZ7X0gDyEfp/gcWAIMAbYCjwY2X5YXhsR6QC8DnzPGLM/1a4JtrWq6+MFQdc2vXEYY7ZEfu8A3sB57NsuIn0AIr93tNwIWwXJrsdh/30yxmw3xoSNMTbwJA22wWF3bUQkgCPmLxpjpkQ2e/a74wVB1za9MYhIexHp6L4GLgAW41yTGyO73Qi82TIjbDUkux5vAdeKSKmIDAKOAj5rgfG1GK5YRbgS5/sDh9m1EREBngKWGWMeinnLs9+dnNvnFgujbXrj6QW84XwX8QMvGWPeFZE5wKsicguwAbi6BcdYVETkZWA80ENENgE/B35DguthjFkiIq8CS4EQcIcxJtwiAy8CSa7NeBEZg2MXrANug8Pv2gCnAzcAi0RkQWTbj/Hwd0crRRVFUdoIXrBcFEVRlAxQQVcURWkjqKAriqK0EVTQFUVR2ggq6IqiKG0EFXRFUZQ2ggq6oihKG0EFXVEUpY3w/wF96uBsijPY8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(calphas.resnums, rmsfer.rmsf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the spikes in residues between 100 and 200, which come from PBC artifacts: the proteins needs to be made whole before the fitting. This should become possible in MDAnalysis 0.20.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}