ideal_chain/ideal_chain.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "sixth-invasion",
   "metadata": {},
   "source": [
    "| <img src=\"randomly_coiled.png\" alt=\"randomly_coiled.png\" width=\"300\"/> | <img src=\"100ps.png\" alt=\"100ps.png\" width=\"300\"/> | \n",
    "|:-:|:-:|\n",
    "|Randomly coiled chains, `0 ns`| Chains after `0.1 ns`|"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "behavioral-constraint",
   "metadata": {},
   "source": [
    "### Example: Ideal chain\n",
    "We setup a system of non-interacting polymer chains, run a short simulation, and examine the resulting trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "assured-fortune",
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "from scipy.spatial.transform import Rotation as R\n",
    "\n",
    "\n",
    "def setup_system(box_size, n_molecules, chain_length, out_path, bond_length=0.5):\n",
    "    \"\"\"Setup a simple system of polymer chains\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    box_size : list of float\n",
    "        Simulation box size in nm.\n",
    "    n_molecules : int\n",
    "        Total number of molecules of `chain_length` monomer units.\n",
    "    chain_length : int\n",
    "        Length of each polymer chain in number of monomers.\n",
    "    out_path : str\n",
    "        Path of the created structure file.\n",
    "    bond_length : float, optional\n",
    "        Length of intramolecular bonds in the chains.\n",
    "    \"\"\"\n",
    "    box_size = np.asarray(box_size, dtype=np.float32)\n",
    "    n_particles = n_molecules * chain_length\n",
    "    \n",
    "    with h5py.File(out_path, \"w\") as out_file:\n",
    "        position_dataset = out_file.create_dataset(\n",
    "            \"coordinates\",\n",
    "            (1, n_particles, box_size.size,),\n",
    "            dtype=\"float32\",\n",
    "        )\n",
    "        types_dataset = out_file.create_dataset(\n",
    "            \"types\",\n",
    "            (n_particles,),\n",
    "            dtype=\"i\",\n",
    "        )\n",
    "        names_dataset = out_file.create_dataset(\n",
    "            \"names\",\n",
    "            (n_particles,),\n",
    "            dtype=\"S10\",\n",
    "        )\n",
    "        indices_dataset = out_file.create_dataset(\n",
    "            \"indices\",\n",
    "            (n_particles,),\n",
    "            dtype=\"i\",\n",
    "        )\n",
    "        molecules_dataset = out_file.create_dataset(\n",
    "            \"molecules\",\n",
    "            (n_particles,),\n",
    "            dtype=\"i\",\n",
    "        )\n",
    "        bonds_dataset = out_file.create_dataset(\n",
    "            \"bonds\",\n",
    "            (n_particles, 3,),\n",
    "            dtype=\"i\",\n",
    "        )\n",
    "    \n",
    "        head_positions = np.random.uniform(\n",
    "            low=np.zeros_like(box_size),\n",
    "            high=box_size,\n",
    "            size=(n_molecules, box_size.size,),\n",
    "        )\n",
    "        \n",
    "        # Randomly coiled chains\n",
    "        for i in range(n_molecules):\n",
    "            position_dataset[0, i*chain_length, :] = head_positions[i, :]\n",
    "            for j in range(1, chain_length):\n",
    "                position_dataset[0, i*chain_length + j, :] = (\n",
    "                    head_positions[i, :] \n",
    "                    + R.random().as_matrix().dot(np.array([\n",
    "                        bond_length, 0.0, 0.0,\n",
    "                    ]))\n",
    "                )\n",
    "                            \n",
    "        types_dataset[...] = 0\n",
    "        names_dataset[...] = np.string_(\"A\")\n",
    "        indices_dataset[...] = np.arange(n_particles)\n",
    "    \n",
    "        # molecules\n",
    "        molecules_dataset[...] = np.repeat(np.arange(n_molecules), chain_length)\n",
    "        \n",
    "        # bonds\n",
    "        bonds_dataset[...] = -1\n",
    "        bonds_dataset[:, 0] = np.arange(n_particles) - 1\n",
    "        bonds_dataset[::chain_length, 0] = -1\n",
    "        \n",
    "        bonds_dataset[:, 1] = np.arange(n_particles) + 1\n",
    "        bonds_dataset[chain_length-1::chain_length, 1] = -1\n",
    "\n",
    "# 15 chains in a 30nm x 30nm x 30nm simulation box.\n",
    "setup_system([30, 30, 30], 15, 10, \"ideal_chain.HDF5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "primary-abuse",
   "metadata": {},
   "source": [
    "### Run the simulation\n",
    "Disable the field interactions with `--disable-field`.\n",
    "\n",
    "```bash\n",
    "mpirun -n 6 python3 ${HYMD}           \\\n",
    "                    ideal_chain.toml  \\\n",
    "                    ideal_chain.HDF5  \\\n",
    "                    --disable-field   \\\n",
    "                    --seed 1          \\\n",
    "                    --verbose 1\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "promotional-lotus",
   "metadata": {},
   "source": [
    "### Examine the trajectory\n",
    "We calculate the average radius of gyration of the chains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "proof-fiber",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABwm0lEQVR4nO3dd3xb5fX48c+RvGM7cWInMXEGhJWQBEJSdtm0QCmUWQoUKC20UH7QAt/S0m9Z3UDp+KaUhlFGyy6FQGmZYQYCDhmQASSQYWfYTmzHjrd9fn/cK0cWHrIt6V7J5/166WXp6lo6VmL56HnOcx5RVYwxxhhjUknA6wCMMcYYY2LNEhxjjDHGpBxLcIwxxhiTcizBMcYYY0zKsQTHGGOMMSknzesAYqGwsFAnTZrkdRjGmG4sWrSoSlWLvI5jMOw9xhh/6u39JSUSnEmTJlFaWup1GMaYbojIOq9jGCx7jzHGn3p7f7EpKmOMMcakHEtwjDHGGJNyLMExxhhjTMqxBMcYY4wxKccSHGOMMcakHEtwjDHGGJNyLMExxhhjTMqxBMcYY4wxKccSHGOMMcakHEtwjDHGGJNyLMExxhhjTMqxBMcYY4wxKccSnBR052trWLCmqsuxBWuquPO1NR5FZIwxxiSWJTgpaEbJcC5/aHFnkrNgTRWXP7SYGSXDPY7MmNS1vamVix8o5YBfvsTFD5SyvanV65CMGdIswUlBh0wuZM45M7n07+/zq+dWcvlDi5lzzkwOmVzodWjGpKwrH17M/FUVVNQ189pHFVzz+FKvQzJmSLMEJ0VNLc6nvqmVua9/ynkHTrDkxpg4W1pWS1uHAtDSrixZX+NtQMYMcZbgpKg/vPQJ7Qq7FQ7j7wvXf64mxxgTW3uPzeu8nhEU9pswwrtgjDGW4KSiBWuqePCddQBkpQeZc87MLjU5xpjYu/aEvQHIz0rjyL1Hc9uZ+3ockTFDW5rXAZjYe2H5Fto7lLysNDZUN3TW5Cwrq7WpKmPiZHhWOgA3nrwPp+1f4nE0xhgbwUlRGWkBvnnQROqa2qhtbOWQyYV874jJXodlTMoKBgSAdrcOxxjjLUtwUkxzWztPLSnny/uMZdo4Z1l4WXWDx1EZk/pCCU6HWoJjjB9YgpNiXlpRQU1DK2fOKqGkIBuAsupGj6MyJvWFEpw2G8ExxhesBifFPFa6gV2GZ3Ho7oXUNjqNxizBMSb+AuKO4FiCY4wv2AhOCtlU28jrn1RyxqwSggGhICednIygTVEZkwBpVoNjjK9YgpNCnny/HFU4Y9Z4AESEkoJsG8ExJgECNkVljK9YgpMiVJXHSjdw0G4jmTAqp/P4+IIcS3CMSQArMjbGXyzBSRHvfraNdVsbOGv2+C7HnREcm6IyJt52TlF5HIgxBrAEJ2U8VlpGbmYaJ0wr7nK8pCDH6YXTYDsbGxNPoSLj9g7LcIzxA0twUkB9cxvPfbCJr+67C9kZwS73hZaKb7BRHGPiKmgjOMb4iiU4KeDfyzbS2NrOmbM/3x6+pMCpx7E6HGPiy81vaLcaHGN8wRKcFPBYaRm7j85l5vgRn7tvZ7M/G8ExJp5EhGBAbIrKGJ+wBCfJra6oZ9G6as6aXYK4NQDhRuSkMywjaCM4xiRAUMSmqIzxiYQmOCKSJSLvishSEVkuIjd1c86FIlIpIkvcy3cSGWOyeWJRGcGA8LWZ47q93+mFY0vFjUmEYEBsmbgxPpHorRqagaNVtV5E0oE3ReQ/qvpOxHmPqurlCY4t6bS1d/DP98s4aq/RjM7L6vG88SNtqbgxiRAMCG3tluAY4wcJHcFRR717M9292LvBAL32cSWVdc2c1U1xcbiSghzKqxtR+2RpTFwFxBr9GeMXCa/BEZGgiCwBKoAXVXVhN6edLiLLROQJERnfzf2IyCUiUioipZWVlfEM2bceK91AYW4GR+09utfzSgqyqWtuY3tjW4IiM2ZoSgsGbC8qY3wi4QmOqrar6n5ACXCAiEyLOOUZYJKqzgBeBO7v4XHmqupsVZ1dVFQU15j9aGt9My+vrODUmeNID/b+z2i9cIxJjICI7UVljE94topKVWuA+cDxEce3qmqze/NuYFaCQ0sK/1pcTluHcubsbge4utjZC8cSHGPiKRiADktwjPGFRK+iKhKREe71bOA4YFXEOeF7DZwMrExYgEkitLHmfuNHsOeYvD7P39kLx1ZSGRNPaYGANfozxicSvYqqGLhfRII4ydVjqvqsiNwMlKrqPOAKETkZaAO2ARcmOEbfW1ZWy8db6vnVqdOjOn94djq5mWmW4BgTZ4EAVoNjjE8kNMFR1WXAzG6OXx92/SfATxIZV7J5fNEGstIDnLRvcd8nE+qFY0vFjYk3p9GfJTjG+IF1Mk4yTa3tPL1kIydMKyY/Kz3q73MSHBvBMSaeggGxKSpjfMISnCTz/PLN1DW1dbuxZm9C3YytF44x8RMMCO3W6M8YX7AEJ8k8VrqB8SOzOWjXUf36vpKCbOqb26htbI1TZMaYgNgIjjF+YQlOEtmwrYEFa7Zyxv7jCQQ+v7Fmb0JLxTdss2kqY+IlLSi2TNwYn7AEJ4n88/0yAE6f1f3Gmr3ZuVTcCo2NiZegNfozxjcswUkSHR3K46VlHLZ7YedoTH+M72z2ZyM4ZmgQkSwReVdElorIchG5yT2+q4gsFJHVIvKoiGTE6jkDtpu4Mb5hCU6SePvTrZTXNEbVubg7+dlp5GWm2QiOGUqagaNVdV9gP+B4ETkI+C3we1XdHagGvh2rJ0wL2DJxY/zCEpwk8XjpBvKz0vjS1DED+n4RYZwtFTdDiDrq3Zvp7kWBo4En3OP3A1+L1XPaXlTG+IclOEmgtrGV/3y4mVP2G0dWenDAjxNaKm7MUCEiQRFZAlTgbN67BqhR1Tb3lDKg26I2EblEREpFpLSysjKq5wsGrMjYGL+wBCcJPLN0I81tHZw1wOmpkPEjnW7G1gvHDBWq2q6q+wElwAHA3v343rmqOltVZxcVFUX1Pdbozxj/sAQnCTxeuoG9x+YxbVz+oB6npCCHHS3t1DRYLxwztKhqDTAfOBgYISKhbWpKgPJYPU/QanCM8Y1Eb7ZponTna2uYUTKcUcMyWVpWy89Omsrbn25lWVkt3zti8oAeM3xX8YJhMVs4Ysygicgt/ThdVfXaKB6zCGhV1RoRyQaOwykwng+cATwCXAA8PYCQu2V7URnjH5bg+NSMkuFc/tBiDtp1JOlBYdyILC5/aDFzzvncXqVRCyU4G6obmF4yPFahGhMLZ/bjXAX6THCAYuB+EQnijFY/pqrPisgK4BER+QWwGLin39H2wEZwjPEPS3B86pDJhcz5xkzOvXshuxUN47p/fcicc2ZyyOTCAT9mSWcvHFsqbvxFVXeNw2MuAz73iUBVP8Wpx4k5S3CM8Q+rwfGxGeNHoMCayh2cd+CEQSU3AMOz08nLSrOVVMbEScCKjI3xDUtwfOzFFZsBOHbKGP6+cD0L1lQN+jFtqbhJBiKym4j8RUQ+EJFy9+sdIrKb17H1Js2WiRvjG5bg+NSCNVXc8PRyAM47aAJzzpnJ5Q8tHnSSU1KQbVNUxtdEZBawBDgdeA94wP16OrBYRPb3Lrre2V5UxviHJTg+tayslm8f5nxYLczNdGpyzpnJsrLaQT3ueHcEx3rhGB+7Daf4d5KqXqSqP1HVi4Bd3eO3eRpdLwI2gmOMb1iC41PfO2Iyo3KdpdxFeZmAU3g80CXiISUF2TS0tFNtvXCMfx0A3KKqXYYa3du3AQd6ElUU0qwGxxjfsATHx6rqmwEYGcOeNTt74dg0lfGtRmBUD/eNBJoSGEu/BGwVlTG+YQmOj1XVN1OQk056MHb/TDuXiluhsfGtfwO/EZHDwg+6t38NPONJVFGwRn/G+If1wfGxqroWCnMzY/qY40LN/rbZCI7xratwugu/JiIVOBtljnYvbwNXexhbr6wPjjH+YQmOj1XVN8c8wRmenU6+9cIxPqaqW4HDROR44As4HYk3AQtV9QVPg+uDJTjG+IclOD5WVd/M9JIRMX9cpxeOjeAYf1PV/wL/9TqO/rDdxI3xD0twfKyqvoXC3NhvillSkM1nVTti/rjGxJKIZAG7AFmR96nqisRH1LdgQOjo8DoKYwxYguNbTa3t1De3xXyKCmD8yBze+KQKVUVEYv74xgyGiJQAc4Evd3c3zmabwYQGFSWn0Z9lOMb4gSU4PlVZ5ywRj9cITmNrO9t2tDAqDgmUMYP0ILAbcDmwGmjxNpzoBQJCh2IfHozxAUtwfCrUAyceIzjhS8UtwTE+NBs4V1XneR1If6UFnKSmQyFo+Y0xnrI+OD5VVe98aI1PghNq9mcrqYwvrQByvA5iIIJugmPTVMZ4zxIcn+ocwcmLfYIzzroZG3/7f8C1InKo14H0V8CdlrL8xhjv2RSVT1W5NTijYrhNQ0h+VjrDs9PZYAmO8aclwLvA6yLSAtRFnqCqoxMdVDRCU1S2VNwY71mC41NV9c3kZaWRlR6fxSIlBdk2RWX86m7gTOAJkrDIGKC93RIcY7xmCY5PVe1ooSiOBcAlBdmsqbReOMaXTgV+qKp3eh1If4UKi20ExxjvWQ2OT1XVxX6bhnDj3W7Gam/Exn8qgfVeBzEQQXdjXNuuwRjvWYLjU1X1zRTmxb7+JqSkIJum1g627kia0X8zdNwMXCMiuV4H0l9Bt8jYEhxjvGdTVD5VVd/CoXGdotrZCyeeI0XGDMBXgD2A9SJSCtRE3K+q+vWERxUFdwDHpqiM8QFLcHyopa2D2sbWuCYeJSN3LhXfb/yIuD2PMQNQiFNcDJAOFHkYS78EA06G02EjOMZ4LqEJjrt53utApvvcT6jqDRHnZAIPALOArcDXVXVtIuP02tYd8etiHDJuhDX7M/6kqkd5HcNAhUZw2izBMcZzia7BaQaOVtV9gf2A40XkoIhzvg1Uq+ruwO+B3yY2RO9V1YW6GMevBicvK50ROenW7M/4goi8IyI/FZGZXscyGAGrwTHGNxKa4Kij3r2Z7l4i3wlOAe53rz8BHCNDbNe6eHYxDldSkM2GbTaCY3zhSeA4YKGIlIvIXSLyNREZ5nVg/ZEWmqKyGhxjPJfwVVQiEhSRJUAF8KKqLow4ZRywAUBV24BaYFQ3j3OJiJSKSGllZWWco06sSjfBiWcfHICSETk2gmN8QVVvUdUjceptrgIygLnAVhF5QUSuFJHdvYwxGp1TVNbozxjPJTzBUdV2Vd0PKAEOEJFpA3ycuao6W1VnFxUlTQ1iVOK5k3i4UDdj64Vj/EJVa1X1UVW9ABgDHAm8A5wPfCwiH4nI7SJytJdx9qRzLyr7nTLGc571wVHVGmA+cHzEXeXAeAARSQOG4xQbDxlVdS0MywiSnRGfbRpCxo/Mobmto3PncmP8xJ3SfkdVr1fVWTgfim4FJgH/8jS4HqQFrQbHGL9IaIIjIkUiMsK9no0z574q4rR5wAXu9TOAV3SIDTFU1TczKgG9aUpsV3GTRFR1o6reraqn4Swl953QCI6tojLGe4nug1MM3C8iQZzk6jFVfVZEbgZKVXUecA/woIisBrYBZyc4Rs9V1TfHdQVVSHizv5kTCuL+fMb0REQui/JUxdl8c52IvKmqTXEMq9+CAZuiMsYvEprgqOoy4HPLQFX1+rDrTTg7CQ9ZVfXNTBoV/8Uj4wqsF47xjTkD+J4qEfmSqi6JdTADFUpwbIrKGO/ZXlQ+VFXfEvcl4gC5mWkUWC8c4wOqGoj2AuQAM4APgP/zNvKubC8qY/zDEhyfaWvvoLqhJWH7Q5UU5LDBRnCMj7gdz3szSlU/BH5DNyPCXrIRHGP8wxIcn9m2owVVKEpADQ6ElorbCI7xlXki0u0vgIhMBt50b76DszGnb3QmOFaDY4znLMHxmdCS7cSN4GRTbr1wjL9MAJ5w20R0EpF9gDeAjwBUtU5VX/Mgvh51JjjW6M8Yz1mC4zOJ2qYhJNQLJ9Q92RgfOAaYCjwkIgEAEZkNvAa8C5zsYWy96tyLyj4wGOM5S3B8JlFdjENKbCWV8RlVLQeOBQ7EaStxJPAy8Dxwuqr6tjNlqNFfh9XgGOM5S3B8ZmeCk6ganJ29cIzxC1VdCxztXl4GHlbVc1W13dPA+hC0Rn/G+EaiG/2ZPlTVt5CZFiA3MzH/NONGWDdj4z0RuaWHuxYBBwPbw85RVb02MZH1T8Aa/RnjG5bg+ExVXTOFuZmI+0kw3oZlpjFyWIaN4Biv9dbcsz7ifgV8meCk2TJxY3zDEhyfqaxvTliBcUhoV3FjvKKqu3odQyzYXlTG+IfV4PhMVX1LwnrghJQUZFO2zaaojBmszr2oLMExxnOW4PiMs9FmokdwciirabQ3ZWMGKc0a/RnjG5bg+EhHh7JtR+K2aQgZX5BNS1tH5wouY1KBiIwXkfkiskJElovIle7xG0WkXESWuJcTY/WcAavBMcY3rAbHR6obWmjvUEYlfIrKWSq+obqR0fl9bQNkTNJoA65W1fdFJA9YJCIvuvf9XlVvi/UT2mabxviHjeD4SKK3aQjZ2ezP6nBM6lDVTar6vnu9DlgJjIvncwaDluAY4xeW4PhIorsYh4yzbsYmxYnIJJydxxe6hy4XkWUicq+IFPTwPZeISKmIlFZWVkb1PDaCY4x/WILjI6EEpygvsVNUORlpjLJeOCYJiMj1IvJdEcntx/fkAv8EfqCq24G/AJOB/YBNwO+6+z5Vnauqs1V1dlFRUVTPZbuJG+MfluD4SGWdNyM4EOqFY1NUxvduxElQNojIrX2dLCLpOMnNP1T1SQBV3aKq7araAdwFHBCr4GyZuDH+YQmOj1TVt5AeFIZnpyf8uUsKcii3ERzjc6oaAIYBpwHVvZ0rTjvwe4CVqnp72PHisNNOBT6MVXy2F5Ux/mGrqHykqr6ZUcMSt01DuJKCbF5csYWODu1c6mqMH6lqIzDfvfTmUOCbwAcissQ9dh3wDRHZD2fLh7XAd2MVW8BGcIzxDUtwfGRrfTOFCa6/CSkZmUNLeweV9c2MsaXixgdEZCowCxgP3Kuqm0Vkd2CLuyqqV6r6JtBdtv5cbCPtKi0gVoNjjA9YguMjVfWJb/IXEr5U3BIc4yW3KPhe4HScXjZpwH+BzcCvgPXANZ4F2IdAQGyKyhgfsBocH/Fim4aQ8bZU3PjH7cAhwLFAHl1HYZ4DjvciqGgFRWyKyhgfsBEcn1BVtno4gjNuhNPN2BIc4wOnAVeq6nwRCUbctw6Y6EFMUUsLCO0dXkdhjLERHJ/Y3thGS3sHhQnepiEkOyNIYW6GLRU3fpANbO3hvjygPYGx9FsgILR3WIZjjNcswfGJys4mf96M4ACMK8ixERzjB+8B5/dw3xnAggTG0m9BKzI2xhdsisonvNqmIVxJQTYrNm737PmNcf0MeFFEXgIex1nOfaKI/BAnwTncy+D6ErQpKmN8wUZwfMIvCU55daMVSBpPqeobwDFAJjAHp8j4JmA34FhVfc/D8PoUFJuiMsYPbATHJ6o6t2nwpgYHYHyB0wunoq6ZscNtqbjxjqq+BXxRRLKBAqBGVZOiQMxGcIzxBxvB8Ymq+hYCAgU53iU44b1wjPEDVW1U1Y3JktyAk+B0WA2OMZ6LagRHRL4K/NvdnM7EQVV9MyOHZXq6TUJJwc6l4rMneRaGGYJE5JZ+nK6qem3cghmkoDX6M8YXop2iegrYIiIPAvep6sr4hTQ0OU3+vBu9ARvBMZ46sx/nKuDbBCcgtheVMX4QbYIzGfgWztLNa0TkXZxW6o+qqi27iYHK+hZPl4gDZKUHKczNtKXiJuFUdVevY4iVtECAdktwjPFcVDU4qrpWVW9w34SOA1YDvwc2iciDInJUPIMcCqrqvNumIVxJQbYlOMYMgu1FZYw/9LvIWFVfUdVvAnsCi4BzgZdE5FMR+aGI2MqsflJVX0xRQSjBsSkq4x0R+aWI/LWH++4UkZ8nOqb+CAawImNjfKDfCY6IHCEi9wEfAdOAPwNfAp7A6VXxQCwDHArqm9tobuvwyQhODuU11gvHeOobwBs93PcGcE4CY+m3oE1RGeML0a6imghc4F4mAa8ClwBPqmqze9rLIvI28PfYh5naqupbAG+b/IWMH5lNa7taLxzjpV2A8h7u2+je71tBwRIcY3wg2hGcT4GLgYeA3VX1GFV9OCy5CVkOvNvTg4jIeBGZLyIrRGS5iFzZzTlHikitiCxxL9dH+8Mkq62hLsYeFxnDzqXiG2yaynhnM7B/D/ftD1QmMJZ+cxr9WYJjjNeirZc5CXi+rz44qvox0FvBcRtwtaq+LyJ5wCIReVFVV0Sc94aqnhRlbElv5zYN/qjBAWep+BcmjfQ4GjNEPQZcLyKrVPXfoYMiciLOPlVzPYssCrbZpjH+EO0IzpnAxO7uEJGJInJvNA+iqptU9X33eh2wEhgXZQwpq9KdoirywRTVuBFugrPNVlIZz1wPLASeEZFKEVkmIpXAM8DbOEmOb9kIjjH+EG2CcwFQ1MN9he79/SIik4CZOG9kkQ4WkaUi8h8R2aeH779EREpFpLSy0tcj1n2qqmtGBEYO834EJys9SFGe9cIx3lHVJlX9EnACcA/Oe8Q9wPGqekI3U+O+EhBLcIzxg2inqASne2h3ptHPOXERyQX+Cfygm0aB7wMTVbXeHZJ+Ctgj8jFUdS7uUPXs2bOT+t2kqr6ZgpwM0oL+2BqspCCbshqrwTHeUtXngee9jqO/0mwvKmN8occExy0ADhUBK/CUiER+csoCxgD3RfuEIpKOk9z8Q1WfjLw/POFR1edE5A4RKVTVqmifI9n4pQdOSElBDsvKarwOwwwhIpIT2lBTRHL6Ot/Pm28GA0JbuyU4xnittxGcFTiJiABXAfOBTRHntACrcIoC+yQigjPUvFJVb+/hnLHAFlVVETkAZxptazSPn6yq6lt8sUQ8pKQgm/9+uIn2DiXo4eafZkipE5GDVfVdoJ6eR4xDggmIaUACYiM4xvhBjwmOqr4IvAggInXA3araU2+KaB0KfBP4QESWuMeuAya4z3kncAZwqYi0AY3A2aqp/W5RVd/MviUjvA6jU0lBqBdOE8XDs70OxwwNFwFrwq4n7e98WtBqcIzxg6hqcFT1plg8maq+iTMi1Ns5c4A5sXi+ZOGXfagA7nxtTef1DdsaKR6ezYI1VSwrq+V7R0z2MDKTylT1/rDr93kYyqBZkbEx/tBbDc5jwE9UdY17vTeqql+PbWhDQ2NLOzta2hnlkxqcGSXDufTv7wNOL5y2jg4uf2gxc86Z6XFkZqgQkU+BU1V1aTf3TQPmqepuiY8sOtYHxxh/6G0EpwhId6+PJomHjP0s1OTPDz1wAA6ZXMgfz96PC//2Ho++t4FPKuqZc85MDplc6HVoZuiYBPT0C5EDlCQulP6zPjjG+ENvNThHhV0/MiHRDEGVnds0+GMEB+DIvUaTl5nGws+2ccXRu1tyY+JORPKBEWGHxorIhIjTsoCz6XmfKl8I2hSVMb4QbR+cbonICFWtiVEsQ1JVXWibBn+M4AAsWFNFU1s7I7LT+fvC9Rw0eZQlOSbefgjcgDNSrMC/ejhPgKsTFdRA2AiOMf4QVWc5EblURH4Udns/ESkDtorIIhHx9ZCxn/lpJ3FwkpvLH1rM1/YbR01jK7eeMYPLH1rMgjUp24bI+MNDwFeBU3CSmP8BTo64HA9MUtXfexVkNILW6M8YX4h2BOf/AX8Ku/0nYCNwDXAt8BvgvNiGNjSEanD8UmS8rKyWOefMpLmtg8cXlZGTkcacc2ayrKzWRnFM3KjqJ8AnIpIJnA8sVdUPPA5rQIIBoc1GcIzxXLQJzgTgIwARKcLpZ3OMqr4qIi0MsWXdsVRV30x+VhqZaf7oWxZaCh5KvD4sr+Xiw3ez5MYkhKo2i8hdOKM1ScmWiRvjD9FuftQMhIYYjgIagDfc29voWhxo+qGqvpnCPH9MT4UrzM2keHgWH5TXeh2KGXo+APb0OoiBSgsIHZbgGOO5aEdw3gW+79bdXAH8V1Xb3ft2w5muMgPgt20awk0fN5wPLcExifdD4D4R2YTzXtPmdUD9YVNUxvhDtCM4VwP74HyyGg/8NOy+rwNvxTiuIaOqvtk3PXAiTR83nE+rdlDX1Op1KGZoeQrYBXgaaBKRShGpCL94G17vAlZkbIwvRLtVwwpgsoiMArZF7A11DbA5HsENBVV1zRTu7o8C40jTSoYDsHzjdg7abZTH0Zgh5M8kcWPRNFsmbowv9KsPjqp+blfvZF3p4AfNbe1sb2rz9RQVOIXGluCYRFHVG72OYTCc3cRBVRHpdes9Y0wcRZ3giMhs4DScNulZkfer6lkxjGtI2BrqgePDImOwQmNjBiIYcJKa9g4lLWgJjjFeiSrBEZFLcZaCbwU+AVriGdRQEVqK7dcRHIBp44ZbgmMSTkQOBr6Ns5qquw9UByQ8qCh1Jjiqg2sVb4wZlGh//64B/gZ8L9lWNPjZzgTHnzU44ExTvbRyC3VNreRlpff9DcYMkogcBzwHvAwcBvwHyMbpv1UGvOZddH0LH8Exxngn2lVUo4GHLbmJrao6f23T0J3p44aj6hQaG5MgNwN/BL7i3v6Zqh6NM5rTCrzqUVxRCYolOMb4QbQJzn+AA+MZyFBUmSRTVID1wzGJNBXnPacDZzXVMABVXQfcSNc2Fb4TGsHp6PA4EGOGuGinqP4MzBWRdOBFoCbyBHcpuemHqvpmhmUEyc7wxzYN3SnKy2RsvhUam4RqAgKqqm6zv8ns7Jy+HWehg2+FEpw2y3CM8VS0Cc589+sNwPUR9wnOpyz//pX2qar6Ft+uoApnhcYmwZYCe+F8mHoZ+ImIlOMsbrgZp+GobwXCioyNMd6JNsE5Kq5RDFFVdc2+np4KmVEynJdXbaG+uY3cTFsXYuLuD8Cu7vXrgGeA593bZcCp0TyIiIwHHgDG4HwIm6uqfxSRkcCjwCRgLXCWqlbHKHbSbIrKGF+ItpOxr1ctJKuq+mZ2KxrmdRh96iw0Lq/lQGv4Z+JMVZ8Lu14uIrOA3XFWUq1S1WjbVLQBV6vq+yKSBywSkReBC4GXVfU3IvJj4MfAtbGKP1RkbFNUxngr2iJjAETkBBH5mYjMFZEJ7rHDRWSX+ISX2qrqk2MEJ1RobNNUxgvq+ERVl/UjuUFVN6nq++71OmAlMA44BbjfPe1+4GuxjDdgIzjG+EJUCY6IjBGRhThDxRfgNOAqdO/+FvCz+ISXulrbO6huaE2KBCdUaGwrqUyiiMh0EXlIRFaLyA736z9EZMYAH28SMBNYCIxR1U3uXZtxprC6+55LRKRUREorKyujfq40q8ExxheiHcH5PyAX2Nu9hPcffwk4JsZxpbxtO/y9TUMkKzQ2iSIiXwMW4SQkT+B8gHoC2B8ode/vz+PlAv8EfqCqXRo6uRsHd5uJqOpcVZ2tqrOLioqifr7OImMbwjHGU9FWjB4PXKCqq0UkcrVUGc6wr+mHyjqnB06Rj7sYh5s+zgqNTcL8Fngap/i3M/kQkZ8Aj7v3PxXNA7mtLf4J/ENVn3QPbxGRYlXdJCLFQEUsg9/Z6C+Wj2qM6a/+1OD01MW4EGiMQSxDytYd/u9iHG56ST6qsMI6Gpv4Gw/cHZ7cQOdoy13u/X0SZyvve4CVqnp72F3zcKbacb8+PeiIw9hWDcb4Q7QJzhvAFRGjN6Hf3ouAV2Ia1RBQVef/LsbhQoXGy8pqvA3EDAWlwD493DcNeD/KxzkU+CZwtIgscS8nAr8BjhORT4Bj3dsxYwmOMf4Q7VzDtcCbwIfAv3CSm4tFZB9gOnBQfMJLXZ0bbSZJDc7ovCzG5GdaobFJhKuAR9zppadwppBG4/S/+Q5wtojkhE5W1YbuHkRV36RrvWC4uNUNBt2PjVZkbIy3ou2D86Hbi+JGnB4S7cBpOF1Gv6Oqn8QrwFRVVd9MVnqAYT7epiHSdCs0Nonxrvv118Cvwo6HkpWFEef76pcoGHAyHBvBMcZbUVeLquoanOFeEwNV9S0U5mYi0tMHTP+ZNm44L6+qsEJjE28X0cPKpmRgu4kb4w/2V8ojydLkL9yMkuGdhcYH7DrS63BMilLV+7yOYTACoSkqS3CM8VSPCY6I9KtwWFWPHnw4Q0dlXTMlBTl9n+gj4R2NLcExpntpbobTYTU4xniqt1VUWyMuewJfBHKAevfrYcAeQFV8w0w9VfUtFOUlRw+cECs0NqZvoSLjNhvBMcZTPY7gqOqZoesi8m1gL+AQVV0fdnwC8CzwYjyDTDXtHcq2Hc2MGpZcU1RghcbG9CUgob2oLMExxkvR9sH5KXB9eHID4N6+EbguxnGltOqGFjoUCpOki3G4aeOGs6aynh3NPfV9NGZoS7NVVMb4QrQJzligp+GGDJweFSZKydYDJ9z0cW6h8SbraGxiT0QyReSnIrKv17EMVMCmqIzxhWgTnFeB34rI7PCDIvIFnH1hXovmQURkvIjMF5EVIrJcRK7s5hwRkT+5uwcvE5H9o4wxaVTVJdc2DeGmhwqNy2yaysSeqjbjjBiP8DiUAQt1MrYiY2O8FW2CcwmwDVgoIhvdlucbgXfc45dE+ThtwNWqOhWn+/H3RWRqxDkn4BQu7+E+7l+ifOyk0TmCk4QJzuj8LEbnZVodjomnhTg7hyelNNuqwRhfiLaTcRmwv7uPyxdwpqw2A++p6nPRPpmqbgI2udfrRGQlzk7kK8JOOwV4wN1Y7x0RGRHa+Tfa5/G7UIJTlIQJDlihsYm7HwEPiUgr8BywhYjGfz1tz+AHAWv0Z4wv9KvRn5vMRJ3Q9EZEJgEz+Xzb9XHAhrDbZe6xlElwKuubyQgGyM9Ozj6L08YN55WPKtjR3MYw62hsYi/0nvAn4I89nOOr7RnC2WabxviDJ3+dRCQX+CfwA1UdULWqiFyCOzU2YcKEGEYXf1V1LYzKzUiqbRrChRcaf2GSNfwzMZfcWzWEEhyrwTHGUwlPcNwdgv8J/ENVn+zmlHJgfNjtEvdYF6o6F5gLMHv27KR6J0nGbRrCTS/ZWWhsCY6JtWTfqsFGcIzxh2iLjGNCnCGLe4CVqnp7D6fNA853V1MdBNSmUv0NwNYdzUnZAydkjFtobB2NTTyJyFQR+aaIXCciY91ju4tIntex9cY22zTGHxI9gnMozo7kH4jIEvfYdcAEAFW9E6fG50RgNdAAfCvBMcZdVV0LU8bmex3GoFihsYkXdwr7XuAMoBXnfeq/OAsbfgWsB67xLMA+2DJxY/whqgRHRKar6geDfTJVfRPotfDEXT31/cE+l1+pqjOCk4RN/sJNGzec+R9V0NDSRk6GFRqbmLodOAQ4BngLaAq77zmc5Mb3CU5buyU4xngp2imqpSLynohcKiIj4hlQqqttbKW1XZO6BgecEZwOhRUbraOxibnTgGtVdT7QHnHfOmBi4kOKXsBGcIzxhWgTnKNxetXcAmwUkYdF5DhJ1mVAHtrZ5C95a3AgrNDYpqlM7GUDW3u4L4/PJz2+Yo3+jPGHqBIcVX1VVS/AafB3OU5fmueBdSLycxGZHMcYU0qlu01Dsjb5CxmTn0VRXqZt2WDi4T3g/B7uOwNYkMBY+i3U6M/2ojLGW/1aRaWqO1T1XlU9HNgLWItTJPyxiLwmIqfGIcaUkswbbUayQmMTJz8DThORl4Dv4PTEOVFEHgTOBG7wMri+dBYZW4JjjKf6vUxcRCaJyI04IzgH4xT9XYLTTv1REfl9TCNMMcm8D1WkaeOGs6aynoaWNq9DMSlEVd/AKTDOBObgLEy4CdgNOFZV3/MwvD51LhO3GhxjPBVVgiMiOSJyvojMx1m+fS5wFzBBVb+qqveo6lnAd4Fvxy/c5FdV30wwIIzITvc6lEGzQmMTL6r6lqp+EcjHafaZp6qHqupbHofWp0BAELEaHGO8Fu363i04ydCTOJ+gXu3hvPfouTjQ4PTAGTkso3OlRTKbPm5nofFs62hs4qMJpxdOo9eB9EdQxBIcYzwW7RTVj4BdVPWbvSQ3qOqHqrprTCJLUcm+TUO4MfmZTqGx1eGYGBORE0VkAU6CsxloEpEFIvIVj0OLSjAgNkVljMeiTXCKgGHd3SEixSJyfexCSm1OgpPcS8RDRITp44bblg0mpkTku8AzQD1wJU5h8ZXu7Xnu/b4WDAjt1ujPGE9Fm+DcgDMP3p1d8PmqBj+pqm9J+iXi4aaNG87qCis0NjF1HfBXVf2Sqt6pqk+6X7+EU/v3U4/j61NQbATHGK9Fm+AIzlLN7pQA1bEJJ7WpKpX1yb9NQ7hQofHKTVZobGJmFPCvHu77J+D7gq9gUGyZuDEe67HIWEQuAC5wbyrwFxGJ/CuWBUwHXohPeKmlrrmNlraOlJmigrBC47JaZk30/d8dkxzmA0cAL3Zz3xHA64kNp/+CItbozxiP9baKqoGdK6IEqAW2RZzTAvwHuCP2oaWeqrrU6YETMiY/k8LcTD4otxEcEzN/Au4WkVHAU0AFMBo4FTgB+I6ITA2drKorvAiyN4GA2F5UxnisxwRHVR8HHgcQkb8BP1fVTxMVWCqqqne2aUilBMcpNM7ng/Iar0MxqeN59+t33YvifMgK+a/7NTR1HkxcaNFJC9gycWO8FlUfHFX9VrwDGQq2plAX43DTxw3ntY8raWhpIycj2tZKxvToKK8DGKyATVEZ47neanBuAf6kqmXu9d6oql4b29BSz859qFKnBgeclVShQmOrwzGDpaqveR3DYAUDVmRsjNd6+7h9JvAPoMy93hsFLMHpQ2V9CyIwMie1EpzpJVZobEy4tIBgbXCM8VZvNTi7dnfdDFxVfTMjczJIC/Z7j1NfG5ufRWFuhhUaG+MKBIT2jg6vwzBmSEutv7Q+V1WXOts0hLOOxsZ0ZXtRGeO93mpwTuzPA6nqc4MPJ7VV1TenXP1NSKjQuLGlnewM3y1qMSahggGh3QZwjPFUbzU4z/L55Zk98eVSTb+pqm9h5oQRXocRF6FC4xWbtjNrYoHX4ZgkJiKjVHVr32f6V9CmqIzxXG8JjtXdxFhVfTOjhqXeFBXsLDT+sLzWEhwzWBtF5Gngb8Dzqpp0mULAioyN8VxvRcbrEhlIqmtoaaOhpT1lp6h2FhpbHY4ZtO8CF+KMIm8WkfuB+1T1Y0+j6oc0WyZujOd6LDIWkZzw631dEhNu8qqqS70uxuFEhGlWaGxiQFXvU9UjgT2Ae4BzgJUi8paIfFtEcj0NMArOXlRJN/BkTErpbRVVnYgc4F6vB+r6uJheVLpN/opSNMEBp9D44y11NLa0ex2KSQGq+qmqXq+qk4DjgHZgLs6ozn0isr+nAfYiEADLb4zxVm81OBcBa8Ku23jrIFSl6DYN4azQ2MSaOzp8Fs6U1WHAcpwNOL8EvCciP1bVWz0LsAdpgQCN7ZboG+Ol3mpw7g+7fl9CoklhqbpNQ7jp46zQ2MSGiBwOfAs4HefD1SPAj1T1XfeUn4nIj4AfA75LcAIB24vKGK/1q9GfiIwQkcNE5Ez364g4xZVyQjU4qbqKCqB4eBajhlmhsRkcEVkDzAd2B64AilX1u2HJTcjLQI+ZtIjcKyIVIvJh2LEbRaRcRJa4l371+4pWULAiY2M8FtXWzyKSBvwS+D4QXlDcICJ3AD9V1dY4xJcyquqbGZ6dTkZa6jaPtkJjEyNPAPeq6ke9naSqi+j9Q9p9wBzggYjjv1fV2wYVYR+CgYB1MjbGY9H+tb0duBL4FTAVKHS//hr4f8Dv4hJdCqmqb6YwN3WnpwDufG0NI4el80lFPU2tTv3BgjVV3Pnamj6+05idVPXavpKbKB/ndWBbDELqt2AAS3CM8VhUIzjAN4HrVPX2sGPbgF+KSBPwvzhDyaYHToKTutNTADNKhjPnldW0dygrNm2nqbWdyx9azJxzZnodmkkiInJ+L3d3ANuBpYPo1XW5+xylwNWqWt1DHJcAlwBMmDChX08QDAjtagmOMV6KNsHpwFm90J0PsRVWfdpa38KUXfK9DiOuDplcyK9Pm87/e3gxf3zpEz4or2XOOTM5ZHKh16GZ5HIfO99TwreKCT+mIvIscK6q1vfjsf8C/Nx9rJ/jjD5f1N2JqjoXZ1k6s2fP7td7XDAQsBocYzwW7RTVg8B3erjvYuDvsQkndVXWN6d0D5yQk2YUk50e4LWPKznvwAmW3JiB2B9YjbNCagrOlPgU4Cc4rSuOBS4Avgj8tj8PrKpbVLXd3f7hLuCAvr5nIIKCraIyxmO97SZ+WdjNtcAZIrIcmAdUAKOBU4A8IK4Fe8muqbWduqa2lK/BAXj70620tiuj8zL5+8L1HDR5lCU5pr9+B9yhqn8MO7YNuEVEWoAbVPUIERkDXI2z+CEqIlKsqpvcm6fijEDHXCAgVoNjjMd6m6Ka082xXXA+SUW6HfhjN8cNsHVHam/TELJgTRWXP7SYL+8zhhdXVnDPBbM7a3AsyTH9cDA9j8ysxFnsALAIGNXTg4jIw8CRQKGIlAE3AEeKyH44U1Rrcfa9irm0gNBhNTjGeKq3Rn+pu545warqUr+LMcCyMqfmZnNtE//+YDNj87OYc85MlpXVWoJj+qMMp3PxC93c9y33fnB64Gzt6UFU9RvdHL5nsMFFI2iN/ozxXLRFxmYQdnYxTu0E53tHTAZgxcbtztdN2zllv3GW3Jj++inwsIhMA54BKoEi4Ks47SnOds87DnjDkwj7EBDbTdwYr/UrwRGREmBPICvyPlV9Lorvvxc4CahQ1Wnd3H8k8DTwmXvoSVW9uT8x+tHOfahSvwYHYPfRuaQHhVWb6zjF62BM0lHVx0XkU5wi43OAscBm4D3gW26DP1T1sp4fxVtptkzcGM9F28k4D3gMZ4M72Ll0M/w3OBjFQ91H951Fw72hqidFE1eyqKofGjU4IRlpASYX5bJy03avQzFJRkQygTOAd1X1TK/jGahAQGhvtwTHGC9FW2fza2ACzrJMwVl9cCTOfPZnwEHRPIiXnUW9VFnXTG5mGlnp0eSAqWFqcb4lOKbfVLUZuBtnQUPSCoqN4BjjtWgTnBNx9qJa6N7eqKqvq+olOFNK/xPDmA4WkaUi8h8R2aenk0TkEhEpFZHSysrKGD597A2FbRoiTSnOZ8v2Zra5K8iM6YcPcKbCk1YwaMvEjfFatAnOGGCDqrYDO4CRYfc9x86pq8F6H5ioqvsC/wc81dOJqjpXVWer6uyioqIYPX1s3fnaGhasqeqyTcNQ2ZtpSrHTtdlGccwA/BD4kYic5G70m3SCYgmOMV6LNsHZgNNNFOATnELhkAOBplgEo6rbQ23X3aLldBFJ2iU4M0qGc/lDi1m/rYHC3MzOPjEzSoZ7HVrc7V2cB1iCYwbkKZwpqqeBJhGpFJGK8Iu34fXN9qIyxnvRfjp6Eac9+r+A3wP3i8gsoBk4nBjtJi4iY4EtqqoicgBOAtZjnwu/O2RyIXPOmcm5dy1k1LCGIdX0rjA3k6K8TFZuqvM6FJN8/kyS728XDAiq0NGhBALS9zcYY2Iu2gTnWiAHQFUfFJF6nJUO2cDlwF+jeZAeOoumu497p/uYl4pIG9AInK2a3B+DZk0sQIEPyrdzxdG7D4nkJmSKFRqbAVDVG72OYbCC4iQ17aoEsATHGC9EleCoagPQEHb7XzijOf3SQ2fR8Pvn0P0WEUnr38ucbW+OnTJmyO3NNKU4j7+t2UprewfpQWuMbfpHRAqAacB44D+qWi0iWUCLu1mmb4VGbdo7lCG0eNIYX+lvo7+9gC8AxcBGoFRVP4pHYKlgwZoqbpi3HIDvfHFXOnTSkJqmmlqcT0t7B2sq69l7bL7X4Zgk4RYW/wpnE81snOmqLwDVwD+BUpzRX99KcxMc24/KGO9E9bFaRPJF5FFgOU6Tvp8BDwLLReQxEbG/Xt1YVlbL2V8YD0BJQXZnTc6yslqPI0sMW0llBuiXwMU409+7QZc5nqdxtmzwtaCb4Nh+VMZ4J9p5gztwloKfDwxT1XxgGHABzn4wd8QnvOT2vSMmk5kWJBgQxuY7u1scMrmwc8+mVLdr4TAyggErNDb9dT7wY1X9G84KznBrcJIeXwu4NTi2H5Ux3ol2iuoU4Ieq+lDogKo2Av8QkRzg9ngElwrKqhsYm59F2hCsQUkPBthjjG3ZYPptBE4i050MotsWxlNpwZ01OMYYb0T7V7ce2NTDfRtxmv+ZbpTXNFJSkO11GJ5xVlLZCI7plw+hx31aT8BpCOproREcS3CM8U60Cc6fgWtEpMtfanf05hpsiqpHZdWNjBviCU5VfTOVdc1eh2KSxy9w2kXcjdN/S4H9ROTnwHdxCpB9LVSDY83+jPFOj1NUInJLxKE9gA0i8iJQAYzGqb9pxFnVYCK0tHWweXsTJQU5XofimSlhHY2L8vy5pYbxF1V9WkTOAW4BLnIP3w2UA99U1ec9Cy5KwYCN4Bjjtd5qcM6MuN3qXsJ3Dg/NPZxObDfcTAmba5tQZUhPUU0NW0l1+J6W4JjoqOpjwGMisifONjHbgI+SpfFn0KaojPFcjwmOqu6ayEBSUVm10xuxZMTQTXBG5GRQPDzLCo3NgKjqx8DHXsfRXzaCY4z3knKn3mRRVt0IMKSnqAD2HpvHqs1WaGyiJyK74GzqWwJkRdytqnpt4qOKXtAa/RnjuagTHBHZDWca6jBgJM6Q8RvAbar6aXzCS25lNY0EBMYOj3x/HlqmFOfzxidVNLe1k5nm+xW+xmMicirwMM5y8AqgJeIUxdkfz7es0Z8x3osqwXF3Dp8PNAHPAluAMTi1N+eKyFGq6vulm4lWVt3AmPwsMtKGXg+ccFOK82nrUFZX1LPPLsO9Dsf436+AF4ALVXWb18EMhC0TN8Z70Y7g3AYsBk5wN94EOpeJP+fef3Tsw0tu5dVDuwdOyM4tG+oswTHRGA/8v2RNbiBsLypfbwlqTGqLdmjhAOCW8OQGOncZvw04MNaBpYKy6sYhX38DzpYNWekBKzQ20VoA7OV1EIOxc4rKMhxjvBLtCE4jMKqH+0biTF2ZMG3tTg+ccUN4BVVIMCDsNSbPEhwTratwtoGpB14EaiJPiPyw5TcBKzI2xnPRJjj/Bn4jIp+q6puhgyJyGPBr4Jl4BJfMNm9vor1DbYrKtffYfF5YsRlVRUT6/gYzlC1zv/4Np6C4O76uVk/rXCbucSDGDGHRJjhXAU8Dr4lIBTs7GY8G3gaujk94ycuWiHc1pTiPR0s3UFHXzJj8ob2qzPTpInpObJJCqMjYpqiM8U5UCY6qbgUOE5HjgS8AxTibby5U1RfiGF/S2png2AgO7Cw0XrFpuyU4pleqep/XMQxW0IqMjfFcnwmOiGThDBlfoar/Bf4b96hSQLmb4BSPsD/mAHuHbdlw1F6jPY7GJAMRmQrMwllVda+qbhaR3YEtqurrzpG22aYx3uszwVHVJhEZAdhnkX5weuBkWmM71/DsdMaNyGblJl//XTI+ICK5wL3AGTj736XhfLDajNMjZz1wjWcBRmHnVg32tmmMV6JdJv4P4FvxDCTV2BLxz5tSnG8rqUw0bgcOAY4B8oDwqvTngOO9CKo/dm626XEgxgxh0RYZrwfOEpH3gP/gdDIOH3tVVf1LrINLZuU1jew3foTXYfjKlOI8Xlm1habWdrLSbWTL9Og04EpVnS8ikf9R1gETPYipX2yzTWO8F22C8zv3azHOnHgkBSzBcbV3KBtrGjlpRrHXofjKlOJ8OhQ+2VLP9BLraGx6lA1s7eG+PKA9gbEMiCU4xngvqikqVQ30cbGP42G2bG+irUNtiirClLBCY2N68R5wfg/3nYHT6djXgu47qxUZG+OdqHcTN9Err3FWUI2zJeJdTByZQ05GkBWW4Jje/Qx4UUReAh7HGSE+UUR+iJPgHO5lcNEIBpwMp8NGcIzxTNQJjohkABfi7EvV2QcHuF9VW+ISXZIqq3a6yFsPnK4CAWGvsbZlg+mdqr4hIscAvwHm4BQZ3wS8Axyrqu95GV80gp2N/izBMcYrUU1RicgU4BPgz8A0nDnwae7t1W6/CuMq2+aO4Ng+VJ8TWkmlNnRveqGqb6nqF4F8oATIU9VDVfUtj0OLijuAYyM4xngo2mXic4FaYLKqHqSqJ6vqQcDuOBvh3Rmn+JJSeU0jhbmZtlKoG1PG5rG9qY1NtbY/q+mbqjaq6ka/b64ZKc3NcKwGxxjvRJvgzAauV9X14Qfd2zfgbN9gXE4PHBu96Y4VGpuhIDSCY1NUxngn2gRnLdDTngNZOH1yjKususESnB7sbQmOGQJCNTg2RWWMd6JNcH4M/EJEDgw/KCIHAT8Hro11YMmqo0PZWNNkK6h6kJuZxoSRObZlg0lpnVNUluAY45loV1H9L06x3wIRqQAqgNHuZStwnYhcFzpZVQ+IdaDJorK+mZb2DuuB04spxbaSyqS20BSVJTjGeCfaBOdD92L6YEvE+zalOJ8XVmyhoaWNnAxrxWRSj+0mboz3ovrroqq20WaUyqqdJeIltkS8R3uPzUcVPt5Sb/t1GQDcfe6izgb8PkpsWzUY4z37+BxjoQTHanB6NjWs0NgSHONaTj8SnGiJyL3ASUCFqk5zj40EHgUm4SygOEtVq2P5vDt3E7cExxivWIITY2XVjYwalmFTL70oKcgmNzPN6nBMJ1W9ME4PfR9ON+QHwo79GHhZVX8jIj92b8d0oYSN4BjjvWhXUZkoldc02uhNHwIBYW/bssEkgKq+DmyLOHwKcL97/X7ga7F+XhEhINBhNTjGeCahwwzdDRdH3C/AH4ETgQbgQlV9P5ExDlZZdQN7j83zOgzfm1Kcz1OLy1FVxB3ONyZERCYB5wF70k0PLlU9axAPP0ZVN7nXNwNjeojhEuASgAkTJvT7SYIBsUZ/xngo0SM49wHH93L/CcAe7uUS4C8JiClmVJXy6kZbIh6FKcX51DW3ddYsGRMiIrNwanLOdS974HRTPwM4CCiM1XOpsylat1mIqs5V1dmqOruoqKjfjx0QsUZ/xngooQlOD8PF4U4BHlDHO8AIESlOTHSDV1XfQnNbh22yGYUpxc4ol01TmW7cCjyOs6GvAN9W1d2Aw3CSkVsG+fhbQu8r7teKQT5et9ICYjU4xnjIbzU444ANYbfL3GOfIyKXiEipiJRWVlYmJLi+WA+c6O01Ng8RrKOx6c5+wMNAh3s7C0BVFwA3Ab8Z5OPPAy5wr18APD3Ix+tWwKaojPGU3xKcqA12+DgeOnvg2BRVn3Iy0pg0apiN4JjuKNDiTh9VABPD7tuAM2UVFRF5GHgb2EtEykTk2zgJ0nEi8glwLINPmLoVDIgVGRvjIb+tZS4HxofdLnGPJYXyGuuB0x9TivNYvtESHPM5K4DJwHyc5OSHIlIKtAA/AtZE+0Cq+o0e7jpmsEH2xaaojPGW30Zw5gHni+MgoDZstYPvlVU3MCInndxMv+WN/jRlbD7rtjZQ39zmdSjGX+YCY93r1wHFwCrgU+BA4BqP4uqXgFiCY4yXEr1M/GHgSKBQRMqAG4B0AFW9E3gOZ4n4apxl4km1RURZdaPV3/TDFLej8UebtzNr4kiPozF+oaoPhl1fKSJTgIOBbOAdVY1LUXCsBW0ExxhPJTTB6WW4OHS/At9PUDgxV17dyG5Fw7wOI2lM2SW0ZUOdJTimR6paD7zodRz9FQyIbbZpjIdsLiVGVJWy6kYO39MfBc/JYJfhWeRn2ZYNBkTkROBNVd3uXu+Vqj6XgLAGxUZwjPGWJTgxsm1HC42t7TZF1Q8iwt7F+ZbgGIBncZr4veteV5weON1RIJiguAYsaDU4xnjKEpwY6VxBZU3++mVqcT6PlW6go0MJBGzLhiFsV2BT2PWkZ8vEjfGWJTgxYj1wBmZKcR4NLe2s39bApEKrXxqqVHVd+E1gk6q2Rp4nImnALgkLbBCCAaGt3RIcY7zit2XiSSvUxdh64PRPaCWVTVOZMJ8BM3u4b1/3ft8LiI3gGOMlS3BipLy6kbysNIZnp3sdSlLZc0weAYGVm23LBtOpt7nKLKA5UYEMRlrQanCM8ZJNUcVIme0iPiBZ6UF2LbQtG4Y6EZmBswdVyIkisnfEaVnAWcDHiYprMAJie1EZ4yVLcGKkrLqRCaMswRmIKcX5LNlQ43UYxlun4jT+BKcG5/oezvsM+G5CIhokKzI2xls2RRUDqkp5TaOtoBqgKcX5lFU3sr3pczWlZuj4FZAH5ONMUR3t3g6/ZKrqZFV9ybMo+8H64BjjLRvBiYHaxlbqm9usB84ATXULjVdtquOAXa2j8VDkrpgKZbgp8cErKEJbR4fXYRgzZFmCEwO2RHxwwldSWYJjQkRkT6AEp/ami2TpZNzcZiM4xnjFEpwY2Jng2AjOQIzJz6QgJ51Vm63Q2ICITAUeAfah+xVVydHJOCBYGxxjvGMJTgyEeuBYgjMwIsLeY/NZscmWihsA/gpkAqcBK4AWb8MZGKcGx6aojPGKJTgxUFbdSG6m9cAZjCnF+Tz07jraO5Sgbdkw1M0EzlbVZ70OZDACIrRbfmOMZ1KimM9roRVUIvaHeaCmFOfR1NrB2q07vA7FeG8N3dTdJJu0gNBhq6iM8YwlODHgNPmz6anBsC0bTJirgetEZDevAxmMYMBWURnjJZuiioGy6gYOmFTgdRhJbY8xuaQFhJWbtnPSjKTYS9HEz6+BccAqEVkL1ESeoKoHJDimfgsEBBvAMcY7luAMUm1jK3VNbbbJ5iBlpgWZXJTLKis0NvChe0lqadbozxhPWYIzSOXWAydmphTn8e5n27wOw3hMVb/ldQyx4BQZW4JjjFesBmeQbIl4bNz52hqyM4JsrG2ipsFZFbxgTRV3vrbG48iMV8QxXkQOEZFhXsfTX8EAluAY4yFLcAapvMYZwbF9qAZnRslwnl22CYCVm+pYsKaKyx9azIyS4R5HZrwgIpcB5cA64A1gL/f4kyLyAw9Di1owEKDdNts0xjOW4AxSWXUj2elBRg7L8DqUpHbI5EJ+e/p0AO54dTWXP7SYOefM5JDJhR5HZhJNRP4HuB24C2fTzfD+C68CX/cgrH6zERxjvGUJziCVVTdQUmA9cGLhxOm7MCwjyBufVHHuARMsuRm6vg9cr6o34IzehPsI2DPxIfVf0GpwjPGUJTiDVF7TaCuoYmTBmqrO6/e89VmX22ZIGQss6uG+DpKkCWAwELBGf8Z4yBKcQbImf7ERqrm56/zZHLP3aFrbOrj07+9bkjM0rQaO6OG+w3H2p/K9YADaLMExxjOW4AxCfXMbNQ2ttkQ8BpaV1To1N7sX8uvTppOTmUZRXgZLNtR4HZpJvD8APxaR/wX2cI+NFpFvA1cBv/cqsP4IBMSKjI3xkPXBGYRQDxxbQTV43zticuf10flZ3HzKPlz5yBICVts05Kjq3SJSAFwP3OQefg5oAG5U1Yc8C64fbC8qY7xlCc4gWA+c+Dl53134zwebuf2Fjzl679HsOSbP65BMAqnqrSJyJ3AwUAhsA95W1VpvI4teUMSmqIzxkE1RDUKZdTGOGxHhF6dOIzcrjasfW0pru21aOFSIyPkiMkpV61T1BVV9SFX/q6q1IjJSRM73OsZoBALO6KON4hjjDUtwBqG8ppHMtACFudYDJx4KczP55dem8UF5LX951ToaDyF/Ayb3cN+u7v2+l+YmOFaHY4w3LMEZhLLqBsZZD5y4OmF6MSfvuwt/evkTlm9MmtkJMzi9/UKNArYnKpDBCI3gWC8cY7xhNTiD4CwRt+mpeLv5lH14+9OtXP3YUuZdfhgZaZaXpxoROQU4JezQz0SkMuK0LOCLwHsJC2wQgmIJjjFesr8Ug1Be3WgrqBJgRE4GvzltOqs21/Gnlz/xOhwTH6OB6e4FnCmq6RGXicALwHe9CLC/gjZFZYynbARngBpa2ti6o8VWUCXIMVPGcMasEv7y2hqOmzqGfceP8DokE0OqehfO3lOIyHzgMlVd6W1Ug9OZ4LRbgmOMF2wEZ4DKO1dQWYKTKNd/dSqj8zK5+vGlNLW2ex2OiRNVPSo8uRGRdC/jGSgbwTHGW5bgDFBZjSU4iZaflc5vT5/B6op6bn/xY6/DMXEkIoeIyH9EpA5oEpE6EXlORA72OrZoBW2ZuDGesgRngKwHjjcO37OIcw6cwF1vfErp2m1eh2PiQESOA14FSoBbgcvcryXAqyJyrHfRRS9UZGzN/ozxRsITHBE5XkQ+EpHVIvLjbu6/UEQqRWSJe/lOomOMRll1AxnBAEW5mV6HMuRcd+IUxo3I5prHl9LQ0uZ1OCb2fgnMA2ao6s2q+lf36wzgWeBX3oYXHVsmboy3EprgiEgQ+DNwAjAV+IaITO3m1EdVdT/3cnciY4xWeXUju4zI6nwTM4mTm5nGrWfsy9qtDdzy34+8DsfE3nTgLtVui1fmsnOlla+FGv11WA2OMZ5I9AjOAcBqVf1UVVuAR+ja+yJpWA8cbx08eRQXHjKJ+xasZcGaKq/DMbFVQ8+djCe79/teqAbHpqiM8UaiE5xxwIaw22XusUini8gyEXlCRMZ390AicomIlIpIaWVlZD+w+HMSHCsw9tK1x+/NpFE5/OiJZdQ321RVCnkc+LWInCciWQAikiUi5+FMTz3maXRRCogVGRvjJT8WGT8DTHLn218E7u/uJFWdq6qzVXV2UVFRQgNsam2nqr7Zmvx5LDsjyG1n7kt5TSO/ei6pW6aYrq7FqbW5H9ghIrXADvf2s+79gyIia0XkA7fOr3Swj9cd24vKGG8lOsEpB8JHZErcY51UdauqNrs37wZmJSi2qJWHloiPtATHa7MnjeQLk0by0ML1vP7xzpG8BWuquPM126AzGalqo6qeC+wDfAun6PhCYB9VPU9Vm2L0VEe5dX6zY/R4XYTq89qs0Z8xnkh0gvMesIeI7CoiGcDZOKslOolIcdjNkwHffTS3JeL+ctmRkwkI/ODRxdQ2trJgTRWXP7SYGSXDvQ7NDIKqrlLVB1T1FlV9UFVXeR1Tf4SWiVuRsTHeSOhWDaraJiKXA88DQeBeVV0uIjcDpao6D7hCRE4G2oBtOJ/cfCXUxdimqPzhyL1Gc9PJ+/Czp5dz8pw32d7Yyp/P3Z9DJhd6HZoZIPcD0IU4CxOKgU3AQuB+d4HCYCnwgogo8FdVndtNDJcAlwBMmDCh308QDNoycWO8lPC9qFT1OeC5iGPXh13/CfCTRMfVH2XVDaQFhDH5WV6HYlzfPHgSr6yqZP5HFRQPz2KfYhu9SVYiMgX4L7ALsAioAKYB5+PsMn68qq4Y5NMcpqrlIjIaeFFEVqnq6+EnuEnPXIDZs2f3O0ux3cSN8ZYfi4x9r6y6kV1GZHcuAzXeW7CmiqVlNZwwbSybaps46f/eYHNtrEo1TILNBWqByap6kKqerKoHAbvjLBG/c7BPoKrl7tcK4F84I0UxFbRGf8Z4yhKcASivabTpKR8J1dzMOWcmfzlvFj89cW82VDdy0p/eYE1lvdfhmf6bDVyvquvDD7q3bwC+MJgHF5FhIpIXug58CfhwMI/ZHdts0xhvWYIzAGXVDdYDx0eWldUy55yZnTU3Fx8+mV+dOo2G1nbOvPNtlm6o8TZA019rgZ7mf7OA9T3cF60xwJsishR4F/i3qv53kI/5OTaCY4y3LMHpp+a2dirqmm0FlY9874jJnysoPufAiTx3xRcZlhnkG3e9wxufJL4ZpBmwHwO/EJEDww+KyEHAzxlkHxy3k/q+7mUfVf3lYB6vJwGrwTHGU5bg9NOmmiZUYZyN4PjepMJh/PN7hzBx1DAuuu895i3d6HVIJjr/C+QDC0Rkk4gsFZFNwFvu8etE5N3QxdNIe2F7URnjrYSvokp2O3vgWIKTDEbnZ/Hodw/iO/eXcuUji9lW38yFh+7qdVimdx8Sh5qYRAtaoz9jPGUJTj+V1zQAluAkk/ysdB646ACueHgxNz6zgq07WrjquD0RsVVwfqSq3/I6hlgIWKM/YzxlU1T9VFbdSDAgjLUeOEklKz3IHefuz9lfGM//vbKa6/71odVGmLhK62z053EgxgxRNoLTT2XVjYzNzyItaLlhskkLBvj1adMpzM1kzvzVbNvRzB/PnklWetDr0EwKCo3gtHVYhmOMF+yvdD+VVzfa9FQSExGu+fJe3PjVqTy/fAvH3f4aL63c0uUc26jTxELQioyN8ZQlOP1UVt1gK6hSwIWH7sofz96PjbWNfPfBRTz3gbPCyjbqNLGSFrApKmO8ZFNU/dDa3sHm7U3WAydFnLLfOApyMrj4gVIuf2gx5x64jX9/sKlL00BjBirQmeBYhmOMF2wEpx821zbRobaCKpUcvmcRj373YDLSAjz4zjr2Gz+Cg3cb5XVYJgXs3GzT40CMGaIswemHDdXuEnHbhyqlNLS0kZUWpKQgm1dWVfC1O95iy3bbqNMMju1FZYy3LMHph51N/myKKlWEam7uOG9/Xv+fo7jwkEks3VDL0be9ap2PzaB0Jjg2hGOMJyzB6Yfy6kYCAmOHWw+cVBG+UWcgINx48j787swZDM9J54qHF/P9h95n244Wr8M0SahzisoGcIzxhBUZ90NZdSNj8rPISLO8MFV874jJnzt2+qzxnLLfOP76+qf84aWPWfjpNn57+nSOmTLGgwhNsgq6jf46rKGkMZ6wv9T9UFbdYAXGQ0RaMMD3j9qdeZcfRmFuBt++v5QfPbGUuqZWr0MzSSLY2ejPEhxjvGAJTj+U1zRa/c0QM6U4n3mXH8b3j5rME4vKOP4Pb7BgTZXXYZkkEHDfXa3RnzHesAQnSm3tHWyqbWKcraAacjLSAvzPl/fmiUsPISMtwDl3LeSmZ5bT1NrudWjGx9LcDCe059n2plYufqCUA375Ehc/UMp2Gw00Jq6sBidKm7c30d6hNkU1hO0/oYDnrvgiv/3vKv721lqeXlLOD47dk/MPntR5zoI1VSwrq+22tscMLfXNTgJz52tr+KC8lqbWdt5es5W2DuXVjyq48pHFpAUCLN1Qw77jR/C7s/YlPyvd46iNSR02ghOlclsiboDsjCA3nrwP//jOgQhw/dPL+cEji2lp67BtHkwX1zy+DICGlnZeWrGFNz6p6qzHaW1XXl1Vycsrt1BR18xrH1VwzeNLvQzXmJRjIzhRCvXAsX2oDMChuxcy/3+O4vKH3uepJRt5e81Wmts6uOO8/W2bBwPA0g01ndcVCIizw3hbhxIQUIVQ/XFLu7JkfU13D2OMGSAbwYlSKMHZZYT1wDGO/Kx0HrjoQE6cPpYtdc00t7V3rpwxZt/xI8hwl4pnBIUj9iri6CmjGZ2XybFTx3DkXkWkB3f+f9ljdK5XoRqTkmwEJ0rlNQ2Myc8kMy3odSjGRxasqeKdT7dx3oETeOjd9Xzjrne44av7cP7BExFLdoa03521L9c8vpQl62vYb8IIbjuza43N9qZWrnl8KYvWVlPX1MZnW3dQWddMUV5m1M+xvamVqx9banU8xnTDEpw+3PnaGmaUDKesurFzBZUVkhrYuc1DqBPykXuP5tK/L+KGectZWlbDr06dTla6JcRDVX5WOnO/OTuq+5eV1XDWnW9zzO9eJTMtyH4Tuk9WwhOaGSXDqW1o5f0NNbR3aGcdT/hzWgJkhjKbourDjJLhXP7QYlZX1FNSkGOFpKZT+DYPAMdOGcN9Fx7AYbsX8uT75Zxx5wLKaxo9jtIkgxklI9irOI/tTW1U1jfzag9Fx5f9/X1ecQuTX1pZwXvrqjuXobe0K4vXVXc5/+rHljJ/VQUVdT0/pjGpykZw+nDI5EKuOm4P/vep5Wysaezyid0Mbd2N4B26RyGH7lHISyu28MNHl/DV/3vT/r+YqGyq2bmDvbPKqoJH3l3P88s38/76aoKBwOf2RUsLCCLO+QB1Ta2cdefbrN26g4mjcliyoabLyq23VlfR3Nbe61S7jfqYVCGaAl02Z8+eraWlpXF57MdLN/DTpz4kMxigrrmNK47enau+tFdcnsukljWV9VzyQClrtzZw3YlTuOjQSUOyLkdEFqlqz3M1SSCe7zEhFz9QymsfVdDS7qyyykgL0NTadSfynIwgre0dtLYrGUHh0D0KSQ8GWLK+homjclhWVktz287vEZzVW+0K4q7cGl+QzajcDDbWNHWbwFz8QCnzV1XQ1uE8x5F7j+51qs0YL/X2/mIjOD1obe/gF8+u4P6317HPLvmU1zRyxdG78/eF6zlo8ij7RG76NLkol6e+fyhXP7aUnz+7gg/Kavj1aTPIzrC6HPN5kUXJt5w+gyNve5Wahp0dj3Myguw/sbDHwuXZv3iR5vqdozwjh6Uza9LIzvNPmzmOqx5bygZ3VegrKyv4/j8WkZWextINNRQPz2L5xu2doz6pvnx9e1MrVz68mA/La9lvQoGNVqUYS3C6UVnXzPf/8T7vrt3GV6aPZcGardxxrtPf5KDJo2yaykQtLyudO8+bxR2vruZ3L37Mx1vq+es3ZzF+pDWMNF11V5T8hUkjO0d1MoLC/hMLeh1NmTmhoMv5syaN/Nz5P3t6OQ0tzjYj7aq88cnWztGdirpmAgLBgHTW9gzP7vsPfl/TWpH333jyVG6ct6LXabBETJVd8fBiXv2oEoCXVmzhtD+/xZj8LD6pqI/Zc25vauWqR5ewZEMNMy2JSiibooqweH01l/79fWoaW/jt6TPYVNvEjJLhXZIZW0VlBmL+qgqueGQxwYAw5xv7c9geQyNBtimqgQstJe9pxGYg54dPhaUFhY4OJXzD88LcDPafWMDi9TUEBTZvb2biqBwaW9p7/KP/7fve47WPK3uc1rr4/lLmf+RMewVFyMkI0tDaTnsv02Dh3zPQqbLekqSODmXvn/2Hlvbu/wamBYSj9h7NXef3/pzbm1r5wSOL+bB8e7evz9f+/CZLNtQCTuJ4zN6jmdvHYw7m50qkUPK2rKw2qsQ2Holsb+8vluCEefS99fzsqeWMGZ7JX8+bzdRd8mMQnTE7ra3awSUPlvLJlnq+ceB4fvm16Z11OamaOFuC4y+RSVBrWwdvra7qHPUJTyTaO5Qjb53fOaUVFGHvsXnkZAZZXVFPUV4mRbmZvLVma5fnyMsMcsCuo1iyoYbh2ems3bqjSxIVqTA3g9L/Pa7zdnlNI0feOr+zeBpgdF4m7/702H79rL0lSX95dQ2//e8qgiK0q3Y2ZQxPeAIC00uGs6mHeqWGlja+8qc3+KyqofP8/caPYNSwTBatryYgUFXftTA8Iy3A2z8+mlG5Pfc76uuP/sUPlPLqRxWdtVjxqJOKJvG44N53ee1jZwQsPegkhF0S2yjivPj+Ul75qKLXZLc3VoPTh5a2Dm56Zjn/WLieL+5RyJ/OnknBsAyvwzIpaFLhMP512aFcdN+7PLRwA59W7uCeC77A0rKazqlPY+Ipciqsu1GfkGBAuhQtt6uyfNP2ztvVDa2UVTcyfmQ2m2qaOmt36prbeXlVBQBbd7SQHhQCCm0dSnpQyMtMZ3tzK21uMlHb0MIpc95kU20TRXmZrKvaQYe7pUWHWyC97/gRXX6Ovv4AN7W289onld3WEy38dCu3vfARX9pnDCgs2fD5ZC8YEFSVpe7oyysrt/Cd+0uZ+81Z/M8Ty3j3s200t7bTFPb6dCi8H1GzlJsZpKWtw3lMEVraOjj416+QnRHkgF1Hdps4XPXoEuZ/VEm7uzFreH+jsuoG3viksjP5a2lX3vy4kvPvWciqzXUxG9EJtRho6yYGcEbA3l5T1Xm7tV1577NtXe5/59OtXeLsrp7rrTVVXVodxLLma8gnOBXbm7j0H++zaF013ztiMv/z5b0IBobeSheTOMMy03jkkoP5yZMf8Mh7G5h+4/MERDh+2ljqm9qoqm+msJdPd8bEUl8NCfcdP6JzSiu0tUT4yEpuZhr/vuKLnUnSvuNH8O5n26ht3FkcPTw7nf0nFnQmUTd8dSo3PbOCJetrmFw0jA/Ka1la5iQSFXXNjMhJ55nvHcIfXv6Yt1ZXsaO5nYmjuu4D+L0HF/H2mq0ozvRv+B/g2sZWLn6glJa2js4kCSAzGOC8uxeyYE0VWelBbjp5H4qH73zcyGTv/XXVnSMw7QrvfraN2b94ifYOJfQKjMhOp6GlrfP1UTeRC8lKD3LI7jsLw6t3tPDe2mpaGjt4aeUWvv+P93nw2wd2nv9heS2vuslN6LV+ZeUWvvz711m3bQfN7sq6zuQPaGjt4PVPnGRj/qoKrnls6aCmwXY0t/FmxOaw70f0WLrnzc86k7Z2dyaouqGVL//+dSrqmmjvUOqa2hDofK0mjupae/j3d9bR0NLeeU4wIOw3YcSA4440pKeoFq3bxqV/f5+6pjZuPXMGJ83YJQ7RGdOzax5fwhOLyhmbn8nWHS2dfzh2LRzG/hMKmDWxgNmTCti9KJeAm3iHumsnS12YTVElt/5MaYWE1/lEM+3whV++RGVdc+ft8OkoVe38MDB1l3wqtjeRk5HG+m0NXR4jPSAcsOtIPtpSR1NrO02t7fzy1Bm8vGoLi9dVk5Ue7JxqA0gLCkf3EVf4z5EeFKaMzefjirouy/dDNUvRvj4H/PIlKsJ+VoDvfHESqyt28O5n22hoae+SKAXESWJCOWVA4KBdR5Kbnd75nO99to3qsNV2aQHh5lP24ZVVFZ31MV0Ku0tG8L8nTeHmZ1awrLyW/cLuf++zbTS1tdPU2jU5zM0M8saPjqZgWAZLN9Rw+l8W8MVQm4INNUwbl8+ystouU3JTxuZRMjKHxeuraWhpJyjwz8sOZc8xebz2cSUX3fceh+w2ksz0IK9+VIkAL151BJMKh/X4bxLJanC68Y+F67hx3nJ2GZHNX785i73HWr2NSaxQV+zzDpzA3xeu5/az9iU3M41F66opXVfN++uq2eo2dsvLSmP/CQXMnlhAVnqAO15dw5/dlX2RW0b4jSU4qSWaQub+Fkf3lRC1tndw4C9fZlvDzj+eWWkB2lVpjRhFCJk1sYB/XnpIl2P73fxCl2X3fdX1dPdzXP3Y0l5j7etnj0yaCnMz2VS7s8ljQODQyaPIzkzrdiSpu7jDHzMoQlZ6gB3uSrmQ8JGUaOw/YQSFeZksWV/D+JE5LCurYXJRLoV5mSxYXUVGMMCLVx3O+JE7k5HI5C08zg3bGvjaHW9R39RGVnqQ+uY2Jhfl8uRlh5CbmcYnW+o48U9vcOL0Yv54dvRT9b6qwRGR44E/AkHgblX9TcT9mcADwCxgK/B1VV07mOcM/8Tb3NbOjfOW8/C7G9jNrYcYnmNL9kxiRSYl4e0HvnvEZL6L88l13daGLgnP7S99jKrzJvjNe95lj9G5rN26g/MOmkhLWwcfb6mjeHgWeT38MRnM6E+yjRyZ+OhrSivac8J1tzFpuPRggEDExkK5WWldRk5K11Z36fS8IWKEBz6/7L6v6ZDufo6+Yu3rZ+/u+4+4ZX7nCEyHwkdb6ntMYLqLu7seSkfc+mqXacLI5KavhKesupEnLzu08/b8jyq46G/vsWpzHeCMLv383yu7/Kzh05mRcY4fmcMeo3N559NtnXVdY/Mzyc100pA9xuRx2ZG788eXP+FrM8dx1F6je4kuOglNcEQkCPwZOA4oA94TkXmquiLstG8D1aq6u4icDfwW+Ppgnje0n9TPT9mHu974jCUbashKD3DzKdMsuTGeiNzH6pDJhcw5ZybLymo7j4kIkwqHMalwGKfPKgGc2oIlG2pYtHYbTy4u73yzufuNz7j7jc86Hz8vM43iEVkUD89mlxFZ7DI8m+IR2QQFLv37+/zx7P04cq/RXRKtvoR+j0Jx9+d7jelNNAlRZI+fyJ5AfSUB0HdyEqtY+/v9s/tIvAaSVB2wa9fHHJ6TQW1DS9S3I2M4aq/R5GendyZNbR2fLwjuK85PK3d0uR16/wq57KjJzFu6kUsfXEReVtqgmy8mdIpKRA4GblTVL7u3fwKgqr8OO+d595y3RSQN2AwUaS+BRjN8PG/pRq58eDFpASEzPcjc82f5cjjfmGiET289+M46bjplGrsMz2JjbRObahrZVNvExppGNtY2sqmmqXOqK1xofj0tIJ31PX3p6FDaOpTxI7PZ0dwe1bSYTVGZWOhr6qe/02J+Eo/YIx8zvLA7mts99VAazPL0aGqzzrxzAe+tdQqao3kO39TgiMgZwPGq+h339jeBA1X18rBzPnTPKXNvr3HPqYp4rEuASwAmTJgwa926db0+t6py2h0LWLyhxvaTMkktcnormhqcptZ2Ntc2dSY8j5du4J3PtjFrQgFf2HVkv57/vc+2sWh9ddS/R5bgGJMaBpuIRfP9vdXxdMdXNTixoqpzgbngvPn0df7bn25l3bYG20/KJL1oprciZaUHO6e7Fqyp4uOK+s7fhau/vGfUvwsL1lTxWOkG+z0yZgiKx/RcpN7qePor0QlOOTA+7HaJe6y7c8rcKarhOMXGA9ZbQae9OZtk011B7yGTC6P6vzyY3wX7PTLGxFss6qRCEp3gvAfsISK74iQyZwPnRJwzD7gAeBs4A3ilt/qbaAzkE68xqWgwvwv2e2SMibfBjhKFS3gfHBE5EfgDzjLxe1X1lyJyM1CqqvNEJAt4EJgJbAPOVtVPe3tMmx83xr+sBscYEy++qsFR1eeA5yKOXR92vQk4M9FxGWOMMSZ1BPo+xRhjjDEmuViCY4wxxpiUYwmOMWbIEZHjReQjEVktIj/2Oh5jTOxZgmOMGVLCtow5AZgKfENEpnoblTEm1izBMcYMNQcAq1X1U1VtAR4BTvE4JmNMjFmCY4wZasYBG8Jul7nHuhCRS0SkVERKKysrExacMSY2LMExxphuqOpcVZ2tqrOLioq8DscY00+W4BhjhppotowxxiQ5S3CMMUNN55YxIpKBs2XMPI9jMsbEWNLuJm6MMQOhqm0icjnwPDu3jFnucVjGmBizBMcYM+R0t2WMMSa1JHyzzXgQkUpgnXuzEKjyMJz+sFjjw2KNj4HGOlFVk7pKN+I9pi/J8m9qccaWxRlb0cbZ4/tLSiQ44USkNFl2LrZY48NijY9kitVLyfI6WZyxZXHGVizitCJjY4wxxqQcS3CMMcYYk3JSMcGZ63UA/WCxxofFGh/JFKuXkuV1sjhjy+KMrUHHmXI1OMYYY4wxqTiCY4wxxpghzhIcY4wxxqScpE1wROR4EflIRFaLyI+7uf8qEVkhIstE5GURmehFnG4svcYadt7pIqIi4tkSvmhiFZGz3Nd2uYg8lOgYw+Lo6//ABBGZLyKL3f8HJ3oU570iUiEiH/Zwv4jIn9yfY5mI7J/oGMNi6SvWc90YPxCRBSKyb6Jj9Ktof88TTUTGu78Hod/ZK93jI0XkRRH5xP1a4HWsACISdH9nn3Vv7yoiC93X9VF3ew2vYxwhIk+IyCoRWSkiB/vx9RSRH7r/5h+KyMMikuWX17O795qeXsMBv0eqatJdcNqrrwF2AzKApcDUiHOOAnLc65cCj/o1Vve8POB14B1gtl9jBfYAFgMF7u3RPo51LnCpe30qsNajWA8H9gc+7OH+E4H/AAIcBCz0Is4oYz0k7N/+BC9j9dMl2t9zj2IrBvZ3r+cBH7u/D7cAP3aP/xj4rdexurFcBTwEPOvefgw4271+Z+h32uMY7we+417PAEb47fUExgGfAdlhr+OFfnk9u3uv6ek1HOh7ZLKO4BwArFbVT1W1BXgEOCX8BFWdr6oN7s13cHYM9kKfsbp+DvwWaEpkcBGiifVi4M+qWg2gqhUJjjEkmlgVyHevDwc2JjC+nUGovg5s6+WUU4AH1PEOMEJEihMTXVd9xaqqC0L/9nj7e+U30f6eJ5yqblLV993rdcBKnD9+p+D8ocb9+jVPAgwjIiXAV4C73dsCHA084Z7ieZwiMhznj/M9AKraoqo1+PD1xNmOKVtE0oAcYBM+eT17eK/p6TUc0HtksiY444ANYbfL3GM9+TZO9ueFPmN1h9vGq+q/ExlYN6J5XfcE9hSRt0TkHRE5PmHRdRVNrDcC54lIGc6+Q/8vMaH1W3//P/uFl79XfpMU/4YiMgmYCSwExqjqJveuzcAYr+IK8wfgR0CHe3sUUKOqbe5tP7yuuwKVwN/cqbS7RWQYPns9VbUcuA1Yj5PY1AKL8N/rGa6n13BAv1/JmuBETUTOA2YDt3odS3dEJADcDlztdSxRSsOZpjoS+AZwl4iM8DKgXnwDuE9VS3CGOB90X28zSCJyFE6Cc63XsZjoiEgu8E/gB6q6Pfw+deYBPO0ZIiInARWqusjLOKKQhjO18hdVnQnswJlO6eST17MAZ+RjV2AXYBjg1QfSfovFa5isb/blwPiw2yXusS5E5Fjgp8DJqtqcoNgi9RVrHjANeFVE1uLML87zqNA4mte1DJinqq2q+hnOfP4eCYovXDSxfhtnvhlVfRvIwtnAzW+i+v/sFyIyA2cK4RRV3ep1PD7h639DEUnHSW7+oapPuoe3hIb53a9eTTeHHAqc7L4PPoIzlfJHnOmINPccP7yuZUCZqi50bz+Bk/D47fU8FvhMVStVtRV4Euc19tvrGa6n13BAv1/JmuC8B+zhVoNnAGcD88JPEJGZwF9xkhsv/6P1Gquq1qpqoapOUtVJOHUNJ6tqqd9idT2FM3qDiBTiTFl9msAYQ6KJdT1wDICITMFJcCoTGmV05gHnuysFDgJqw4ZpfUVEJuC8UX5TVT/2Oh4fieb/oyfcOpZ7gJWqenvYXfOAC9zrFwBPJzq2cKr6E1Utcd8HzwZeUdVzgfnAGe5pfohzM7BBRPZyDx0DrMBnryfO+99BIpLj/h8Ixemr1zNCT6/hwN4jvaiejsUFZ8rhY5yVCz91j92MkxwAvARsAZa4l3l+jTXi3FfxaBVVlK+r4EyprQA+wK3G92msU4G3cFa0LAG+5FGcD+PMgbfifPr7NvA94Hthr+mf3Z/jA4///fuK9W6gOuz3qtSrWP126e7/ox8uwGE4Q/3Lwv7dTsSpb3kZ+MR9vxzpdaxhMR/JzlVUuwHvAquBx4FMH8S3H1DqvqZPAQV+fD2Bm4BVwIfAg0CmX17PHt5run0NB/oeaVs1GGOMMSblJOsUlTHGGGNMjyzBMcYYY0zKsQTHGGOMMSnHEhxjjDHGpBxLcIwxxhiTcizB8TF3t9pXw27fKCJVHobUhYgME5FHRGSrOLugX+h1TL0RkUtE5GvdHF8rIrd5EJIxQ5qInNXd+4aIvCoiT3TzLcZELa3vU4yP3A0843UQYS4Fvgqcj9NVco234fTpEpx+EE9FHD8VsI68xiTeWTjdxe+LOH4ZTn8UYwbMEpw4EZEgEFRnZ+GYUNUynIZIfrE38JGq/tOrAEQkW1UbB/MYqro4VvEYYwZPVVd4HYNJfjZFFSMicp+IlIrI10RkOdAEHCgixSJyr4h8KiKNIvKxiPzCbece/v3jReQ595y1IvKdbp6jyxSViFzoTg3lRpzXZcpFRA4TkTdEZLt7WSIiZ/bx8xSKyP3u9FODO2Q8O+z+tTidJ2e6MfTYMdJtr/1zEalwn/9eETnb/b5J7jnvish9Pbyui93rR7rf82URmSci9cAc976rReQ9EakVkS0i8oyI7B72OK8Cs4ALQvGGhsa7m6Jyh84/EJFmEdkgIr+Unfu3hL/200XkRRHZISKrROS03l7XHl6f0M91pIg8LiL17v+Xy7p5LUpF5CsissL9d/m3iIwUkd1FZL4bR6k4+0UZ41vu7/vpwBFhv5M3uvd1maIKvfeJyIHu/+9GEXlTnK0xRovIU+7vzUoRObqb5/qOiCx3f5/XiciPBhDvJDfGc0TkQRGpc9/Tbog4r0REHnPvaxSRNSLy8/4+nxk8S3BiaxJwC/Br4ATgM5zh123AVTg7ud4KfAv4v9A3iYjg7LkxDSdpuAq4Ejh4sAGJSD7wLM5+Uafj7EHyIDCij299CvgycA3wdZz/K/PDkoZTgedw2oAf3EesPwCuA+50n78R53UKdw9wRniy5l4/A7i3m3OXAie718HZfG0Ozu65FwNBYIGIDHfvv8yN9bmweP/dXbAi8iXgUeB99/H+z30d5nRz+kM4+6ScitNe/BERKen2VejbXe7PdSrOlh1/FpEDIs6ZgLMdxf/iTLkdAszF2ZzwEZzXK82NQwYYhzGJ8HOcfZEWs/N38u5ezs/B+b/+e+AbOL8LD+K0/H8TOA1nqvxxEckJfZOI/A/wF5z3tJPc6z8XkcsHGPetQAPO79pdwA0i8v2w+x/A2RjyEpy/A7/E2SLBJJrXe2WkygVnDlmB/fo4Lw04B2eEJ8M9dqL7vQeGnTcRaANeDTt2I1AVdvtC9/tyI55jLXCbe322e05eP36W493vOSLs2DCcjSr/GvEz97oXEU6isQn4c8Tx59znmOTezgd2AN8KO+cioBkY5d4+0v2e30fxnNlAHXB+2PFS4L5uzu98vdzb7wDzI875EdAOlES89heFnTPK/Tf7Xj//74R+rpvDjqW7r/dvIl7vNmBy2LFb3O8N/zlD/5+meP17YRe79HbB2Yn71W6Ovwo8EXb7xm7eky5zj10fdmyqe+wE93Y+UA/cEPH4NwObccoIoo11kvvYL0QcvwsnsQq4t+uBr3r92tpFbQQnxspVdUn4AXd65gfulEIjTuHcP3Ay+gnuaQcAW1R1Yej7VHUdsCgGMa3B+YV7SEROEZERUXzPAUCFqr4WFs8OnJGgw/r5/OOBsXx+Z+Uut1V1O86b3YVhhy/E2SQ1sgD4cyMvInKQO1W0FScJaABycXY7j5o4tVP742xCF+5RnFGsyJGqF8J+hq1ABc5o0kCEP1YrzohQ5GOtVdXwYu7V7tdXujk2boBxGONHLcAbYbej+b9/MM6Hs8dFJC10cb9nDAP7Xf1XxO0ngV3CHmsJ8Gt3GnsCxjOW4MTWlm6O/QC4DeeX4hSc5CE0nJnlfh2L84cxUnfH+kVVq4HjcEYEHgMq3bqN3Xr5tuIennsLMLKfIYx1v1ZGHI+8Dc500xdFZDcRmQx8kc9PT4Xi6OS+ibyAs+Psd4FDgS/g/AxZn/vu3hXivFaR/5ah25E/f03E7ZYBPGd/Hqu7cyKPh44NNA5j/KhOVTvCbn/u/77uXNQR+r9f6H5djvPhMnSZ7x4fP4A4It8bQ7eL3a9fxxkt/j2wTpyax2MG8DxmkGwVVWx1V2h7Js5Q609DB0RkasQ5m4HR3XzvaJx6lZ40uV8zIo4XdAlK9R3geBHJBo4FbsepHTmoh8fd1EM8Y3Dqifpjs/u1KOJ45G1U9XUR+QRn5EaAjYSNaoSfGnH7eJz5+VPckSbcT2n9TcYAqnDeACN//jHu1/7+/MYY74R+X0+i+w+gHw3gMSPfG0K3NwGoajlwoYgEcD7Q3gjME5EJ3YxGmziyEZz4y8apIwl3bsTt94AxInJg6IA7KrF/H48dWjI+Jez7DsSZd/4cVW1U1WdwRkUik6xwC4HRInJ42OPmAF/BKebrjw04Sc4pEcdP7uH8e4ELcHrrPKCq7VE8RzbQgTM1FXIWn0/g+xxdcZ9vEU5iGu4s9znejiIeY0x0BjPiGY23cT4k7qKqpd1c6gbwmKdG3D4NJ7np0sJDVTvcD5c34XwAmziA5zKDYCM48fcicIWILMSphzkX2D3inOdwVs88LiLX4iREN9H3FNW7OMVtfxKRn+GMWPwI2B46QUS+glOs+xSwHmdu+rt0nbfuQlWfF5EFwKMi8mOcJnjX4CQSt/b9I3d5rHYRuRW4VUQqgbdwkpvp7ikdEd9yP/ALnP+bf4vyaV7BKSz+m4jcA+zjxlsTcd4q4Msi8mX3Z/qsh09UNwDPi8jfcFYmTcdZ8XGXOr2IoibOUvS/Abuq6tr+fK8xQ8Aq4BRxOoyXARtVdWOsHlxVa9yl538UkYnA6zgf7PcEjlLVU8Fp1YAzbXWUqr7ax8PuIyJ/Bf4JHI6z8vVKVe1wV20+j7OS6mOcWsurcT7krYzVz2WiYyM48XczzjLGX7hfW4Arwk9QVcX5o78CZwTj9zhLknsdLXDnm0/FSRKewPlFuhSoDjttNc6Uzq9wpntuAf6Lk/T05ms4ydkfcApuBThaVVf38j09+T3O0vnLcN4UCtx4ICwZc3+mzTgjSG+p6sfRPLiqfoAzrXUgTiH0OTgjMLURp/4C503mMZxRs6/28HgvAGfjrEB7BqeO6nfAQJaV5uD8m9cM4HuNSXV34Lwv3YvzO3lJrJ9AVW9h55Ltp3Heh8+la8FyaFl5NHWPP8IZJf8nzofFn7OzhUQT8AFOm495OB/YGoAv6SAbkpr+E+dvqzGJJSJ3A8ep6sSI4yNxRqUuV9V7uv3mJCIi9wMdqvotr2MxxnRPRG4CDlfVo3o5ZxJOb7OvquqziYrNDJxNUZm4E5FpOCsLFuCMNp2A0+zw2rBz8nDqgq7E6V/zcOIjjYuDcRosGmP86xCcxRcmhViCYxJhB07/nMtxelKsw0lufhd2ziycOfB1OE3rGhIdZDyoar/68BhjEk9Vj/M6BhN7NkVljDHGmJRjRcbGGGOMSTmW4BhjjDEm5ViCY4wxxpiUYwmOMcYYY1KOJTjGGGOMSTn/Hz2ZNC6XY7TNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "with h5py.File(\"sim.h5\", \"r\") as in_file:\n",
    "    time = in_file[\"/particles/all/position/time\"][...]\n",
    "    positions = in_file[\"/particles/all/position/value\"][...]\n",
    "    potential_energy = in_file[\"/observables/potential_energy/value\"][...]\n",
    "    molecules = in_file[\"/parameters/vmd_structure/resid\"][...]\n",
    "    box_size = in_file[\"/particles/all/box/edges\"][...]\n",
    "    n_particles = positions.shape[1]\n",
    "    n_steps = time.shape[0]\n",
    "\n",
    "    \n",
    "unique_molecules = np.unique(molecules)\n",
    "n_molecules = len(unique_molecules)\n",
    "G = np.zeros(shape=(n_steps, n_molecules,), dtype=np.float64)\n",
    "\n",
    "for i in range(n_molecules):\n",
    "    ind = np.where(molecules == unique_molecules[i])\n",
    "    chain_length = len(ind[0])    \n",
    "    \n",
    "    # make molecule whole\n",
    "    r = np.squeeze(positions[:, ind, :])\n",
    "    for j in range(1, chain_length):\n",
    "        for k in range(n_steps):\n",
    "            for dim in range(box_size.size):\n",
    "                dr = r[k, j, dim] - r[k, 0, dim]\n",
    "                if dr > 0.5 * box_size[dim]:\n",
    "                    r[k, j, dim] -= box_size[dim]\n",
    "                if dr <= -0.5 * box_size[dim]:\n",
    "                    r[k, j, dim] += box_size[dim]\n",
    "    \n",
    "    r_mean = np.squeeze(np.mean(r, axis=1))\n",
    "    dr = r - r_mean[:, None, :]\n",
    "    dr_squared = np.squeeze(np.apply_over_axes(np.sum, dr**2, [1, 2,]))\n",
    "    G[:, i] = np.sqrt(dr_squared / chain_length)\n",
    "\n",
    "hist, bins = np.histogram(G, bins=15, range=(0.2, 1.25), density=True)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2)\n",
    "fig.set_figwidth(8)\n",
    "fig.set_figheight(6)\n",
    "font_opts = {\"fontsize\": 15}\n",
    "\n",
    "axs[0].plot(bins[:-1] + 0.5 * np.diff(bins), hist, \"-x\")\n",
    "axs[0].set_xlabel(\"radius of gyration, nm\", **font_opts)\n",
    "axs[0].set_ylabel(\"probability density\", **font_opts)\n",
    "\n",
    "axs[1].plot(time, potential_energy / n_particles, \"-\", marker=\".\", markersize=7.5)\n",
    "axs[1].set_xlabel(\"time, ps\", **font_opts)\n",
    "axs[1].set_ylabel(\"potential energy per particle, kJ/mol\", **font_opts)\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}