Skip to content

Instantly share code, notes, and snippets.

@dwhswenson
Created October 2, 2018 21:24
Show Gist options
  • Save dwhswenson/9830820dbd45aa44c362d66d8c3b7a13 to your computer and use it in GitHub Desktop.
Save dwhswenson/9830820dbd45aa44c362d66d8c3b7a13 to your computer and use it in GitHub Desktop.
Simple OPS LAMMPS example
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# LJ gas $\\to$ droplet transition\n",
"\n",
"This is based on a model used in [Bolhuis and Csányi, PRL **120**, 250601 (2018)](https://doi.org/10.1103/PhysRevLett.120.250601). It consists of 16 Lennard-Jones particles in a box $12~\\sigma$ to a side. Condensation into a droplet is a rare event in this system (due to the low density).\n",
"\n",
"We'll use a time step of $0.005~\\tau$ and an inverse temperature $\\beta \\approx 2.4~\\epsilon^{-1}$ ($T = 0.42~\\epsilon/k_\\text{B}$).\n",
"\n",
"The stable states can be distinguished by their potential energy: the gas phase has $V \\ge -5 \\epsilon$ and the droplet phase has $V < -30 \\epsilon$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"\n",
"import openpathsampling as paths\n",
"from openpathsampling.engines import lammps as ops_lammps\n",
"\n",
"import nglview as nv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Set up the engine"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"lammps_script = \"\"\"\n",
"# Lennard-Jones droplet-gas system\n",
"# Based on Bolhuis and Csanyi: https://doi.org/10.1103/PhysRevLett.120.250601\n",
"\n",
"units lj\n",
"atom_style atomic\n",
"atom_modify map array\n",
"\n",
"region box block 0 12 0 12 0 12\n",
"create_box 1 box\n",
"create_atoms 1 random 16 830606 box\n",
"mass 1 1.0\n",
"\n",
"velocity all create 0.42 87287 loop geom\n",
"\n",
"pair_style lj/cut 2.5\n",
"pair_coeff 1 1 1.0 1.0 2.5\n",
"\n",
"neighbor 0.3 bin\n",
"neigh_modify delay 0 every 20 check no\n",
"\n",
"timestep 0.005\n",
"fix 1 all langevin 0.42 0.42 5.0 8464867\n",
"fix 2 all nve\n",
"\n",
"variable fx atom fx\n",
"\n",
"run 100 # mini-equilibrate\n",
"\"\"\""
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# the PDB file is optional, but will enable interactions with MDTraj/NGLView\n",
"pdb_file = \"\"\"\n",
"CRYST1 12.000 12.000 12.000 90.00 90.00 90.00 P 1 1\n",
"HETATM 1 C1 XXX A 1 0.250 0.750 0.000 1.00 0.00 C\n",
"HETATM 2 C2 XXX A 2 0.250 1.250 0.000 1.00 0.00 C\n",
"HETATM 3 C1 XXX A 3 2.250 2.750 0.000 1.00 0.00 C\n",
"HETATM 4 C2 XXX A 4 2.750 2.750 0.000 1.00 0.00 C\n",
"HETATM 5 C1 XXX A 5 0.750 1.750 0.000 1.00 0.00 C\n",
"HETATM 6 C2 XXX A 6 1.250 1.750 0.000 1.00 0.00 C\n",
"HETATM 7 C1 XXX A 7 1.200 2.200 0.000 1.00 0.00 C\n",
"HETATM 8 C2 XXX A 8 1.250 2.750 0.000 1.00 0.00 C\n",
"HETATM 9 C1 XXX A 9 2.250 0.250 0.000 1.00 0.00 C\n",
"HETATM 10 C2 XXX A 10 2.250 0.750 0.000 1.00 0.00 C\n",
"HETATM 11 C2 XXX A 11 2.250 0.750 0.000 1.00 0.00 C\n",
"HETATM 12 C2 XXX A 12 2.250 0.750 0.000 1.00 0.00 C\n",
"HETATM 13 C2 XXX A 13 2.250 0.750 0.000 1.00 0.00 C\n",
"HETATM 14 C2 XXX A 14 2.250 0.750 0.000 1.00 0.00 C\n",
"HETATM 15 C2 XXX A 15 2.250 0.750 0.000 1.00 0.00 C\n",
"HETATM 16 C2 XXX A 16 2.250 0.750 0.000 1.00 0.00 C\n",
"END\n",
"\"\"\"\n",
"with open(\"topology.pdb\", \"w\") as pdb_out:\n",
" pdb_out.write(pdb_file)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"engine = ops_lammps.Engine(\n",
" inputs=lammps_script,\n",
" options={'n_steps_per_frame': 200,\n",
" 'n_frames_max': 500000},\n",
" topology=\"topology.pdb\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x10e43e650>"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXmQK3d5NXy6tUuz7+udfebO3H2/dgiYQAiQxP4CxHyJKYdAwh9vYkxwknLK1FflIoSLgWCbSsgLJgZsXiiCvyqIQ1Jg53OceMe+3u69s0mjGY00m2a0b61evj+GX9+fNFpaUqs1o9unijI1d9RqadRHT5/nPOdhJEmCDh06dOioPdhan4AOHTp06NiFTsg6dOjQsU+gE7IOHTp07BPohKxDhw4d+wQ6IevQoUPHPoFOyDp06NCxT6ATsg4dOnTsE+iErEOHDh37BDoh69ChQ8c+gbHE39fH+nTo0KGjdDBKfkmvkHXo0KFjn0AnZB06dOjYJ9AJWYcOHTr2CXRC1qFDh459Ap2QdejQoWOfQCdkHTp06Ngn0AlZhw4dOvYJdELWoUOHjn0CnZB16NChY59AJ2QdOnTo2CfQCVmHDh069gl0QtahQ4eOfYJSw4V06CgISZIgCAIAwGAwgGEUZaro0KEDOiHrUAmiKEIQBPA8j1QqJf+cYRgYDAb5fyzLgmVZMAyjk7UOHVnQCVlHRRBFETzPy1UxwzAy4UrSblorIerFxUX09/fDZrPJv2cwGGA0GnWi1qEDOiHrKAOSJEGSJKTTaYiiCAAykRISJj+j/0tLGZIkyVU1x3EZj9GJWseNCp2QdSgGIVGe5/cQsRLQhJ3vcfmIGkCG9EHkD52oddQTdELWURTZRExIsBwipCvoXFBC1JIkZfwOTdDZOrUOHQcJOiHryAvimOB5XibBSirSSghSKVHPzs5ienoaAGSCztVQ1KFjP0InZB17IEkSUqlURjXMspVb1rM1ZjWQTdSJRELWqInOzXFcxu/oRK1jv0InZB0yJEmSHRNzc3Nob29HZ2dnSccIh8Pwer2wWCxwOBxwOBwwm805m37VRKGKGgB4nkc6nc74N52oddQaOiHrkImY53kA173DpZBnIBCA0+kEy7Lo6ekBx3HY3t7GysoKOI6DwWCQpQ9RFDOIWktkOz8IChE1eT9o14c+9KKjGtAJ+QYGPcwBZFaVSqpZSZLg9/uxtLQEs9mMqakpNDQ07CE0YJfo5ubmwLIstre34fF4kEqlYDAY5Eo6u6LWEsWImtbSV1dX0dnZCZvNpg+96FAVOiHfgMg1zJFNICzLyta2bEiShM3NTSwtLaGhoQFHjhyBw+GQ/y3bBQEARqMRFosFra2taG9vl3/O8zzi8ThisRh2dnYyiNput2cQtcVi2RdEHY1G0dnZCYZhMr7Q6MewLAuj0agTtY6SoBPyDYJCwxy5wLLsngpZFEWsr6/D7XajpaUFJ06cgM1mU3wOuZ7LaDSiqakJTU1NGT+niToQCGB1dRWpVAosy+6pqLUm6kKOE/KeiaKYMUIO6EMvOopDJ+Q6R7nDHAzDyL8viiK8Xi9WVlbQ0dGBM2fOwGKxFH3e7OcopamXj6gFQUAsFstL1ES7ttvtsFqtVSO6fMfNJ30A+tCLjuLQCblOUekwB8uy4HkebrcbXq8X3d3dOH/+PEwmU8HHFXJTqOGyMBgMeYk6Ho8jEAggGAzC6/UimUyCZdk90kelRJ3ry0YJShl64XkeOzs76O3t3aNR686P+oVOyHUGepjj8uXLOHXqVMmVVjqdxtbWFkKhEIaHh3HhwgUYjZV/VKppezMYDGhsbITJZMLY2Jj8c0LUsVgMoVAIPp+vYqIul5DzIdffJ5VKwe/3o7u7O+d0om7Rq0/ohFwnyDVVl0wmASifkOM4Dm63G1tbW2hqasKhQ4cwMjKi2jlq6UMmIETd2NiY8fNqELWakCRJJtlc/6YPvdQndEI+4KCHOQgRk4uYEGCxizGZTGJpaQmBQABDQ0MYHx/H5uYmYrGYqudaC0LOh0JEnUgkEI1GEQ6Hsba2hmQyCYZhYLPZ4HA4kEqlkEwmq2rPE0Ux73SkPvRSv9AJ+YAi1zBH9gVMrGv5Lux4PA6Xy4VIJIKRkREcPny4JB9yuee9n2EwGNDQ0ICGhoaMn4uiKFfURFsnjbnsiprkPVeCcmSRcoZeyF1Vc3OzTtT7ADohHzAUGubIRj4vcSQSgcvlQjKZxOjoKI4cOVKSD7lcHOSLnGVZmai9Xi+OHDkCo9GYQdSRSATr6+uyVJRL+lCaCVLoi7RUFCLqUCgEv98Pi8UifwkQuUQfetEeOiEfECgZ5siGwWCQfx8AQqEQXC4XBEHA6OgoWltbS/IhV4r9JFmoBZqoaYiiiEQigVgshmg0io2NDSQSCQCQpQ+6os4mX7Ubh7lArI1kgIV+bgD60EsNoBPyPkapwxzZIFXuzs4OXC4XWJbF2NgYmpubiz6W9iGrhXohZCVkSQ+w0KCJOhaLYXNzMydRp9NpTQguVyVeTPrQh16qB52Q9yGIL9Xv98NkMsmaZCkfbEmSwHEc3njjDTgcDkxNTe1pYBVCJZJFIR+y2iRfC1RSvSol6p2dHSQSCQQCAUUVdbkQBEHxsfShl+pDJ+R9hOxhDmI/y754ix1jY2MDS0tL4HkeExMT6OnpKflcbtSmnhJU4zVkE7XD4UA8HsehQ4eQSCRknXprawvxeBwAYLVaM4jabreXTNSiKGbIFeWi2NAL8cWfPn1a/kLLpVHf6A1FnZD3AfJt5jAYDIorSlEUsba2huXlZbS2tuLkyZNYWloqOuKcD3pTrzCq/VrIdCVN1HQ2tSiKSCaTckXt9/vLImqiIVcL9OQmkTSAwiu5bmSLnk7INUQuIqYvnOymXC4IggCv1wuPx4POzs6MnIlKSLVahFwPFTJQfUImBJYPZIDFbrdnELUkSRnSByFqSZL2SB92ux2CIMBsNlf1tQC7n1O6Ei9UUd/IQy86IdcAhYY5aBQiZJ7n4fF44PP50NPTkzNnopQKOxvVIE+tCFkLh0K1UUleRiGijsfjiEaj2N7eRjweRzKZhNVqRTQazSBqNWQMGjzPK6rEb/ShF52QNYSSYQ4aLMvu+fCl02msrKxgfX0d/f39BXMmWJYtWmEXem69Qq4d1PQhA5lE3dHRIf98YWEBDQ0NMJlMcjMxHo9DFMWc0ke5RJ1dIZdz/vR/CYoR9UMPPYR77723aCjWfoFOyBqglGEOGnSFnEql4Ha74ff7MTg4iJtuuqnoBVupZFEqeUqShLW1NbjdbgDIuJiJT1cnZGXQqsonUkZLS0sGUUuSlKFRV0rUlRJyPhQj6ieeeAL33Xef6s9bLeiEXEWUM8xBw2AwIJVK4dq1a3LOxMTEhOLKqVLrWikNRRJc39bWhuPHj4NlWTkTguQWx2Ix2Y6XvbJJRybUrpDzIR9RkuwOm82Wk6jp5QGxWEwm6uzpRHLsahFyPtB3YwdJutAJWWXQTQmn04mOjg40NTWV/KGIxWJYXl5GKBTCzMxMRs6EUlS7qUc7O9ra2uSGIvkSyg7v8fv9CAQC6OrqQjQaxdbWFtxuN9LpNEwmk1xJk4u5mt3//Q6tCLnU56GJml7FJUkSUqmUXFF7vV7E43EIgiBvdGEYBuFwGHa7XdO/rU7INyBybeZIp9Pgeb6kDwSdM9Hd3Q2z2VyWjxjIrUErRSG9VxRF+Hw+LC8vo6OjA2fPnlVU5RLNvLm5ec+0IMdx8sW8traGWCwmX8y07FGNhtN+hFaSRSmDIYXAMAysViusVmtOol5ZWUEymdxD1NnSh5pEzfP8gfus6IRcIQpt5jAajXuyAPIhFArB6XRCFEWMjo6ira0N0WgUoVCo7HMjkkc5yEXINBF3dnbi3LlzJckNhUjebDbDbDajtbVV/ll21eXxeGQdk7ZwNTQ0qDq9th+gZYVcTdIiRG2z2dDQ0IC+vj4A1ydJo9Eo4vE4fD5fxpdwtvRRDlGHQiFFMQH7CTohl4l8wxx0VVPMRyxJEgKBAFwuFwwGw56ciUpsa0BlLgv6ddA79ZQScaU79cjv56u6cnltASCRSMDtdsvShxpRmLWAVhWyllo1PaTEMAwsFgssFsuev22+uyWz2bxnwW0hog6Hw2hpaanq61IbOiGXiGLDHDTyEbIkSfD7/XC5XLBarXlzJiohVPL4SghdkiR4PB6srKygq6ur5Io4G2rZ3vJ5bUVRxMsvvwy73S5HYSYSCXmIgpB0Q0NDVcPl1UCxwRC1oJZkoeR5lFTiNFG3tbXJPy+HqIPBoF4h1yuUDnPQMBgMGRounTPR2NiIY8eOwW63F3x8LQhZFEXZFZFMJhUtN1WCavuQyYBAV1dXBtmSdU204yOVSsFgMOxpJO4XxweRv7TAfiLkfChE1Ol0WpY+1tfXEYvF8K1vfQtvv/02jEYjvvWtb+HIkSM4efJkwestG5/4xCfw5JNPoqurC2+//TYAYGdnBx/96EfhdrsxPDyMH/3oRxkyW6XQCbkISh3moGE0GpFIJPbkTJw6dQpWq7Xo47UmZEEQsLq6itXVVXR3d8PhcGBiYqLs598vyLeuied5ueLKdnysJEz4+VIKEmPA/3WyF++b6db0nLWSErRCNX3IZrMZbW1tGUT9j//4j3jsscfw8ssvg+M4PP7442hpacHMzIziY3/84x/Hn//5n+POO++Uf3bp0iW85z3vwb333otLly7h0qVL+NKXvqTa69EJOQ8q9RCTxwQCAbzwwgvo7OxU7EagH18JlBIyIWKPx4Pe3l65It7c3Kzo+bOx3yb1jEZjTsfH5eUd/ONP5yFKIiBJ+OJ/zMPpdOLmQ9eziiORSFUdH/Uw/k2jVj7kU6dO4c/+7M/KOsY73/lOeciJ4Cc/+QmeeeYZAMAf/dEf4ZZbbtEJuZogRJxKpXD58mWcP3++5AuD5Ex4PB6YTCbVbvlLRTFCFgQBHo8Hq6ur6O3txcWLF1WxHeV7v7QgZPIclZDZfy4GAYZBi223CRXnBFxLOvCJw+Py9g8i6WQ7PsqNwcxGPVbIWvvKg8EgBgcHVT3mxsYGent7AQC9vb2qFy06ISP3Zg6j0QhBEEq6sNPpNJaXl7GxsYH+/n4cP34cKysrNZujz0fINBH39fWpRsTFsN8q5Hwwsgzo05Sk3Z8RDdNsNmN6evpX/7Y7uRaNRvc4PojVq5zlp/VWIdfCE6y7LA4Ycg1zlCNN0DkThw4dknMmiAG+VsjWoEnl7vV60d/frxkRExwUQr71eDeemvMjlODBMADLMPjDc305f5eeXMt2fBBrXiHHh8PhkCfZaGhRIWv5t6hFxR8Oh1V3WXR3d2NtbQ29vb1YW1tDV1eXqse/IQmZWNcEQdgzzFEKEokElpaWEAwGMTw8vCdnotKmXKUgFXK5RKzG7X82DgIhj7Tb8dBHjuBf39oAx4t4/0wnjvU3lXQMOlievmiJ44PkQGQ7Pojro9S7s3KgNUlqXfGHQiHVK+Rbb70V3/3ud3Hvvffiu9/9Lm677TZVj39DEbKSYQ4liMVicLlciMViGBkZwfT0dM5jqEHIJOSnnAuHrKh/8cUXMTAwgJtuuqmk20by3Grdah6UChkARjvsuPvdI6ofV4njw+/3IxaL4bXXXoPJZMrwT6uZ8aElIdfi7x4KhSqypP3BH/wBnnnmGfj9fgwMDOD+++/Hvffei9tvvx3f/va3cejQIfzLv/yLimd8gxByKcMcNLLJMBKJwOl0guM4jIyMoKOjoyCZq0HIpMot5cLheR7Ly8tYW1uDJEklE3H2c6tJyAcd1SKWbMdHKBTC6dOnIQiCrE+vr68jGo3mzIGgk9WUot4ah9molJB/8IMf5Pz5008/XfYxi6GuCbkSDzEAOYsiHo/D5XJl5EwogRoEREhdSVVEh9eTivill14qm1DVrmgPUoWcD1rmFDMMA5PJhNbW1j0ZHyQHgiSr0RGYtD5dyPGhlRWtVg3KWCxW0oLg/YC6JGTimKjEQ0yq6suXL8NsNu/JmdAKSqps2t0xODiIixcvqnKhlTvpV0vbWz2h0PuYLweChMpHo9E9jg9a9rDZbDXPXK4myOfsoN0B1BUh05s5XnjhBdx8881lEfHW1haWlpaQSqVw+PBhdHdrO6VFoxAhp9NpuN1ubG5uKt4iUgrUXuNUD4S8n+1o+ULl8zk+RFGEJElYXl6WiTqX46NS1IKQgfIKsVqjLgg531RdKX8MSZLkrRdNTU04duwY3G63KtkGlVzEuUiR4zgsLy9jc3Mzw2anNsohUEmS4PP5EA6H0djYmBGNqRNybZDP8bG9vY319XVYLBYEg8Gcjg86jKlc1IKQtcwCURMHmpCJlpbPQ6zk4iEZvysrK3tyJkrJM86HShtjdIXMcRzcbje2trYwNDRUNSImKKVCJl9oS0tLaGtrQ1NTE+Lx+J7b5kQige3t7QORuJYLB/0LJRtWq3XPAgTSNyHbqZeXl5FOp2E0GjP0aYfDoWjoqRaETAqCg4YDTcj0gsPsC7sYEQqCAK/XC4/Hg66urpw5E2oQMiHUSgg5lUphbm4Ofr9fEyImULLolE6wa2lpwZkzZ2AymcBxXMY5iqKIcDiM2dnZDP8tfZGrbeuqFg7al0g+5NOQjUYjmpqa0NSU6b1Op9OyPr2xsYFYLAae5+X4S5qs6c97LQg5GAweuCk94IATMpC/ijOZTEin03s+CGRIwufzoaenp2DOhBq2tUqOwXEcdnZ2sLa2hvHx8ZIWnNIo9za70KJTorU7nU40Nzdn3FnkegzLsmhsbITJZML4+Lj8c/oipzNuabfAftoIchAli3wolShNJhNaWloyiC47pzjb8eFwOOTejpY2u2oMhWiBA0/I+ZBd3XIch5WVFTlnQokTwWg0lr2Tjj5GqYScSqWwtLSEnZ0d2O12DA4Oor+/v6znr0QyyfVlR8L1nU4nGhoacPLkSdhsNkXHy6Uh57vIabfA1tYWEomEHExfy6D5eiJkNQiyUE4x+Ruur68jHo/j1VdfBZDb8aH2e3oQ1zcBNwAhk5yJ7e3tkp0IJM+4EhgMBsWyB03Ew8PDmJqagsfjqcjpQNZAlUPI2QS6vb2NxcVF2Gw2HD9+vGDYd64LTGlTr5BbgFRigUAAHo8HHMfJskdDQ4PsO98vQfP7GdWUEui/IcdxaGlpwcDAQIbjIxqNYn19HclkEgD2NBIrcXwcxGAhoA4IuZBX0+l0IpVK5cyZUAI1NeRCSCaTWFpaQiAQkImYvC6WZSs6B7IGqpzEOVIhBwIBLC4uwmw248iRI2hoaCj7fCoBkT2ymzVkY0QsFkM6ncabb76ZMSRBV2KVVoRaVMhaNQ5FUdQkiZAebKIdH9nnQhqJoVAIPp8PyWRS/n36rshkMhX9G+gV8j4ByZnY3t5GT08PTp8+XfYFVG1Cpol4ZGQEhw8f3nOuBoMBHMeV/fyVeIk5jsPCwgJsNhsOHz5ccde6WrY3epptY2MDp06dAsMwGbGYm5ubqsgeWhByvWycJshecJoLLMvKX540sh0fKysr8l1RdiOR/nIJBoM1nR8oFweekMnFEQ6H4XK5wHEcRkdH5W/HSi4eNQg51zGSySRcLhdCoVBeIiaoxV69cDiMxcVFxGIx9Pf3Y3R0tOTnraXWShaE5orFzN6vl0v2IBd5LrLSqkLW4r3bbwtOc6GY4yMWi+1xfDzxxBNYW1vD6dOnVRmf/trXvoZHHnkEDMPg2LFjePTRRxWtYCsHB56QI5EIrl69CgAYHR2VZ/7X1tZU0X/VdFkkEgm4XC6Ew+GCKXFqnkMphByJRLC4uAhBEDA+Po5AIFC0ssmFQkM5tW6I5UtbI7JHNBqFz+fb4/YgJK0F6rFCVtvKWMjxcfHiRXzve9/Ds88+i5/+9KcQRRHPP/98WZ89r9eLhx9+GFevXoXNZsPtt9+OH/7wh/j4xz+u4qu5jgNPyAaDAePj43v0IrUcEmpIFvF4HFeuXEE4HMbo6ChmZmYUfzgqHV9W8vhoNCqn2I2Pj8tfasFgUNXR6f2MfCE+RPaIRqPY3NxENBoFx3G4evVqRkWtZhORVPjVhlbEr9W2EOL4+O3f/m3867/+K+655x6cPn264jsOnueRSCRgMpkQj8fR15d7WYEaOPCEnG91u8lk0qQhVwjxeFyutmZmZkoiYrXOoRAhx+NxLC4uIpFIYHx8PCOohjxWbc33IE265doGEo1Gsby8jMHBQUSjUezs7Mi6pslkknXNQrJHMWg19nsQJItyQUdvVvJe9vf34y//8i9x6NAh2Gw2vO9978P73vc+tU5zDw48IeeDGhVyuX9IEtcZjUbR0dGBxsbGsle9qEHI2Y9PJBJwOp2IRqMYGxvLm+uc67FKUU9+3WwwDFNU9vB6vfIKL9p3S4ZcCr039ShZ1GJ0Wg2XRSAQwE9+8hN5EvX3f//38fjjj+NjH/uYCme5FweekPN9sNWokEtFPB6H0+lEPB7H6Ogojhw5gmAwiLW1tbKPWalkYTAY5KqUbiaOjY3hyJEjBYmhEsvdftWQK0WhL5p8sgftu6XdHrTnlpY9tPoyq+f4zUgksqcRWA6eeuopjIyMyHdIH/rQh/D888/rhFwq1NB/lSIWi8HpdCKRSGB0dDSj4qy0wlWjQk4mk7h27RoCgQBGR0cVNROBwqPTNypKJUtis7Pb7XndHnSAj8lkkrNAwuFw2bKHEmglWdRiM4kkSaq8b4cOHcKLL76IeDwOm82Gp59+GmfPnlXhDHPjwBNyIbuYGmRSaKddNBqFy+VCIpHA2NgY2tvbc/qIa0XIHMfJaWtTU1MF7XW5UA0N+aBDreo1n9uD4zhsbm5iY2MjIxfCZrPtyfao9Dy0kiy0lq/U/MxeuHABH/nIR3D69GkYjUacOnUKn/rUp1Q7fjYOPCFXG6TSphuHxJWQSqUwNjaGtra2feUjpoPrGxsb0dPTg97e3pKfW6+QtQdJTmtqasLExASAvbLHxsYGEolEzim2Utwe9arz05vk1cD999+P+++/X5VjFUNdEHI1g89pQo5Go1hcXATHcXJFrPTx5aKUDxXP83C73djY2JCD630+X9lfCGpvDKkH1GIwpJDsQYYjsmWP7EjTfJWwFoSsNelHo9EDmYUM1Akh50MhuUEpjEYjQqEQ5ufnkU6n5YpYKdQYLikGnuexsrKCtbW1PQFKLMuW7TbRK+S92E+j0waDIecUW74FqET2ICSthRxVC8krFAqp0tCrBeqCkPNVyLnkhlIQiUQQCAQQiUQwPT1d1krxal68giDA4/HA6/XmjRStREsvV0MOBALw+XxyZkRDQ8O+D51XCq1IrJLPjdlsRltb2544zEQiIdvySCTmL3/5y4pkj2KolQf5IAYLAXVCyPlAQupL/YCFw2E4nU4IgoDW1lZ0d3eXRcbVgiiK8Hg88Hg86Ovrw4ULF/ISXiVe4lIli3A4jIWFBbAsi56eHiSTSWxsbMjvpdVqRTKZxObmpmqNqVpgv1TIpYCWPQheeeUVnDp1Stamc8ke5H92u70sYtW3hZSGuiDkfBdIqfotCdURRRFjY2NobW3F0tKS5n7mXJAkCZIkwev1YmVlBT09Pbh48WLRyrMSHVhphRyLxbCwsACe5zExMYGmpiZwHLdnv2EymcTly5flQBjSmKIv/P1eTddLuBD5uyqRPVZXV/PKHsW+VGtVIeuEvA9BKuRiCIVCcDqdkCQJY2NjGX9MLTRgJVhdXcXKygq6uroKrp3KRiWEXExDTiaTciocPXqdi8TJGLLJZMLw8LB8EfM8n+EeoKtpmqT3SzW9nzTkSlDsdRSTPSKRCNbX1/e4PQhRk7tSnZBLQ10TcrEKORQKYXFxEQByBhSRY6RSqYrOo9zmoiRJ8p65aDSKc+fOlSy/VKNC5jgOLpcLOzs7iib+aBC9n/y+0WhEc3NzxnufHepDqmmynr5QNb0fSLtS7FfSp2UPOgqAuD2i0Sj8fj/cbrcsFZJrMBKJlC17lIpQKKTIAbUfUReEXGiUNVeFHAwG4XQ6wTBMXiImMBqNiMViFZ1fqc1FssnZ5XKhra0NLS0tGBkZKavZomaFTNvqsjeblIJiMkiuUB/y/MWqabJMs1oXvlZkWW3ZJlflGkqkkeJFtDvMMLDKX2Mh2cPr9SIYDMLj8SAej0OSpD3ZHlarVdX3NBKJlJXhvR9QF4ScD9nVbSAQgNPphMFgkLVOJcfQKjUue5Pz6dOnYbVa8frrr9fES0weS9wcq6ursr+5WHVVKMuiXKeCkmo6mUzitddey9gooaY2XS8aMl0hS5KEH722hv+4ugWWAQZabfiLd4+gxV7Zeiez2Qy73Q6WZTE0NCQ/Vz7Zg3Z6kFVN5UCXLGqMQk09shBzcXERRqMRk5OTJXkUtSDkYpucKyHVSjRwcvG8+OKL6O3tVdRELAa1iSa7mg4EAjh+/DgAHChtmoYWGjKdY/GGN4x/v7KJzkYLWAbwBBJ47GUv7rplWJXnoStxJbLH1tZWhuyRPeRS7L3RCXmfIpFIyMb4cnfCGY1GVbeGZINscrbb7Xk3OVdCquWQOZFMCIldvHix5GolX5VXzalKGmpr0/Qx6q1C9gZTAMPIMkWzzQSXvzKZjkDptpBcsgfZAEKIOlv2oImalj3oLOSDhrog5OwP7/b2NpxOp7yl+OTJk2Uf22AwqFIhZx+D3uR89OjRguuBtCJkulJvamrC6dOn8frrr6u6mVgrQs733OVq06Sa3q8Nt3Keg1SuXY1mSJIEUZLAMgwiyTSmu9UZPeZ5vqw1YMD1DSAWi2WP2yMejyMWiyESiWBtbQ3JZBKJRALf/va3sbGxgStXrqC1tbXi5l4wGMSf/Mmf4O233wbDMPjnf/5n3HTTTRUdsxDqgpCB3T/Szs4OnE4nLBYLZmZmYDQaceXKlYqOq9aiU0KowWBQlk+mp6f3bNnNBTV04GIIBAJYWFiA1WqVK3XifVYTtSTkfCilmhYEASaTSf6ydzgcqjfgtCB9WrI4c6gZvz7WhuddATAs0G43448uDqj2PGq/P3SWNC17JBIJSJKEz33uc3jqqadYHVOqAAAgAElEQVTw4IMP4tixY3j44YfLfq67774b73//+/HjH/8YHMchHo+r8RLyoi4IWRAEvPLKK7BarZiZmZFJjuf5islUjYAdg8GASCQCj8cDhmFK1rErqZCLEWAkEsH8/DwYhsH09HSGrFMpKeQilv1IyLmQr5r2eDxywPz6+jqi0WhGNU1IuhJtWusKmWUYfPLmQXzwaBdSvIi+ZissRnWeX0sfss1mwy233AKWZfHwww9X/PkNh8N49tln8Z3vfAfAbpNSzbHyXKgLQjYajTh+/Pie1dxqDHVU+kclt1SSJOHYsWNlzdhXSsi5EI/HsbCwAI7jMDExoWoThEQf5iPeg0DI+WAwGGC329Hf3y//jK6maeeAwWDIkDyUVtNaLDnNJn2GYdDXrP5qe60HQ9S8q3O5XOjs7MQf//Ef44033sCZM2fw0EMPVXX7eF0QMgBZ36NRy+45iepMp9Po7OyE2WwuO/CkkjyKbCSTSTidTkQiEYyPj6Ojo0OV4yrFfnM0lIp8VX8hbZqQNKmmiQ+3sbExpw9XiyWn9brglHCAGu8fz/N47bXX8PWvfx0XLlzA3XffjUuXLuHzn/98xcfOh7oh5GpDqa5HNjknk0mMj4+jra0N6+vrFQ2XGAyGihe2chyHpaUlbG9vY3R0NGMDtjeYwH8vbINlGNwy1YGuxvKaMEpwUCSLfChF382nTZOwebohRVfTqVSq6u9RvS44jcfjqlWwAwMDGBgYwIULFwAAH/nIR3Dp0iVVjp0PdUPIxW6RK/nGJDpyoQ8W2eQci8X2rHNSY41TMpks67E8zyOVSuGVV17B8PAwJicnM94L51YMf/r4ZcQ5AQwDfOt/3PjOx0+jv8VW4Kjl46ATMlBZ9ZUvbJ6upuPxOGZnZwEgI8xHzak2LaYBAe0JORgMqha92dPTg8HBQczNzWFqagpPP/00ZmZmVDl2PtQNIecDcUlUYt0ix8j1wSISQDgcxtjYGDo7O6uyV6/UxiKJ6FxdXQXDMHmn6/73fy8hwQmwm3dfWyTF49HnV/C5D06Vfb6FcNAJuVoOCLqaDgQCGB8fh9VqzZhqy1VNl6JN06jXBadqD4V8/etfxx133AGO4zA6OopHH31UtWPnQt0TMsmzUIOQaT9lKpWCy+VCMBjcIwHkenyle/WUPp5EdC4vL6OnpwcXLlzAyy+/nPfcgvE06NgChgECca7sc1V6jgcVWvqQ80218Twv2/FI+BTRppVW01pJFoC2fQO1w+lPnjyJX/7yl6odrxjqhpDVykQudgxaix0ZGVG0ybnS4RIlFTYdSNTe3p6RDJedsEbjPYc78bYvAl7YrcBZhsFvTHXu+b1yIUkSXvOEsR3jcLSv8cA39bRAsarSaDSipaUloxLMzogoVk1rXblqhYO8LQSoI0LOB6WZyIVgNBrl7Rd+vz+nFlsIakgWhR7v9/uxuLiIxsZGOZCIBtHAc12Av3+6H+Ekj3951QuGYXDnxUF88Gh3xu+UUhWKooiVlRV4PB4AwLeuiHhjiwfLMBAl4J4Ljejvr16FXA9jzeU8R6nVNPHo8zxflcQ1Aq3vhg5yjgVQR4RcrQqZ53mEQiFsbGxgbGwMFy9eLLmyqJSQ8w2nBINBLCwswGw2583BKPT43X9j8KfvGMafvmM4578Xqq5pkOzmpaUl9Pb24ty5c/ivhW286Z9DkpcA7F6YX3sphPGGebS0tMi2L7vdfmAq54M2Op2vmn7rrbfgcDhU16ZpaGHfy4ZOyPscRqOxrAqZ3uRMBgEGBsobJ1VbsohEIlhYWAAARaFJamRhFCIIv9+PhYUFtLS0yFJJOp3GZjSN7AIpzgMDg4dgMu5OL25tbSEej8vxi4SkGxoaNN80oQT7tUIuBQzDgGVZdHR0ZHyJq6FN09BSpyYIhUIYHBzU9DnVRN0Qcqkh9fkgCAJWVlbg8/nkTc4+n6+i8elKqx1CqLTHeWJiQnGiVbWyMEKhEObn52E2m3HixIk9FfrR3r05HR02Bk+7ouhobsCZQ30YHr6+6ocmg2g0ClEU5c3VhKjNZrOiav2goxZVuBraNI1q5FgUg14h73OYTCZFgSC0Tayvrw8XL16Uv92NRiMSiUS1TzUvSOXyxhtvyNN1pVywaq9xImPX6XQaU1NTeXM5Znob8Ve/OYZLP18EA6DJasBvDFthNzHgRQnPLGzjt6Y70Wg1wmAw5ByiiMfjiEajCAQC8Hg84DhOzsglRG2z2TRrUGlRIWsBpf5gpdo0+QKlq2mWZTVvHIbDYZ2Q9zOKaciiKGJ1dRUej0e2iWV/q6vh1CgH6XQaS0tL8Pv9MBgMuHjxYllkoNYaJ47j4HQ6EQwGMTExoWjs+vbTvfi9E92IpQT855wf8aAfRoaBw2xENJVCMJFGozX3x5BO9eruvt5oTKVSMhn4/f4MySOZTCIcDqO5ubkq1Vm9EHKlOrWSajoYDCIajeLVV1/ds1m8WlIG+dsfVNQNIZfa1BNFET6fD8vLy0U3OWtNyLR+PTQ0hIsXL+LFF18smwgqrZDT6TRcK6vYWN/A4YlRRVY/GiYDixY7C5vZgKB4vdoWJQnGEna3EZCMXDrrlkgewWAQGxsbcLvdEARBljxINa1E8iiEepBDgOoMbGRX04FAAH6/HyMjI0WrabWcHgc5nB6oI0LOh2wNmbgB3G43Ojo6FG1yVoOQlWyeJrKJx+PBwMBAhmxSCcolZFEUEY3Fcff/+SVe3dytDN+9sY0HPtQDs7H0C+fsoWb8y+o6NqMcolIS/S021XIziORhsVgwMTEBk8kkV2yRSAShUAherxepVAomkymjeUj2vilFPVTIQPVfB5FFKtGmS62mdQ15n6BYhZy9yfns2bOKs03ViPEkx8h14UuSBJ/PB7fbje7ublV21+V6bqWQJAmbm5twOp349yUOr2/v2uMA4NkFP77xXy7c/Z7xgsfI9ffobLTg14cdEAxWdHe2o6vRUtJ241JBV2y05MFxHCKRCKLRKLa3txGPx2V5hCbqXH+DciQLSZIgSCjrbuAgI1/cAJBfm06n0/L2lnKq6VK2u+9H1A0hA7lzEgwGg7yos6WlJefgRDGoueiUlkVo4mtra1NUrZeDUirkQCCA+fl5OBwOnD59Gg+/8TJ4QYDRsPtFwosSXnQHcHeZ59JoMcBuN6GrQPbutfUI/scZAMsA75pox3inuvmzZrMZ7e3teyQPEu5Dr3EiRECIulRvrcsfwzMLO0ilRQy12fDuqXbYTLW382khvZQTLGQymcqqpg0GQ8Zi4IOKuiJkGvR+OJ7nceHChZKJmEAtQqaPQZabEuJTcm7lNpSUEHI0GsX8/DwA4MiRI/LWlZ4GI675r0s+DAMMtlb2wS9EBnMbUXzzOQ9spt0vgCtrUfz5u4Yx2pF76EUt5FuySYiASB7hcFju5NODLbnufHZiHH5+zY8WmxFtdhM8wQT+e3EH75tWbzR9P0MQBFUKDCXV9PPPP48HHngAa2truP3223H8+HHcdttt8gbycs//7Nmz6O/vx5NPPlnx61CCuiNkSZLkJadkk/Prr79eNhkDua1fpYIEDIVCISwsLMBoNBZdbpp9DuUa7QsRcjKZxOLiImKxGCYnJ/c0RP7weDPe9gsIJXfjORssJtzz3omSz4GgWNrbc64ArEYWrfbdOwl/lMNL7kDVCTkXchHB3Nwc2tvbwbIsotEolpeXMyQPupoOxHcHY6y/qog7HGYs7xS2T9ZL0xCofvQmXU0PDAzgAx/4AG699VZ84QtfwJtvvgmOqywk66GHHsL09DTC4bBKZ1wcdUXI5HbbarXuIbta25VEUcTs7CyMRmPJO/WA65KHWoTM8zxcLhf8fj/GxsZw5MiRnO9Ps9WIR//vSSyGWUgAzg+1oiGPTS0b+bZrFIKRASTkdmJIkoRfzPrxnCuAdocJf3i2Hx0N2uqFkiTBbDajqakpYxMykTzozdUb0TTWNhkYOTtsNhs4GNBkKXy+Wo1ma3Et1CoLeXJyEpOTkxUda3V1Ff/2b/+G++67D3//93+v0hkWR10Rcjwez1hySkDITOupIWA3uH5xcRE7Ozs4dOgQRkdHyzpOpePPxGlCwn+8Xi8OHTpUNJuDZVlYDAzerVICXLEK+Z0T7XjTF8FmJAVJ2k2fu2l0t2p/7GUv/vkFD0RJgiQB/zm3je/eeQLNtvKjVUtFPsLMJXmIogjDlXW8sRqAPx6FkE7hXCfw2msbchXd2NiYIXloveC0mtCakNV0WHzmM5/BAw88gEgkosrxlKKuCHlwcDAnaRHrW6WEXEr1Quclj4+Py7e+5aISLzFpbPp8PiwtLeUdgFH7eXOhGCEPt9tx97tH8JI7CJZhcGG4Bf0tu3LTYy97YWQZGH5FWMFEGv/jDOC3j3blPV4twbIs3ne0FyeH2sEJItrsJjgsRnAcJ/tyieQBAA0NDbDZbOB5vuIM70LQKnqzFoSsxlDIk08+ia6uLpw5cwbPPPNM5SdWAuqKkPNBTZdEMRIj03VbW1sYHb0+RBGPx6ueiZwP0WhUnkQs1cmh9oYPhmGwEkjila11sCyDk/2NGGzL/KIabLXlbBwKogQD9X0o/epnWqJUSYFhGHQ3ZXqtzWYz2tra9kge8XgcgUAAPM/jzTffhCAIsFqtGVY8NYYn6nXBaTAYVKVCfu655/DTn/4UP/vZz+TJz4997GN4/PHHVTjLwrghCFmtTGSe5/MSsiAIWF5ehs/nw9DQ0J6VSdXORM6FcDiM+fl5iKKIzs5OTE9Pl/y85VbI+UjDG+LwX0tRDHZbIErAk1e28HsnutHTVLzp+ttHOvHTtzbBiru+XruJxcWRzAvwoCaxGQwGNDY2wmQyIRQK4dixY5AkCclkEpFIJMPuZTQaM0ja4XCURLBaSRaFrpdqQC3J4otf/CK++MUvAgCeeeYZfOUrX9GEjIE6I2SttobQoLMw+vr6cNNNN+X8sBsMhoq6vqUQI0mFS6VSmJychCAI2NjYKPt5S62QJUmSQ9ANBoMc9wgAS8E07CYGDZbdj15aELGwGVNEyHf/xiha7GY859pBm92E//XO4apuyK4F6IYbwzCw2Wyw2Wx77F5ksMXj8cgbzbMHW/JJHvUsWRzkKT2gzgg5H9QiZLpCpQPZu7q6imqylZ6DkgqZ4zi4XC4EAgFMTEzIm6+DwaAq4ULFQIhbEARIkgSDwQBJkiCKonzuLESkRQmStHtMQZRgMigjByPL4JM3D+KTN9cu77baLgglxzeZTHskD1EU5cGWra0tLC0tged5WK3WDCue1WrVjJAlSdI07S0cDqOnp0fVY95yyy245ZZbVD1mIdQVIauViZwL9Aj21tYWnE5nRiB7MVRTsiByydraGkZGRjA1NZXxXqgRLlQMhHhJhUcyDADIhLy8vIxWMYRNcyvWgklIDGAxMBhrs4DnefnirfQirqaXt9qEXC5ZsiyLxsbGjGUFRPLInnAjr8Hn86GxsbFkyWO/Qq+QDwiMRiNSqVRFxzAYDAgEAvJ03cmTJ0sa1VRjjVP24+kN0/39/Xt060KPVYpiTT2aiMnv04RFMkTcbjf6+vrwm79+EReSAlxbUTAsg6FWKxothowqmvyXHIuWPGqNag9uqFlV0pJHZ+d12+L6+jq2t7chCMIeyYOupqvl8qgWdELeZ6hWhRwKhbC2tgaj0Yhjx47t8TkrgRoVMiE9UqUvLi6ivb29YHQoUNmkYb7quhgRA7uDOgsLC2hubsaZM2fkO4kWO4vTQ7kjEkVRlI9N/gtcl0FYlq05SVe7Qq52Y5JlWTgcjoxVR0TyiEajOSUPQtTVWoaqBnRCPiAoV7+NRqNYXFwEz/Po6emRN1VoeQ4EBoMBqVQKwWAQ8/PzsNlsijMw1NipR0CTJbn1zb5AyUYRSZJw5MgRxePh5PnIORMUImn6cVoQxX6VLEpBrmYbLXn09vYCyJQ8otEo1tfX5VCfbJdH9vFqteD0IGchAzcIIZdaIScSCTidTsRiMUxMTKCtrQ3r6+vyrV05qLRCTqfTWF1dRSAQwPT0dNHFpjQq3RhCyBe4XqnmqlDT6TRcLhdCoRDGx8czmk6VQAlJS5KEaDSKRCIBQRDAcZzs8KCPUSn2Q1OvUigl/XySRzqdlkl6dXVVvi7sdntG81DrBacHfX0TUGeEXKntjV5RNDY2hs7OTvmYWrgkciGVSmFxcRGBQABNTU04ceJEycdQY8kp3bDLropJsD7xYE9OTmpy2w3svq9kGCcUCmFmZkZ2EuSTPCppHmpByFqMTlfiDzaZTGhtbc2oRkVRRDweRyQSwfb2NkKhEGKxGN54442Matpms1Xt/eM4rqIQsf2AuiJkIHcTqhiZ8jyPpaUlbG5uYmRkJOeKIjWacqXouOScyMRfT08PNjc3FT8+kuThCyVhZBkMtFjLImRiXdvY2ADP83JWg8PhkN/nra0tuFwueQ2WllURaWp6PJ6CXwTky4RU0vu5eajFrb4gCLBY1PVvk52GRNKLxWJYWlrC+Pi47Jne2NhAIpGQc4wJUeeSPEpFvaTk1R0h50I+MhQEASsrK/D5fBgcHMzrUgC026tHb78eHByUw39CoZDiLwR/NIUfvrKKRFqEKEnob7FisIQRY7ph53A4cP78edk65Xa7EYvFIEkSOI6D3W6X5QktyWxnZweLi4tobW3FuXPnClZ8ubYfl9s8rAcNWSud2mg0wmq1wmq1ZkgeZGN1JBKB1+uVP0+05NHQ0FBWlvJ+bTgqRd0RspLsBVEU4fV6sbKygt7eXkW766pNyJIkYX19Pe+gSSmyw38vbEOUgL5fhfKsBhKQhOKPzdewo7MXSH5yKpXC0NCQPAXodDoBQK58mpqa0NjYqHrFTBqGAHD06NGyA5vKbR7Seno1sJ805EpQaEov1469bMljeXkZ6XQaFoslo5rOJ3mkUqkDvbqJoO4IuRBEUcTGxgaWlpbQ0dFR1C5Go5qEvL29jYWFBTQ1NeHMmTM5bydLkUwiSR428/ULzmRgkOQKkwi5lc/XsON5Hm63G9vb2xgbG5OnAGmIoohoNIpwOIy1tTU5R8PhcMgETfIaSgWRcHZ2duRGq9ooRNIcx8HtdsuETJrERNZRq3lYa7JUC4X26eVCtuQB7L63qVRKrqazJQ9C1A6HQ7Wkt1rjhiBkQjYvvfSS7IktVUNTi5DpCigSiWB+fh4GgwHHjh0raA8rhZAnux14Zt4PczMLXpTA8RI6bbkv8mJ+YlEU4fP54PF4MDg4iHPnzuUlDJZlc2YCx+NxhMNhecKR53nY7XaZpJuamvJWN2REfXl5GYODgzh//rymt6UMw8iDLYODgzh8+LBcIdPadHY1Xa4urVVTTyvJohIwDCNLHh0dHfLPieQRjUbh9XrxxBNP4IknngAAPPDAAzh58iRuvvnmsiyqHo8Hd955J9bX18GyLD71qU/h7rvL3SBZOuqOkPMNJ/A8j6NHj5btU1SDBEgeRjqdxsLCghz+o+SbnR4MKYZzw21I8RIue0IwGxjcerwHQbc/43eUDHaQnYRkAWs5F1i+yofcngYCASwvL4PjONhsNrmKbmpqQiKRkAdLzp49q/nkGEnLa2hoyHj+fNVwMclDCUlX6oBQgoNehWdLHp/73Ofw3ve+F//0T/+E/v5+/OIXv8Dg4GBZ6YZGoxFf/epXcfr0aUQiEZw5cwa/+Zu/iZmZGbVfRu7n1+RZaoBwOIyFhQWwLIuZmRk4nc6aj4IyDIO5uTmEw2GMj4+jo6NDMdGXMv5sYBm8a7ID75q8XlU87979rxIiJgtPTSYTjh8/rvo2X7J/zuFwyGEwdNTkzs4Orl27BkEQ0NjYCJZlsbOzg8bGxqrapgiI/TEWi2Fqakqx57uYLk2/7/mah1poyFrkIau14FQpotEohoeHcccdd+COO+4o+zi9vb3yYExjYyOmp6fh9Xp1Qi4X8Xgc165dQzqdzqg+1QgYKhckWCcUCqG9vR0zMzMlX3RqLFolkku+CbtUKgWn04l4PI6JiQlNNTnSPIxEIgiFQjh69Cja29uRSqUQDocRiUSwvr6OeDwOk8mU0TgkNrxKIUkSVldXsbq6iuHh4Zz2x1JRavMwkUjAbDZDEISq2fC0yEOuh+hNt9uNy5cv48KFC6oetxDqjpA5jsPg4CDa29szfl4NDVjJ75Lwn76+PnR3d5dUFasFSZJgNBpx7do1NDc3Z3iJgev2v42NDYyMjGB6elrTcyQOE7fbjf7+/gyd2mKxoLOzM8M2xXGcHNq+tbWFeDwuj/MSoi41wYxIW0psdJUiF0nHYjG5n9DW1rYntlTN5qFWFfJBJuRoNIoPf/jDePDBB0teSFwJ6o6Q29raclbCJpNJtZD6YtKHJEnw+/1YXFxEW1ub7Oa4du2aJl5m+jzIbfKxY8f2eIlZloXBYEAsFkN3dzfOnj2r+SLYUCiE+fl5NDY2ZgQQFYLZbEZ7e3vGly7P8zJJLy8vIxqNyvo1IemGhoY9JJFKpbCwsIB0Ol1y7oYaEEURbrcbW1tbmJyc3DP9Rv6rZvPwoDT1SkEoFCp7gXA20uk0PvzhD+OOO+7Ahz70IVWOqRR1R8j5YDQaVctELkTIdPhPdkRnpdN+SpFLJ87e4xYIBDA/Pw+j0YiBgQHE43H88pe/lEmMuCUaGhqqcvESPzPHcZieni47tInAaDTuGecVBEG24Xm9XkSjUUiSJG/WiMfj8hJaugLXCsTuSHYdZr/P+Ua8K20eauHkKNX2VinUqpAlScInP/lJTE9P47Of/awKZ1YabihCTiQSFR8jX4Ubi8WwsLAAQRDyhv9Um5CVNOzoJLajR4/uqQgFQUAkEkE4HIbH40E0GgWADDkgV6WpFERP39zcxNjYWFUlHIPBgObm5gwtnKzcWl5ehsVigcFggNPpxNraWoYNr5oN4GQyibm5OTAMg5MnT5acv1Bq8xDIJGkt5CitJQu1goWee+45PPbYYzh27BhOnjwJAPi7v/s7fPCDH6z42EpQd4RcKBM5HA5XdOzsNU7A9UZYOByW1yblgxqEnEvDVhKJqTSJzWAw7JmioivN1dVVmaRJJU3IutAFKEkSNjc34XK50NfXh/Pnz2ueG5FIJDA3NweWZXHu3DmZCCVJktcf+f1+LC0tIZ1Oy6O8hKQrzX8QRRErKytYX18v+lkpFUqbh+FwWLZeksepnYgH1EZDVqMJ/Y53vKOmuRh1R8j5oEZTz2AwyMcgk2ubm5sYHR1V1AhTKzGO1uaKTdiRitDr9ZadxJav0iSarc/nQzQahSiKGXIHIWni57Xb7Yp1YjUhCAKWlpawvb2dc8qPYRjZK93T0wMJAINdAo9EIggGg/B4PEilUrBarRkOD6WB7UQi6uzs1OzLKFciXjwex4kTJ+R9h9VIxCPH0pqQD3oWMlCHhFwoglMNDTmdTmNlZUWeXCPhP0pQ6eZpmpCVrE6qZhIby7I5SZoenZ6dnUUymQTDMOjt7UVnZ6emVTFdlWe7N3LhBdcO/t/X18DxIs4Nt+L2033o7raju7tbPh5tw/N6vUgmkzCbzRl3Cna7Xf5b0E3Dani6i4Gs0FpaWspr5cvVPKS/5IHSm4da6NQ0dELex8gVMFSpy4JMl3k8HvT39xfdMp0Lau3VEwShoE5MKlKbzYZTp06pHrVY6PyIxryysoJQKITDhw/DbrfLWQQLCwtyvgVdaebTbBNpAWuhJAwsg/4WK4wKL/JoNIq5uTlYrVZFVfnCZgyPv+xFZ4MJLTYTXnDtwG5i8aFTffLv0KO8XV1d8s9TqZR8t7CxsYF4PA6j0QiGYRCPxzE8PIyBgYGaSDSzs7Mwm80F34NqNQ+1RDKZ1PzLrhqoS0LOhUoq5J2dHSwsLIBhGAwMDGB8fLys49CSR6kgFcf6+jo6Ojpgs9n2XAS0c6GUCTO1QFfl3d3dOHfunFyV58u32NzchNPphCAIe/ItYmngey95EEzwgCRhpN2Oj57rh9mQ/+InWnk4HFY8lg4ATn8MBhawmHbPt73BjLd9EXzoVPHHWiwWWCwWOW8hFAphdnYWdrsd7e3tCIVC8Pl8GWuSqulgEUURy8vL2NjYwNTUVNmVYyXNQy11WPJcBz16E7iBCJloZqWAhP+wLIujR4/KgSblIldTsBjoCmV4eBibm5tYWFiQJ7rIEEQoFJItXLUYPiHj1haLpWhVTudb9PXtVqDZjTWXy4X/bzmJ9YQBva12WCwWzG9E8NZqGGeG9nbT6RCiQ4cOlayVN1uN4IXrTdE4J2CorbSKi+M4LC4uIplM5nSwkFCcbAdLKc3RYggGg5ibm6uaVq2EpFdWVjImY9VOxMsHnZD3KSr9hk4mkzLpTU5Oyo6DZDKp6RqnbC2PXLD0eS4tLWFhYUG+HfV4PAiFQiU3ncoFyX2IRqMlVaTZoBtrJEvgcsoNQywFVkrvhhEF43jp9R0Yg9erzKamJtlG1tTUVHYI0ZlDLXjJHcTiVhQsGNjMBnzoZK+ix0qSBJ/Ph5WVFYyMjKC7uzvne54rB5g4WHI1R8nfW0lkKQmsSiaTOHbsWNk50eWAECyJLWhvb8fJkydlia2azUNg94tO64GmaqE+XoVKILe7Ozs7e3bqAZVJDuTxSghZiZ94e3tbngT8tV/7NfkDmUql8NbyFu7/iRMbUQ4jjQw+dsyB3vZmmcAsFkvFJE27N/KtvaoU450OPBtMYqClATa7AwmDA+8914+Bxl3nht/vx9WrV+X1UiaTCcFgsCyLmtnI4n+9axiLmzGkBRFD7Xa02IoTeyQSwezsLJqamsoauc7nYCF3C9mSDq27m83mjLHz4eFh9PT0aF4pCoIAl8uFYDCImZmZjCEfmmjp5mH2aHglujQpQOoBdUnI+T6QDMPkHBulVzkVsoZVe9Gp0iS2hYUFGI3GnF37uMDg//nFKk/s/uMAACAASURBVGIpEWaDCdfCEv6P04D7R5oynAEWi0Um6FJJmoyFk9viatmb3jHejkiKx5urYTAMg9+a6cJkdyNEUYTf70cgEMDU1BQ6Oztl90O2RY224BV7jWYDi5leZbp7Op2W7wwOHz6sql5Pa820pEN09+3tbbjdbiSTSaTTaVitVgwPD9fEZbCzs4P5+Xn09fXh7NmzBd/fcpqHhdZpEVQjWKhWqEtCzgeia5Hqidxqut1uRaucKiXkfI9XQsR0JGShJLara1FwvAi7efd1GFgJV9ZjsDe3yc4A2r5Fhj1yEVj2BBkJwDEajWVNmJUKs4HFbcd78cEj3WAZBgaWkf287e3tGV8GZF09bVEjcZ6hUCinj7icuwW6Ih0aGsLU1JQmFSkdWdrd3S037UiDORKJYG1tDalUSv6yJa+zGrIVnel94sSJihwO5a7TIkRdL9tCgDol5EJeZJ7nYTab5SqPpHspGVaolJCzIzSVTNhlJ7EVkwYsRhaidH2iT5R2hxxoZ0Iu+xZN0tkE1tDQgFgshlQqhampKc2rEZOBRTKZxJVfrYRSopEyDCOTdL4vouy7hWIERqx0dru9JoH5wPWmHfGWEzLL9kqT8fe1tTW5AUw7PGivdKkgMko1JRIlJE3+/89+9jN4vV7Vz6EWqEtCzgej0YhgMIirV6/CYrHsCf8pBjWzKIpN2NGVWG9vr+KO+cnBJkz3NOBtXwSCKMFoYHDnxUGYjYUfm4ukSe6E1+tFQ0ODHOFJNnvQVWa1QFu4iIOkXOTzEZNKOhwOw+fzZQx7EPJaW1tDMBjE1NRUTaoxpU07+jXmiiwlq7ToyFI6V7rQZyyVSmF2dhYsy9Zk4jKbpDc3N3HPPfeAZVk89NBDmp5LtcCU6Eao3ZB3CRAEYU8lG4/H8eqrr8JgMODo0aNlNwGef/553HzzzWWf2/PPP48LFy4UlCdINm9TUxNGR0dL/uBzvIh/v7KJ9XASR3ob8WtjbSVXMaRp2N7ejuHhYblZRaQAUmWGw2F5/RKtSatxsZIdfN3d3RgaGtJ0+IBU0mtra/D7/TAajbJPmhBYJVWmUtBfzIUcHOUgnU7LAy2RSCQjspS8RuKVJi6SiYmJir4U1YAkSXjiiSfw5S9/Gffffz9+7/d+7yBY3hSdYF0SsiiKsgeSeEPD4TAcDgc6OzvltUHloFxCJrdXb731FmKxmFyZNDc3yx96OoltYmJC82xeYPeLa35+HgzDYHJyUtEdhCRJcu4DTdK0K6AUko7H45ibm4PRaMTExETVtepC52A2mzExMQGz2QyO4+Sx6XA4vGd7SaVSQK5zmJ2dhc1mw/j4uCYSCUn7o/+WiUQCFosF/f39aGlpke+WaoGNjQ189rOfhcPhwIMPPljzL4cScGMTciKRgNvtxsbGBkZHR9HT0wOPxwOGYTA4OFj2sUsl5FwNOzrzgSYw4PpOr1I3XlQKnuextLSEQCBQMA1OKQhJk9cXiURkkqarTJqk6RCg7LB2rSAIAtxuN/x+v6JzoKWASCQiSwG0Jl3qiik6tL4Wmj05B9K7mJiYgNFolF9jJBIpafxdrfP58Y9/jK9+9av4/Oc/j9tuu+0gVMU0blxCjsfjePHFFzEwMIDBwUGZ2Hw+H1KpFEZGRso+9gsvvIALFy4UJUu66ZCvYUd7eQcHB+FwOGQCi8VissZH1i5V4xaZHmoYHBxEf39/1T7o2SQdDoflmEvSLR8cHMShQ4dqko1AGr29vb0Zn5tSQaQA8hpL0WsDgQDm5uZqItMQhMNhzM7Oor29HSMjIznPgR5/J6+V5/k94+9qSFfr6+v4i7/4CzQ1NeHBBx9UNbZUQ9y4hEwq5Ozbqq2tLQQCAUxOTpZ97FdeeQUnTpwo+EHLbtjlyi+mk9iGhoZy2u1yXdhGozFDq61kCzPRqltaWjAyMlIT1wC5+FmWhcPhQCwWQzqdhsPh0CwwPpFIZMg01ZBIaL2WfOESvzH5O/p8PqTTaRw+fLgmQTn0gEc5W1yIV5p+nRzHyXZKOldayWdWFEX86Ec/wte+9jV84QtfwO/+7u8etKqYxo1LyJIk5Yy5DAQCWFtbq2il9+XLlzE1NZWzy63ETxwOh7GwsACr1Yrx8fGSHQrkFjkUCsn6HnEEkEq62Ac+kUjIqWu10qrJYEUkEsHU1FRGk5UegiD/o8OH1LpFph0cagfGKwHP87L9zu/3w2QywWw2a7JCKxtkwKO/vx8DAwOqER/tByfVNLEa0ncM2YXF2toaPvOZz6CtrQ1f+9rXKpbQ9gFuXEIGdrvk2YhGo3C5XDh+/HjZx33zzTcxMjKSMZmlhIjpJLaJiQlVJ7tob204HEYymcwY8iAkTWu04+PjNbn1oyWSoaEh9Pb2Krr46fChfCTd1NSkuNlECKiW0kAsFpNT4UjTjl6hRZwPgHortLJBD3hoWZnTudKRSASJRAKJRAI/+MEPYLfb8cwzz+DLX/4ybr311oNcFdO4sQmZ47g9AUPJZBJXrlzBmTNnyj7u1atX0dvbi9bWVkVETDaL+P3+qu+QIyDDAaSKJrfIPM+jtbUVAwMDaG5u1txHGgqFMDc3h5aWFoyOjlbcqSckTTcOBUGQ5Q5SfdHPk0qlMD8/D0EQMDU1VRNpoNSmHb1CKxwO70mJK5ekNzY24HK5apaBkY3FxUXce++9iMVi6OjowNLSEj784Q/jvvvuq+l5qQSdkLNfG8/zePXVV3HhwoWyjzs/P4+WlhZ0dHQUbNhJkgSv1wuPx4OBgQH09/fXpAoLBoNYWFhAQ0MD+vv7M6QAWqsttcIsBRzHyVXY5ORkxRumC4FuNmWTNCG28fHxiqyPlYBU5j09PRU1L+mUOPI6AcgpcfQKrWyQdDyDwYDJyUnNv5izIYoivv/97+Mf/uEfcOnSJXzgAx+Qr6dc2TMHFDc2IafTablyJZAkCS+88EJFgx1OpxOSJGFgYEDeCpENOolteHi4Js0yWiLJR4L5Ksxce/HKAXGR+Hw+jI6O7knP0wqBQACzs7OwWq2wWCx79v8RAqvmDjjypcRxXNWkAdpOSdvT6CjPaDQKr9e7LwY8AMDr9eLTn/40+vv78ZWvfKVuQoJyQCfkbEIGKh/siEajWFlZkSsSesBDkiQ4nU4YDAZMTEzUrFO+vLyMzc1NOUK0FJDoR5qkRVHMGH5obGwsWrWQLSsdHR0YHh7WdOElAV2ZT01NZTQv871Otb6MCOjg/NHRUXR1dWn6pURe59bWFlZXVyFJUkbjkPxdtR70EEURjz32GL7xjW/ggQcewG/91m/VXDKpMnRCVoOQC+nEpAETCARkj7PNZkNbW1tVJreKnSdZZtnX11eRjzYbZMN0toZJfxkRXy1xcEiSpHjST23QclEpJEjIKxQK7akwyyFp0rRzOBwYGxuryZ0ScZJsbm7i8OHDaG5uzshbpu+M1Hax5MPq6iruuusuDA8P44EHHqibpLYiuLEJmef5nEFAzz//vKJN0UoadtlJbF1dXeB5PsOWFo/HS7allQqy1JRc+FpogrQbgJB0Op2GIAjo6+tDX19fyRNqaiAcDmNubg7Nzc2qNA6zpyqj0SgkSSrYUBNFUXaz1CqMCLju8SZ3KYU+87R8RQ96ZE/jVfLZEkUR3/3ud/HNb34TX/nKV/De97633qtiGjoh5yLkYoMdSibsspPYijVnaFtaKBSSYy0JQZc70ZRKpeQdbpOTk5ovNQV23wu/3w+n04nOzk60tLRkZD2QMWLyhVTJIEshZAfGV7txSJM03VAzGo3w+/3o6+urmZ1OEAQ4nU6Ew+GK3otsP3gkEpEnK0udxltZWcFdd92FsbExfPnLX67JZ7XGuLEJOVfiG1B8sKPQhB1w3bXQ2NhYVhIbeZ5kMplhS6MdD4So890a0wMNtWyWkcB6k8mEiYmJnEMu6XQ6wyNN3zGQ/1USoE5rtLW0byWTSVy9ehXJZBINDQ1IJBIAMv3DSrT3SrG9vY2FhQXVBzwI6Gk8OqOEpP3R03jA7mf10UcfxSOPPIKvfvWreM973nMjVcU0dELORchvvfUWhoaG9kyGFZMnqp3ERt8yEg2TbqYRnXZ7exsul6ti21Q58AYT+OJ/LMK9E0OPTcJHxhhcPH645M54rkEWEhJPyzrFEI1GMTs7i4aGhppptPSgS7Zenc8/TDdI1ZrES6fTmJ+fB8dxmJ6e1jQhL1fa31NPPYWf//zniEQiGBoawoMPPoiJiQnNzmkf4sYmZDqCk8bs7Cw6OzvR3t6uiIjT6TSWlpYQDAZVSUErBXQzbXt7Gzs7O2BZFp2dnWhtbS0rSaxcJNMC7vzOZWyGEzCABw8j+lvt+M4fncrYRlIOsjd5hMNhuUGaK2OZ53m4XC6EQqE9Y9dagm7ajY+PK9Krc5E02bhdDklLkoTNzU24XC7V85LLhSiKeOSRR/DDH/4QH/jAB8BxHF577TXcfffdeP/731/Tc6shFP1RbqiNIcDu1hDSfCpExHQS29DQECYmJjT/oLMsK4fO8DyPs2fPwm63y01Dl8uFWCwGk8mUUV1WY4fatVU/NkNRWE0GWMy7zpHtGAdfMInh9spWzudbKUWC8Hd2duB2u5FOp2EwGJBIJNDT01M05KlaIBGd5TTtcm2ZphukHo9HJmm6ks6VDpdMJjE7Owuj0ViTDR65sLS0hLvuugtHjhzB008/XZOclIOMG6pCliQJy8vLiMViGBgYgN1uz7k6iSSxdXZ21sxDK4oiPB4PfD5f0cqHBKcTuSM7y6KSMWnSLHNthvHV19JwWIxgGQaiKCGeFvCDT5xBZ2P1VjgRkGrUaDSira1NbjYVG5VWG0SjrTSisxhyuVjobdTJZBJbW1uYnJzcF3GUgiDgkUcewfe+9z08+OCDeOc731m1AuYTn/gEnnzySXR1deHtt98GsOt7/+hHPwq3243h4WH86Ec/qkmedgHc2JJFduIbadilUin4fL6MTQiksmRZFm63G1arFWNjYzXZVEG7FgpFcxY7RnaWBQmHp50dhYiL9vIODw+ju7sbX37KiV9c2wIvSjCyLG493o1Pv3u00pdcECQQaWdnB5OTk3v0ajIqTV4r7R0mr1WNMB6O4zA/Pw+e52uWgSEIgrxglGEYGAwGsCyb8WWk9WIDAHC5XLjrrrtw4sQJfOELX6h6Vfzss8+ioaEBd955p0zIf/3Xf422tjbce++9uHTpEgKBAL70pS9V9TxKhE7IJM+ikE6cTCbh9/uxsrKCVColb+dV4nZQG9FoFPPz8zCbzRgfH1f1C4G2MJGmIT0mTVZJGQwGBINBzM/Po7W1FSMjIzJxi5KE/17YxmowiaE2W1m7+krB1tYWFhcXZceAUqLJZ0srp5lGN+3GxsYylqNqiVwDHgBk33uuxQa03FGNv5MgCPjmN7+J73//+3JVrBXcbjd+53d+RybkqakpPPPMM+jt7cXa2hpuueUWzM3NaXY+CnBjEzIhnpaWFpmElSSxAdjjdpAkKedUmlqgs4EnJiY0m+fPtUoqmUyCZVkMDAygs7MTDQ0NmmvniUQiI/xGja3WhSQA8sWbTVzExdHY2IixsbGa7ZELhUKYnZ2V75iKffa0IOnFxUXcddddOHPmDP72b/827xbsaiGbkFtaWhAMBuV/b21tRSAQ0PSciuDGJuSXX34Z99xzD0KhEA4fPowzZ87g3LlzOHHiBEwmE5577jmYzWZFSWz0xRwKheQPOLmQm5uby2qkiaIIr9eL1dXVmnpoiV69traG4eFhWCyWPRczfcdQreEOEku5ubmJycnJqjtaChFXMplEIpHAzMxMzSbtBEHA4uIiIpEIpqenK5IC8q2VypY7iv1dBUHAN77xDfzwhz/Eww8/jHe84x1ln1Ml0Al5FweGkAnS6TSuXLmCF198Ea+88gqeffZZhMNhnDlzBrfddhvOnj2LycnJkmUJMvAQCoUyGmk0SRfyxZJEuPb2dgwPD9es+iJNqmKrpOimIdHe6aZhpVUsOY9a+KtpbGxsYGFhQb4LIpulae29Gi6WbFR7wAMoTtLZWSzz8/P49Kc/jfPnz+Pzn/98TXR0Al2y2MWBI2QaDz/8MF544QXcd9998Pv9MknPz8+jo6MDZ8+exZkzZ3D+/PmS/Zy0TYuQNM/zGc0lUnktLCyAZdmaJcIBmXvkyjkP8lqzfcM0cSkZ1Egmk5ifn69pGBGw27Sbm5uDKIqYmprK0O9pFwvdDM7eyKIGaabTaczNzYHneRw+fFjzxjL58qWJ+m/+5m/Q1taGxcVFfPGLX8Ttt99e84zibEL+q7/6K7S3t8tNvZ2dHTzwwAM1Pccs6IScDUEQclaAZPz25Zdflkl6c3MT4+PjOHPmDM6ePYtTp06VrKfSEY+BQADb29sQBAFtbW3o7OyU9WgtZQp6zf3ExMT/397ZB0V1Xn/8swiLi4DiG6K8CBGoLwiyUM2MMdqMpqNOkhFjmzEN1rHOWEuJHdtqaCppXrCT1GjTmsRYnRqtUdNOTI3FSU2ZMVEWzIIUEYog0QU0EnlZ3tnl/v4w9/52cUWBu9wVns+MM8Ife88F9uy55znf71GtLeC4UVqeeHAcSRs9erSTU5pjm2TatGmaefM6TpP05dDO1QfSQEYNHd36tLDpvBulpaVs2rRJmcUvKioiICCAffv2aRbTM888Q05ODnV1dQQHB/PSSy/x1FNPsWrVKq5evUp4eDjHjh3ztD18IiEPBLvdTllZGSaTCZPJREFBAV1dXcyePVtJ0jNmzLhnFej4hg8PDyc4OJjm5malipaFHY6tDrXd4OQ45Plqd8/Qyjh+IDU2NipOaXq9nubmZsaNG0d0dLQmkmdQ99DO8Qmp56ihYyXt6l7b29u5dOkSer2emJgYzX4ejthsNv70pz/xj3/8gz//+c8D2rIjAERCVp/W1lYKCgrIy8sjLy+PkpISAgIClASdnJzsNJ5148YNvvrqqzvGx3rS2dmpVJWyG1x/Hv/vRktLC2VlZfj6+vZr07VayG2BtrY2xo0bR1tbG83NzcpBmny/7vaQltfd19fX853vfMdt0uv72Z7d1NTE9evXe916/eXVBj6/fAv/kSN4YvYkJvi79/d36dIl0tLSePTRR9m2bdugtk3efPNN9u7di06nIy4ujv3792uiB3ADIiG7G0mS+Oabb8jLy8NkMpGXl8e1a9cYP348VquVhIQEXnjhhT5PT9zt8b9nP/p+xp8qKytpaGjQ1JdXkiQsFgsWi0XZYuL487DZbE6Hhj0d4dR8aqirq+Py5cuKif9gtwVkE6mbN29y7do1AHx8fO7YyCK3dj4ru8n2U5exd99+643x8+G91fGMHaW+TNpms7Fr1y4+/vhjdu/eTXJysurX6I3q6mrmz59PSUkJBoOBVatWsXTpUtasWTOocbgJ4WXhbnQ6HePHj2fp0qUsXboUgD/84Q8cOHCAZcuW0dzczI9+9CNaWlqYMWMGSUlJJCUlMXv27F6rVJ1Oh5+fH35+fsoyTsdtFhaLBavV6qTSGj16tFJZOlpShoeHa+LDISNvmg4KCuK73/2uyx6+LId27Pk5PjXU1NTcMcXSVw9pedt0d3c3CQkJmlVdshlQXV0dCQkJBAYGOv1ua2pqnPyV3z7ThBcSowzegI761i5Ol93k6cQpqsZVUlJCWloa3/ve9/j88881e4qy2Wy0tbXh4+NDa2srkydP1iQOrRAVssqUlpbe4fzV2dlJUVGR0o/+73//i16vZ86cOUqSnjZtWp97uq4qSy8vLzo6OggICCA6OtqtRu290dXVRXl5OW1tbcTGxg44jp5TLK48pF35WDj28KdNm9bnHYNq0heBh+wK99zBi7R02hhBN6Cj1QY/jB/Hcw9HqCJQ6urqYufOnXzyySfs3r2bpKSkAb3eQNm1axcZGRkYDAaWLFnCoUOHNI1HRUTLwlORJImmpiby8/OVVkdFRQUhISFKPzopKalPxvOdnZ1cvnyZlpYWgoODlVGtjo6OPnlYqHFvstTY3WKXnluzm5qanDykvb29uXbtGmPGjFFlnVN/GYjAY/+5qxzKs+Azwgu7/bb8f9uiCQTq2p2ELP3pvxcXF/Pzn/+cJUuWkJGRoVlVLFNfX09KSgpHjhxhzJgxPP3006xcuZJnn31W07hUQiTkBwm5z5qbm6scGsqGOnKCTkhIuOMN59ifdeUK19PDomfSkj0s1Ji4sFqtlJaWEhgYSFRUlGZLPWVr0ubmZnx8fPD29naadBjMUcO6ujrKy8sJCwtjypQpfVdzShKH86v5rKwOP/0I1s+PIG7K/x9C9hR3ONqxyv96Kiu7urrYsWMH2dnZvP322yQmJqp2vwPh2LFjZGdn85e//AWAAwcOkJuby+7duzWOTBVEQn7QsdlsXLp0SZmNLigoQJIk4uPjSUpKQpIkLly4wMaNG4mMjLxvtaGj8b08jjaQ3Xc2m81ph5uW+9LkQztHhZurxbNyknaXh7TsDme32+8QmrgbV0IWvV7P0aNHCQ4O5u9//ztPPPEEGRkZHuGhLGMymVi7di35+fkYDAbWrFlDUlISaWlpWoemBiIhDzXkR/RTp07x6quvUl9fz8SJExk5cqTT6N3kyZP7nFzkSssxaTlak7oSOjguew0PD+/XddWio6NDkcrGxMTcMwG68pC+1/3eD44/k6ioKIKDg/t1P2ojK+4KCwvx9/ensbERo9HI/v37tQ7NiW3btnHkyBG8vb2ZM2cOe/fu1byVohIiIQ9Vjh07hp+fH8uWLVNO7eUDw/z8fGpra4mMjFQMlebMmUNgYGCfk2XPQ7TOzk7lEE2v11NTU4Ofn5+m4g7Hls1ADu0c10g53q88Mywn6t760G1tbZSWlnqUwAPgwoULpKens3z5crZs2YJer0eSJG7duuUR5vbDBJGQhyvd3d2Ul5cr/Wiz2Ux7ezuzZs1SkvTMmTP7XAFKkoTValXaEz4+Pk62ju6wJu0NuWctH9qp7Vt9N2FHTw9pLy8vZbvLYLjU3S8dHR28/vrr/Oc//+Hdd99l9uzZWoc0nBmeCTk7O5v09HTsdjvr1q1jy5YtWofkEXR0dFBYWKj0o4uLi/Hz8yMxMVE5NJw6depdk6ksva6oqCA0NFTpz8rjWXJVKfejB2pN2ht2u52KigrFWnUwe9aOM8NNTU00NDTQ3t6OwWAgNDSUMWPGaLK1oyeFhYWkp6fz1FNP8atf/WrQq/WGhgbWrVtHcXExOp2Offv28fDDDw9qDB7G8EvIdrudmJgYPv30U0JDQ0lOTubw4cPMmDFD69A8DkmSqK+vJz8/X0nSVVVVhIaGKgnaaDQyduxYSktLsVqtSnviXj29u1mTygn6XtakvSFvEXH8UNCC7u5urly5wjfffKMIb3r6Zff3kHQgdHR08Pvf/54zZ87wzjvvEBcX5/ZruiI1NZVHHnmEdevW0dnZSWtr66AtXvBQhl9CPnfuHJmZmZw6dQqArKwsALZu3aplWA8M8poguR9tMpmoqqrC19eXNWvWsGDBAmbPnt1ni0xXog6bzcaoUaOcpOC9tRxkm0647X2r5UFPQ0MDZWVlBAcH39W72dFD2vGQ1HEcTe3JC7PZzPPPP09KSgqbN2/WrIfd1NREfHw8lZWVHuFY5yEMP+l0dXU1YWFhytehoaGYTCYNI3qw8PLyIjIyksjISFJSUliwYAGbNm1iwYIFmM1mDhw4QFFRESNGjGDOnDkkJiaSnJxMdHR0r8lUp9NhMBgwGAzK1IE8MdJTLuxYVcoCCovFQnV1taY2nXB7vE8W38yaNatXgYePjw/jxo1zOjRzXDxrsVhUM5Fqb28nKyuLc+fO8de//pWZM2f26/7UQt7Y/uMf/5gLFy5gNBrZtWuX25efDgWGVEJ2Ve2LT+j+4ePjw+nTp5VdaXPnzmXDhg3Kwd6XX35Jbm4ur7zyCuXl5UyYMMFp9O5eBv86nQ5/f3/8/f2ZMuW2L4O8KksWdlitVmXSISIiAn9/fyRJ0uR36ijwiI2N7VcMvr6+TJw4UfFcdjSRqquro7KyslcPaVecP3+eTZs28YMf/ICcnBzN1IiO2Gw2zGYzb731FnPnziU9PZ3t27fz8ssvax2ax6P9b09FQkNDFQctuF1ZDTdzEjVxtbhSp9MRGBjIokWLWLRoEfD/cmnZ4P/dd9/l5s2bREdHYzQaMRqNJCYm3lMhN2LECMaMGYO/vz8dHR10dHQwffp0RX13/fp15QBNLWvSe+G4SSQxMVHVVokrEynHJ4fa2lplm4o8yWIwGJQPptdeew2TycTBgweZPn26anENFLm/L3sor1y5ku3bt2sc1YPBkOoh22w2YmJiOH36NFOmTCE5OZm//e1vqj7CXbt2jeeee47r16/j5eXF+vXrSU9PV+31hwp2u53S0lLFq8NsNmO32+8w+O9Z0cmHdneTGsv9aPnA0NUo2v1Yk94LR4FHXzaJuAN5kqWpqYm8vDxeeuklZUP5unXrWLhwIVOnTtUsPlc88sgj7N27l9jYWDIzM2lpaeH111/XOiwtGX6HegAnT57k+eefx263s3btWjIyMlR9/draWmpra0lMTMRqtWI0Gvnoo4/EJMd90NraitlsVrw6Ll26RGBgIEajkaioKE6cOMHmzZtJTk7uUyXqOIomS8HlSr6nNen9IAs8fH19NRW99KStrY1XXnkFs9nMG2+8QXNzM/n5+UyePNnjDHgKCwuVCYuoqCj2799PUFCQ1mFpyfBMyIPNk08+yc9+9jMWL16sdSgPHPIeuczMTI4fP058fDwWi4WIiAin0bvRo0f3uWfb07+ipaUFvV5/x6qsnvF4osADIDc3l82bN/Pss8+Snp6uughG4HaG35TFYFNVVUVBQYHYN9ZPdDodAQEBREREUFlZicFgoLu7m4qKCkwmE59++ilZWVm0trY6GLSCPgAABzZJREFUGfzHxcXds4L29vYmKCjIqSpzlEbLUw6yNFqv11NdXU1QUBDJyckek/BaW1t5+eWXKSws5IMPPiAmJmbQY7Db7SQlJTFlyhROnDgx6NcfTogKuZ80Nzfz6KOPkpGRwYoVK7QOZ0jT2dnJhQsXlNno4uJifH19nQz+H3rooT73jeUDtMuXL9PY2Iher1c+JNS2Ju0PZ8+e5Ze//CWpqamkpaVp9iGxY8cOzp8/T1NTk0jI/Ue0LNxFV1cXy5cv5/HHH+cXv/iF264jKhPXSJJEY2Ojk8F/ZWUlkydPVmajk5KSGD9+fK+tjoaGBkpLS5k0aZIi8Oju7naSglut1kFX3bW0tPC73/2O4uJi9uzZQ3R0tNuudS8sFgupqalkZGSwY8cO8XfYf0RCdgeSJJGamsrYsWPZuXOnW68lKpP7R5Ikrl69qiTo/Px86uvr7zD4NxgMNDQ0cOXKFSRJYvr06S7H+xyRV2X1xZq0v/fwxRdf8Otf/5q1a9fy05/+VPPWycqVK9m6dStWq5U33nhD/B32H9FDdgdffPEF77//PnFxcSQkJADw2muvKUtO1cJisfDJJ58olYmgd3Q6HREREURERLBq1SrgdiK9ePEiJpOJI0eOsGXLFsVbIzU1lSeffPK+pjlcLWGVpeANDQ1cvXrVyZpUTtR9SaYtLS1kZmZSWlrKhx9+yEMPPdT3H4LKnDhxgokTJ2I0GsnJydE6nGGBqJA9FFGZqE9aWhrV1dWsXr2aiooK8vPzKSsrY+zYsU4qw5CQkD63JGSrTrmKbmpqchJ03M2aVJIkzpw5w5YtW/jJT37Chg0bNHeKk9m6dSvvv/8+3t7eygfQihUrOHjwoNahPYiIlsWDyokTJzh58iS7d+8mJydHJGSVqKmpuUO5KY/eORr8X79+naioKCeD/4CAgL7vw/t2VZYra1Kz2UxsbCwHDhygoqKC9957j8jISDVvV1XE3+GAEQn5QUVUJtrS3d3N//73PyeD/87OzjsM/vsjGOnq6qKxsZEXX3yRs2fP0tbWRnx8PPPmzeM3v/mNx3qviIQ8YERCHgq4+40gjMTvj/b2dieD/4sXLzJq1Cgng/+IiIh7thusVisvvvgiVVVV7Nmzh4iICCwWC0VFRSxbtmyQ7kagAeJQT3Bv0tPT+f73v8+HH36oGIkL7mTkyJHMmzePefPmASg76WSD/6NHj/LVV18RFhbmpDIMCgpCp9MhSRI5OTm88MILbNy4kXfeeUdJ3mFhYU62sYLhi6iQhzHCSFxduru7qaqqUlod58+fx2q1EhMTw9dff43BYGDPnj2Eh4cPalzCEMsjEC0LQe8UFhayfv16ZsyYIYzE3URXVxdFRUX885//5Le//a0mExTCEMsjuK+E7BnzNQJNkI3EN2zYQEFBAaNGjRK+tSrj4+OD0WgkMzNTs3G2kJAQEhMTAQgICGD69OlUV1drEougd0RCHsa4MhI3m81uudabb77JzJkzmTVrFs888wzt7e1uuY6gd4QhlmcjEvIwZtKkSYSFhVFWVgbA6dOn3fIYW11dzR//+EfOnz9PcXExdrudDz74QPXrCHqnubmZlJQUdu7cSWBgoNbhCFwgpiyGOW+99RarV692MhJ3Bzabjba2Nnx8fGhtbRWrtQaZrq4uUlJSWL16tXAn9GDEoZ5gUNi1axcZGRkYDAaWLFnCoUOHtA5p2DCYhliCuyIO9QSeQX19PcePH+fKlSvU1NTQ0tIiVIeDiGyI9dlnn5GQkEBCQgInT57UOiyBC0TLQuB2/v3vfxMZGcmECRMAWLFiBWfPnvW4PXBDlfnz59PHJ2GBRogKWeB2wsPDyc3NpbW1FUmSOH36tGpr69euXcvEiROZNWuW8r1bt26xePFioqOjWbx4MfX19apcyxPIzs4mNjaWadOmiRHFIYhIyAK3M3fuXFauXEliYiJxcXF0d3ezfv16VV57zZo1ZGdnO31v+/btPPbYY5SXl/PYY48NmcRlt9vZuHEj//rXvygpKeHw4cOUlJRoHZZARcShnuCBp6qqiuXLl1NcXAxAbGwsOTk5hISEUFtby8KFC5XRvgeZc+fOkZmZyalTpwDIysoCbrsDCjwecagnGJ7cuHGDkJAQ4LZK7euvv9Y4InWorq52MiEKDQ0VirshhkjIAsEDgqunWWEKNbToa8tCIPA4dDrdVOCEJEmzvv26DFgoSVKtTqcLAXIkSYrVMERV0Ol0DwOZkiQ9/u3XWwEkScrSNDCBaogKWTAU+RhI/fb/qcBxDWNRk3wgWqfTRep0Oj3wQ27fq2CIICpkwQONTqc7DCwExgM3gG3AR8BRIBy4CjwtSdItrWJUE51OtxTYCYwA9kmS9KrGIQlURCRkgUAg8BBEy0IgEAg8BJGQBQKBwEMQCVkgEAg8BJGQBQKBwEMQCVkgEAg8BJGQBQKBwEMQCVkgEAg8hP8DUe2izYmL4cQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# simple view of the system\n",
"snap = engine.current_snapshot\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"ax.scatter(*zip(*snap.xyz))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "0922f2fdacbc4230b8c43e756844539e",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"NGLWidget()"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# better 3D view\n",
"traj = paths.Trajectory([engine.current_snapshot])\n",
"mdt = traj.to_mdtraj()\n",
"view = nv.show_mdtraj(mdt.center_coordinates())\n",
"view.add_spacefill(\"all\", color=\"blue\")\n",
"view"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Set up collective variables"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# create a CV for potential energy (from a given lammps engine)\n",
"potential_energy = ops_lammps.LAMMPSComputeCV(\n",
" name=\"ops_cv_pe\", \n",
" groupid_style_args=\"all pe\", \n",
" extract_style=0,\n",
" extract_type=0,\n",
" engine=engine,\n",
" cv_time_reversible=True\n",
")\n",
"\n",
"potential_energy.enable_diskcache();"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.396679316834829"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"potential_energy(engine.current_snapshot)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define states"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"gas = paths.CVDefinedVolume(potential_energy, \n",
" lambda_min=-5.0,\n",
" lambda_max=float(\"inf\")).named('gas')\n",
"\n",
"droplet = paths.CVDefinedVolume(potential_energy,\n",
" lambda_min=float(\"-inf\"),\n",
" lambda_max=-30.0).named('droplet')"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gas(engine.current_snapshot)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"droplet(engine.current_snapshot)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define transition network and move scheme"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"network = paths.TPSNetwork(gas, droplet)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"scheme = paths.OneWayShootingMoveScheme(network, engine=engine)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create initial conditions"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"visit_all_states_ens = paths.join_ensembles([paths.AllOutXEnsemble(state) \n",
" for state in [gas, droplet]])\n",
"visit_all = engine.generate(engine.current_snapshot,\n",
" running=[visit_all_states_ens.can_append])"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1138fae50>]"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnXd4XMXVxt+zVdWyLcsFy71hcLcA2xhjwHSCIQkECKGGTkIPEEhCEiCQBEgIIeCEhA9Cr3GwIdjGdNx7b7jIVbYkS1bbNt8fe+fu3Lt3i7S72tXu+T2PHu3ts3fvfefMmTNnSAgBhmEYJvuxpbsADMMwTPvAgs8wDJMjsOAzDMPkCCz4DMMwOQILPsMwTI7Ags8wDJMjsOAzDMPkCCz4DMMwOQILPsMwTI7gSHcBVLp16yb69++f7mIwDMN0KJYuXXpQCFEWa7+MEvz+/ftjyZIl6S4GwzBMh4KIdsSzH7t0GIZhcgQWfIZhmByBBZ9hGCZHYMFnGIbJEVjwGYZhcoSUCz4RnUVEG4loCxHdl+rrMQzDMNakVPCJyA7grwDOBnAMgEuJ6JhUXpNhGIaxJtUW/vEAtgghtgkhPABeBzA9xddMO7NX78WBuuZ0F4NhGMZAqgW/N4BdynKltk6HiK4noiVEtKSqqirFxUk9C7Ydws2vLMPxj85Ds9ef7uIwTERe/OpbfL3lYLqLwbQjqRZ8slhnmDVdCDFDCFEhhKgoK4s5MjijaWjx4ZIZC/Tl9Xvr0lia1rFpfz1afFxB5Qob9tXhof+uw2X/WJjuooTxs7dX4qGZa9NdjKwk1YJfCaCPslwOYE+Kr5k21u4xCvye2sx16zR7/fj3gh0IBAS+2FyFM576HC9/E9fo7Kzg1Cc+xfS/fpXuYqScbVVHsP1gQ9j6A3UtaShNfLy5pBIvfr093cXISlKdS2cxgCFENADAbgCXALgsxddsd3bXNuHaFxdjXL8uhvV7apvSVKLY/HX+Fvzlky3olO9EQ4sPQHiFlc1sqwoXwWzjSIsPpz7xGQaVFWLeXVMN25oUd2OLzw+3w97OpWPSQUotfCGED8CtAP4HYD2AN4UQGdVWm/S7ebjh5cQStr21ZBc27KvHqwt3GtYfavAkdN5UUq2VrbbRA48vAAAp63MQQmDBtkMQQsTeWaGyphEnPvYJXlm4A4/MWodT//hp0srTGnbXNuHut1am1eV1+T8W4rK/L4i9o8LqysMAgK0WlZv6Wze0tP/3qm7wxLyfrf2dmNikPA5fCDFbCDFUCDFICPFIKq91oL71LpQ9h5vxv7X7sfdw263x0kJX2LriPAcON2Wu4Dtswe4Vr1/ogt+UZMFfvL0a26qO4N1lu3HJjAV4b/nuVh2/cFs1dtc24dn5W/H3L77FtoMNaPIkXsbqVlbE972zCm8vrcSCbdUJX7utfLnlIL7eeiiufW/691Kc+sSnWLP7sL7uF++vMeyj3scjzb7kFLIVjPvtHFz9r8Vh61WR353BLeSOStaMtF2+swbHPzIP/1kRv6ioD9fE333SZtE/YmEh9eyUh5oGb5vOFwl/QCRF8ADAaQ/+9B5fAB5/UPA/3VgFr/Y5GVz03Dc49YnPsKumEQDw2qKdMY4wctdbKwEEXROSSY/NS9jy27i/Pu59g30cwUgWp80qBiG1fLG5ChMenacvn/bEp5i/4QBmr96LFp/fsvL6cM0+bKtqwJo9IcF/eYGxf0at3Otb2vacfrapCre+ugwzV7atW86qAlPL9fqiXeh/3yx8vqnjR+9lClkj+AePBB/8215fEfcxvzJFArQ1quajNXvD1nUpcKGmMdhs/fH/Lca6JPjH73lrJYb/8qNWHVNV34JPNx4IWy8l80iLFy2+kMh/a9HBlygrd9UCABZvr8GdbwZ/ny83H0RNFEv7cFNIhPwBgXF9OwMAahq92HzgSKvLsOjbav16CzVLvXOBM+ZxtUo5mtPg0nlo5lrsU8Z0bK1qwNUvLsbNryzDlf9chHG/nRPx2P2msSDqPVWF9ZLnF+BwY+tF/8p/LsIHq/bip68tx98/3wYAWLP7MJbtrIl6XLQKe+/hUJlna+/VSzkUTJBqskLwV+yqxXUvhfzw8VqA5gepvo1N25WVh8PWdcp3oL7Zh79+sgVz1x/AOU9/0aZzq7yruUQCgfgt3PvfXYWr/rXY8BLOXLkHL3z5LQCgxRvQXTpAajqa528MWWjvLtuN/vfNwuUvLMQlMxZg4TZrN8W2qqCo9yrJQ4PHZ+hUNAtZLAIBgYuf/wbf+9vXAILiDwBFbge2VR2BP8r9rG4IRbOoz8eCbYfwf19vT2qLyAq1dWMmlouprsl47PwNoYpfbSnWt/jw1tJdSIRHZq8HAJz3ly/x3We/xpLtkcum3m/5LB880oL9dc2GMsroojW7DxsqK6btZIXglxW7DcvRXpJo1CXpoXrzhonIdzmwbm8dnv5ki74+WaNvVYs8FvLdkhYYAPz0teX6Z4/fKPiqhZVqNu6vxw9mLMDT8zaHbatpDFrjo8pLIARQ1+xFr5I8ANF98EdafPjJa8tRVR8S6kbNmt2mCUijJ/h8VNY04dQnPsPzn2+NeL5DR0LXkoIvhMAlMxbgVzPX4i/K77u7tklvzSSLePzrkQycA/XG0Mvb3wi1fq1cg0IIPDRzLZbuiG6hAwirJM2tpe8/903EY33KsbL1UvHwXJzw6Dw8PGs93A4bxvbtrD+7++qacc6fEzeYmCwR/F6d8gzL9727uk3nqWuDhe+xEN/jB3RFvjP81ppfwLYSKZrm759vC4u57tEpuki+9M0OLNp+CHlaeZNV6UWzms08OWcTpvx+vt5Jd6C+Gde8GGyx9SrJBwDUNoYEP5rb7oOVe/DflXvwx/9t1NdJgZeYO6c37I3s05cVDxASfPV3VO/3iY99gul//Qofrg538UXicKMXWw7U478R/OANcfTZrNhVa6jgJAePWD9vzV4//qG18AxlafLixa+36y2haJjv4dDuxQCA3p3z9XWRDJxn54cqye2Hwl2ILb4AuhUZjTjuwE0OWSH4NlNn2qxV8b9wKou3V+OjNftadUykTtQCV/gQh31Jsp6tfMnVDR48Mns9rvjnIsN6KXaqcJlZs7sOnfKCFlqyInVaG+K5s7oR92idtEu3hyxMKSA1jZ6wlpwV+a6g66dBEXnzb9RoWjaLiyqUaoe8/E5qa6mLRT9APBYyEAybHP2bjzHtyc/xk9eWG9xDQgj844ttYcfYCFjz6zMN6y589mtMe/Izy2uoAgwE3WEnKJ3AKlaVRiQatVZ0z055OO3o7jjU0ILqBo9BmM1umC0H6tH/vlmGVu+VpucVCN7TbkXhkW9fb+U0EImSFYLfGvbXNaP/fbPw9tLKsG2fbqzCjf9e2qrzNXqN1uPiB6YBAPKc4QNZWhsOGIlmbwBfbK7Cs5+GXhzpa5ZN5OoGD2av3qtb7NUxIobynHa4HbakRQG1JaZf3p88V+je9dSs+kZPcHDQRePLdUvfzBMfb9Stf/V7mOPMzWXrlB+qnF9esAMVD8/VY9jV1oGsaFXruKQgXJhUq9zrD6D/fbPwysLwjscVu4wVw6lPfIparWKes24/Hp61PuyYkeWdUeR24J4zhxnWH27yYld1Y9j+x/Xvgt9MPxYnDw2mLTnh0XmW/vCtVUda1WH/ieZrv/fsYSgtcmFrVUNYB3Jds/E6/10Zboh5/eEtQbfDrlfCR/cs1tdf9vfMSwPR0cg6wb9qUn+4HJG/1lYtwuPdZSHBnzK0DCX5sSM2rNi4z+gOkFZogSZaEweW6pVAIgN3fIr1t2JXDX70wiL8/qOQ20JGKUlXynUvLcHNryzTO0xrGj1RO7N3VjeixRfA859vC3OBtIVmC1fXUz8YjXNG9tSX1/3mTPzl0rH68oZ99fD5A/rI314leSjKC4mxy2FDvssesTJ59tOQL95g4SuVsj8gwiz8Zm+wrDUNHj1e/a/ztxj2ddoJLd5AWIf5DM3/r97bBqUPqVGrbB6bvSGsvB6T2O2qbsJri4Kdp5U11i6MH08eAAC45ZTBGG8a2S3j7tUGr8thwxUT++PEwaVh5/rpqYP1z68t2oXrX47f2JFuU7fDjkK3sTV7cUU5gHAL36qDmyi8D8JuI13wS4tcOG9UL33bCY/OjbuMTDhZJ/id8hzw+AKtitXu17UA/UoLou4z9Q/zLTuOrtIGjzx96VhsfPgsfX2+ZuF3LnDCrfnHW9PZKvH6g99FHbR0xxsrw/aTnYvyZd9kqoj8AYGzo3R89VT6QW5oxYsfCSnKD5wzXF933qij8Myl4wAAJw3phgKXI6wl9OScTbqv/O2bJuGokpBLYmvVEeQ77ahp9OIdixZaodIyaPKG7nWtEnK4ZHt1mOCv3FWLmgaPYQzHR2v34dWFO9DY4gNRUNhe/Ho73llmvG6zN4Amj9/w26oVu197Dr2B0PbZq/filleXYWtVeHipWzNWrFqDb1w/Ad8ZfZS+LIU1/Byh+yDHW3QtDHeH3T5tqOXxraG6wQNSciQWux24aWqwIjEL/qb9xu87pHsRhAhVuBKioNADgI0IU4aEkiruz+AcQB2BrBN8tyYgHgtrosXnxyZt0A0pVpCAQLFiSf76v2vDrPHthxqxLkqcfovXmI9EtjK6Frr0l7i1gn+grhmjf/0xHv9oo2XnFhCymqRVKcug+vnlqNoN++pR32zt2nnxmuP0z19sPogvNyfmL5UulW7FIZeHw0aw2Qhz75yC5y4fDwB6Z7FkyY4a3K9Zj0UuBwZ0K9S3Ld9Zq1cQclCWSnFeqJWmVvi3vLpM//yDGeHpCb7ZdghjfzsHs1cb+2/eWlqJRo8f+U67Hvl1z9urwo6vN41lUDus5e/j06z5jfvqcfMryzBr1V58YNFRK42DQw3hwlZi6i8wt2Q/XLMPry3aabi+jLqyGg1u7vtSiRVuKl0tpw3vbnhX3rtlEnp0csNlt2HtbuP7sqrSGMEkW9VPzgm2VGWLxW4j/Tm22wgXH9fH0AJvbT8bEyJrBP/+s4/GwG6Furiu2R0uzg+8twYP/Xdd2HohgGJ36IH611fbDf7G31gcA4QidNwOG84fc5Rh2+Qh3XDy0DJcXNEHLs3KammlX3vu+gNo9Pjx3GdbsavauokvrVXZ2SpfYekb7dkpD/+4skLf/1EL1wIADOlebFi+/IWFuOetlfib4iZpDbICKisKtRxIq2UHdy/W3QBmC191mRS47QZR+87oo/ROWSv6dA21BlZVHkZtowdefyDMgozEIlPs+KrKw2jw+C074KcOC1md9c0+g+iprhpd8LXvpXboWkWFSaGTo7QLXXb9mTa7HW1kFOyZK/fg/ndXG4ydQ1oHdFcLwY+GfK4+31RlcKEJIbDvcDMONXhwdM9i9CrJ1+/vr75zDAZ3L0aBy4FR5SX4aush/Py91XrfgrllJY2sv38R7BMZXFakf6++XYMt7smDuwEAJgzsqh/X2n42JkTWCP4NJw/CJ3dP1QXCKrTsi83WQ7QFgpaEihpu+c+vwkPYgGByLwB49MKRYdkGB5UV4f+uOR6j+3QGEcHtsEW18Cc//klY8rVZq0MWYKTh69KSli+lSQNQ6LYbXvZI6Q3M3x8IWriPf7ShTTmKpA9d9cFbkWe6b0sUQZTuCPmdnrx4tO4qA4K/0aYoaRKe/3yb7h66aHzI/fG9ceX4xXnxzbTZ5PHp/TEqP1dcVct21BgqZLUVpXZKBgICtTHyKzntwS/b4PFhVHkJVv7qDMy982Tce9bRBreb+dyRkFZ8906xI5xUGj0+bD/YgCv+uUhvcQHAnW+uxITfzUNVfYv+rsnKTn3OSotcWL+3Dq8u3IkP1+yFEAKNHh9+NKGfvk9RnrECk/1fNgKG9SzGZ/dMxbVan4W5n8AqHJqJTdYIviTaKFT1IVH9jkIAN00dZHAf+CP0AXywao/u7pARMb06W0eNqEQTfI8vgMqaJvz8PeP4AXXQDwD07VqAe8862rBOCqu0ssik+AUuR9zW3S/POwa3TxsStv5Pc8MHRkVDCIH3lwcrqCK3wzJ0UeKwW7sVZOckAHx576l47+ZJcNptBsEf+uCHOOOpz/UKyRyNExBCH7g0tm+og3N4r2JcO3kAtj92bszv8v6KPShyh1daXZTonHveXmUYzKaOyFVdIwN/PjsshcHwXp0wXWkdStdPk8ePIrcDDrsNfboW4Kapg8J+23hG+cpn3lxZTBoU3omrsnh7jX7+95bv1sMt1b4kaSTcMGUQuhW5dGscMN6fuiYfWnwBBEQw6urS4/sCQNh9lc9pqdbf0K+0UP/OhaZWluqmY+In6wTfKjpEEtkiEhjRuwTXnTRQX1PX5MX8DQfCkk7d+upyXP5CMDxMvryd82MLqttpjyj4P/xH0K8sXT+SBlO0TJdCF26aOsiwrr7Zh+teWoKn5m4CEIrmkc3/qyb1j1vwr5k8ALdPGxp2DXO5ovH11oN4a2ml3iIpznPgq/tOxZf3nhL3OQCjRde7c74u2OVd8sP2XbCtGg/NXBs2wppAemhgaZFLd8OYrUUVq5ZOsUUrpXOBE+/ePElfloLYpcCJLQeO6JFOZku0sqYJRCGxmzCwK359/rH6dimyjR6/ZctCJZLgFxg6r2XLj3D5hL76+n9fe4LhmIGKsRMsZ6PBNXTiY59gh6kfSX63keUlWPLg6ShVxjPsUcac7Ktr1luiBS673qfUyXRfB3Uvwi/POwbP/HAszBS4jfdizrr9YfswsUn1BCjtjupv9PoDulsAiNwMnDw4KATq4KQtB47g5W92GBJXmZFRCPEk4Qpa+NY+/MXaQCNzB+aRZh9cDptebpmtce6dU7B2Tx1ue30FLpnxjcFH7fEH8J8Vu3G4yYsLxhyF740vt4xY+v74cvzwhL5Ysr0Gw3oa/ffmDj4rCzcS5ljpIrcDBS6HpR8cCFmCI3p3Qr7Trt8Lc6UjmTioFKWFLsNcAw+8t9oyDxJRyNoudjv0ZyGakBa5HWho8RmG/xfnOdGrJM+QdsJpt2Fsn85hx5cWuVHT6MVV/1qMN2+YGCbKXn8ATptNr5wGdCs09FPI/Zu8fuRHuGehfYNlPG9UL3ygDDYsK3Zjx6Ggu3FAaUjIH75gJP69IOjSk66e/946GVurjmD6mKPwxMebcKTFhxe/3g6vT+CpOcaW3cl/+NSwHG2shRoxVd3gwXNa+GqBy65Xqubn6tijOunjBcwUaffimF6d4A8IbNxfH/Z+M7HJurulip85/7pqsQgtX+R1Jw3AuVqcr9oJtnJXLaqjjE4FQoIfTwx/JJeOWgmZOyQbWvzoofheZafX4O7F6KN1apk7JL1+oQ8+kkJidgUAwB8vGo2xfbvguikDMcX0kplbBPFmBbYa0BPLSu1Zkoc5d0zBuzediLdunITvjSvH3WcMtRy4BgS/i3QJSNTRpCN7l+ifhQhFu5QUOHV3n3ruYk10pDuvyO3A5kfOxjOXhazM4jwH5t89Naxit7qv/bXwXpmgzdyqbPD4DG6syYO7GVpQcv9Gjw8FEe6B5PzRR2F4r05hbr7iPAdunjoIfbrm46lLxhi2vXj1cfjgJ5P15ZHlJbhgbG8QEe4+cxgeOv9Y2G2EmkYP5q6PbkVH6wx/9MKR+udPNhzA858FRw3nuxy47IS+yHfaceG43vqz9fAFI8JGPKsUaL+Txx/A5ROD/QA3v7IMB+qb8ae5m1qVUDCXyToL//vje+O5z4LWRKPSxDf7TqWvtLPia7z6xP7Ic9qw41Aj3lyyC10KnBHjfldXHsbaPXVw2CimqAHB6IsWixdETWVg7pD0+AMoLXTrHYJqhRDPNVXL8e0bJ6JbkRsLth3CUJNFb0YK/qCyQhyoa7HM92/F6F9/bFh+9ccnWIqimSE9QuV54uLRMfc3R3uoFdQkbYDR6t2H0ez1Y8XOWjhshEFlRQhoLR2HUoO9fsMEzFyxB8t31uLbgw0ocjtARIYQz+K84HiBwWVFhk5lK8q7GMdz+EwW/pFmH5x2G371nWHoWujCQC0yJfTdfAhoA76iRSQBQUv+w9tOCltf4HLgZ2cdjZ+ZKgIAmDqse9RzAsGO4/fjmFcimoXfJYIb0WW3YWiPYqz/bXDMyri+XbBkR01Mt2GR5tJp8fnRTTv3nHX7caC+BSt31eKckb2Q57Cjd5d8S7ccEyTrBH+wEl6oNsuPmPzhCzULTH048px2XH3iAMz4fCsaPf4wYVH5zjNfAgha9/GImttp7dJRUwBIy3PdnjrM13LYnzCwK1ZoGRjVSBNzJ5YVanO3on8wrK2/yVdrhax4ivOcaPT4caQNE2T07pyPSUonXjIxjyVQJ9Jo8Qbw6nUnYORDH+sTYU8aVIo8p13Pvqi25I49qgTHHlWi53SRUUWqu0G6o565bBx+/t5q/EQZoapS0a+LofVg9oMDwVTETjvh6hMHmA8HAPzx403YWhWc2SueSt2K/Bgtg1g47TbDYLVImFMnxMMpRxtbk9K15HREf4fkb9DiDaCbklNJZifdcagR1720BD+ePAAPxhmBlYtkneCrBBTfdST/vVVWx+7FsaNuJPG+XJFcOgYL32XH4SavIXf+4LIibH/sXBxu9BoG3sSy/oDwgTnxMqxnMZx2wp2nD8VvPljXpnTTI3p3atO14yHavAWNHp/BOgeAMZqvXT4PVgOOpLjKDt1iU0oHIOh++udVx4UdK7lt2hB9diwAmPz4fD2GXnKk2QeHLfrvIl2Rbf392nqcfnycfvF4wkLNmMOX5U8RyxdfqFv4ActBZHI+jLeWVrLgRyHrfPgqqshHEnyrwVCtGaQSyzKRuBzWUTpq3pp8pz1suL20NM2jLK2sP7MPtDXRNSqdC1zY/Mg5mDK0DEVuR5hL56w/fY43FofH86uuktIo/thEcVuknpY0WbjNpJicOzLYV2OOSAFCFqT06asWvitC6KiZfKfdkPcFCB9d3dDiswxFVQeNSRytcE2oPvNEBT+W+N53drirqK3I1lasZ1W2aFt8foOFb4YnSolOVgu+aglGFHyL9V0sMiDeeLJ11Ei82SXdDpteuRxu8uKKfy7C/rpmgx80z2kPC32LNHDJPGAJQFhK2bYKvkpxngNHlKa7xxfAhn31uPed8DkHuhW5dYstlS/eQ985Fr+0sOIq+nWxHEcgBfAHx/XBht+epXd4q8gKVFas5qRt8ZDntGNUeXjkjkqDx2/5u3zxs1PD1tljtARULjuhLx48NzgYLNHf3WzEPPbdkYblsiI3Lq4ox7M/HBf1PP+62tgaUvMASaTgO2PcY9ny8vgCeqUciV/+Zw3+YjGpDpNCwSeih4hoNxGt0P7OSdW1IqH64D1+a2G2Enw1GuPn5xyNP18yJqI/Nd5JU9xKeOW7yyrx+aYqPDt/i6GMASHCOnYjxYxbuSU6maKFYr1E8VDochhcOrVRIpdafH5MG94DAHDJcX0SvnYkuhS6cM3kcB/42zdNwiBTJygQSkhGRBGjf2RfjpxMR+0jiTf0z+y+iUSkwWZh+7Wy81GW0xnn+WOdBwAePHd42DwEbqcNv//+aJwzspf5UAOnDOtuiKBSM6NKZHeKOU2EGfkeBETwd9z+2Ll6+gUzL32zA0/M2RT1fLlKqi38p4QQY7S/2Sm+ls77t5wIwCj4VhEygHXKYjXC4JoTB2D6mN74Zqv13KvxDvF2O+xo9PhR0+DRxcUvhKGF4PMLQ0czED0G3pzy1jyQJRkWflGewzCCVQ1VnfDoPJz0+0/05RZfAP1KC7D9sXNx0hDreOp0EI9gy3sn4+DVzvx4LXyzfzoSkXz455rcQdGSm1khn51EO23lc9OrJA8/Pmkg+pUaXWDxBAxIYjVSZMBDIEZ220ILg0udx0ASb6Wbq2Tl3RnTpzOGdC8y+MdbtBf5/rOPxtkjQjnZrSqCQpcdR/csxpShZXBoD/+tESIz4sXttGFfXTPG/naOnkzKHzB22voCgbBO5GijQof2MIZXmjsrk2HhF7kdhqgYNZHcvrpmQw4Zjy+QsP84Wbx940T9czxlunbyQJx5bA9cosT4f3jbSeiU59BbLbGQA+deuLICN00dhGuUSJwXlAR2kSxwp0ngW2vhN2nPe6wBW7GQFaRsDQ3uXoSlD07D98YFo8TG9+8S8Vgz8v16S/k9VORXjJXOvMDiPTCPfzllWFmbUpDnEqmO0rmViK4AsATAXUKI+OZ+SwKbDxzB5gNHUFXfgrJit26Jnzi4G7YcCHWMWrkGiAgf3T7FIL6xco/EQrU83lkWjMIQioU/sKzQ0sK3ikjQywmjIHQ3Nb3NFn9bKHTbUdfswx1vrEBVfYshBFLF5w/AFxBxW7nJYN5dJ+O0J6yn9qvo31UfpRxPS6ekwInnf1RhWDe8VyeseujMCEeEI7/7acN74LThPfDSN9v1bacN74G+XQuws7oxYovDYVrf2nhy2aK1soZbg6yQ1Ge2tMiNhy8YgdunDdGnw4wHKcDmqRYl0pUTiKHTVoPQZIbbfqUF6FdaiImDSvUJfxhrEjLHiGguEa2x+JsO4G8ABgEYA2AvgCcinON6IlpCREuqqpL/Y8k4XSn4eU6bHuVxXP8uGKGMzDSjvnDmWPvWNpv3KHN9yskt1BmVivOc8AYE/NqT/8xlYzHjR+Mj+pwB4KenDcalx/fFkgen4VffOQZnKS0XINyn3xZk9Mp7y3fjyy2Rc+TLePP2tPCt/PUqeVpZktHSicT3lbER5ugh8xzGMoFZJB++2fJvq+BbWcOtQVZIbtOzl++yW3Z4R0O6TK3yEQFAZ+0ZjdWvYeXesmvHXDS+HC9dc3yrwqlzlYSeDCHEtHj2I6K/A/ggwjlmAJgBABUVFUkfHy07YKXgu+x2uOzBBzkRa7TF58f1UwZiVHnkCkNFzpGq4hdCd+l0ynOgusGjW/jThveIKvZAMHzyd1oExdUnDsB60wQtrcmBE4lYZQCCLRV1boB08LaFyyDPGWydJKMvIxKPf2+UPj+y+btfMLa3YdpFOT9vRAvflpiFryYoSwRZaSfjt7zttKF4/KMNEf3+vzr/WAzpUWyY1SoSQ3sUYfoWPMoaAAAcCklEQVSY3vqydAPJZ9QcupxJLsZMIWUuHSLqJYSQGZ0uBLAmVdeKhnTLqBao7CCKlhgtFgFhzIkei9tPH4qfmWZL8geE3uwvcNlxoE7A7w8f/h8v5hc00RcfiK8l0+IL6E339nTpAMC7N0+Cw0aW4ZBycFoqKyG7jTB1WBk+3VgV1goc2qMYC39+mu7HlpOvR6qIrztpIBZvr8bG/fUQArDHMYJb5Y7Th6Ku2WsZ/tgazD78RLhp6qCIifCAoB8+2naVj+842bAs3f7ydzb79OubvSkdD9IRSWX193siWk1EqwCcAuCOFF4rDBkC5jFNAZjvsutpFcw5wmPx+T2nGJJqtYaLK/rgZtOD/Z8Ve/RcPQ67DV7NDw603roDjO4UIhjy+7eVfFfsR6Shxadbl+1tUY3r2yVi7Lv0Nac6t8pzl4/HogdOs9zWo1Me+moJ1Xpoz1ukiJS+pQX46PYpMV0/kehZkoe/XT4+4ZadlQ8/E5H3sSCi4Ld+hHi2kzILXwjxo1SdOx6kr1GG2cn0sZ3yHLhj2hA88P4avHBVRcTjrehbWmCIqmktPzvraFQ3ePD64l1h2/KddjR7/fAHBOw2iis/jxnpuigtdGHpL05vczlVrAZ4mfndhxsiujXSyeDuRVi9+7DlhODJJM9pj8sali4dq9HAKlLo05UELJkWfiqR8Q3yGe3XtQBdCpyo0fIAtSXXT7aTOW9nkpHi5/EJHG706gnIiAhnHNsTix+Y1ib3Qzy576NhZX09c9lYFLjs2HO4GV9srmrziy6ta3V2p0SRL1WkTjcAutirZcgEHjx3OL43rhzTjokvrDLVyA7K5hijs52aL7+1Lp1kId+dTKq8rZANJdmh67Db8LvvjtK3s4UfTmb/ogng0oaHe/0BPDp7fdLOa5V2oVXHW4RZnjfqKD3AcmXl4TbP19m5wIV3bpqIP5tyoCeCTP1gTvsLWIeqZpJIlBa58cTFo5PSeZ0MpGtn0uDoIb6ZY+Fnzm9pjZYMT6kY1VtmzqrKZLHgy4fW6w/gjSXhLpS2kqgFe5Rp/tt5dwU7orxJmsBhfL+uUQdrtZbB3YOhj5cdH54qQZ0SUtLenbYdifIuBfjy3lPwk1PD8/2oyBw6rfXhJwuZSyfTf8tQumvr7Z42ZPPMdrJe8FVr+RcZkDa1V4lxAEo3bcJmf4Y+nKP7dMaSB6fh8gnBWYacdsLdZwzFtOHdLSu/THLpZCLlXQpiWu6y0zRWfplU0VEsfNmaVO+nes94FqxwMqOtmwLkQ6tOc2iVFrctvH3jxLimNbTC7BKSkzN7Yw01TCMy7fK/rjoO/bsV6tE/i7dXh+2bSS6djooMyY2VNz9VyBHcmW7h/2b6CPTpWoDJykQ7ah0ZKz9PLpK1gi87nmQIJhA9j3prkLNHtQVzbLysmFI5OChZnHK0cXo8qwFELPiJI4U+TXqvC2WmW/hlxe6wsTCq4FtNbpTrZPYvmgBWroVMCDNTZ6qaMDBUcSRzUon2wqqSynSrsCMgfffpMlDl2JVMeF8SgQ38cLJW8K0yEmbC1Maqha9Ol9e5wIWfJpiRs72xqlTT1dGYTTiUgIN0EErR0PEcAGpCQT8rfhhZK/hWHWOZ0MJTBzKZy3jnGcPauzgJobpvZIbGaPH6THzINMm+NHXky1HpiWbdTAfHDeiqj5VhH344WSv4ViNV+3SxTtHanqhZ/9I1sCZZqBb+LacOxrZHzwnLyc+0HtlK8qWpI19m3czvgIJf5HZgjpZzJxMMvEwjJ8yxz+6Zih6d8jLOJ5mugTXJQvXh5znsrZ6hibFGJqxLl2BJoU9Geu10IB9DDssMJycE3zxFW6Zg1Qp55rKxqG3sGCMEVQs/WRFQTHC8SOcCF05PU0qI3313JE4Y0BVj+0SfkD1TkYYUu3TCyQnB70icNyqx1LbtidoBzdE5yaO0yI2Hzj82bdfvVuTGjy1GUXcUpCHFYZnhZLVZ9tI1x+OPF41OdzHC+PHkAW3Kd59pEJE+CCvTY7aZ3EFa+Gzgh5PVFv6UobFn0UkHD553DB7MgDQPyUC+XGzhM5mCtKU4LDMcNsuYhJCRRjzClskU9InRWfDD4LeUSQgZmZNpEVBM7qILvubDr3h4Lu58Y0U6i5QxsOAzCSEjM61GNjNMOtDDMjUD/+CRFryrJFHMZVjwmYSw683nNBeEYTQ4LDMyLPhMQvDLxWQaMixz+c5aPPnxxjSXJrNgwWcS4q4zhqE4z4FhPYvTXRSGMfDZpio8/cmWdBcjo0hI8InoIiJaS0QBIqowbbufiLYQ0UYiOjOxYjKZyomDu2H1Q2eiE+fQYTIcX5qyj2YSiVr4awB8F8Dn6koiOgbAJQCOBXAWgGeJiMM4GIZJGx4W/MQEXwixXghh5SSbDuB1IUSLEOJbAFsAHJ/ItRiGYRJBnd86V0mVD783gF3KcqW2jmEYJi2w4MeRWoGI5gLoabHpASHEfyIdZrHOMoyDiK4HcD0A9O3bN1ZxGIZh2kQLC35swRdCTGvDeSsB9FGWywHsiXD+GQBmAEBFRQXH9jEMkxLYh586l85MAJcQkZuIBgAYAmBRiq7FMAwTE3bpJB6WeSERVQKYCGAWEf0PAIQQawG8CWAdgI8A3CKE8CdaWIZhmLYiBf/peZtx15sr01ya9JBQemQhxHsA3ouw7REAjyRyfoZhmGQhXTpPztkEAHji4sybKyPV8EhbhmGyljtPH6p/ZpcOCz7DMDlCi8+fs64cCQs+wzBZixof3uIN4J1llfrygfrm9i9QmmHBZxgmJ6hv9hmWr3gh9wIHWfAZhslaSDHxH/1wvWHbtwcb2rk06YcFn2GYrIUUxa9t9ALI7TkcEgrLZBiGyUSevnQsyorcWLazJmyb3UbwBwS8fhZ8hmGYDs/5o48CAEvBd9oInvYuUIbALh2GYbIWskjj6LDnruzl7jdnGCYnOWVYGQCga6ErzSVpf1jwGYbJWi4a3wejykv05ZOGdMOvp4/A4O5F6FdakMaSpQcWfIZhspayYjdm3jpZX3752hNQku/EwG6FaPIkJ5/jntom9L9vFt5fvjsp50sl3GnLMEzW848rKjCwrFBfznfZ0eRNjuBv3F8PAHh3+W5cMDazJ/ZjwWcYJuuZdkwPw3K+095mC9/jC2DYLz6EDON/+tKxAKyn+cs02KXDMEzOke9qu+DXNnmgjtm6+61gQjZbB1B8FnyGYXKOfGfyXDoy7bLNKgY0w2DBZxgm5yhw2eELCHgt5rlt8YVXBA++vxrTnvwMACKO0CUCKmsaUdOQucO6WPAZhsk58px2AECjya2zteoIhj34EWau3GNY/+8FO7HlwBEAgDfCRCpEhMmPz8fkxz9JQYmTAws+wzA5R74rKPjNJrfOxn3BiJvZq/ZGPNaqVQCEfPgNSQr3TAUs+AzD5Bx5jqDgmztuXVraBU8EUQ8ERMRt1AHidFjwGYbJOWSKZL8pRbLToQm+5rbZcqAetY0hn/zgB2ZH9OHbOoCaJhSHT0QXAXgIwHAAxwshlmjr+wNYD2CjtusCIcSNiVyLYRgmWdhkTvyAUbzNFv60Jz83bA+IyC6djmDhJzrwag2A7wJ43mLbViHEmATPzzAMk3TsZG3hSzy+AHwRhN0TodO2A+h9YoIvhFgPGGeVYRiGyXRkhmS/ycKXy15/ADXaDFlmovn3M51Uep0GENFyIvqMiE5K4XUYhmFahRwkFTBpt09bsXZPHfbUNlkeGykssyPMoBXTwieiuQB6Wmx6QAjxnwiH7QXQVwhxiIjGA3ifiI4VQtRZnP96ANcDQN++feMvOcMwTBtx2K1dOqrFL1MmmHn8ow2W6zfsC5O3jCOm4AshprX2pEKIFgAt2uelRLQVwFAASyz2nQFgBgBUVFRkfhXJMEyHR1r4ZpeOT1nerA20MrO1qsFyfWVNqEXgDwg9EiiTSIlLh4jKiMiufR4IYAiAbam4FsMwTGuRYhyIYuEnQsSO3TSTkOAT0YVEVAlgIoBZRPQ/bdMUAKuIaCWAtwHcKISoTqyoDMMwycEeh4WfCFb5eDKBRKN03gPwnsX6dwC8k8i5GYZhUkWkOPyfvrY8KedvyUYLn2EYpiMSaaRtsmjxsuAzDMNkBJE6bZNFprp0WPAZhsk5pIW/ctfhlJyfXToMwzAZguy0fWruJnx70DrMMhHYwmcYhskQ1MyWajZMyU9PG5LQ+dmHzzAMkyGog6KsBkiV5DtjnuOyEyJnBmCXDsMwTIZgVxI+nv/MV7jgr18ZtudrUyBG49ELR0bcxi4dhmGYDMFmsupX7Ko1LOe7rKVxwsCucZ2/mV06DMMwmYE9Rkp3t8Pawj97RK+4zs8WPsMwTIYQKbGZ22HDDScPhNNuLY2R6olR5SWGZfbhMwzDZAhml47EHxCwE8Fpt94eqV3wt8vHG5Z/+Z+1+PPczYkUMSWw4DMMk3NEcun4hYAjWlpjIrxz00S8dM3xhtXFeQ50K3IZ1j01d1PC5Uw2ic5pyzAM0+GwWZi6H6/dByGC1n91Q3hsPhC08Mf3C3Xc3jFtKGoaPeiU58R/fzIZK3bW4qZXlqWo1InDgs8wTM5hZeFf//JSAIDDRhjZO+iTP2dkT8xevS/ieW6bFhqg1askH71G5ie5pMmFXToMw+Qc0WajsttsGNKjGNsfO9dgzQORO207Ciz4DMPkHJE6bQFADdBp9hrDKylit23HgAWfYZicwxbFVLcrDn5zeGU8Fn7frgVtLleqYcFnGCbnKHI78PAFIyy3qVE6104egO+PL8e5o+IbcAUA8++eCgCYPLhbQmVMBSz4DMPkJJdP6IcXrz4ubL3q7inJd+KPF41GkSsY3xKPQ8duI4zv1yVZxUwqLPgMw+QsLke4BFrF4Y/o3QkA0L9bYVznddgIXn/mjbblsEyGYXIWh0VAvlUEz+UT+uGEgaUY2qM4vvPaKSNz4idk4RPRH4hoAxGtIqL3iKizsu1+ItpCRBuJ6MzEi8owDJNcHBYpFKxi9IkobrEHghWJN0Xz5SZCoi6dOQBGCCFGAdgE4H4AIKJjAFwC4FgAZwF4lohiJ5hmGIZpR6zcN1aVQFvO6w9kmYUvhPhYCOHTFhcAKNc+TwfwuhCiRQjxLYAtAI63OgfDMEy6sHLflBW7Ez6vw07w+bPPwle5BsCH2ufeAHYp2yq1dQzDMBmD1UCq3p0TT49QWdOEDfvqcdFzX6O+2Zvw+ZJFTMEnorlEtMbib7qyzwMAfABekassTmVZ3RHR9US0hIiWVFVVteU7MAzDtAmrgVQ9S/ISPu/aPXUAgMXba/Dygh36+poGD1748lsIkR7rP2aUjhBiWrTtRHQlgPMAnCZC36ISQB9lt3IAeyKcfwaAGQBQUVGReW0ghmGyFivBjzTbVVsJKJ2397y9EnPXH8D4fl0wpk/nKEelhkSjdM4CcC+A84UQjcqmmQAuISI3EQ0AMATAokSuxTAMk2p6dErcf2/Gpwh+bWPQvdPiTc8UiIn68J8BUAxgDhGtIKLnAEAIsRbAmwDWAfgIwC1CiMyc5JFhmJzF7MO3GoiVKH5F8GUOHwGg0ePDaU98iqU7apJ+zUgkGqUzWAjRRwgxRvu7Udn2iBBikBBimBDiw2jnYRiGSQeDuxehT9dQJ22syc3jpbxL6Jyq4Mv6JSAE1u2pw9aqBjw8a11SrhkPnFqBYZicxW4jPHrhSH05Wtrk1vDGDRP1z6rg62cXQIGWn6fJ037ODxZ8hmFyGtWqjzqfbSsoyXfqn6UP3+sPYOG31QCAgABcjuC1Gjy+8BOkCM6lwzBMTkOK4EfLk98a3EpfgLTwG1pCwn6kxYdpTy4EwBY+wzBMu6Ea9dGmPmwNaktBToiuht5v2Fenf25kwWcYhmkfVL99sgRfbTXMXBkcguRVcuuoaRdY8BmGYdqJVFj4Vqgi701TYjUWfIZhchzFwk+SD98KNVrH60tPUgEWfIZhcpqje4by3CcrLFMl3xlM1aDOgOVjC59hGKb9KXQ7sPqhM1Ba6MJdpw9N2XXUFAveNKVO5rBMhmFynuI8J5b+4vSUXkO18NM13y1b+AzDMCkkoMVj/m/tfn2dx8eCzzAMk3X4AwI+fwBPz9usr2vuoNkyGYZhmCj4AsLgvweAZrbwGYZhshOzz54tfIZhmCylxWTRd9QJUBiGYZgYmDtpm73s0mEYhska7j4jFNNvtvCbfWzhMwzDZA23njoEvzzvGADA7pomwzb24TMMw2QZDnswVcPlLyzU1zntxC4dhmGYbMM8ocpVk/qjb9cCtvAZhmGyDfOUiWP7dobTbgvz6QvRPrl1WPAZhmFShDm/vo0ILke47LZXMrWEBJ+I/kBEG4hoFRG9R0SdtfX9iaiJiFZof88lp7gMwzAdB+nD15dtBKc9XHbbK11yohb+HAAjhBCjAGwCcL+ybasQYoz2d2OC12EYhulwmCc6sdsITnt4zn1z6oVUkZDgCyE+FkLIqdgXAChPvEgMwzDZQV2z17BstxEKXOFZ6X0dwaVj4hoAHyrLA4hoORF9RkQnRTqIiK4noiVEtKSqqiqJxWEYhkkvdc0+w7LdRijOsxL8DHHpENFcIlpj8Tdd2ecBAD4Ar2ir9gLoK4QYC+BOAK8SUSer8wshZgghKoQQFWVlZYl/I4ZhmAzBnFLBYbNZCv53nvmyXcoTc8YrIcS0aNuJ6EoA5wE4TWixRUKIFgAt2uelRLQVwFAASxIuMcMwTAfhpqmDMHf9fmw5cASAtPCd+vZzR/bCrNV7sb+upV3Kk2iUzlkA7gVwvhCiUVlfRkR27fNAAEMAbEvkWgzDMB2NknwnHr5ghL5stxE6KYI/rl+Xdi1PonPaPgPADWAOBUeULdAicqYA+A0R+QD4AdwohKhO8FoMwzAdDjUM0+zDt4rJTyUJCb4QYnCE9e8AeCeRczMMw2QDLkXwHSbBL3TZ27UsPNKWYRgmhTgdobh7s0unc4HT6pCUwYLPMAyTQlwRXDp9uubDYWtfCWbBZxiGSSHhPvygVW8nCsu1k2pY8BmGYVKI2jHbu3O+nl+nOM+JdkqSqZNolA7DMAwTBdXCL3Q7MLBbIW49ZTB+cFwfbD/U0K5lYQufYRgmhZhDL4kId585DH26FsDfTknTJCz4DMMwKcQqO6ZEdem0h/iz4DMMw6QQZ5RIHNXd0+DxRdwvWbAPn2EYJoXYbITrThqAs0b0DNs2aVApBnYrxLaDDahv9hli9FNSlpSenWEYhsED5x6D8f26hq232Qh3nTEMAFBvyp2fCljwGYZh0ogciFXfnHqXDgs+wzBMGgkJPlv4DMMwWU2hOyj4jR5/yq/Fgs8wDJNGbMHU8hyWyTAMk+3IfDqBdsizwILPMAyTRuy6hZ/6a7HgMwzDpBFN79nCZxiGyXZ0lw778BmGYbIbKfh+tvAZhmGyGxmlwxY+wzBMlqNb+B1B8Inot0S0iohWENHHRHSUtp6I6Gki2qJtH5d4cRmGYbILPUqnHVLjJ8PC/4MQYpQQYgyADwD8Ult/NoAh2t/1AP6WhGsxDMNkFTJ7codw6Qgh6pTFQgCy1NMBvCSCLADQmYh6JXo9hmGYbKI9O22Tkg+fiB4BcAWAwwBO0Vb3BrBL2a1SW7fXdOz1CLYA0Ldv32QUh2EYpsOgd9pmSpQOEc0lojUWf9MBQAjxgBCiD4BXANwqD7M4Vdg3EkLMEEJUCCEqysrK2vo9GIZhOiTtGaUTl4UvhJgW5/leBTALwK8QtOj7KNvKAexpVekYhmGynFCUTuqvlYwonSHK4vkANmifZwK4QovWmQDgsBBib9gJGIZhchhN7zuMD/8xIhoGIABgB4AbtfWzAZwDYAuARgBXJ+FaDMMwWQURwUYZ5NKJhhDiexHWCwC3JHp+hmGYbMduI06twDAMkwvYiDpGHD7DMAyTGHYbdYzUCgzDMExi2IldOgzDMDkBEdAOes+CzzAMk27YpcMwDJMjcJQOwzBMjsBROgzDMDkCu3QYhmFyBBtH6TAMw+QGdhu7dBiGYXKCYKdt6q+TlAlQGIZhmLZz8tAy9OlakPLrsOAzDMOkmYfOP7ZdrsMuHYZhmByBBZ9hGCZHYMFnGIbJEVjwGYZhcgQWfIZhmByBBZ9hGCZHYMFnGIbJEVjwGYZhcgQS7THNSpwQURWAHQmcohuAg0kqTkeF7wHfA4DvAZBb96CfEKIs1k4ZJfiJQkRLhBAV6S5HOuF7wPcA4HsA8D2wgl06DMMwOQILPsMwTI6QbYI/I90FyAD4HvA9APgeAHwPwsgqHz7DMAwTmWyz8BmGYZgIZIXgE9FZRLSRiLYQ0X3pLk+qIKI+RDSfiNYT0Voiuk1b35WI5hDRZu1/F209EdHT2n1ZRUTj0vsNkgcR2YloORF9oC0PIKKF2j14g4hc2nq3trxF294/neVOFkTUmYjeJqIN2vMwMdeeAyK6Q3sP1hDRa0SUl2vPQWvp8IJPRHYAfwVwNoBjAFxKRMekt1QpwwfgLiHEcAATANyifdf7AMwTQgwBME9bBoL3ZIj2dz2Av7V/kVPGbQDWK8uPA3hKuwc1AK7V1l8LoEYIMRjAU9p+2cCfAXwkhDgawGgE70XOPAdE1BvATwFUCCFGALADuAS59xy0DiFEh/4DMBHA/5Tl+wHcn+5ytdN3/w+A0wFsBNBLW9cLwEbt8/MALlX21/fryH8AyhEUtFMBfACAEBxg4zA/EwD+B2Ci9tmh7Ufp/g4Jfv9OAL41f49ceg4A9AawC0BX7Xf9AMCZufQctOWvw1v4CP3wkkptXVajNUnHAlgIoIcQYi8AaP+7a7tl6735E4CfAQhoy6UAaoUQPm1Z/Z76PdC2H9b278gMBFAF4F+aW+sfRFSIHHoOhBC7AfwRwE4AexH8XZcit56DVpMNgk8W67I69IiIigC8A+B2IURdtF0t1nXoe0NE5wE4IIRYqq622FXEsa2j4gAwDsDfhBBjATQg5L6xIuvugdY/MR3AAABHAShE0HVlJpufg1aTDYJfCaCPslwOYE+aypJyiMiJoNi/IoR4V1u9n4h6adt7ATigrc/Ge3MigPOJaDuA1xF06/wJQGcicmj7qN9Tvwfa9hIA1e1Z4BRQCaBSCLFQW34bwQogl56DaQC+FUJUCSG8AN4FMAm59Ry0mmwQ/MUAhmi98y4EO25mprlMKYGICMALANYLIZ5UNs0EcKX2+UoEffty/RValMYEAIdlk7+jIoS4XwhRLoToj+Bv/YkQ4ocA5gP4vrab+R7Ie/N9bf8ObdkJIfYB2EVEw7RVpwFYhxx6DhB05UwgogLtvZD3IGeegzaR7k6EZPwBOAfAJgBbATyQ7vKk8HtORrAZugrACu3vHAR9kfMAbNb+d9X2JwQjmLYCWI1gREPav0cS78dUAB9onwcCWARgC4C3ALi19Xna8hZt+8B0lztJ330MgCXas/A+gC659hwA+DWADQDWAHgZgDvXnoPW/vFIW4ZhmBwhG1w6DMMwTByw4DMMw+QILPgMwzA5Ags+wzBMjsCCzzAMkyOw4DMMw+QILPgMwzA5Ags+wzBMjvD/h+f9pCvAt/sAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x112ab8e10>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(potential_energy(visit_all))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"No missing ensembles.\n",
"No extra ensembles.\n"
]
}
],
"source": [
"initial_conditions = scheme.initial_conditions_from_trajectories([visit_all])"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Trajectory[139]"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"initial_conditions[0].trajectory"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create the path simulator and run it!"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"storage = paths.Storage(\"lammps_lj.nc\", 'w')"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Working on Monte Carlo cycle number 1000\n",
"Running for 14 minutes 30 seconds - 0.87 seconds per step\n",
"Estimated time remaining: 0 seconds\n",
"DONE! Completed 1000 Monte Carlo cycles.\n"
]
}
],
"source": [
"simulator = paths.PathSampling(\n",
" storage=storage,\n",
" move_scheme=scheme,\n",
" sample_set=initial_conditions\n",
")\n",
"\n",
"simulator.run(1000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"storage.close()"
]
}
],
"metadata": {
"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.15"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment