readi-bonds.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## reaction diffusion simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import uuid\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from IPython.display import HTML\n",
    "\n",
    "from copy import deepcopy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### define particles and trajectory classes\n",
    "only property of particle-types is diffusion and a geometry object\n",
    "\n",
    "**to clean up**\n",
    "particle types contain their uni-molecular reactions but this should be separated to mirror the case of bi-molecular reactions. i.e define inducts (uni and bi-), each induc class has a method to get their instances in the simulation.\n",
    "traj has uni- and b- reactors, that call the method, loop over the reactions associated with the inducts, and effects them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class energyterm(object):\n",
    "    pass\n",
    "\n",
    "class bondterm(energyterm):\n",
    "    def __init__(self, k, d, part1=None):\n",
    "        self.elastic=k\n",
    "        self.distance=d\n",
    "        self.partner = part1\n",
    "        self.name=uuid.uuid4()\n",
    "        \n",
    "    def assign_partner(self,part1):\n",
    "        self.partner = part1\n",
    "    \n",
    "    def force(self, position):\n",
    "        partner = self.partner\n",
    "        pos1 = partner.pos\n",
    "        ## we use geometry from partner, ... what happend if they are different?\n",
    "        if(partner.ptype.geometry.bc == 'pbc'):\n",
    "            box = partner.ptype.geometry.size\n",
    "            pos1 = pos1 - np.floor(pos1/box) * box\n",
    "            delta = (pos1 - position)\n",
    "            delta = delta - np.rint(delta/box)* box\n",
    "        else:\n",
    "            delta = (pos1 - position)\n",
    "        dist = np.sqrt((delta*delta).sum())\n",
    "\n",
    "\n",
    "        factor = self.elastic*(dist-self.distance)\n",
    "        force = factor*delta/dist\n",
    "        return force \n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class hardcore(energyterm):\n",
    "    def __init__(self, k, d, part1):\n",
    "        self.elastic=k\n",
    "        self.distance=d\n",
    "        self.partner = part1\n",
    "        self.name=uuid.uuid4()\n",
    "\n",
    "    def force(self, position):\n",
    "        partner = self.partner\n",
    "\n",
    "        pos1 = partner.pos\n",
    "        ## we use geometry from partner, ... what happend if they are different?\n",
    "        if(partner.ptype.geometry.bc == 'pbc'):\n",
    "            box = partner.ptype.geometry.size\n",
    "            pos1 = pos1 - np.floor(pos1/box) * box\n",
    "            delta = (pos1 - position)\n",
    "            delta = delta - np.rint(delta/box)* box\n",
    "        else:\n",
    "            delta = (pos1 - position)\n",
    "        dist = np.sqrt((delta*delta).sum())\n",
    "\n",
    "        if dist<self.distance:\n",
    "            factor = self.elastic*(dist-self.distance)\n",
    "            force = factor*delta/dist\n",
    "        else:\n",
    "            force = np.array([0.,0.])\n",
    "            \n",
    "        return force "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#p1 = ptype(1e-2,'p1', geometry(2,20.))\n",
    "#part1 = particle(ptype=p1,age=0,pos=np.array([0.,0.]),energyterms=[])\n",
    "#part2 = particle(ptype=p1,age=0,pos=np.array([1.,0.]),energyterms=[])\n",
    "#bond1 = bondterm(1,1.5,part2)\n",
    "#bond2 = bondterm(1,1.5,part1)\n",
    "\n",
    "\n",
    "#part1.add_energyterm(bond1)\n",
    "#part2.add_energyterm(bond2)\n",
    "\n",
    "#part2b = part2.step(10)\n",
    "#part2b.pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class ptype:\n",
    "    \"\"\"basic particle type, with name and physical properties\n",
    "    it also has a list of unimolecular reactions A -> products\n",
    "    \"\"\"\n",
    "    def __init__(self,diffusion,name,geometry):\n",
    "        self.diffusion=diffusion\n",
    "        self.name=name\n",
    "        self.unimolecular_reactions = []\n",
    "        self.unimolecular_probs = []\n",
    "        self.geometry= geometry\n",
    "        \n",
    "    def add_reaction(self,reaction):\n",
    "        \"\"\"adds a unimolecular_reaction to the list for this particle type\n",
    "        and updates the probabilities function\"\"\"\n",
    "        (self.unimolecular_reactions).append(reaction)\n",
    "        rates = np.array([rea.rate for rea in self.unimolecular_reactions])\n",
    "        self.unimolecular_probs = rates\n",
    "        \n",
    "    def choose_reaction(self,dt):\n",
    "        \"\"\"chooses one reaction (or none) for a given time step\n",
    "        that reaction can then be effected on a particle with effect_reaction\n",
    "        \"\"\"\n",
    "        if len(self.unimolecular_probs)==0:\n",
    "            return None\n",
    "        else:\n",
    "            rates = self.unimolecular_probs\n",
    "            probs = (rates/rates.sum())*(1 - np.exp( - dt* rates.sum()))\n",
    "            #print(\"Cumulative probs: \",np.cumsum(probs))\n",
    "            alpha=np.random.uniform(size=1)\n",
    "            conf = np.cumsum(probs)>alpha\n",
    "            if any(conf):\n",
    "                ii = np.argmax(conf)\n",
    "                #print(\"Chosen rection %s\" % self.unimolecular_reactions[ii].name)\n",
    "                return self.unimolecular_reactions[ii]\n",
    "            else:\n",
    "                return None\n",
    "        \n",
    "\n",
    "class geometry(object):\n",
    "    \"\"\"a basic basic geometry specification\n",
    "    ideally it would include boundary conditions, dimensionality\n",
    "    external potentials, etc. \n",
    "    It's used to evolve the dynamics of particles, \n",
    "    so each particle type is attached to a geometry object at construction\n",
    "    \"\"\"\n",
    "    def __init__(self, size, dim=2):\n",
    "        self.dim=dim\n",
    "        self.size=size\n",
    "        self.bc='pbc'\n",
    "        \n",
    "\n",
    "class particle:\n",
    "    \"\"\"particle instance, it has an id, position, age, types, etc.\"\"\"\n",
    "    def __init__(self, ptype, pos, age=0, id=None, energyterms=None):\n",
    "        if id==None:\n",
    "            self.id = uuid.uuid4()\n",
    "        else:\n",
    "            self.id = id\n",
    "        self.ptype = ptype\n",
    "        self.age = age\n",
    "        if energyterms==None:\n",
    "            self.energyterms = []\n",
    "        else:\n",
    "            print(\"creating particle with %d energy terms\" % len(energyterms))\n",
    "            self.energyterms=energyterms\n",
    "        self.pos = pos\n",
    "        \n",
    "    def add_energyterm(self, energyterm):\n",
    "        (self.energyterms).append(energyterm)\n",
    "            \n",
    "    def step(self,dt,kT=None):\n",
    "        if len(self.energyterms)>0 and not kT==None:\n",
    "            #print(\"force!\")\n",
    "            forces = [ eterm.force(self.pos) for eterm in self.energyterms]\n",
    "            force = np.vstack(forces).sum(axis=0)\n",
    "            dx = force*self.ptype.diffusion*dt/kT\n",
    "        else:\n",
    "            dx = 0\n",
    "        pos1 = self.pos + dx + np.sqrt(2*self.ptype.diffusion*dt)*np.random.normal(0, 1, 2)\n",
    "        \n",
    "                                        \n",
    "        if(self.ptype.geometry.bc == 'pbc'):\n",
    "            box = self.ptype.geometry.size\n",
    "            pos1=pos1 - np.floor(pos1/box) * box\n",
    "        #print(\"type %s, age %d\"% (self.ptype.name, self.age))\n",
    "        #print(\"pos\",pos1)\n",
    "        self.pos = pos1\n",
    "        self.age += 1\n",
    "        #pnew = particle(self.ptype, pos1,\n",
    "        #                age = self.age + 1, id = self.id)\n",
    "        return self\n",
    "\n",
    "    def deghost(self, swap=None):\n",
    "        \"\"\" removes energyterms referring to this particle \"\"\"\n",
    "        for eterm in self.energyterms:\n",
    "            part = eterm.partner\n",
    "            for partner_eterm in part.energyterms:\n",
    "                if (partner_eterm.partner)==self:\n",
    "                    print(\"deghosted %s\" % part.id)\n",
    "                    if swap==None:\n",
    "                        (part.energyterms).remove(partner_eterm)\n",
    "                    else:\n",
    "                        partner_eterm.partner=swap\n",
    "    \n",
    "    \n",
    "        \n",
    "class traj:\n",
    "    \"\"\"a traj is a list (frames) of lists of particles\"\"\"\n",
    "    \n",
    "    def __init__(self, particles=[]):\n",
    "        self.ptraj = []\n",
    "        self.nframes = 0\n",
    "        (self.ptraj).append(particles)\n",
    "\n",
    "    def __getitem__(self, key):\n",
    "        if isinstance(key, slice):\n",
    "            print(\"sliced trajectory\")\n",
    "        outt = traj()\n",
    "        outt.ptraj=self.ptraj[key]\n",
    "        outt.nframes = len(outt.ptraj)\n",
    "        return outt\n",
    "        \n",
    "    def describe(self):\n",
    "        n = len(self.ptraj)\n",
    "        print(\"Length %d (%d)\" % (n,self.nframes))\n",
    "        \n",
    "    def step(self,dt,kT=None):\n",
    "        \"\"\"updates the latest frame\"\"\"\n",
    "        newvals=[]\n",
    "        for pa in (self.ptraj)[self.nframes]:\n",
    "            pa1 = pa.step(dt,kT)\n",
    "            newvals.append(pa1)\n",
    "        (self.ptraj)[self.nframes] = newvals\n",
    "        #(self.ptraj).append(newvals)\n",
    "        #self.nframes = (self.nframes)+1\n",
    "        \n",
    "    def store(self):\n",
    "        \"\"\"hard-duplicates the latest frame\"\"\"\n",
    "        newvals = []\n",
    "        for pa in (self.ptraj)[self.nframes]:\n",
    "            pa1 = particle(ptype=pa.ptype, pos=pa.pos,age=pa.age,id=pa.id)\n",
    "            newvals.append(pa1)\n",
    "            \n",
    "        (self.ptraj).append((self.ptraj)[self.nframes])\n",
    "        (self.ptraj)[self.nframes]=newvals\n",
    "        self.nframes = (self.nframes)+1\n",
    "                \n",
    "\n",
    "    def react_unimolecular(self, dt):\n",
    "        ###### unimolecular reactions would be defined \n",
    "        ###### on unm_inducts to mimick other orders\n",
    "        ###### then here we can pass the set of unim_inducts \n",
    "        \"\"\"loops over particles and reacts them with their unimolecular reactions\n",
    "         unimolecular reactions are associated to the incoming pair\n",
    "        \"\"\"\n",
    "        newvals=[]\n",
    "        for pa in (self.ptraj)[self.nframes]:\n",
    "            react = pa.ptype.choose_reaction(dt)\n",
    "            if(not react==None):\n",
    "                pa1 = react.effect_reaction(pa)\n",
    "                if not pa1==None:\n",
    "                    for pa11 in pa1:\n",
    "                        newvals.append(pa11)\n",
    "            else:\n",
    "                newvals.append(pa)\n",
    "        \n",
    "        (self.ptraj)[self.nframes] = newvals\n",
    "        \n",
    "    def react_bimolecular(self, dt, bim_inducts, verbose=False):\n",
    "        \"\"\" bimolecular reactions        \n",
    "        \"\"\"\n",
    "        frame = (self.ptraj)[self.nframes]\n",
    "        if verbose:\n",
    "            print(\"reacting frame %d\" % self.nframes)\n",
    "        #loop over inducts, get instances, choose reaction, effect reaction\n",
    "        # remember to add the particles that are not in any instance\n",
    "        for induct in bim_inducts:\n",
    "            #print(\"induct name: %s\"%induct.name)\n",
    "            induct_instances = induct.get_ec(frame)\n",
    "            ##if len(induct_instances)>0:\n",
    "            ##    print(induct_instances)\n",
    "            for EC_instance in induct_instances:\n",
    "                react = induct.choose_reaction(dt)\n",
    "                if(not react==None):\n",
    "                    pa1 = react.effect_reaction(EC_instance)\n",
    "                    ## remove incoming particles from frame\n",
    "                    for part in EC_instance:\n",
    "                        #### EC_instance is a list of two (index,particle)\n",
    "                        #print(\"removing part %d\" % part.id)\n",
    "                        frame.remove(part)\n",
    "                    ## add educts\n",
    "                    if not pa1==None:\n",
    "                        for pa11 in pa1:\n",
    "                            frame.append(pa11)\n",
    "                            \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#col_dict = {'type 1':'red','type 1b':'blue','type 2':'orange'}\n",
    "#col_dict = {'type 1':1,'type 1b':2,'type 2':3}\n",
    "#size_dict = {'type 1':10,'type 1b':20,'type 2':30}\n",
    "\n",
    "### http://stackoverflow.com/questions/9401658/how-to-animate-a-scatter-plot\n",
    "\n",
    "def plt_traj2(traj, col_dict, size_dict):\n",
    "    \"\"\"plots a trajectory\"\"\"\n",
    "    ptraj = (traj.ptraj)\n",
    "    \n",
    "    # the first particle determines the size (should be the max here)\n",
    "    size = (ptraj[0][0]).ptype.geometry.size\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,8))\n",
    "    ax = plt.axes(xlim=(0, size), ylim=(0, size))\n",
    "    scat = ax.scatter([],[], c=[], marker='o', lw=2, vmin = 0, vmax=20)\n",
    "        \n",
    "    def init():\n",
    "        scat.set_offsets(np.hstack(([], [])))\n",
    "        return scat,\n",
    "        \n",
    "    def animate(i):\n",
    "        x = [p.pos[0] for p in ptraj[i]]\n",
    "        y = [p.pos[1] for p in ptraj[i]]\n",
    "        cols = [1.*col_dict[p.ptype.name] for p in ptraj[i]]\n",
    "        sizs = [size_dict[p.ptype.name] for p in ptraj[i]]\n",
    "        scat.set_offsets(np.vstack((x, y)).T)\n",
    "        scat._sizes = np.array(sizs)*35\n",
    "        scat.set_array(np.array(cols))\n",
    "        return scat,\n",
    "\n",
    "    anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                               frames=traj.nframes, interval=20, blit=True)\n",
    "    return anim\n",
    "\n",
    "\n",
    "### reshape trajectory for analysis/postprocessing\n",
    "def reshape(traj):\n",
    "    \"\"\"converts a frame-based trajectory into a id-based data-frame\n",
    "    for lifetimes and concentrations\n",
    "    \n",
    "    \"\"\"\n",
    "    pp=[]\n",
    "    for i,frame in enumerate(traj.ptraj):\n",
    "        for part in frame:\n",
    "            pp.append([part.id,part.age,part.pos[0],part.pos[1],i,part.ptype.name])\n",
    "    gg = pd.DataFrame(pp,columns=['ID','age','x','y','frame','type_name'])\n",
    "    return gg\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## define reactions classes\n",
    "\n",
    "reaction class is the used for uni- and bi- molecular reactions\n",
    "each has a method to effect it (i.e. apply it to the inducts and yield the products)\n",
    "and one single rate\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class reaction:\n",
    "    \"\"\"base class for reactions, just a name, a rate and a counter\"\"\"\n",
    "    def __init__(self,name,rate):\n",
    "        self.rate = rate\n",
    "        self.name = name\n",
    "        self.counter = 0\n",
    "        #print(self)\n",
    "        \n",
    "    def report(self):\n",
    "        reaction_type = type(self).__name__\n",
    "        print(\"   %s %s: counter %d\" % (reaction_type, self.name, self.counter))\n",
    "        print(\"        rate: %d\"% self.rate)\n",
    "        \n",
    "    def reset(self):\n",
    "        self.counter = 0\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### uni-molecular reactions\n",
    "ptypes are the inducts (this should be changed)\n",
    "the can be\n",
    "* decay A->0\n",
    "* convert A->B \n",
    "* fission A-> B+C (e.g. unbinding, or dimerization breaking)\n",
    "\n",
    "once could add multiple fission types A->B+C+D, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### ============== UNI_MOLECULAR REACTIONS ===============================\n",
    "        \n",
    "class decay_reaction(reaction):\n",
    "    \"\"\" unimolecular decay with given rate\n",
    "        A -> 0\n",
    "    \"\"\"\n",
    "    def __init__(self, name, rate):\n",
    "        super(decay_reaction,self).__init__(name,rate)\n",
    "        \n",
    "    def effect_reaction(self,particle_in):\n",
    "        self.counter += 1\n",
    "        \n",
    "        ## since we destroy a particle, we need to deghost\n",
    "        particle_in.deghost()\n",
    "        return None\n",
    "    \n",
    "class convert_reaction(reaction):\n",
    "    \"\"\" unimolecular conversion with given rate to given type \n",
    "        A -> B\n",
    "        \n",
    "        the new particle is actually the old one, with new age and type\n",
    "        to allow keeping the energy terms in case this is needed.\n",
    "        default is to keep_energy, use keepenergy=None not to\n",
    "    \"\"\"\n",
    "    def __init__(self,  name, rate, outtype,  keepenergy = 0):\n",
    "        self.outtype=outtype\n",
    "        self.keepenergy = keepenergy\n",
    "        super(convert_reaction,self).__init__( name, rate )\n",
    "        \n",
    "    def effect_reaction(self,particle_in):\n",
    "        self.counter += 1\n",
    "        if self.keepenergy==None:\n",
    "            particle_in.energyterms=[]\n",
    "        ## otherwise energyterms are kept\n",
    "        particle_in.ptype= self.outtype\n",
    "        particle_in.age= 0 \n",
    "        particle_in.id = uuid.uuid4()\n",
    "        ## return a singleton list for consistence with other reactions\n",
    "        return [particle_in]\n",
    "        #return [particle(self.outtype, particle_in.pos )]\n",
    "\n",
    "class fission_reaction(reaction):\n",
    "    \"\"\" unimolecular fission with given rate to given type \n",
    "        A -> B+C\n",
    "        default is to loose all energyterms\n",
    "        use keepenergy=0 to keep them to the first educt, \n",
    "        keepenergy=1 to pass it to the second\n",
    "    \"\"\"\n",
    "    def __init__(self,  name, rate, outtype1, outtype2, keepenergy=None):\n",
    "        self.outtype1=outtype1\n",
    "        self.outtype2=outtype2\n",
    "        self.keepenergy = keepenergy\n",
    "        super(fission_reaction,self).__init__( name, rate)\n",
    "        \n",
    "    def effect_reaction(self,particle_in):\n",
    "        keepenergy = self.keepenergy\n",
    "        self.counter += 1\n",
    "        pos = particle_in.pos\n",
    "        types=[self.outtype1,self.outtype2]\n",
    "        if not keepenergy==None and keepenergy<2:\n",
    "            particle_in.ptype= types[keepenergy]\n",
    "            particle_in.age= 0 \n",
    "            particle_in.id = uuid.uuid4()\n",
    "            particle_out1 = particle_in\n",
    "            particle_out2 = particle(types[1-keepenergy], pos )\n",
    "        else:\n",
    "            particle_out1 = particle(types[0], pos )\n",
    "            particle_out2 = particle(types[1], pos )\n",
    "        return [particle_out1 , particle_out2  ]  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99c96c7f-cab5-4f2d-96d5-67e6fd818426\n",
      "998dabce-a43e-40dc-89d8-615efb9f964b\n"
     ]
    }
   ],
   "source": [
    "membrane=geometry(dim=2,size=20)\n",
    "p1 = ptype(1,\"A1\",membrane)\n",
    "p2 = ptype(1,\"B1\",membrane)\n",
    "\n",
    "#conv = convert_reaction(name=\"con1\",outtype=p2,rate=10,keepenergy=0)\n",
    "#fiss = fission_reaction(name='fission',outtype1=p1,outtype2=p2,rate=2,keepenergy=0)\n",
    "#p1.add_reaction(conv)\n",
    "\n",
    "part1 = particle(p1,np.array([0,0]))\n",
    "part2 = particle(p1,np.array([0,0]))\n",
    "\n",
    "part1.add_energyterm(bondterm(1,1,part2))\n",
    "part2.add_energyterm(bondterm(1,1,part1))\n",
    "\n",
    "print(part1.id)\n",
    "\n",
    "print(part2.id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "part1.deghost()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<__main__.bondterm at 0x7ff11592e908>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "part2.energyterms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#part1.add_energyterm(bond1)\n",
    "#print(part1.id)\n",
    "#part1.energyterms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#part1 = particle(p1,np.array([0,0]))\n",
    "#print(part1.id)\n",
    "#part1.energyterms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#pout = conv.effect_reaction(part1)\n",
    "#print(pout[0].id)\n",
    "#pout[0].energyterms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#pout = fiss.effect_reaction(part1)\n",
    "#pout[1].energyterms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### bi-molecular reactions\n",
    "encounter-complexes are the inducts of bi-molecular reactions\n",
    "they contain their reactions as ptypes contain the uni-molecular reactions\n",
    "they can be:\n",
    "* annihilation A+B -> 0\n",
    "* fusion A+B -> C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#### abstract types of EC\n",
    "### zz=bim_encunter_complex(\"bim1\",p1,p3,.9)\n",
    "#### add reactions to the EC_type\n",
    "### r1 = annihilation_reaction(\"anni1\", 1.1)\n",
    "### r2 = ...\n",
    "### zz.add_reaction(r1)\n",
    "### zz.add_reaction(r2)\n",
    "\n",
    "#### don't confuse with realizations of ECs in the simulation that can be obtained with\n",
    "### zz.get_ec(t.ptraj[100])\n",
    "\n",
    "### for a timestep, cycle over ectypes, get all ECs, \n",
    "### choose reaction for each ec, apply reaction\n",
    "\n",
    "class encunter_complex_type:\n",
    "    def __init__(self, name):\n",
    "        self.name=name\n",
    "        \n",
    "class bim_encunter_complex(encunter_complex_type):\n",
    "    \n",
    "    def __init__(self, name, typ1,typ2, dist):\n",
    "        self.type1 = typ1\n",
    "        self.type2 = typ2\n",
    "        self.dist= dist\n",
    "        self.bim_reactions = []\n",
    "        self.bim_probs = []\n",
    "        super(bim_encunter_complex,self).__init__(name)\n",
    "        \n",
    "    def describe(self):\n",
    "        print(\"Name: %s\" % self.name)\n",
    "        print(\"   Reactions:\")\n",
    "        for reac in self.bim_reactions:\n",
    "            reac.report()\n",
    "        \n",
    "    def add_reaction(self,reaction):\n",
    "        \"\"\"adds a unimolecular_reaction to the list for this particle type\n",
    "        and updates the probabilities function\"\"\"\n",
    "        \n",
    "        ## should check if reaction type is ok ...\n",
    "        (self.bim_reactions).append(reaction)\n",
    "        rates = np.array([rea.rate for rea in self.bim_reactions])\n",
    "        self.bim_probs = rates\n",
    "\n",
    "    def choose_reaction(self,dt):\n",
    "        \"\"\"chooses one reaction (or none) for a given time step\n",
    "        that reaction can then be effected on a particle with effect_reaction\n",
    "        \"\"\"\n",
    "        if len(self.bim_reactions)==0:\n",
    "            return None\n",
    "        else:\n",
    "            rates = self.bim_probs\n",
    "            probs = (rates/rates.sum())*(1 - np.exp( - dt* rates.sum()))\n",
    "            #print(\"Cumulative probs: \",np.cumsum(probs))\n",
    "            alpha=np.random.uniform(size=1)\n",
    "            conf = np.cumsum(probs)>alpha\n",
    "            if any(conf):\n",
    "                ii = np.argmax(conf)\n",
    "                #print(\"Chosen rection %s\" % self.bim_reactions[ii].name)\n",
    "                return self.bim_reactions[ii]\n",
    "            else:\n",
    "                return None\n",
    "            \n",
    "    def get_ec(self,particles_):\n",
    "        \"\"\"returns a list of instances of hte ECtype\n",
    "        each of these are a vector of (two) (index, particle)\n",
    "        bimolecular reactions effect_reaction functions take one of the instances\n",
    "        and modify it\n",
    "        \n",
    "        the bimol react function of the traj loops over the instances,\n",
    "        picks a reaction and reacts it\n",
    "        \"\"\"\n",
    "        p1 = [i for i,p in enumerate(particles_) if p.ptype==self.type1]\n",
    "        p2 = [i for i,p in enumerate(particles_) if p.ptype==self.type2]\n",
    "        ec=[]\n",
    "        #report_pairs=[]\n",
    "        #print(p1)\n",
    "        #print(p2)\n",
    "        used_ones=[]\n",
    "        for q1 in p1:\n",
    "            for q2 in p2:\n",
    "                if not q1 == q2 and not ((q1 in used_ones) or (q2 in used_ones)):\n",
    "                    vec1= particles_[q1].pos - particles_[q2].pos\n",
    "                    #print(vec1)\n",
    "                    dist1= np.sqrt((vec1 * vec1).sum())\n",
    "                    if dist1 < self.dist:\n",
    "                        #report_pairs.append([q1,q2])\n",
    "                        ec.append([particles_[q1],particles_[q2]])\n",
    "                        ## we want particles to be in only one instance\n",
    "                        ## otherwise we can react safely\n",
    "                        used_ones.append(q1)\n",
    "                        used_ones.append(q2)\n",
    "        \n",
    "        #print(report_pairs)\n",
    "        return ec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class annihilation_reaction(reaction):\n",
    "    \"\"\" annihilation reaction\n",
    "        A + B -> 0\n",
    "    \"\"\"\n",
    "    def __init__(self, name, rate):\n",
    "        super(annihilation_reaction,self).__init__(name,rate)\n",
    "        \n",
    "    def effect_reaction(self,encunter_complex_in):\n",
    "        \n",
    "        for part in encunter_complex_in:\n",
    "            part.deghost()\n",
    "        self.counter += 1\n",
    "        return None\n",
    "\n",
    "class fusion_reaction(reaction):\n",
    "    \"\"\" fusion reaction\n",
    "        A + B -> C\n",
    "        default is to keep all energy (e.g. keep bonds to both inducts)\n",
    "    \"\"\"\n",
    "    def __init__(self, name, rate, outtype,keepenergy=[0,1]):\n",
    "        self.outtype = outtype\n",
    "        self.keepenergy=keepenergy\n",
    "        super(fusion_reaction,self).__init__(name,rate)\n",
    "        \n",
    "    def effect_reaction(self,encunter_complex_in):\n",
    "        self.counter += 1\n",
    "        pos = None ## get particle pos in encounter complex in and average\n",
    "        poss = [part.pos for part in encunter_complex_in]\n",
    "        #print([part[0] for part in encunter_complex_in])    \n",
    "        if not self.keepenergy == None:\n",
    "            ## if only one, copy it and deghost the other\n",
    "            ## if both, copy one, deghost the other assigning to the kept particle\n",
    "            raise\n",
    "        return [particle(self.outtype, poss[0])]\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class ensemble(object):\n",
    "    def __init__(self, trajs=[]):\n",
    "        if not isinstance(trajs, list):\n",
    "            trajs=[trajs]\n",
    "        self.trajs=trajs\n",
    "        self.params=[]\n",
    "    \n",
    "    def describe(self):\n",
    "        print(\"Ensemble has %d trajectories\" % len(self.trajs))\n",
    "        for t in self.trajs:\n",
    "            t.describe()\n",
    "    \n",
    "    def add_traj(self,traj):\n",
    "        (self.trajs).append(traj)\n",
    "        \n",
    "    def lifetimes(self,doplot=True):\n",
    "        \"\"\"plots the lifetime distributions of the different species in the ensemble\n",
    "        returns a grouped set by species with the individual lifetimes\n",
    "        \"\"\"\n",
    "        res_list=[]\n",
    "        trj_index = np.array(range(len(self.trajs))).astype('str')\n",
    "        for t in self.trajs:\n",
    "            ## get lifetime for all particles, flatten and drop ID\n",
    "            tr = reshape(t).groupby(by=['ID','type_name']).agg({'age':'max'})\n",
    "            tr2 = tr.rename(index=str).reset_index().drop('ID',axis=1)\n",
    "            res_list.append(tr2)\n",
    "        lifes = pd.concat(res_list,ignore_index=True)\n",
    "        if doplot:\n",
    "            plt.figure()\n",
    "            axes = lifes.hist(by='type_name', bins=30)\n",
    "            #for ax in np.ravel(axes):\n",
    "            #    ax.set_xlabel(\"lifetime (frames)\")\n",
    "            #    ax.set_ylabel(\"freq\")\n",
    "            plt.tight_layout()\n",
    "            plt.show()\n",
    "        lifes2 = lifes.groupby('type_name')\n",
    "        return lifes2\n",
    "    \n",
    "    def concentrations(self,species,init_val=0,doplot=True):\n",
    "        ### get concentrations as a function of frame\n",
    "        ### group by frame and type, then count particles of each type\n",
    "        res_list=[]\n",
    "        for i,t in enumerate(self.trajs):\n",
    "            ## get lifetime for all particles, flatten and drop ID\n",
    "            tr = reshape(t).groupby(by=['frame','type_name']).agg({'ID':'count'})\n",
    "            tr['trj_in']=str(i)\n",
    "            tr2 =tr.rename(index=str).reset_index()\n",
    "            res_list.append(tr2)\n",
    "        counts = pd.concat(res_list,ignore_index=True)\n",
    "        \n",
    "        \n",
    "        ## group_by species and select the one we want\n",
    "        rr=counts.groupby('type_name').get_group(species)\n",
    "        ## clean up the pandas mess and sort\n",
    "        rr=rr.fillna(0).drop('type_name',1)\n",
    "        rr3=rr.groupby('frame').agg({'ID':{'mea':'mean', 'sd':'std'}}).rename(index=str).reset_index()\n",
    "        rr4=rr3.set_index(rr3.frame)\n",
    "        rr4.index = rr4.index.astype('int')\n",
    "        rr5=rr4.sort_index()\n",
    "\n",
    "        ## fit with specified initial value\n",
    "        import scipy\n",
    "        vals1=scipy.optimize.curve_fit(\n",
    "        lambda t,a,b: a*(1-np.exp(-t/b)),\n",
    "            np.array(rr5.ID.mea.index),  np.array(rr5.ID.mea)-init_val )\n",
    "\n",
    "        if doplot:\n",
    "            fun = lambda t: vals1[0][0]*(1-np.exp(-t/vals1[0][1]))\n",
    "\n",
    "            plt.plot(np.array(rr5.ID.mea.index),  np.array(rr5.ID.mea), label=species)\n",
    "            plt.plot(np.array(rr5.ID.mea.index), fun(np.array(rr5.ID.mea.index)))\n",
    "            plt.legend()\n",
    "            plt.show()\n",
    "            \n",
    "        return vals1,rr5\n",
    "\n",
    "        \n",
    "#    def save(self,filename):        \n",
    "#    def load(self,filename):\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## define a class to run the simulation\n",
    "wraps in progress bar class, calls the integrator for diffusive motion, then reacts uni-molecular and bi-molecular (tri-molecular, etc.) reactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyemma._base.progress import ProgressReporter\n",
    "\n",
    "class runner(ProgressReporter):\n",
    "    def __init__(self, nframes, timestep, stride=1):\n",
    "        self.nframes = nframes\n",
    "        self.timestep = timestep\n",
    "        self.stride= stride\n",
    "        self._progress_register(nframes,\"Running simulation\")\n",
    "    \n",
    "    def run(self, t2, reactions_mono= None, reactions_bim=None,verbose=True ,kT=None):\n",
    "        if verbose:\n",
    "            t2.describe()\n",
    "        for i in range(self.nframes):\n",
    "            if (i % self.stride)==0:\n",
    "                #print(\"storing: %d\" % i)\n",
    "                t2.store()\n",
    "            t2.step(self.timestep,kT)\n",
    "            t2.react_unimolecular(self.timestep)\n",
    "            if reactions_bim!=None:\n",
    "                t2.react_bimolecular(self.timestep, reactions_bim)   \n",
    "            self._progress_update(1)\n",
    "        if verbose:\n",
    "            t2.describe()\n",
    "        ### report bimolecular complexes reactions\n",
    "        if verbose:\n",
    "            if reactions_bim!=None:\n",
    "                for zz in reactions_bim:\n",
    "                    zz.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def random_init(n, types_list, fractions=[], cellsize=20):\n",
    "    \"\"\"aux function to create a randomized initial state \n",
    "    with a concentrations or the \"\"\"\n",
    "    if not isinstance(types_list,list):\n",
    "        types_list=[types_list]\n",
    "    if len(fractions)==0:\n",
    "        fractions=np.ones(len(types_list))\n",
    "        fractions = fractions/fractions.sum()\n",
    "    counts = (n*fractions).astype('int')\n",
    "    init_pos=[]\n",
    "    for i, parttype in enumerate(types_list):\n",
    "        for j in range(counts[i]):\n",
    "            p = particle(parttype, \n",
    "                         pos=np.random.uniform(high=cellsize,low=0,size=2),age=0)\n",
    "            init_pos.append(p)\n",
    "            \n",
    "    return traj(init_pos)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SIMULATION "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### initialize the system\n",
    "define which particles we have, and how they react."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# analysis and post processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### do a movie of the trajectory (can be slowish to generate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "col_dict = {'mono':0,'dimer':7, 'mono_mem':6, 'dimer2':15, 'mono_mem2':14, 'dimer4':20}\n",
    "size_dict = {'mono':2,'dimer':5,'mono_mem':2, 'dimer2':5, 'mono_mem2':2, 'dimer4':5}\n",
    "aa = plt_traj2(sett[1],col_dict, size_dict)\n",
    "HTML(aa.to_html5_video())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "aa.save('im.mp4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### liftimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## analyse life-time distribution and concentration of all species as a functin of time\n",
    "## remember final semi column to prevent double plot\n",
    "zz=em1.lifetimes(doplot=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t=zz.get_group('dimer2')\n",
    "uu=t.drop('type_name',1).dropna()\n",
    "hist, bin_edges = np.histogram(uu.age,bins=50)\n",
    "plt.plot(hist,marker='o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### concentration and dimeric fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ff = (em1.concentrations())\n",
    "gg=ff.groupby('type_name')\n",
    "rr=gg.get_group('dimer4')\n",
    "rr=rr.fillna(0).drop('type_name',1)\n",
    "#rr2 = rr.pivot(index='frame', columns='trj_in', values='ID')\n",
    "rr3=rr.groupby('frame').agg({'ID':{'mea':'mean', 'sd':'std'}}).rename(index=str).reset_index()\n",
    "#rr3.index = rr3.index.astype('int')\n",
    "#rr3.head()\n",
    "rr4=rr3.set_index(rr3.frame)\n",
    "rr4.index = rr4.index.astype('int')\n",
    "rr5=rr4.sort_index()\n",
    "#rr5.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rr5.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### concentration fluctuations at equilibrium\n",
    "calculate the full distribution of concentrations at steady state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### do for different concentrations and diffusions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def setup_system(copynumber, diffusion, rates):\n",
    "    \"\"\"\n",
    "    setup a trajectory for different system parameters\n",
    "    \"\"\"\n",
    "\n",
    "    ## create particle types and their unimolecular reactions\n",
    "    membrane = geometry(dim=2,size=20)\n",
    "    q1 = ptype(diffusion['mono'],\"mono\",membrane)\n",
    "    q2 = ptype(diffusion['dimer'],\"dimer\",membrane)\n",
    "    \n",
    "    ## dimer disruption into states that remember their orientation\n",
    "    q2.add_reaction( \n",
    "        fission_reaction(\"fission 1\", \n",
    "                         rate=rates['on'], \n",
    "                         outtype1=q1, outtype2=q1, keepenergy=None) )\n",
    "    ## Dimerization 1: monomers bind to form interface, from any angular state\n",
    "    zz1=bim_encunter_complex(\"dimerization_EC\",q1,q1,dist=0.3)\n",
    "    r23 = fusion_reaction(\"dimerization 1\", rate=rates['off'], outtype=q2, keepenergy=None)\n",
    "    zz1.add_reaction(r23)\n",
    "\n",
    "    ## collect all the bimolecular reactions\n",
    "    bim_inducts = [zz1]\n",
    "\n",
    "    t2 = random_init(copynumber,[q1])\n",
    "    t2.describe()\n",
    "    return (t2 , bim_inducts)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "n_replicas=20\n",
    "totvals=[]\n",
    "\n",
    "params = itertools.product([1e-1,5e-2,1e-2],[0.05,0.1,0.3])\n",
    "for iset, (diff,rateon) in enumerate(params):\n",
    "    sett = []\n",
    "    for i in range(n_replicas):\n",
    "        t2, reacts = setup_system(30,\n",
    "            diffusion = {'mono':diff, 'dimer':1.5e-2} , \n",
    "            rates = {'on':rateon ,'off':10. })\n",
    "        run_1000steps = runner(20000, timestep = .1, stride=10)\n",
    "        run_1000steps.run(t2, reactions_bim = reacts ,verbose=False)\n",
    "        sett.append(t2)\n",
    "    \n",
    "    em1 = ensemble(sett)\n",
    "    em1.describe()\n",
    "\n",
    "    val, rr5 = em1.concentrations(\"dimer\",init_val=0)\n",
    "    totvals.append((diff,rateon, val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res1=pd.DataFrame([[va[0],va[1],va[2][0][0],va[2][0][1]] for va in totvals],\n",
    "            columns=['Diffusion','RateOn','fraction','time']\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res2=res1.groupby('Diffusion')\n",
    "res2.indices\n",
    "res2.get_group(0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "species='dimer'\n",
    "init_val = 0\n",
    "## group_by species and select the one we want\n",
    "rr=counts.groupby('type_name').get_group(species)\n",
    "## clean up the pandas mess and sort\n",
    "rr=rr.fillna(0).drop('type_name',1)\n",
    "rr3=rr.groupby('frame').agg({'ID':{'mea':'mean', 'sd':'std'}}).rename(index=str).reset_index()\n",
    "rr4=rr3.set_index(rr3.frame)\n",
    "rr4.index = rr4.index.astype('int')\n",
    "rr5=rr4.sort_index()\n",
    "\n",
    "## fit with specified initial value\n",
    "import scipy\n",
    "vals1=scipy.optimize.curve_fit(\n",
    "    lambda t,a,b: a*(1-np.exp(-b*t)),  \n",
    "                     np.array(rr5.ID.mea.index),  np.array(rr5.ID.mea)-init_val )\n",
    "\n",
    "return vals1,rr5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zz=em1.lifetimes(doplot=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "col_dict = {'mono':0,'dimer':7}\n",
    "size_dict = {'mono':2,'dimer':5}\n",
    "aa = plt_traj2(em1.trajs[0],col_dict, size_dict)\n",
    "HTML(aa.to_html5_video())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test for bonds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n",
      "[1, 2, 1]\n"
     ]
    }
   ],
   "source": [
    "membrane = geometry(dim=2,size=20)\n",
    "q0 = ptype(0,\"anchor\",membrane)\n",
    "q1 = ptype(5e-2,\"mono\",membrane)\n",
    "\n",
    "\n",
    "rel=decay_reaction(name='release',rate=.005)\n",
    "q0.add_reaction(rel)\n",
    "\n",
    "q = particle(q0,pos=np.array([10,10]),age=0)\n",
    "init_pos=[q]\n",
    "bonds=[bondterm(10,.1,q)]\n",
    "\n",
    "for i in range(1,3):\n",
    "    p = particle(q1,pos=np.array([10,10+i]),age=0)\n",
    "    init_pos.append(p)\n",
    "    b = bondterm(10,.1,p)\n",
    "    bonds.append(b)\n",
    "\n",
    "print(len(bonds))    \n",
    "\n",
    "for i,p in enumerate(init_pos):\n",
    "    if i<len(bonds)-1:\n",
    "        p.add_energyterm(bonds[i+1])\n",
    "    if i>0:\n",
    "        p.add_energyterm(bonds[i-1])\n",
    "\n",
    "print([len(p.energyterms) for p in init_pos])\n",
    "t = traj(init_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "runnr = runner(2000, timestep = .1, stride=10)\n",
    "runnr.run(t,verbose=True, kT=.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "col_dict = {'mono':0,'anchor':7}\n",
    "size_dict = {'mono':2,'anchor':5}\n",
    "aa = plt_traj2(t,col_dict, size_dict)\n",
    "HTML(aa.to_html5_video())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "membrane = geometry(dim=2,size=20)\n",
    "q0 = ptype(0,\"anchor\",membrane)\n",
    "q1 = ptype(5e-2,\"mono\",membrane)\n",
    "\n",
    "\n",
    "#rel=decay_reaction(name='release',rate=.005)\n",
    "#q0.add_reaction(rel)\n",
    "\n",
    "q = particle(q0,pos=np.array([10,10]),age=0)\n",
    "init_pos=[q]\n",
    "bonds=[bondterm(10,.1,q)]\n",
    "\n",
    "for i in range(1,5):\n",
    "    p = particle(q1,pos=np.array([10,10+i]),age=0)\n",
    "    init_pos.append(p)\n",
    "    b = bondterm(10,.1,p)\n",
    "    bonds.append(b)\n",
    "\n",
    "print(len(bonds))    \n",
    "\n",
    "for i,p in enumerate(init_pos):\n",
    "    if i<len(bonds)-1:\n",
    "        p.add_energyterm(bonds[i+1])\n",
    "    if i>0:\n",
    "        p.add_energyterm(bonds[i-1])\n",
    "\n",
    "print([len(p.energyterms) for p in init_pos])\n",
    "t = traj(init_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "runnr2 = runner(10000, timestep = .1, stride=10)\n",
    "runnr2.run(t,verbose=True, kT=.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pos4 = [np.vstack([frame[i].pos for frame in t.ptraj]) for i in range(4)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[plt.plot(pos4[i][:,0],pos4[i][:,1]) for i in range(4)]\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "col_dict = {'mono':0,'anchor':7}\n",
    "size_dict = {'mono':2,'anchor':5}\n",
    "aa = plt_traj2(t,col_dict, size_dict)\n",
    "HTML(aa.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### hard spheres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "membrane = geometry(dim=2,size=5)\n",
    "q1 = ptype(5e-2,\"mono\",membrane)\n",
    "a = particle(q1,pos=np.array([2,2]),age=0)\n",
    "b = particle(q1,pos=np.array([0,0]),age=0)\n",
    "init_pos=[a,b]\n",
    "hca = hardcore(50.,.5,a)\n",
    "hcb = hardcore(50.,.5,b)\n",
    "b.add_energyterm(hca)\n",
    "a.add_energyterm(hcb)\n",
    "cov1 = bondterm(k=10.,d=.8)\n",
    "\n",
    "ec1=bim_encunter_complex(\"grab\",q1,q1,dist=1)\n",
    "rea1 = bond_formation_reaction(\"dimerization 1\", rate=.1, bonds=cov1)\n",
    "ec1.add_reaction(rea1)\n",
    "rea2 = bond_breaking_reaction(\"disruption 1\", rate=.01, bonds=cov1)\n",
    "ec1.add_reaction(rea2)\n",
    "t = traj(init_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length 1 (0)\n"
     ]
    }
   ],
   "source": [
    "runnr2 = runner(200000, timestep = .1, stride=20)\n",
    "runnr2.run(t,reactions_bim=[ec1], verbose=True, kT=.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "col_dict = {'mono':0,'anchor':7}\n",
    "size_dict = {'mono':2,'anchor':5}\n",
    "aa = plt_traj2(t,col_dict, size_dict)\n",
    "HTML(aa.to_html5_video())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXn8FVX5xz8PoKCgLCYoEOKWioaIhZipXykX+CmQUlaW\nIWrkkrgrtoCmZWqG/oTUMJes9BcqUgJiAiqoYMBXEVAhRfCLiAsgKIvg8/vj3OnOnTvLmZkzc87c\ne96v133N3DNnzj7nOetziJlhsVgslvqkme4AWCwWi0UfVghYLBZLHWOFgMVisdQxVghYLBZLHWOF\ngMVisdQxVghYLBZLHSMtBIioGRHNJ6JJAc9vJ6KlRNRIRL3UBdFisVgsWRGnJzACwGK/B0TUH8C+\nzLw/gOEA7lQQNovFYrFkjJQQIKKuAAYAGB9gZRCABwCAmecAaEtEnZSE0GKxWCyZIdsT+D2AKwAE\nbS/uAmCl639TycxisVgsBhMpBIjofwC8x8yNAKj0s1gsFksN0ELCzlEABhLRAAA7AdiFiB5g5jNd\ndpoAfNH1v2vJrAIisoqKLBaLJQHMnEkDPLInwMzXMHM3Zt4HwHcBTPcIAACYBOBMACCivgDWMfN7\nAe7V7G/UqFHaw2DjZ+NXb3Grh/hliUxPwBciGg6AmfluZp5MRAOIaBmATwCcpSyEFovFYsmMWEKA\nmZ8B8Ezp/i7PswsVhstisVgsOWB3DCukoaFBdxAyxcavuNRy3IDaj1+WUNbjTRWeEXGe/lksFkst\nQERgXRPDFovFYqldrBCwWCyWOsYKAYvFYqljrBCwWCyWOsYKAYvFYqljrBCwWCyWOsYKAYvFYqlj\nrBCwWCyWOsYKAYvFYqljrBCwWCyWOsYKAYvFYqljrBCwWCyWOsYKAYvFYqljrBCwWCyWOsYKAYvF\nYqljrBCwWCyWOsYKAYvFYqljIoUAEbUkojlEtICIFhLRKB87xxLROiKaX/r9PJvgWiwWi0UlkQfN\nM/MWIjqOmT8louYAZhPRFGae67H6LDMPzCaYFovFYskCqeEgZv60dNsSQnD4HRScyfmXFovFYskO\nKSFARM2IaAGA1QCeYuaXfKwdSUSNRPQEEfVQGso6ZsoU4JxzdIdCH++/D3z6abQ9i8WSDNmewOfM\nfBiArgCO8Knk5wHoxsy9ANwBYKLaYNYv48cD99yjOxT66NgRGDpUbxguuAD4znf0hqFojBsHvPhi\nsnfvvx+YPVtteCzBRM4JuGHmj4loBoCTACx2mW903U8honFE1IGZP/K6MXr06P/eNzQ0oKGhIUGw\nk9G8ObB6NbD77rl5mRvbtgEtYuVmcXjnHb3+P/gg8PHHesNQNC64AOjXD3j66fjvDh0K9O4NzJun\nPFiFYebMmZg5c2Y+njFz6A/AFwC0Ld3vBOBZAAM8djq57vsAWB7gFusEYL7tNq1BiM2pp4pwRwEw\nr16dfXjyBmA+8ki9Ydh1V7k8yIqtW5mXLIn3zq9/zbx8eTbhkQFg7tcv+PnSpcwrVwa/27t3NuGK\nw8SJzJMnx3vnxReZ+/dXH5ZS3RlZXyf5yQwH7QlgBhE1ApgD4ElmnkxEw4noxyU7Q4jo1dK8wRgA\np6sSUqq56SbdIZDja18D3nor3jsbNmQTFlVMnw789a/x32O/ZQga2LAB+Owz9e6+9ZboyQUxbhxw\n0EHx3LzmGuC++1IFSzk33wxs3y7u999flPEgovKcCJg2TV3YHD78sHw/eDBw6qnx3p80SczjFYlI\nIcDMC5m5NzP3YuaezHxDyfwuZr67dD+WmQ9h5sOY+WvMPCfrgCelqUl3COR44YX4Y6ovv5xNWFRx\n9tnAGWfoDkV8nInpXXcVQxybNql1f599gD/8Ifj5xo3Bz8LQLTy9/l95JfDBB+X/adPx3/+OZ1+m\nkfSFLwAv+S17qWHsjmFD6NTJvzVIMRbePvecuvDoZtUqYMECca+rMluzRvjtzpdZs7LpTa5bp95N\n3UIgirDwxSn3MqxdK4S4DN68eP994JFH1IbHJKwQyJAJE4Bbb5Wzu2YNsGVLpVncDyFsSMEEli+X\nt3v66WJyEKisLD7+GPjkE6XBCqRTJ+Avf6k2V90TUE2cdDYNpwyrFmCbN4vr3LnAiBHx3r31VmDI\nELXhCYMZ+Pzz/PyreSGw006iVamDb38buOwyeftpWz9Jhw1MJEig9egBNDTk96GsWVNtZnoLu1cv\ncXXCSZTNXEYS3GXcLx3/9Kds/b/nHuD227P1Iy5NTaIsL1wo/p93HtCtW37+17wQ2LwZePNN3aFI\nzuuvy9t1Jt10MmWKmkrSnWdu95qagDfeEEsQO3dO708S8hYCcRsH69eLqzucpgiBKKLG7d1lnEiU\nh/feU+O3M6eWJn+TNOS6dgV+8QugZ0/x/4UX8p27rHkhAOhruXXpkt6NV1/1N7/zTuCWWyrN2rVL\n719aBgwAFi8G5s8HVq6M/77zjl8L3M2cOeo+/jD8yo5jtmVLucI1HW881q9P/13ceGNlI2XqVGBg\nDO1hYRWmX0/wtdfEXphZs8pm++wDHHmkGj+ffTb6/e3bs+mBvv9++V71fEgUhRYCw4fHl5idO+cn\nFOL64838sPfPOw+44grRtXXGO02BGTj8cOC00+K919SUbzc4KcxiiPHoo/0F78iRwNat4W4sWQIs\nWlTpZhSXXALMmBEvrEG0a1dernvbbcANN8R3Y+RI0RhxeOQR4B//CH/HXcbdyzG9OEMjbpz0WrGi\nbLZ1a7g7fv7Gxf1u167AT35S+Xz58vJmQicf16wxf47OodBC4O675dYKuz+wVauAp57KLkxJSCOU\nRowwd2dl3HgFCTM/d4I+6kWLKpchpiXInwED/JcSMosWsl8v6N57gWXLxH3PnsAhh8QLy5gxYs+A\nLO50c+6ff14MPQDAD34gelPXXAP8XIHy9x12ENff/z69W6ayenV1vu+9NzBsWKVZp07VPXU3t9wC\nXHed/zN3mZs0KVk441BoISCLd7WEaV145wP1VnayrRfTJyrz5JBDgB/+MFs/mINXCIWplxg2DPjN\nb7IJUxROa/nii4Hrry+b//3v6vxw1JZceqm4xhHecch6zD6J+2vXVru/enWw/V/8AhhVdTJLNYMG\nZd/Trwsh4N1UkteYm1OYXngB2HPPaHu1LAQaG5O/q6pHkQS/PJAJj9fO9OlqwuMwdy7QurV8GLp3\nF1dvK3bs2HThcPuR93dlkrt+786fL65btsTThOtNx6zTtfBCIEnG5T0n8Pzz4a0CVeH5/HOxBG7M\nGDXupcVdeA87LLk7WbUokxKWX0HPgobs4vZKnXj/+9/RFYussHLcTJum3qWXcRX/bdoUPZcChMfr\nsceAXXYJfh43jkTx5urcOJs3W7WKFtg6yV3v5Cef6E8QlWPGYTgVf1TBS9sTcOOcPXDxxfHfVU1W\nwvbjj8NVBjhdcxWo6gkEubvXXvHC4CyhVJm2WQ03xA1jjx7xNro56eL259lnxX6ZdevUrZZTkd9Z\nv5+G3HsCOiLbsmXl/z32yNd/WSHw4YfAR1XKt6MpwnBQGuLGL4kOpa1by5usoli5MnmaO+856/bj\n9gScdfIq8zzP/SUjRwY/U7nT2bv73uHdd+O5s2mT/wSv6uXJYcNqWX/fhRwOat06XsJ885vA229n\nF54oZIXAQQcBRxwh/57DP/+ZLFyWMuvXC+HxjW9El61HHw1+FtSriyLNh/7ee6LceMtBkhZsknCs\nWCFfVsePj+++g7fl71zdG8z8egdu84khx10FvevdaLd5c/4NySwppBCIe9wgUXnFgok4hW7TpmQH\nqDibZ0zrEegct0/K9OmVm4EuuURfWMJw5/W0aWKi+JRTKu288EJ8d91Lronklh/HbV0nJY6ATbKf\nwnE3qmeU5fr/7dvVDAvHoZBCAIhf4cVVO6sS2Z5A3Pe875smBMJgli/sRYqXG79hqXvvBR5/PB//\n3ZvRZPHqz/eqXJkzRz4/ss63sO/Du08jTlhOPlnEMw5JN4Z6zwyJoyZGFVqFQFNTcOJt3So+GC/e\nMVTZxHfvMtTNzJmV44x+m3oAeSHgtFzz1DyYluOPFz8TkVk4EDTm7OTfpZeKOQY/5YUqGiSqJivj\nDgf17StWu+nECWOY6vW0Auhb3wp+pnKvw9Sp1eY1rzbC/YF17Qo8/LC/veefL+/Ce+ih8k5LZyWD\nW0Oi6XjDOGqUUPngkHattVP5F6nF/PTT8ufP5h0vGU2WUY2KJ54QV7+hBRXxUZUmScqb7HBI0Kqj\ntI0VJ+5ufTtBdoiEwI4bz6RDXHH98TYS/HrINTcx7B3zdi/nGzKk3MJyJ+b3vgf88pfZhy0r4gwH\nue3KFiinZWlKT0C20Jq6Ga6pSTQ88iQsLfzKRJJxaZlVQCr33QSpcE97HkNYA9CvJyCz9yAOjh9J\nh24mTiwLSPfubUCE2zlMKS9yFwJhBfGRR8raI2UrTmdXXpFRtevSFCHgoPrjy4uxY0XDIytUCLWr\nror/Tpgum6Qw6zvEJiwdr7lGrXt+9ty6n555xv8AIj/8lOOFhSHrlY2RQoCIWhLRHCJaQEQLichX\n4wUR3U5ES4mokYgCV1zLJnQzSfGU1ylTafBW7F6VtV4h4NcbksE0IRA0bh6Fd5goj56A7mHFuBVQ\nkrX9jtZQN1Fr0mXClaXAdJN0vD+rcXbv9/aDH8i959c4ClvcceCB8cIVF5mD5rcAOI6ZDwPQC0B/\nIurjtkNE/QHsy8z7AxgO4M5ql+IRpD/Dm0j33ZfWp+yJuzrIKVxxC65pcwJJw+PWF1+L+KVLVmfY\nus9leOWV6ufeMpZUcEeRRneUF90C2/E/aaPLOwTkRscqRqn2NjM7K/NbQqia8BbjQQAeKNmdA6At\nEXVKEzBvRjvSs4jLIZMuEXXj9wEncSdPZFenOPs+nG5vUGvUu5zOS5qdr6oqv6Blr1FlIIsjUJmF\nSuM0qKpwk+iO8uand7d1nLDJlEXdwgUQ6vHzRkoIEFEzIloAYDWAp5jZq0m9CwD3ytymklkV3sw4\n/3z/j9ebIRMmyITUTKI23MhMDP/f/0X7E1c/fdZEtZSIgD//uaxLylk5FvTBHnxwuHuLF/v7IbNe\n/vDDo+0kwfRGS9YVXxq17UGT395DXdwk1fUUx16W6BjSlVIgx8yfAziMiHYFMJGIejCzzycXzX33\njcYzzzj/GgA04PPPgebN47kTJ8N0Z26UPqCglm85neTmPrIeO4xLVLp//rl/69572IZzxGbSVSXL\nlkULkDS6YNavB9q2lbevuzzmyfz5wHHHVZqpqpRlBVi7dsCvfiVnV4Ys8s9ZkloWAjMxevRM9R75\nEEuLKDN/TEQzAJwEwC0EmgB80fW/a8msih/9aDT69RP3117r54e4ejO4Q4c4ITWLsEJz222VxzC6\n4+0+oakIE+Be4m5WSvtx6erOt2sX3OKXDVOYKhS30jWi6CMcHYoubKLC/7e/ybslu3InLqrS2BkN\nKQuBBowe3eCqI30qS0XIrA76AhG1Ld3vBOB4AK95rE0CcGbJTl8A65hZum3lTsiNG/3tDB4c/I4M\nBxwQz35eXHxxWe+4g1/FITNmbdpHb1p4dOCXBn5mQb2cf/2r8j8RMHmynN8qVWq7USm4FyyQGwLp\n3bt6Kaq393bFFcHLYFWeJuj+Pr0K6ZI2RHRu+JSZE9gTwAwiagQwB8CTzDyZiIYT0Y8BgJknA3iL\niJYBuAvA+UGORX0UJ5/s/17aVp4OnRwOUbpiVO0QNKXS1TUObsLEXhBJJ36zPj866zSLcr93b7me\nzYIF0StnstgHEYVbg2kaHGGfp1pvh8jhIGZeCKC3j/ldnv8XynjYsaOfH+X7qI0Rpk+0JcHdEkrz\nUZqWJqr028ii+7QxouDhIPdpb6blUxpUrHwL2lSYdzol8c/bi0kbZmf+K09y3zG8887VZn7dy1rT\nKBlG3rpC8qJW4iGDe/5GJX6buXSnaxHPEs4qzXTnhQq0nSwWd1wx79N28kSV2gjT0iTv5W46ewKP\nPSZvN04+6R6uiQqrqhU8eROmdygKlXMisv5kiTYh4K4g4iRqLQ4HveaZZpedTExiJ0/qaTgor13O\nquOY1j3TylwehNVXpgq9MHI/aN6PsIIUNFFSS4XvN7+p/F+0ieGvf11sVHOW/iYNTxE/oCSk6QmY\nNtkus4kxqR/euKrwSzW1UA/lLgRkh4O8uoK8BSXNUIPJGVfE4aDZs8VJTs7GL5NXB336KdCyZfzN\niWlIozdHt2CMykO3bqKs/VY9vEgkeuFplo+bXJfIou1ksaTjaioqGJMzbt06s8Mng66KS8bf1q2B\n0aOzC4Nf3iXRmxPXjzxxDnjKGieeQXuHVLB6dbr3vYJJ5c7kvDBmYnj27HD7Di1a+Js7LFoUrT5B\nd+sqCr+4vfNO9IYxnZWDn99xW255reB4441s/ElD377+5ronhv2eqz70hLn6mEU3plWseU0MXyi1\n6D492oSA1+yEE+TtA8HqlmfPjt4YZooQWLpU3u6MGcDPfhZuR3cL0asxs54mhuPgF+fddpN714Q4\npl3ddMMNlf+XLAH69w92Q4WG16yWm6bRORXFuHHZue3G2OEg2aVrSSbOTPiQgOBdpEFxiBp/1SkE\n3GnqbNHXLZQc8sxvb5zT6tHX3RPYsEH9PNWNN0bbcZvLHjAVhkyvNEk8o5RDFgEjhoOWLatWoPXx\nx9X2/N5v1kxonbzjjmg/TSOocJsaXlluvVV3CMLJUyj4HYaepNHjJuqQewe/jZkyYXGzbp283aRE\nrQ5SkV8ybiSJX9G/VcCQnoCfTpCxY6vtuf+7hcCoUeLnZ98PU3oCfvrvw0i7cSdLku5tqHVky1qc\nMpmVMjSHqVOTHWQv634SO2lR0eAKOxayyBjRE/DrqnnHlYOWizn2nG6Z296JJ6YLa9YEHY7hxGGX\nXSrNn3kGuDPk4E7TCqRMF1zVbmlATJ7LkGWlo1oFSNAxq1nSvz/w5JPi3oSGhYo4q4iHc6JZlFnR\nMGZiWNa+d8ext4C4N5dNm+bvpik9gSCCCv/KlcB550W/pwOidBubtm1LH/6BA9O9nxdpJ1ZV9i66\n+J7/ly9RrfS//jWfcCRZKbVkSTZhyRMjhoPiqE/1CoGjj658PmVKunBlyebNukOQHWnHudMKgDff\nzHY9eVJk0yBIk2YaZNL0O99R729cotJI1fnPafFb9fbss3rCohIjhoP8TlXy2vMWaEdwHH98pbnT\nPTOxQpDVPa5CbYTO3k7Syj+NILjLpdj8W98Kt5tn2sieeRs0Bq9rdVCcZb5ZhcHx+4wz0vuhcrmy\nacOuaTGiJ+A9WQsot/S9yyK9PQE369eXn/utyigKjrBwVkjJortwpqmwjjlGXThMQvWcQNqhJBXk\nVc5Uqh8JWl5tQnrqxoiegDcjDjusbOZUiN6Wid9xbG+9Vf4fdiav7soyin32Sfae6fHywwnziy8W\nVxOnG5k8cB8wE+e9uKg+lzpJGIukCyvJnEAtYOTEcGNj2czZ/Ru0OsjNli3hPQWHHXaQC6tqmNPr\nKgGCVyQUUQi4lweH7fWoV9JMDMswdGj487RlSub9u+8Of25KuchziCxPjBgO8quwnedeZVVhlXwR\nxvT69EnvxocfpncjCYMGAWefHfzcCVecj+Wpp9KHy+2nH1lMusqS5KASFe5liervx3uAvC5k9+HI\nbtYrCpFCgIi6EtF0IlpERAuJ6CIfO8cS0Toiml/6/TzIPb/hIL+KYNgwcQ06w9NvOMj9f+bMwChp\nRcU28xdf9DfPWrhNmgRMmBD8fO5ccc3zI/nVr6IrRq8iLt3DQX7kfQpbXhRpOKhekekJbANwKTMf\nDOBIABcQkZ+ezmeZuXfpd32Uo1EZ3Lq1uAZV9lEfzYgRyf02HZ3nLzMDt9wSPqz10kvimker9Ze/\njLaTpc57VUSt0skbVUMfv/td+rBkiTt9Tehl6SBSCDDzamZuLN1vBLAEgN8WE6kklJ3xd9THPv+8\n//vOElG3O507+7vb1CQTsnxQMVn329/6m+cl3K64Avjb36Lt6RK2RfyYg8Lsp05Yd/zi5OukSdmF\nIwqZXnfRG4QqiDUnQETdAfQCMMfn8ZFE1EhETxBRjyA3nnlGzi9npYijZ9wrPPy0MwbtPLzsssr/\nRc/4F17QHYJiY+JwkO6KPQ55LxFNyqZNasJR60gfL0lEbQBMADCi1CNwMw9AN2b+lIj6A5gI4Et+\n7lx88WisXetsEGso/eLjjIu7C4p7w5CJZP3x+Lk/aZKY0FXt9xNPAJdcEu53kSq2orBoEbDXXvn5\nd//94rphQ3mIVvWhMhY/ZpZ+2SMlBIioBYQA+DMzP+597hYKzDyFiMYRUQdm9umQjcbo0cAHHwA3\n3ywf0KgdxABw/fXAscfKu5UlssrMVBIkBLLg6aej7ejemWsCqsO1ciXQvbtaN/0YMEBcHU23mzeX\nVVM7cz6WLGlAZQP52sx8kh0O+hOAxcx8m99DIurkuu8DgPwFQJm0FbGzrEzGHR0VxIMP5u+nX1oE\nHVyjGln1CPWGTYPk5LFHwSLREyCiowCcAWAhES0AwACuAbAXAGbmuwEMIaLzAHwGYBOA06PcjZtB\n3h5A0OqUww+vnndQcTJRXNq0yd9PL6tXlxXqEYn18lltlIurHbZeePRROXum9mTc5J2fq1YBe+yR\nr59h1Gp5jhQCzDwbQPMIO2MBjJX1NE1iBukScvA7ZL65J/R5ZKbXzzzwxsu7SWr7dn27pU0izwo3\n7AB1N3HKpI7KSIeQuvXW7HvUtVqxx0HLjuE5c5L3BP75T3EN2nk7fny1mY6egA68aeqdH1GxIclU\nhVumrPdWoQE2iiL0GlRQqxvoTENL9ZjkABHvcFC3bpX/HZxdq26OOMLfrVrDGy/vdvxajTcQrpBQ\nJyaFJQ1E+cfFzgnkgxYhMHw4MG9evHecDHUOvu7UKdiulyj98nmhY4lonv4XhTxb0lkcP2h7AhaV\naBECixcDf/lLsnedvQDOcjWZii1Iv1CW6PhQX301/Hktf1SmVIz1sFksL0xTJV2rGKFFNIl926qt\nJkobY9G681n5W68fe1qKVn5sPsuhTQjEbZUGKZCTKSgPPxzuVl7Y4aBw9t47+bsym9fywPQ0LhJ5\nzAmYOpeUJ9qEQNqjE53/Mkqinnwynl9FJaoQmz4clKbl5rcgQAe1OtSoA9PKa60KCW1CQHb9tEOQ\nELjNdw9zPLfyIuuPt+hCIK98MVGBXBryyFfm/OOWR7xMWVqsk8KsoC/anIBfgdJdCaso5EHpblJ+\nRIXlvvuAU0/NJSi5sHJl+T6rfPAOqVpqh8IKgTQVqq4Kq4vfKQwKMW1OIK7QybMl9thj2bibdKe4\n6a3QjRvL5ad9e71hkcWkhkkUjoZWHdSVEMjzQ9NxmpVOIZD1gehFIWmcTa+w3OE76KD8/SyyHzIc\ncog+vwsnBA4+uPK/l6FD5d3KEpljD2sJFQrkVAkNUz7srAhLp1qKe9q4LFqkJhy1TmGEgMPPS0fY\nBxWQ/fbLLyymYdpwUC1VSLLUapx1TAynJewc7Dw47DC9/stSGCHgFEBnzLWIcwJZY5oQ0EW9xNOP\nrOJez2nqEDcN4mjs1Zm+hRMC555b+d9LWFe51seo8yhI27f7m6uYE6jniiZOWulOp7xWudWTAjkr\nBCRwEmn9+sr/ljJ59AQ2b87ObVXoDItJ6aCSIg4HxWWj9+T0HLFCQALvkr40iVarhbnow0G13lPL\ng1oaDsp7ddCoUdm5bTKFEQJedG+8SouzykkVl16q1r0iU0RhUpQKI0/yThPVar/jhN/2BBKQZp+A\nCR+c6opq27bgoRqHelkdVMThIFVhzjLupm829GLCd14EIoUAEXUloulEtIiIFhLRRQH2bieipUTU\nSES91AcVOPLI8n3RM1j1kZe33w68+GK4naKnmSy1Hk8dPR0d2jaLno9FmeyPPGgewDYAlzJzIxG1\nATCPiKYx82uOBSLqD2BfZt6fiI4AcCeAvqoD667kghLthhui3TGhcOn+kC3FpV56WEUvrzUzHMTM\nq5m5sXS/EcASAF4tOIMAPFCyMwdAWyKKcQBkfIKGgz79NEtf1VHEceu41EMcvehWdVD0ilMlqtMi\nrns1IwTcEFF3AL0AzPE86gLApcsQTagWFEpJkmgmzgncd19+ftZLS86E/K017HBQtpg+HAQAKA0F\nTQAwotQjSMho131D6Refoq8OcoRAmzb5+Wn6R6UqfE88IdLX9Piqwt3jymNiuF7SVS8zS7/skRIC\nRNQCQgD8mZkf97HSBOCLrv9dS2Y+jI4VwCCKvk9gxx3FtR6HTIIoumAH6mM9fV6Y8J2mId1wUAMq\nG8jXpg1OILLDQX8CsJiZg87xmgTgTAAgor4A1jHzewrCF0iSAmJSoWrZUlzzFAImxd9SLNxlJ68y\nW0/l1eg5ASI6CsAZAPoR0QIimk9EJxHRcCL6MQAw82QAbxHRMgB3ATg/01AjWatx3Digc2czCpfz\nIdWyEKhH3UE6WuV5DAe51UbUQj5ZykQOBzHzbACR5yUx84VKQiRJkoL49tvAccepD0sSdAgB04nK\n02OPBZ55Jp+wFIl6mfAvGjW5OihrBgyQt5ukJzBvnriaUIh1rFTK8lzdPOLRQnoZQ/2SVT5cd102\n7oaRhxZRE+oC3RglBHr2lLfrHC4TByJzWt7Tp+fv54IF+fsZh1r4IJPGIc6RjV4/8jqa0PHXlG8o\nirh54dZIoMK9ODhzhDowSghkXbhM2icgS1hYTWoZq8i7IuWLag48UN6uN52SHm5vOnmXh8GD1boX\nJ/wDB6r1Ow5ahYC35Z+HEChKK0aGTz4BzjlHdygEKj7YnXZK70ZRiZN+YZV+LSmQy/tQGZ07jFXr\nEouDViHwyiv5+meiAEgTph13BPbcs9p8jnc/tybiflQ6u8SqyKOi9JaZvMt1XsIgbbyKtDrNTgyX\nqOfhoI4dgf32i/+eX5rtskv68JiIifmmA286FK2Fnpc/MvWJLVN1JgTy8iMOrVuL6ymnAG3b+tt5\n803gwQe1JN1eAAAe/klEQVT9n6k427co1PIHa88TUO+f7vKi239ZjBICWeOMu5mUOd27i2vYmODe\newPt2/s/86vwdY4vuokrjGpBeOkoW7WQblmge06gKBhSXeSDiRPDrVqJq0qFZ6bFUZao+Bc1Xlmj\nQ8NnHhSt51FUv7ULgYULy/f1PCegkk2bdIdAUOvprJIi7C6txfzM8vSvIuQpYIAQ2GEHOXtxdhMH\nYWJPIIsVEKZo46zFSiMK3bqDaglm4Kyzkr/fKdNjrWoHrULg5JOjD0d36No1vX9vvCGueXyoxx4r\nZ082LEEf+le+Ivd+EahHoaGaWlIrzZzu0CW/5dN+ftQ7WoUAEbBlS+X/IFRk1oYNZreaZAWHm699\nrdrMlInhuGzdGv7cfrACXWW41rWI6twsVrfDQd7hmaxP2dq2TVxNbC01bw7stpsav/OqJCZPDp9n\niSuMli5NHybd6K4gTSzbutAdTt3+y6JdCDgnbAHAiBHBdv/4x2j3nJU2QfTta15PwCkoSfW/6Fwi\n6mhlDcK0tDaZqLLrpl4mhmvdP1PQLgTcStDcAiELvv99cTUps52wRE1ixalQdVS+efhpUr6ppksX\nebvudFiyJHoYTTVWuKunroeDVEY+yq1jjjH3aLydd1YXNvuR6iNpeU76Xo8ewK23pnfHRGopLiaj\nXQjowKTClUVY8lItXKublEzHm9bvZXqad7W/RdEdFPdQmSj7tVrGtQqBZs3yTdg89wnIxsuZDE4a\nLpPVRliywVu2li0LfpaH/1mR9lvVXWnXzOogIrqHiN4jIl/Fz0R0LBGtKx1AP5+IpM/8Ulkhjxwp\nn5C6C4cbZ0WUn4CaMSP/8CTFLy91LrnTRS1Vwrop+sRwzQgBAPcCODHCzrPM3Lv0u17Wc5VCYNQo\nf3P3JjMTewIOREC/fpVmDQ2Vz4PekzFzUHmGgx0OMo8lS7JzO+/8LroW0aIQKQSYeRaAtRHWElWt\ntaw6OkkBPOKI4GdBq4fiCoGrr44XJpOwH3U0TzyRjz/vv5+PP2mw5UUOVaPHRxJRIxE9QUQ9ZF+S\nXR3kLO2Mws+tPIYpZMMSRpSw6tULmDmz2rxNG2D27Eqz3XcPdidPvUKmrFLaeef8/KqHiqcoE8Om\n+2cKKo4qnwegGzN/SkT9AUwE8KVg66P/e7dmTQOAhkgPgiqTZs3KlVqQQHFPkpo8HCRD0IYir+qI\nsNVBKoVAVFqqToOkeXe99AClPlSNH+d1qIwpAl4lZs0JzCz9sie1EGDmja77KUQ0jog6MPNH/m+M\n/u/dHnuk81umQitKT0AlYR/o6tXq/ClKSy3ObtwiYEKLNQ8h8Omn6d6Pe8yq6nT97LM0bzegsoF8\nbRrHQpEdDiIEjPsTUSfXfR8AFCwAPJ77LBE98EA/P+QCWeSeQF7noWY5cVhkDjhAjTsmVNBZ4Y5b\nHsuQGxvTvR/3W1+zJvx53Lz9z3/i2ddFZE+AiP4KIZJ2I6IVAEYB2BEAM/PdAIYQ0XkAPgOwCcDp\nsp67M8kZ737uueoxbZnMJPLvGdRST0DGzc6dw9MrqzmBoqZpEdE1HORGZWOqXTtg3Tp17iVFZ/nS\n6XekEGDm0GlZZh4LYGwSz+fNK0feUaOs+uB0r6rqIvcEVGDKgTN5Umvj17qFwCuvAJMmqXNv2LBK\n1ReqiJsWOr8N0/cJZMbixdVmaRLD711HfbQqP9KEJQ83dQk6EytaHR+W7rKV18TwrFnZ+ZMntndp\nwPGSXmSXeUbZufTSanMTj5dUTa3HLytsZRAPq5qkdjAuK1VVYo6K6qLMCaisvHUIAr/46lYbUaSK\nvUhhBYpRXotUXup2OAiQi3yaVkdTU/k+z55A0rMRHnwQ6N3b/5lpFYXtdVSjO4/yGg7SHc8sqMU4\nyaBdCHjVISSZGO7YMV6FlEdm77CD/w7fIJzwn3FGumMmTZoTUB2OuO7FUX1cpApANqxDhqj11338\np8r0yirtFy5UG44ilZE4aBcCnTunT9w4+tTzrCBV+2VyITRxiK1W0TUxPGFCNu5mFeaswltraBcC\nXrJsUYYdiq6aLISADCYN0ahO57hL+Gp1dZAsWc7JmBTPIOL2AHXGyW8VY14YJwT8xv9ll0bKkHdP\n4MIL5ey7w5W2MJq6F8JkaikuDron5mXZtCkbd+MSNY+3fXt2ft9yS3ZuR2GcEGjXLtrOySdH2xk2\nDPjFLyrN8mwlO0LgppuSvRvH3A2RWBk1fXp8f1WjezhMRaXVvXs+fmahQE53+suyeXM27sYlT22z\nJmGcEACqDz7xFuaNGxHJAQcA111Xbb5uXT5nsuoaDnL7r5uitERN8zMKWSGQ10qhKLJsQaskarjR\nxLKgAiOFwEEHhT9/881qM9kNZQ8+CNxwQ7JwxcERAiqHqeIMi5lcYJNqhyzSuu8skY2XKXMCLSKU\n05iST3Z1kEF4K0Rv4jvP27f3f/9nP1Mfprik6QkkKWyDBiXzSwfr1yd7L6lulzyXiCZ1R1Xv1JTW\nvwnuyqByLi6ItCrzs8ZIIeAlKHPWRh166SHN8Awz8Oqr8f3Lqydw1VWV/9PsNfBj6lS17uVBlGpg\noHIzoR95VVDnnSdv14QWqUmr0NIQp3eTNN1l5jl1ouJkMeVEFbDDDwfeflvcjxhR/U4WBXTWLOCY\nY/Lpike927FjdW/HibNz7dUruf9+9O9fbeZOZ1lVwDvumM/kKQCsWlW+v/Za4OOPK583NQFdu5pR\nqfqRpBxnOSdgck8g7HzuKD76CNiwAbjjjnB7ppaTtBSiJ+Bl333L92PGVD8PUjPh/qjOPbdcKaxa\nBdx/P/D73ws748YB55xT+a6zguF73wOmTRP3CxaUj7V0Jr+2bAHuvBN46aXkPYGgwnb44eJ6443A\nRRdFuyMLs386Xnkl8Mkn0a1lQH6FR9jRl1Ek/QhXrhSrtLzqit1LE3UPBxWBrASMCrfmzk3u7mmn\nya0Cq9W8NVIIpG3Jy+gaGj8eePllcd+vHzB0aFnz6OWXA/fc4//eQw8BJ54o7pcvF9f77itPfrVp\nU+7aq+6RtGkjrkk1rQaxYQNwySXV5jffLFrQXbv6vxflpylLFG+80X8yuqgfdVi4f/97/WGIy+TJ\n6tyKCzPw4YfZ+mH60JmRQsBL3AJ3yin+5n6ZsX07sGxZpZnf5pUpU4LdGzZMXF97rXLnX1BPIGpo\n5brr4h+ykVVBizvv4sbJN/dcyqZNyfcwxC0H//ynWveyxt2bCqvMdYV7wYLy/b33qnP3I6nDaOMj\nszT1yivldQwl2c/wwQfmHzNppBCIWh0kM2Yu689ddwUXlhkzykMhfh+lN5zepa1BFfODD1abucN8\n9NH+LXMZ/JbPRjFjhrh26AA8+mjls7DJ1agW1Nq1Ynjsy1+uNP/hD+OHEYhf+f3qV+HPnRU5s2cD\nK1ak95MondDcvFm4sXVreOtY1xLRRYvUupc1Mi38OHHae+/4YfjWt0R+AmYcoemHkUIgLXFaxWGV\nXL9+QN++/s+8k4xB4ZDRhfTBB8DAgdHuRfkly6ZNwr5zsttpp4nr2rXi/rbbynbDjhH83/+N9qtV\nK/lwhfHhh+kPHvdyzDHi+vWvq9Pd8re/JX/XaYy0bAn861/B9twT3mGoFgKmqHeQRac+Hgd3L3jg\nQDNPZIsUAkR0DxG9R0SvhNi5nYiWElEjESlZl5LmvE+ZiWHnf9SQwTvvVE8SA0DbtsDgweHvBlXM\nXvPddks/nBPnfWd7vNMK8vaELr7Yf6JNJ1/4gu4QyDF+fPJ3ZRoWgHyL8oknkoelFpg/X3cIKvPq\nuedEL980ZHoC9wI4MeghEfUHsC8z7w9gOIA7VQQszSYO2Qpx3Tpx2H0UQZPEUbgPuXeTtsJX1cK7\n/HKxismPNEvuaok8x9/32Uete0VR12DRS6QQYOZZAMJGOgcBeKBkdw6AtkTUKcR+bPxa8GHI9gRk\nFNGlYcWKbCZsVanbXrEC6NMnfXgsFktxUTEn0AXAStf/ppKZMrytsT33DLdvyg49Z0LIi+lLxixl\nbGvaUuto2DE82nXfUPrF48gjw58nPd9XNVu35nPsYlZuWoB339UdAkt9MrP0yx4VQqAJwBdd/7uW\nzAIYHduDuMNBphC0O9YKAYvFEk4DKhvI12bmk+xwEJV+fkwCcCYAEFFfAOuYOVAn4i67xApfas46\nq3yf9yabvn3z6wmEraZKsr7ZYrHUBzJLRP8K4HkAXyKiFUR0FhENJ6IfAwAzTwbwFhEtA3AXgPPD\n3Iu7EzYt7nXqebeWZZeIqiBsSaftJVgsliAih4OY+fsSdiRP0k23/t8hSaW2fHn+cwVxFMjFddeL\naSoQLBZLMch9x3ASIfD445X/41Ssjt299orvb1ry7Ak47L9/5X9TVkpZLBYzyV0IDB4cfyPShg3Z\nhEUXRGVVDapw9kZ490i0aCH0HsmoeLDUN1//uu4QWHSQuxDYY49oxV5esjygJUtkN62pwHuojNt8\n4EDgQukBO4vFUk8UQoGcdwipKBOdOoaD/v737Ny26CGNUro4FOW7sqilEELAS5I5AR1kJQT8Jrgd\nN736Z95/P51fFv3IHJJksSRFS/HKs2KutVUzc+cCp58e/DwsbU3ZSR3FI4/oDoFZ2Ba6JUuMPGg+\niqJ8FFn0BL761XA3w9xWsTw3D049VXcIzCKv8l6U78qiFqM7mrffXj73101RCmuecwK1JAQslRSl\nvFuKidFC4Kc/BX73u0qz3XcPPvjcNJyP1zsklcUmMisELGmxwqY+MVoI+LFgAbDTTuK+Q4do+zrn\nBIL8zuJjcyYP7YecHbpON7O6nyxZUjgh0L59uaL74hfD7ZqKruEgSzreC1SLmC177KHHX0t9ULjV\nQbaSC8emjyUptuzUJ4UQAocdVvmuLazVyAwHffvb+YQlS3beWXcILJbaohDDQV/5Svk+rgDQOScQ\nFta8JoYPOqh8f/jhav3UgbMZTodiPF2ND7tEtL4ZODBb97UKgX794r/DnG78u1Yn2RxhF3YKWy2s\nDnLGx1euDLdXBHbZBejVS3coLPWO1uGgqVPl7KtcZXPuufHfKTLuNKqF3dNjxgCrVgFt2mTv10sv\nZev+qFFitdv48eH26rGF7lWJXjRuukmdW1l/t1p7AjvsEP8d9wfRuXP896+4Iv47WaD6w5b5aJyl\ntUWmVStgzz3z8cs9DJkF7dtn635c8hA28+bJ2bvuumzDkSXMQM+e0fZM+R4LMSfgxj0xLLNE1CtF\nWyhSlLF0qbj+9Kfy73zjG2r89iNIILo/7Kw32T31VLbuA9EV1de+ln0YVBG0mTBLwnrfeeiWkp3L\n+epXgSefVO+/ijPOO3ZM7wZgTs9cSggQ0UlE9BoRvUFEV/k8P5aI1hHR/NLv5+HuJQ1uZUHNQhun\nLPvtJzLx9tuBtm1Fq+6hh8Lfad06uX9RBI33u9Ood+9qLaMy/Oxnlf/HjgWGDhX399/v71cQl10W\n3/843Htvtu7rYvfd073fu7e4nnhi2cydX9/8JvDAA8Dw4fJutmxZ+f8f/0gePqCsF2vTJmDffYET\nTkju1n77+ZsnGWZau7by/zHHVP7v37/6HfeKRkCkr0NUPeB9V/twEBE1A3AHgBMBHAzge0R0oI/V\nZ5m5d+l3veJwesIkrs2bJ3vfqQjdvYLvR56kHMzq1cCKFcCQIeH2ssxMGbf33Rf4z3/K/wcPTuZ2\nx47AUUeJ+06dyubuSqVtW3+3ooZy/D4oN1GCZv/9y70B73JSd0U6YEC4OzJ4K4OkROUdUXVZj7vS\na9iwajP3qpOePUX6yJSjt94SV6/APfjg8PeinjvDw61aRYehRYvw9Hd66l6efRbo27fSbOTI6jAA\nwHHHiWtQ7+VPfwJmzQIefbT6WceO4ntz+J//Kd87y7m9ae0slLn++nxPApTpCfQBsJSZ32bmzwA8\nBGCQjz3pdrmqsdCketaditBdgUWt0hg5UmS6H61aVU5Udu8eHQbV468yPQGHtm1Fa+Sxx+TcdrdM\nDj1UHI3puBvUWtu+3d88qpI55xxxXbAAOP98cS/bhd9rLxGuJ58EGhuBJUuAt98Wz84+G1izprw6\nzB2+734X+Owz/zkAt1Aiqgy/u1f1y1/KhdEtmIKOc7zxxmozb7r9+9+V/7t0Edeg4cnzzwc++aTS\nzC9dr7sOmDixsofnpXt34Ne/Bk46yT+MF1zg/16bNv75f6Bfk9KD05Nx2HHH6DnFn/uMR7RuDdxy\ni7h3jrl1N94eeqjcOBrkV8uhHIezzhKNIUdonXpqZfxeew24+WZx7250ON/Od75TNhs3rhzHAQPy\nHSqSqUa7AHAvyHunZOblSCJqJKIniKhHmIOHHhrvsJNLLqn8H6cnEJaYHTuWWwGXX1793F2Btmol\nMh0QlUsYTkspjLyEgB/r1gEbNwY//+53xQllZ58t/g8ZUk7HSy4RYe/fX3zs7njILEf1yw93S8t5\n3qtXuet+zDGie3/ooZWC2+GII4DXXweee078b9NG2O3WTfzGji0vCHDc37at/P5++4mWpdNSdX+w\nkyeL6267VfrptBId3PsxwjjhBDGECJTjd/rpYphmv/1EJXTllcC0aeV3vGXFr5f1pS+Jq3eI88tf\nLrshs9GuUydR+Z15JnDKKWVz72qskSNFY+7kk8tmTtp68/iDD8T15pv989+JXxe/WqXE9Onl+0mT\nRAPGr4wde2x57i3qGNtrrxVXRycUs6jIH3tM3Lvn0NyNnagejUOLFuU6ysnrn/4UuOYace9uDBMB\nI0YA99wj/rvjFtU7TouqieF5ALoxcy+IoaOJUS/EUcbVI0CkpD1xaeedRddr0yb/Svmjj4DlyyvN\nmINbLs2aARddFOxfhw7l1rfz0apCxTJa947iIUPEAfXuow2ZgR/9SNx37gzccUe1G6+/LtLtlFNE\nxeZt1YZVAu3blystoPwhTJwILFokWvd+Fdkuu4j0DFoocP75wAEHiPujjhJ23T0Bp+d2553lCsuL\nO1yAqBTcCgzd8Xr99coKxBFcl1wiDsw54wzg7rvLz3fdVUzYLl0qBAIRcPzxlWPRjvubNgkh7uC0\nxpOcZ73rrsHPAOC++4CXXwZeeQX4zW/87Tz8cPneyS9vJekIUPc37zeM6B4+8dK2bXlu6pRTRPr7\nra6ZOhX417+C3XFz4okiXbt1A5qawu0eemj5Pmh/U1iD04njmDHAVVeJXqf73AwiUWacYTu3W06P\nOCuII/odRNQXwGhmPqn0/2oAzMy/DXnnLQCHM/NHHnMGRrlMGko/i8VisZSZWfo5XAtmzmYRLzOH\n/gA0B7AMwF4AdgTQCOAgj51Orvs+AJYHuMVJEbJR3G/aJO6vvDLYjsO551abBdGnD/Nll1W7M306\n84YNiYP+37ANHZrOjTAWL5aPpxuA+dZbxXXYMHE944z4bgDMM2ZE27n7bnFduZJ55Ehxf/TR/mGf\nOjU8TjfdJJ5/85vxwsvM/OqrzI8/LmcXYG5oCH4GML/7LvNZZ5XNu3cvP+vZU1wvvzx+OAHmDz4Q\n6faTn1Q+mzKF+T//EXaGDxfXK64Q18mTw92cMkXcr13LvH498/bt0WHZvp152zZ/9wDmVauiy+AH\nH1TaAZjfeUfcT5zI3KpVtduTJon766/3d//ll5l/97twf3/wg/K7s2bJfSsTJpTtXX65f/3iDufg\nwdXmzrfl2Pn88+r3AOa77qo0/+MfK/0r1Z2R9XWSX+SACjNvB3AhgGkAFgF4iJmXENFwIvpxydoQ\nInqViBYAGAMg5BTc9DjjbDJDHXEmWObMEZNG3tUJxx2Xzw7VNBx4YLJ1+pMni2WBI0cC550nzJLO\nV8i8N2yYmKzt2lVMLjKLOQg/nO56EGk2/h18sBqdLG++KXYW77FH5cKBQw4Rwx8jRwJXX53c/cZG\nMZxy7rnAH/5Q+eykk8Tk9Lp15XkGQKSZ7Dhyu3ZiWEhmaLVZs/B5uD33BLZuDXdjt92C83TQIDHc\n5SWqXPXs6X8CoUrcK4iCiBpe85ujcOaWvHE866zgiWnVSG2dYuapAA7wmN3luh8LYKzaoFVy0UXl\ngu6sCvArHLLjgWH06FE57loEiCrXIsviVBa//nWlW3HYsEFUUu5xUz+OOEJUIt26VZrnuRIiKUFp\nsvfe/vqoJkwQ8w7OHMb48ZWTqLJEpSlQOVEsu5Epq93BSbQAhDFxYuXeBhUkCWOHDsCyZcFqZ5Yu\n9V+04MZvtdL06f6akZs3D58jUUlhDpq/8kpg/fpKM7+C7N0MkmR333PP1YaytSQMGRJ/z0SbNpUT\nyEH88Y/+5ml2cV9zjbr1+kHMnRu+csUP70aqp59WF54g3n1X3+lnSXn66XD1L1m0hr/61eoltjLs\nu2/lKiU3QZvTTjstetIZSL7nSQWFEQJduojVCm789HN4W5WjRgWvWw4iK531RdCF//e/Z+d20Mc+\ndKj8sjsvN9yQODjSODtZTSfOCWSmKKVLokk4CX37it4ZIOKel1r1bt3K+xKCePrp8uZLN2eeGT28\npoLCCAEvskMIO+6YTNGcahYtqh4GqSfC8qtly+CNUxaLCi64IH5jMC+CBOGhh+azc7iwQqBoBO11\nsFh0YEpPwGS8Q3q1SuG0iEah48Qpi6VIjBgB9Omj1s0gXT1FZsCAZHMHRaOmegJFWGVisehmzBj1\nbpow5KqaZs1q40jWKGpKCFgslvyxja9iY4WAxWIpBP37AzNm6A5F7RGpO0ipZ0Scp38Wi8VSCxBR\nZrqDam5i2GKxWCzyWCFgsVgsdYwVAhaLxVLHWCFgsVgsdYwVAhaLxVLHWCFgsVgsdYwVAhaLxVLH\nWCFgsVgsdYwVAhaLxVLHWCFgsVgsdYyUECCik4joNSJ6g4iuCrBzOxEtJaJGIuqlNpgWi8ViyYJI\nIUBEzQDcAeBEAAcD+B4RHeix0x/Avsy8P4DhAO7MIKzGM3PmTN1ByBQbv+JSy3EDaj9+WSLTE+gD\nYCkzv83MnwF4CID3+OdBAB4AAGaeA6AtEXVSGtICUOsF0cavuNRy3IDaj1+WyAiBLgBWuv6/UzIL\ns9PkY8disVgshmEnhi0Wi6WOiTxPgIj6AhjNzCeV/l8NgJn5ty47dwKYwcwPl/6/BuBYZn7P45Y9\nTMBisVgSkNV5AjIni70EYD8i2gvAuwC+C+B7HjuTAFwA4OGS0FjnFQBAdpGwWCwWSzIihQAzbyei\nCwFMgxg+uoeZlxDRcPGY72bmyUQ0gIiWAfgEwFnZBttisVgsKsj1eEmLxWKxmEVuE8MyG85MhIiW\nE9HLRLSAiOaWzNoT0TQiep2IniSiti77I0ub5pYQ0Qku895E9Eop/mN0xKUUjnuI6D0iesVlpiw+\nRLQjET1UeucFIuqWX+wC4zeKiN4hovml30muZ4WJHxF1JaLpRLSIiBYS0UUl85rIP5/4/bRkXiv5\n15KI5pTqkoVENKpkrjf/mDnzH4SwWQZgLwA7AGgEcGAefisI+5sA2nvMfgvgytL9VQBuLN33ALAA\nYpiteynOTm9rDoCvlu4nAzhRU3y+DqAXgFeyiA+A8wCMK92fDuAhA+I3CsClPnYPKlL8AOwBoFfp\nvg2A1wEcWCv5FxK/msi/kp87l67NAbwIsQ9La/7l1ROQ2XBmKoTqHtMgAPeX7u8HMLh0PxAi0bcx\n83IASwH0IaI9AOzCzC+V7D3geidXmHkWgLUeY5Xxcbs1AcA3lEcihID4ASIfvQxCgeLHzKuZubF0\nvxHAEgBdUSP5FxA/Z79R4fMPAJj509JtS4jKnaE5//ISAjIbzkyFATxFRC8R0Tkls05cWv3EzKsB\ndCyZB22a6wIRZwfT4t9RYXz++w4zbwewjog6ZBd0aS4koddqvKu7Xdj4EVF3iB7Pi1BbHk2L35yS\nUU3kHxE1I6IFAFYDeKpUkWvNP7tZLJqjmLk3gAEALiCioyEEg5tam11XGR8TlgWPA7APM/eC+Ph+\np9Dt3ONHRG0gWnkjSi3mLMujCfGrmfxj5s+Z+TCIHlwfIjoYmvMvLyHQBMA9QdG1ZGY8zPxu6fo+\ngIkQQ1vvUUk3UqlrtqZkvQnAF12vO/EMMjcFlfH57zMiag5gV2b+KLugR8PM73NpkBTAHyHyEChg\n/IioBUQF+WdmfrxkXDP55xe/Wso/B2b+GMBMACdBc/7lJQT+u+GMiHaE2HA2KSe/E0NEO5daJSCi\n1gBOALAQIuxDS9Z+BMD5GCcB+G5phn5vAPsBmFvq4q0noj5ERADOdL2jA0JlC0FlfCaV3ACAbwOY\nnlksgqmIX+nDcjgVwKul+yLG708AFjPzbS6zWsq/qvjVSv4R0RecoSwi2gnA8RDzHnrzL8dZ8ZMg\nZvuXArg6L39ThnlviJVMCyAq/6tL5h0A/KsUn2kA2rneGQkxi78EwAku88NLbiwFcJvGOP0VwCoA\nWwCsgNjY115VfCAmvP6vZP4igO4GxO8BAK+U8nIixBhs4eIH4CgA211lcn7pu1JWHg2NX63k35dL\ncWosxednJXOt+Wc3i1ksFksdYyeGLRaLpY6xQsBisVjqGCsELBaLpY6xQsBisVjqGCsELBaLpY6x\nQsBisVjqGCsELBaLpY6xQsBisVjqmP8HXIZQypbLqLAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff0ffb97400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dis = np.array([pbcdist(frame[0],frame[1]) for frame in t.ptraj])\n",
    "plt.plot(dis)\n",
    "plt.axhline(y=.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8VXWd//HXGxHRMEBKEPCCtwRzQksma5pOF0XN1GxS\nMlPTJictnalpkvoV0PSYstvoNOnM5A2dfAgzzYgWKpKeaZwKSTExEClFDpeDiEoSqFw+vz+++3i2\nx3PZ55y991rr7Pfz8dgPNt+z9trfszas9/7e1lJEYGZmjWtQ1hUwM7NsOQjMzBqcg8DMrME5CMzM\nGpyDwMyswTkIzMwaXI9BIGkPSYskLZG0VNKMUvkMSWskPVR6nFj2mumSVkpaLumEsvJjJD0i6XFJ\nV9bmVzIzs95QJesIJO0VEVsl7Qb8H3ApcBLwQkR8r8O2E4FbgGOB8cBC4LCICEmLgM9ExGJJ84Gr\nIuLu6v5KZmbWGxV1DUXE1tLTPYDBQFt6qJPNTwNujYgdEbEKWAlMkTQG2DsiFpe2uwk4va8VNzOz\n6qgoCCQNkrQEaAXuKTuZf0bSw5KulTS8VDYOaCl7+dpS2ThgTVn5mlKZmZllqNIWwa6IOJrU1TNF\n0iTgauDgiJhMCojv1q6aZmZWK4N7s3FE/EFSM3Bih7GBHwJ3lJ6vBfYv+9n4UllX5a8hyRdAMjPr\ng4jorMu+W5XMGnpDW7ePpD2B44HHSn3+bc4AHi09vx2YJmmIpAnAocADEdEKbJY0RZKAc4F53fwy\nhX3MmDEj8zo0Yt1d/+wfrn+2j76qpEWwHzBb0iBScMyJiPmSbpI0GdgFrAIuKp3Al0maCywDtgMX\nR3sNLwFuBIYC8yPirj7X3MzMqqLHIIiIpcAxnZSf281rvgF8o5PyB4GjellHMzOrIa8sroGmpqas\nq9BnRa47uP5Zc/2LqaIFZfUmKfJYLzOzPJNE1GKw2MzMBjYHgZlZg3MQmJk1OAeBmVmDcxCYmTU4\nB4GZWYNzEJiZNTgHgZlZg3MQmJk1OAeBmVmDcxCYmTU4B4GZWYNzEJiZNTgHgZlZg3MQmJk1OAeB\nmVmDcxAMIBFw3XVw6aVZ18TMisR3KBsgfvc7+NSnYNMmWLMm/WlmjcV3KGtQO3bAt74Fb387nHIK\nPPggvPwybN6cdc3MrCgGZ10B67slS+DCC2HUKHjgATj44FQ+YQI8+SRMnpxt/cysGNwiKKj//V+Y\nOhUuuwwWLGgPAUjPn3giu7qZWbH0GASS9pC0SNISSUslzSiVj5S0QNIKSXdLGl72mumSVkpaLumE\nsvJjJD0i6XFJV9bmV2oMixfD2WfDeeeBOvQIOgjMrDd6DIKIeAl4T0QcDUwGTpI0BbgcWBgRbwLu\nBaYDSJoEnAlMBE4CrpZeOVVdA1wYEYcDh0uaWu1fqFGsXg0HHND5z9q6hszMKlFR11BEbC093YM0\nrhDAacDsUvls4PTS81OBWyNiR0SsAlYCUySNAfaOiMWl7W4qe431UndB4BaBmfVGRUEgaZCkJUAr\ncE/pZD46IjYAREQrsG9p83FAS9nL15bKxgFrysrXlMqsD9wiMLNqqWjWUETsAo6W9HrgvyUdSWoV\nvGqzalZs5syZrzxvamqiqampmrsvvO6C4KCDYNUq2LULBnk6gNmA1dzcTHNzc7/30+sFZZK+AmwF\nPgk0RcSGUrfPfRExUdLlQETEFaXt7wJmAE+1bVMqnwa8OyI+3cl7eEFZN7Ztg5EjYevWrk/0Y8em\nKaXjx9e3bmaWnZotKJP0hrYZQZL2BI4HlgO3A+eXNjsPmFd6fjswTdIQSROAQ4EHSt1HmyVNKQ0e\nn1v2GuuFNWvSCb67b/vuHjKzSlXSNbQfMFvSIFJwzImI+ZJ+BcyVdAHp2/6ZABGxTNJcYBmwHbi4\n7Ov9JcCNwFBgfkTcVdXfpkF01y3Upm3A+F3vqk+dzKy4egyCiFgKHNNJ+bPA+7t4zTeAb3RS/iBw\nVO+raeVWr4b99+9+G7cIzKxSHkosoN60CMzMeuIgKKBKgmDCBAeBmVXGQVBAlbYI3DVkZpVwEBRQ\nJUEwdmy6J8G2bfWpk5kVl4OgYCKgpaXnweLddkth8dRT9amXmRWXg6BgNm2CPfeEYcN63tYDxmZW\nCQdBwVTSLdTGA8ZmVgkHQcH0Jgg8YGxmlXAQFEwli8nauGvIzCrhICiY3nYNuUVgZj1xEBRMb7uG\nnngizTQyM+uKg6BgehMEI0akaaSbNtW2TmZWbA6CgmlpqTwIwAPGZtYzB0GBvPwyPPMM7Ldf5a/x\ngLGZ9cRBUCBr16YQ2G23yl/jAWMz64mDoEB6Mz7Qxi0CM+uJg6BA+hIEbhGYWU8cBAXSm8Vkbdwi\nMLOeOAgKpC8tggMOSDe737GjNnUys+JzEBRIX4Jgjz1g9OgUBmZmnXEQFEhfggDcPWRm3XMQFERE\n34PAA8Zm1h0HQUFs3gyDBsHw4b1/rVsEZtadHoNA0nhJ90r6raSlkj5bKp8haY2kh0qPE8teM13S\nSknLJZ1QVn6MpEckPS7pytr8SgNTX1sD4BvUmFn3BlewzQ7gcxHxsKRhwIOS7in97HsR8b3yjSVN\nBM4EJgLjgYWSDouIAK4BLoyIxZLmS5oaEXdX79cZuPoTBL7ekJl1p8cWQUS0RsTDpedbgOXAuNKP\n1clLTgNujYgdEbEKWAlMkTQG2DsiFpe2uwk4vZ/1bxj9DQK3CMysK70aI5B0EDAZWFQq+oykhyVd\nK6mt93oc0FL2srWlsnFA+STGNbQHivWgL4vJ2oweDVu2pIeZWUeVdA0BUOoW+k/gsojYIulq4GsR\nEZK+DnwX+GS1KjZz5sxXnjc1NdHU1FStXRfS6tVw8sl9e63UPnPoqKOqWy8zy05zczPNzc393k9F\nQSBpMCkEbo6IeQARsbFskx8Cd5SerwXKv7uOL5V1Vd6p8iCw/nUNQfuAsYPAbODo+CV51qxZfdpP\npV1D1wPLIuKqtoJSn3+bM4BHS89vB6ZJGiJpAnAo8EBEtAKbJU2RJOBcYF6fat2AentDmo48YGxm\nXemxRSDpncDHgKWSlgABfAk4W9JkYBewCrgIICKWSZoLLAO2AxeXZgwBXALcCAwF5kfEXVX9bQao\nHTugtRXG9WNExVNIzawrihze2VxS5LFeWWlpgeOO69/1gubNg2uvhTvu6HlbMysmSUREZ7M5u+WV\nxQXQ3/EB8BRSM+uag6AAqhEEbbOG3NAys44cBAVQjSAYNgyGDIHnnqtOncxs4HAQFEB/FpOVGzsW\n1q3r/37MbGBxEBRANVoEkGYdOQjMrCMHQQFUKwjGjoW1XS7hM7NG5SAogP4uJmvjFoGZdcZBkHMv\nvAAvvQT77NP/fblFYGadcRDkXFtrQL1eIvJabhGYWWccBDlXrfEBcIvAzDrnIMi5agaBWwRm1hkH\nQc61tMD48dXZ1+jRsHFjuoidmVkbB0HOrV+funSqYffdYdQoePrp6uzPzAYGB0HOrV8P++1Xvf15\nnMDMOnIQ5FwtgsDjBGZWzkGQc9UOgnHj3CIws1dzEOTYzp2waRPsu2/19ukWgZl15CDIsaefTiuK\nB/d4Q9HKeQqpmXXkIMixancLgQeLzey1HAQ5tn49jBlT3X26RWBmHTkIcswtAjOrBwdBjtUiCEaN\ngq1bYdu26u7XzIrLQZBjtQgCKe3T3UNm1qbHIJA0XtK9kn4raamkS0vlIyUtkLRC0t2Shpe9Zrqk\nlZKWSzqhrPwYSY9IelzSlbX5lQaOWgQBeJzAzF6tkhbBDuBzEXEkcBxwiaQjgMuBhRHxJuBeYDqA\npEnAmcBE4CTgaumVq+lfA1wYEYcDh0uaWtXfZoCpVRB4nMDMyvUYBBHRGhEPl55vAZYD44HTgNml\nzWYDp5eenwrcGhE7ImIVsBKYImkMsHdELC5td1PZa6wTbhGYWT30aoxA0kHAZOBXwOiI2AApLIC2\n9a/jgJayl60tlY0D1pSVrymVWScioLXVLQIzq72K16xKGgb8J3BZRGyRFB026fj3fpk5c+Yrz5ua\nmmhqaqrm7nPv2Wdhr71g6NDq73vsWHjooerv18zqq7m5mebm5n7vp6IgkDSYFAI3R8S8UvEGSaMj\nYkOp26ftKvdrgf3LXj6+VNZVeafKg6AR1apbCNw1ZDZQdPySPGvWrD7tp9KuoeuBZRFxVVnZ7cD5\npefnAfPKyqdJGiJpAnAo8ECp+2izpCmlweNzy15jHdQyCNw1ZGblemwRSHon8DFgqaQlpC6gLwFX\nAHMlXQA8RZopREQskzQXWAZsBy6OiLZuo0uAG4GhwPyIuKu6v87AUesgWLcujUO8Mp/LzBqW2s/R\n+SEp8liverriCnjmGfj2t2uz/+HDYdUqGDmyNvs3s/qTRET0+uudVxbnVC1bBOBxAjNr5yDIqVoH\ngccJzKyNgyCn3CIws3pxEORUrRaTtXGLwMzaOAhyyi0CM6sXB0EObdkCu3bB3nvX7j3cIjCzNg6C\nHGprDdRyjr9bBGbWxkGQQ7XuFgK3CMysnYMgh+oRBGPGwMaNsHNnbd/HzPLPQZBD9QiC3XeHffaB\np5/ueVszG9gcBDlUjyAAdw+ZWeIgyKF6BYEHjM0MHAS55BaBmdWTgyCH1q9Pg7m15haBmYGDIJfc\nIjCzenIQ5MxLL8ELL8CoUbV/L7cIzAwcBLnT2gqjR8OgOnwybhGYGTgIcqde3ULgFoGZJQ6CnKln\nEIwalS5wt21bfd7PzPLJQZAz9QwCKb3X+vX1eT8zyycHQc7UMwggdQ95nMCssTkIcqbeQTB2rMcJ\nzBqdgyBnsmgROAjMGluPQSDpOkkbJD1SVjZD0hpJD5UeJ5b9bLqklZKWSzqhrPwYSY9IelzSldX/\nVQaGLFoE7hoya2yVtAhuAKZ2Uv69iDim9LgLQNJE4ExgInAScLX0yn22rgEujIjDgcMldbbPhucW\ngZnVW49BEBH3A8918qPObqR4GnBrROyIiFXASmCKpDHA3hGxuLTdTcDpfavywLVzJzzzTFpQVi9u\nEZhZf8YIPiPpYUnXShpeKhsHtJRts7ZUNg5YU1a+plRmZZ5+Ot0sZvDg+r2nWwRm1tdTztXA1yIi\nJH0d+C7wyepVC2bOnPnK86amJpqamqq5+1yqd7cQtLcIItK6AjMrjubmZpqbm/u9nz4FQURsLPvr\nD4E7Ss/XAvuX/Wx8qayr8i6VB0GjyCIIhg1Lt63cvBlGjKjve5tZ/3T8kjxr1qw+7afSriFRNiZQ\n6vNvcwbwaOn57cA0SUMkTQAOBR6IiFZgs6QppcHjc4F5farxAJZFEADsvz889VT939fM8qHHFoGk\nW4AmYJSk1cAM4D2SJgO7gFXARQARsUzSXGAZsB24OCKitKtLgBuBocD8tplG1q61NZsgOPRQ+P3v\n4S1vqf97m1n2egyCiDi7k+Ibutn+G8A3Oil/EDiqV7VrMOvXw8SJ9X/fQw+F3/2u/u9rZvnglcU5\nklXX0CGHOAjMGpmDIEeyCgK3CMwam4MgRxwEZpYFtY/l5oekyGO9aikC9twTnn8ehg6t73vv2JGm\nkT73XKqDmRWTJCKi1yuC3CLIibaTcL1DANJK5gMPhCefrP97m1n2HAQ5kVW3UBt3D5k1LgdBTuQh\nCH7/++ze38yy4yDIifXrYcyYnrerFbcIzBqXgyAn8tAicBCYNSYHQU44CMwsKw6CnFi3Ll0SOisH\nHpguR/3yy9nVwcyy4SDIiVWr4KCDsnv/IUNSEPkqpGaNx0GQE08+CRMmZFsHdw+ZNSYHQQ5s2ZIe\n9bxXcWccBGaNyUGQA23dQlnfKtJBYNaYHAQ5kIduIXAQmDUqB0EOOAjMLEsOghzISxAcfHCaNbRj\nR9Y1MbN6chDkQF6CYOhQ2HdfaGnJuiZmVk8OghzISxCAu4fMGpGDIGMRKQiyXExWzkFg1ngcBBl7\n9lkYNAhGjsy6JomDwKzxOAgylqduIYBDDvF9CcwaTY9BIOk6SRskPVJWNlLSAkkrJN0taXjZz6ZL\nWilpuaQTysqPkfSIpMclXVn9X6WY8hYEbhGYNZ5KWgQ3AFM7lF0OLIyINwH3AtMBJE0CzgQmAicB\nV0uvrJe9BrgwIg4HDpfUcZ8NKW9BcMgh8MQTsGtX1jUxs3rpMQgi4n7guQ7FpwGzS89nA6eXnp8K\n3BoROyJiFbASmCJpDLB3RCwubXdT2WsaWt6CYNgwGD48XRbbzBpDX8cI9o2IDQAR0QrsWyofB5TP\nQl9bKhsHrCkrX1Mqa3h5CwJw95BZoxlcpf1ElfbzipkzZ77yvKmpiaampmq/RS7kOQgG6CE3GzCa\nm5tpbm7u9376GgQbJI2OiA2lbp+nS+Vrgf3LthtfKuuqvEvlQTBQ7dqVLumQlzUEbdwiMCuGjl+S\nZ82a1af9VNo1pNKjze3A+aXn5wHzysqnSRoiaQJwKPBAqftos6QppcHjc8te07DWr4cRI2CvvbKu\nyas5CMwaS48tAkm3AE3AKEmrgRnAN4H/kHQB8BRpphARsUzSXGAZsB24OCLauo0uAW4EhgLzI+Ku\n6v4qxZPHbiFwEJg1GrWfp/NDUuSxXtV2881w551wyy1Z1+TVnn8exo+HF17I/mY5ZlY5SUREr//X\nemVxhvLaIhgxAvbcEzZsyLomZlYPDoIM5TUIIC0sc/eQWWNwEGQoz0HgcQKzxuEgyJCDwMzywEGQ\nke3bobUV9t+/522z4CAwaxwOgoy0tMB++8Huu2ddk845CMwah4MgI3nuFoL2IGiAWbxmDc9BkJG8\nB8GoUenPZ5/Nth5mVnsOgozkPQik1CpYuTLrmphZrTkIMpL3IAB461vh/vuzroWZ1ZqDICNFCILT\nToN5DX9pQLOBz9caysiYMfDQQzB2bNY16dqLL8Lo0WnQ+I1vzLo2ZtYTX2uoQLZuhc2bUxjk2dCh\ncPzx8JOfZF0TM6slB0EGVq2CAw+EQQU4+u4eMhv4CnAqGniKMD7Q5gMfgHvvTa0YMxuYqnXPYuuF\nIgXBPvvA294G99yTWgfW2DZsgAcfTONbgwbB5z8Pe+yRda2sv9wiyECRggDcPTQQRKTPcOJEeO97\n4d//vedW3ssvw89+BjNnwgc/COPGpdd/73vppkWLFsFxx8GKFXX5FayGPGsoA2ecAWefDX/xF1nX\npDKrVsGUKekey7vtlnVtrLdWrIDLLoPVq9tP4tdfn07kZ54JF16YWn0SrFsH8+enx733whFHpOB4\n61vT48AD2+9aFwH/+q/wla/At74F55/vO9plra+zhhwEGTj6aLj22vQfqyje8hb4wQ/gz/4s65pY\npV54Af7+79NJ/0tfgs9+9tUXOWxpgdmz4YYbYK+9YPBgeOopmDoVTj45/bnvvj2/z6OPwkc/Cm9+\nM/zLv8Dw4bX7nax7DoICGTECfv/79uv5FMFXvwrbtsG3v511TQaup59O13Y64oi+7+O559K/rV//\nOoXA8cfDN7/Z/VTlXbvSCnIpdfUM7sPI4bZt8Ld/m1oSt9yS9mP15yAoiOeeS83rzZuL1Yx+6CE4\n6yx4/PFi1bsIImDOnNR9s3Mn/PVfw/TpPXfDbdiQWmmPPQZPPJECYOfOdJvRww6Dv/mb+p+Qb7sN\nPvUpuPXW1KVk9dXXIPCsoTprGygu2sn06KPhpZdg+XKYNCnr2hTDz3+evl0fd1zXn/fGjfDpT8Oy\nZWnh3n77wcc/DgsXpgHd8eNf+5rt2+H734d/+IfUJfOhD6WT/8EHp1Zmlv+2Tj89tXg/8hG48840\n9mD551lDdVa0GUNtJDj1VM8eqsT27fCFL8A558AFF6Rv5zNnvvZGPz/+MfzJn6QT+EMPwbHHphP/\nwoVwwglpDOm22179moUL03jN3Xen7pzvfz+FwZQp8IY35OMLRlMT/PCHaaaRZxQVQ79aBJJWAZuB\nXcD2iJgiaSQwBzgQWAWcGRGbS9tPBy4AdgCXRcSC/rx/ERU1CCBNI/3KV1K3hXVu3TqYNg1e9zpY\nsiStw3jwQbj5ZnjnO9NJ/5xz4P/+L5X/+Mfwjne8eh+77ZYGd9/znjS7bMECuPRS+PKXU2D84z+m\nzyIPJ/2unH56Gu+YOjUFVmctG8uP/rYIdgFNEXF0REwplV0OLIyINwH3AtMBJE0CzgQmAicBV0t5\n/qdcG0UOgne/O40RrF+fdU3y6b77UlfICSfAT3/a3k3ztrfBVVfBmjUpSH/xi3RiXLLktSFQ7rjj\n4OGH0wn16KNTS2DZsnSSLcL/nAsugEsuScdj06asa2Pd6ddgsaQngbdFxKaysseAd0fEBkljgOaI\nOELS5UBExBWl7e4EZkbEok72O2AHi08+OfUJf/CDWdekb84+OwXCRRdlXZP82LUrzcz5/vfTN//3\nv7+6+49I3U1DhlR3v/Xyd3+XxksWLoRhw7KuzcCW1dVHA7hH0mJJnyyVjY6IDQAR0Qq0zUQeB7SU\nvXZtqaxhRMDSpf2bHpi1Rl9lvGtXWpj1s5+lOfOf/3xaWzF/fpqyWe0QgPTtv6ghAHDFFXDkkWkh\n5YsvZl0b60x/Zw29MyLWS3ojsEDSClI4lOvTV/uZM2e+8rypqYmmpqa+1jE3liyBPfdMt4AsqpNO\ngr/8y7RYae+9s65Ndb30UjrJr14Nra2vfaxbl7r2RoyAww9Pg8CHHQZf/GJq6ZUv1rJ2UlqBfM45\n6TjNmzfw/u1kpbm5mebm5n7vp2rrCCTNALYAnySNG7R1Dd0XERM76Rq6C5jRSF1DA2VR1tSpcN55\nqZuoqFasSKtqV61Kj6eegmeeSX33BxyQpnGOGfPax8EHu3ujr3buhIsvTl+I7ryzWAsqi6LuC8ok\n7QUMiogtkl4HLABmAe8Dno2IKyR9ERgZEZeXBot/BPwpqUvoHuCwzs74AzUI3vIWuPrqNHukyO67\nL01ZvP/+4rVuduyA73wHvvvd1LKZNAkOOig99tvP11KqtYg06+yOO9JsqHEN1Tlce1ksKBsN/Lek\nKO3nRxGxQNKvgbmSLgCeIs0UIiKWSZoLLAO2AxcPyLN9F558MnUvvP3tWdek/97zHpg1C045BX75\nSxg5MusaVeY3v0kzWUaNgsWL08nf6ktKA+sjR8K73pUub37IIVnXynyJiTq58sp0ca5rr826JtXz\nuc+lk+tdd+W7f/yll+DrX0/91Fdc4atk5sW//Vv6QnHnnWlhnfWfrzWUc01NaYZJUaeNdmbnznR5\ng9Gj03/qWpxcI1If/uLF6bFkSfo2eeSRqVvnyCPTgG3brJoIeP75dGXNlpY08PvP/5y2ufpqGDu2\n+nW0vpszJy2W+4//gD//86xrU3wOghzbtCkNMra2pllDA8mWLWn65DnnpKtPVsPKlemiZb/8ZTr5\n7757uvzCscfCMcekE/2yZfDb36Y/V69uX6TX0pLunLX//u2PqVPT1EW3AvJp4UL42MdSl9EnPpF1\nbYrNQZBjs2enKXP/9V9Z16Q2WlrS2McPfpBWvfbF88/D3Llw443pKprTpqWxiGOP7XlA8cUXU3hI\n6cTv6+EXz2OPpdby6aenQPCgfd84CHLsjDPSP/Bzz826JrWzeHGaI37jjenSCPvs0/32Een6+23X\n4Zk/P107/7zz4MQT8z3mYLWxaVO6a9/rX5+uvOq1Br3nIMipbdvS/PMnnhj486bvuAO+9rU0R3/o\n0LSCuu0xalT6pr9yZftjt91SH/+0aekx0I+P9ezll9P1iR54AG6/Pd27wyrnIMipO+5I94m9776s\na1I/EenCdI89lh4rVrSPk7Stxj3sMJ/4rXMR6Qqr3/lOaklv3Qp//OOr/5wwIY0nHH+8u5HKOQhy\n6sIL4aij0l2nzKxy//M/acr1Xnuly3qX/7lkSboX87p1aTrwJz7h9QjgIMilnTvTatVFi4p76Wmz\nPFu6FG64IY0pTJqUFgx++MMpMBpRVlcftW788pcpCBwCZrVx1FGp63XNmrQeYc6cNHPsoovSF7AB\n8H2yLhwENTRvXt+nU5pZ5YYMSbPzfvrT1EqYMCHd+/nNb07XlXr66axrmG/uGqqRiHSp4jlz0iIo\nM6uviHRL0OuvT/d+/vjH4fLLUyt9oHLXUM4sX56ucXP00VnXxKwxSWnV+/XXp/+PbdOVP/952LCh\n69ft3JkWSRb8u2ivOAhq5LbbinNvWbOBbvToNJbw6KPptp+TJqUbCrW2prLZs+Gyy1JwjBiRAuPC\nC9O6hkbgIKiRefPSbR3NLD/GjoV/+qd01dwtW9KlyD/84XQF3fHj04LI1avTtNRnnkl35Hv++axr\nXXseI6iB225Lzc/HHvOlEszybNeudJHCzuzcmf4fL1iQBqGLMPvPYwQ50doKf/VX6fo5DgGzfOsq\nBCCNKVx5JXz60+mugotec1PdgcMtgiqKgA98IM0S+vrXs66NmVXLT36SFqtdc03qSsorryzOgWuu\ngeuuSwvJ3BowG1iWLIFTT01XR33jG2HffdOfbY8jj4R3vKP9JklZcBBkbMWK1Hy8//50tU0zG3j+\n+Md0//GNG9MitY0b258/+GA6DzQ1pUupT52aLrRYTw6CDG3fnr4JnH9+uoSumTWmjRvhnnvSLKQF\nC1Lr4ZRT4KyzYMqU2k8ndxBk6CtfgV//Ot1cxesGzAzSjKTf/CbNIpwzJy0wPeus9Jg8uTbnCgdB\nRn7xi3SNkyVLBvbSdTPru4gUCnPmpPtx7757WnD65jfDm96UHiNG9P99ChMEkk4EriRNXb0uIq7o\nZJtCBMHSpenD/M534EMfyro2ZlYEEenWrnfd1X7jphUr0q05jzgiDTqffDK8973pTn+9UYh1BJIG\nAf8MTAWOBD4qqXBDq4sWpVXDxx+fbjjTMQSam5szqVc1FLnu4PpnzfXvmZTGC776VbjlljTI/Ic/\npNtzfvnL6fac3/xmuizGRz4CP/oRPPdcbetU7wVlU4CVEfFURGwHbgUKcSGGCPjZz+B974Mzz4QT\nTkizBz772dduW+T/DEWuO7j+WXP9+2bQoHQfhfe/H77wBfj5z9N9vU86CebOTeHw3vemq6f+6Efw\nyCPVvQ7kaRtyAAAD8ElEQVTS4OrtqiLjgJayv68hhUMuRMDmzWl1cMfHffeln02fDmef7XUCZlZb\n++6bFrFdcEGatnrffWksct68tGB11ap0e86jjkqti4kT+/5e9Q6CXLvqKpgxA8aMaX+MHp3+/PKX\n06ph3yjbzOrtda9L01BPOaW97MUX0+W1ly6F4cP7t/+6DhZLejswMyJOLP39ciA6DhhLyv9IsZlZ\nDuV+1pCk3YAVwPuA9cADwEcjYnndKmFmZq9S166hiNgp6TPAAtqnjzoEzMwylMsFZWZmVj+Z3Y9A\n0omSHpP0uKQvdrHNP0laKelhSZPrXcfu9FR/Se+W9Lykh0qP/5dFPTsj6TpJGyQ90s02eT723dY/\nz8ceQNJ4SfdK+q2kpZIu7WK7XH4GldQ/r5+BpD0kLZK0pFT3GV1sl9dj32P9+3TsI6LuD1IA/Q44\nENgdeBg4osM2JwE/LT3/U+BXWdS1H/V/N3B71nXtov5/BkwGHuni57k99hXWP7fHvlS/McDk0vNh\npHGzIv37r6T+uf0MgL1Kf+4G/AqYUpRjX2H9e33ss2oRVLKw7DTgJoCIWAQMlzS6vtXsUqUL43J5\nCbqIuB/obq1ino99JfWHnB57gIhojYiHS8+3AMtJa2zK5fYzqLD+kNPPICK2lp7uQRon7dg/nttj\nDxXVH3p57LMKgs4WlnX8h9Rxm7WdbJOVSuoPcFypaflTSZPqU7WqyPOxr1Qhjr2kg0itm443QizE\nZ9BN/SGnn4GkQZKWAK3APRGxuMMmuT72FdQfennsvaCsdh4EDoiIrZJOAm4DDs+4To2iEMde0jDg\nP4HLSt+sC6WH+uf2M4iIXcDRkl4P3CZpUkQsy7pelaqg/r0+9lm1CNYCB5T9fXyprOM2+/ewTVZ6\nrH9EbGlrwkXEncDukvapXxX7Jc/HvkdFOPaSBpNOojdHxLxONsn1Z9BT/YvwGUTEH4D7gBM7/CjX\nx75NV/Xvy7HPKggWA4dKOlDSEGAacHuHbW4HzoVXViQ/HxEb6lvNLvVY//I+RUlTSFN1n61vNbsl\nuu5HzPOxb9Nl/Qtw7AGuB5ZFxFVd/Dzvn0G39c/rZyDpDZKGl57vCRwPPNZhs9we+0rq35djn0nX\nUHSxsEzSRenH8W8RMV/SyZJ+B/wR+EQWde1MJfUH/kLSp4HtwDbgrOxq/GqSbgGagFGSVgMzgCEU\n4NhDz/Unx8ceQNI7gY8BS0t9vQF8iTQLLfefQSX1J7+fwX7AbKVL4g8C5pSOdSHOPVRQf/pw7L2g\nzMyswWW2oMzMzPLBQWBm1uAcBGZmDc5BYGbW4BwEZmYNzkFgZtbgHARmZg3OQWBm1uD+Pwt7Yog8\nUhinAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff0fee51908>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "val,bi = np.histogram(dis,bins=50)\n",
    "plt.plot((bi[:-1]+bi[1:])/2.,val)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class bond_formation_reaction(reaction):\n",
    "    \"\"\"\n",
    "    creates a bond of given type between two particles\n",
    "    \"\"\"    \n",
    "    def __init__(self, name, rate, bonds, ):\n",
    "        self.bonds = bonds\n",
    "        super(bond_formation_reaction,self).__init__(name,rate)\n",
    "        \n",
    "    def effect_reaction(self,encunter_complex_in):\n",
    "        p0 = encunter_complex_in[0]\n",
    "        p1 = encunter_complex_in[1]\n",
    "        if not(any([term.name==self.bonds.name for term in p0.energyterms])):\n",
    "            #print(\"%s: effect it\"%self.name)\n",
    "            self.counter += 1\n",
    "            bonds = [deepcopy(self.bonds),deepcopy(self.bonds)]\n",
    "            bonds[0].assign_partner(p0)\n",
    "            bonds[1].assign_partner(p1)\n",
    "            p0.add_energyterm(bonds[1])\n",
    "            p1.add_energyterm(bonds[0])\n",
    "        return [p0,p1]\n",
    "        \n",
    "    \n",
    "\n",
    "class bond_breaking_reaction(reaction):\n",
    "    \"\"\"\n",
    "    breaks abond (type??)\n",
    "    \"\"\"\n",
    "    def __init__(self, name, rate, bonds):\n",
    "        self.bonds = bonds\n",
    "        super(bond_breaking_reaction,self).__init__(name,rate)\n",
    "    \n",
    "    def effect_reaction(self,encunter_complex_in):\n",
    "        p0 = encunter_complex_in[0]\n",
    "        p1 = encunter_complex_in[1]\n",
    "        #look for the term with the right name\n",
    "        n0 = [term.name==self.bonds.name and term.partner==p1 for term in p0.energyterms]\n",
    "        n1 = [term.name==self.bonds.name and term.partner==p0 for term in p1.energyterms]\n",
    "\n",
    "        if any(n0) and any(n1):\n",
    "            m0 = list(np.where(n0)[0])\n",
    "            m1 = list(np.where(n1)[0])\n",
    "            #print(m0)\n",
    "            #print(m1)\n",
    "            #print(\"%s effect breaking ==== \"%self.name )\n",
    "            self.counter += 1\n",
    "            del(p0.energyterms[m0[0]])\n",
    "            del(p1.energyterms[m1[0]])\n",
    "        return [p0,p1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pbcdist(p1,p2):\n",
    "    box = p1.ptype.geometry.size\n",
    "    pos1 = p1.pos - np.floor(p1.pos/box) * box\n",
    "    pos2 = p2.pos - np.floor(p2.pos/box) * box\n",
    "    delta = (pos1 - pos2)\n",
    "    delta = delta - np.rint(delta/box)* box\n",
    "    return np.sqrt((delta*delta).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ensemble has 5 trajectories\n",
      "Length 5001 (5000)\n",
      "Length 5001 (5000)\n",
      "Length 5001 (5000)\n",
      "Length 5001 (5000)\n",
      "Length 5001 (5000)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def setup_dim():\n",
    "    membrane = geometry(dim=2,size=5)\n",
    "    q1 = ptype(5e-2,\"mono\",membrane)\n",
    "    a = particle(q1,pos=np.array([2,2]),age=0)\n",
    "    b = particle(q1,pos=np.array([0,0]),age=0)\n",
    "    init_pos=[a,b]\n",
    "    hca = hardcore(50.,.5,a)\n",
    "    hcb = hardcore(50.,.5,b)\n",
    "    b.add_energyterm(hca)\n",
    "    a.add_energyterm(hcb)\n",
    "    cov1 = bondterm(k=10.,d=.8)\n",
    "\n",
    "    ec1=bim_encunter_complex(\"grab\",q1,q1,dist=1)\n",
    "    rea1 = bond_formation_reaction(\"dimerization 1\", rate=.1, bonds=cov1)\n",
    "    ec1.add_reaction(rea1)\n",
    "    rea2 = bond_breaking_reaction(\"disruption 1\", rate=.01, bonds=cov1)\n",
    "    ec1.add_reaction(rea2)\n",
    "    t = traj(init_pos)\n",
    "    return (t, [ec1])\n",
    "\n",
    "\n",
    "n_replicas=5\n",
    "sett = []\n",
    "for i in range(n_replicas):\n",
    "    t2, reacts = setup_dim()\n",
    "    run1 = runner(100000, timestep = .1,stride=20)\n",
    "    run1.run(t2, reactions_bim = reacts ,verbose=False)\n",
    "    sett.append(t2)\n",
    "    \n",
    "em1 = ensemble(sett)\n",
    "em1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<__main__.traj at 0x7ff0fda8ee48>"
      ]
     },
     "execution_count": 246,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "em1.trajs[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "distances=[]\n",
    "for tr in em1.trajs:\n",
    "    dis = np.array([pbcdist(frame[0],frame[1]) for frame in tr.ptraj])\n",
    "    distances.append(dis)\n",
    "\n",
    "dis=np.concatenate(distances)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYk+X1//H3AUVwYakVUEDFBVlqVaqAK2NRFPf6dUEp\nFWu/RXFB6U8FtWWsVYG2KkotdUNUBBH9igsgog5uIItQkE3cULGAqGyiCM79++Nk6jDMMDNJJs+T\n5PO6rlxknjxJDs9ATu7t3BZCQERE8letqAMQEZFoKRGIiOQ5JQIRkTynRCAikueUCERE8pwSgYhI\nnqs0EZjZg2a20szmlTo2xMwWmdlcM3vKzOqXemyAmS1NPN611PH2ZjbPzN4zs7vS/1cREZFkVKVF\nMAI4qcyxyUC7EMKhwFJgAICZtQXOA9oA3YB7zcwSz/kncEkIoRXQyszKvqaIiESg0kQQQngD+LrM\nsSkhhOLEj9OB5on7ZwBjQghbQggf40mig5k1BXYLIcxMnPcIcFYa4hcRkRSlY4zgt8CExP1mwKel\nHlueONYM+KzU8c8Sx0REJGIpJQIzuxHYHEIYnaZ4REQkw3ZI9olm1gs4BfhlqcPLgRalfm6eOFbR\n8YpeWwWQRESSEEKwys/aWlVbBJa4+Q9mJwPXAmeEEDaVOu9ZoLuZ1TGzlsABwIwQwgpgrZl1SAwe\n/wYYv703DCFk7W3gwIGRx5CPsSv+6G+KP9pbsiptEZjZ40ABsLuZfQIMBG4A6gAvJSYFTQ8h9Akh\nLDSzscBCYDPQJ/wY3eXAw0BdYEIIYVLSUYuISNpUmghCCBeWc3jEds6/Hbi9nOOzgYOrFZ2IiNQ4\nrSyuAQUFBVGHkLRsjh0Uf9QUf3ayVPqVaoqZhTjGJSISZ2ZGqMHBYhERyVFKBCIieU6JQEQkzykR\niIjkOSUCEZE8p0QgIpLnlAhERPKcEoGISJ5TIhARyXNKBCIieU6JQESqbfVq2LIl6igkXZQIRKRa\nvv4a2reHYcOijkTSRUXnRKTKQoALL4SPPoLiYpgxI+qIpDQVnRORGjdqFMybBy+9BMuWwfvvRx2R\npIMSgYhUyUcfwTXXeDLYbTc491x44omoo5J0UCIQyQOp9rRu2QI9e8L118Ohh/qxCy6A0aNTj02i\np0Qgkge6doWnnkr++YMGwU47Qb9+Px478khYtw7mz089PomWBotFctz69dC4MdSvD3PmwF57Ve/5\nM2bA6afD7NnQvPnWj113HeywA9x2W/rileRpsFhEyjVrlk/3vPRS+O1vq9dNtGED9OjhU0XLJgHw\n7qExY1LvepJoKRGI5Lhp07wb56ab4MsvYfjwqj/3mmvgmGN8YLg8hx4KdepoGmm2UyIQyXEliWDH\nHeHRR+GPf4T33qv8ecOHwyuvwN13V3yOGXTvrkHjbKcxApEcFoKPD8ydC82a+bFhwzwhvPmm9++X\n9f333hJ4+WUYPx4OOmj777FkCRx/PHz6KdSunf6/g1SdxghEZBsffAD16v2YBAD69IGGDeH227c9\nf+VK6NLFP9TffrvyJAB+TtOmMHVq+uKWzFIiEMlhJd1CpdWqBQ895C2DWbN+PD5zJhxxBPzyl/DM\nM9CgQdXfR2sKspsSgUgOmz4dOnXa9nizZt7337MnbNwIjzwCp57qx26+2ZNFdZx/Pjz9tHcrSfYp\np4dQRHLFtGn+YV+e88/3MYAOHfwDvKgI2rZN7n323hvatIHJk+G005IOVyJSad43swfNbKWZzSt1\nrJGZTTazJWb2opk1KPXYADNbamaLzKxrqePtzWyemb1nZnel/68iIqV9840P5B52WMXn/OMf/sE9\nY0bySaCEuoeyV1UagCOAk8oc6w9MCSEcBLwCDAAws7bAeUAboBtwr5mVjGD/E7gkhNAKaGVmZV9T\nRNJo1iz4+c+9NERFGjXy8hENG6b+fueeCy+84F1Nkl0qTQQhhDeAr8scPhMYmbg/Ejgrcf8MYEwI\nYUsI4WNgKdDBzJoCu4UQZibOe6TUc0SkBpQ3UFyTGjf2bqbnnsvce0p6JDtY3DiEsBIghLACaJw4\n3gz4tNR5yxPHmgGflTr+WeKYiNSQadPKHyiuSeoeyk7pmjWk1V8iMRKCzxjKZIsA4Fe/8oVoGzZk\n9n0lNcnOGlppZk1CCCsT3T6rEseXAy1Kndc8cayi4xUqLCz87/2CggIKCgqSDFUk/3z0kZeUaNGi\n8nPTqWFDaNfOK5V27pzZ985HRUVFFBUVpfw6VSoxYWb7As+FEA5O/DwY+CqEMNjMrgcahRD6JwaL\nRwEd8a6fl4ADQwjBzKYDVwEzgReAu0MIkyp4P5WYEEnB44/7vP5x4zL/3n37+jqF667L/Hvnuxor\nMWFmjwNv4TN9PjGzi4FBwIlmtgTokviZEMJCYCywEJgA9Cn1iX458CDwHrC0oiQgIqnL9EBxaR07\nenkKyR4qOieSgw4/3FcJH3VU5t/7ww/h2GNh+XY7f6UmJNsiUCIQyTEbN8Iee/jeA3XrZv79Syqe\nzplT/mY2UnNUfVREAB+o/dnPokkC4HsUqHsouygRiOSYKMcHSnTqpESQTZQIRHJMRRVHM6ljR49D\nsoPGCERySAiw117+IbzPPtHFsWaNjw+sWVP+LmhSMzRGIJIHtmzxiqIVWbbM++j33jtzMZWnYUNf\nzPbuu9HGIVWjRCCSRYYO9XLRd9/t3/7LKukWsmp/J0y/Tp3UPZQtlAhEssQXX3jJ6Oefh/vug0sv\nhc2btz4nDgPFJTRzKHsoEYhkiZtvhu7doVs3eOst+Pxz6NrV1wuUiKLiaEWUCLKHBotFssCiRXDc\ncf7nT3/qx374AQYM8JpCzz4LLVv6Y6tXQ7160cYLPp7RsCF89ll6Nr6RyiU7WKzxfJEscO21/qFf\nkgQAateGIUN8zKCgAC67zO/HIQmAzxZq3x5mzoQTT4w6GtkedQ2JxNxLL8HixXD55eU/3qsX/N//\nwT//GZ/xgRLqHsoOahGIxNgPP8Af/uDf/Le39/DRR8P8+fGbs9+pEzz8cNRRSGXUIhCJsYce8g3m\nf/Wrys9t0gR2373mY6qOkhXGGvKLNw0Wi8TU+vXQqpVPF/3FL6KOJnnNmsHrr8N++0UdSe7TymKR\nHDNokE8PzeYkACpAlw2UCERiaNkyGD4cbr016khSpwJ08adEIBJDAwbAlVfmxsYumjkUfxojEImZ\nmTPhrLPgvfdgl12ijiZ1Gzb4QPZXX21/5pOkTmMEIjkgBF88VliYG0kAYNdd4cADYe7cqCORiigR\niMTIhAmwahVcfHHUkaSXuofiTYlAJCa2bIHrroPBg+O3MCxVSgTxpkQgEhMjR3otodNOizqS9NMU\n0njTYLFIDHzzDRx0kFcS7dAh6mjSr7jYV0i//z7ssUfU0eQuDRaLZLG77vJ6QbmYBABq1YIjjoAZ\nM6KORMqjRCBSg776qvLFVKtWwZ13wm23ZSamqKh7KL6UCERqUP/+0Lmzl4ouvZNYabfcAj16wP77\nZzS0jOvY0XdQk/hRIhCpIe+/D0895QvDGjSAn/0MRo/euhLn0qV+7KabooszU445xltHmzZFHYmU\npUQgUkP+/Ge46irYZx8YOhSeeQZuvx1OPdVrCQHceKPvN5APA6iNGkGbNmoVxFFKicDMrjGzd81s\nnpmNMrM6ZtbIzCab2RIze9HMGpQ6f4CZLTWzRWbWNfXwReJp4UKYNAmuvvrHYx07wuzZ/s34F7+A\nfv18E/q+faOLM9NOOAFefjnqKKSspBOBme0FXAm0DyH8HN/t7AKgPzAlhHAQ8AowIHF+W+A8oA3Q\nDbjXzKo9zUkkGxQW+jf9Bg22Pr7jjnDDDZ4AFi2Cv/4Vdt45khAj0aWLEkEcJb2OIJEIpgGHAuuB\np4G7gWFA5xDCSjNrChSFEFqbWX8ghBAGJ54/ESgMIWwzj0DrCCSbzZ0L3br5GEGu1AtKl2+/9W6w\nzz+H+vWjjib3ZHwdQQjhc+DvwCfAcmBtCGEK0CSEsDJxzgqgceIpzYBPS73E8sQxkaywZIl/kFXm\nT3/y2UJKAtuqV8+nkU6dGnUkUlrSFU3MrCFwJrAPsBZ40sx6AGW/yif11b6wsPC/9wsKCigoKEgq\nTpF0WLcOjjrKt4587jkvBVGeGTNgzhwYOzaz8WWTLl1gyhQ4/fSoI8l+RUVFFBUVpfw6qXQNnQOc\nFEL438TPPYFOwC+BglJdQ6+GENqU0zU0CRioriHJBrfdBgsWwN57exmIiRPL34P3pJPg7LOhd+/M\nx5gtZs706qrvvht1JLkn2a6hVGocfgJ0MrO6wCagCzAT2AD0AgYDFwHjE+c/C4wyszvxLqEDAC04\nl9hbt85LQLz2GrRu7buGHXOMtwxK7yf82mu+LiDXSkinW/v2PkawYgU0bRp1NAKpjRHMAMYBc4B/\nAwbchyeAE81sCZ4cBiXOXwiMBRYCE4A++tov2WDYMN9EvnVr//nyy+Hee31AeNIkPxaCLwobOBDq\n1Iku1mxQu7avttbsofhQ9VGR7Vi3Dg44AF5/3auDlvbWW94NdPvt3kq48krv7si1vQRqwj/+4Wsq\nHnoo6khySxRdQyI57557vN+/bBIAHzyeOtVbBhs2+LlKAlXTpYtvwBMCaDVR9NQiEKnAunVeCO6N\nN8pPBCVWrIDhw33aaC0VbamSEKBFC3j1Vd/PWNJD+xGIpNk998DJJ28/CYAPeBYWKglUh5lWGceJ\n/umKlGPtWp8p9Mc/Rh1J7ipZTyDRUyIQKcc993jff6tWUUeSu7p08a6h4uKoIxENbYmUsXatl41+\n882oI8ltzZpB48Zem6l9+6ijyW9qEYiUodZA5pxwgrqH4kCJQKSUktaAxgYyQwPG8aBEIIL3U0+Y\n4LuHnXaapjRmSkGBL8zT9pXRUiKQvLZxo68BaNvWS0T07g3/+lfUUeWPhg392mv7ymgpEUhe+s9/\n/IN/n328kujw4V7yoGdP1QrKNHUPRU+JQPLKl1/6PsLt2sGaNd4tMX68d1Go1EE0tJ4gekoEkhe+\n+w7+9jevILp5Myxe7FVFNRYQvaOP9mJ969ZFHUn+UiKQnFZcDKNHQ5s2XkH09de98mXjxpU/VzKj\nbl3o2BHSsNGWJEkLyiRnvfWWdwMVF8OIEd79I/F0+ukwZgyccUbUkeQnVR+VnPTNN75HwD33wIUX\nqiBc3K1ZAy1b+nage+0VdTTZS9VHRUoZPx6OPBJ+/WslgWzQsKEn7OHDo44kP+m/iOSkxx/3DxbJ\nHldcAffdp8VlUVAikJyzerVvJnPWWVFHItXRpg0ccgg88UTUkeQfJQLJOU8+6UXjdt016kikuq66\nCu6+23cwk8xRIpCcM2oU9OgRdRSSjG7dfOBYJScyS4lAcsrHH8OSJdC1a9SRSDJq1YIrr/RWgWSO\nEoHklNGj4ZxzVC8om/XqBZMnw/LlUUeSP5QIJKdotlD2a9BAU0kzTQvKJGfMm+crVD/6SGsHst3i\nxdC5Myxb5iUopGq0oEzy3uOPwwUXKAnkgtat4bDDNJU0U/RfRnJCcbEnAs0Wyh1XXeXbhqpzoOYp\nEUhOePNN71s++OCoI5F0OflkWL/eiwdKzUopEZhZAzN70swWmdkCM+toZo3MbLKZLTGzF82sQanz\nB5jZ0sT5muAnaaO1A7lHU0kzJ6XBYjN7GJgaQhhhZjsAuwA3AF+GEIaY2fVAoxBCfzNrC4wCjgCa\nA1OAA8sbFdZgsVTH9997xcrZs33rSckd69bBvvvCO+/4n7J9GR8sNrP6wLEhhBEAIYQtIYS1wJnA\nyMRpI4GSii9nAGMS530MLAU6JPv+IiVefNHr1CgJ5J769X2soH//qCPJbal0DbUEVpvZCDN7x8zu\nM7OdgSYhhJUAIYQVQMleUM2AT0s9f3nimEhKtHYgt117rY8Bvflm1JHkrlR2KNsBaA9cHkKYZWZ3\nAv2Bsn06SfXxFBYW/vd+QUEBBdpeSsqxfj1MnOgb0Ehu2mUXGDQI+vaFGTM0Pbi0oqIiitKwx2fS\nYwRm1gSYFkLYL/HzMXgi2B8oCCGsNLOmwKshhDZm1h8IIYTBifMnAQNDCG+X89oaI8gza9bA+edD\np05w7rnQrh1YJT2dxcW++nTiRHjuuczEKdEIAY46Cnr39hIUUr6MjxEkun8+NbNWiUNdgAXAs0Cv\nxLGLgPGJ+88C3c2sjpm1BA4AZiT7/pJbbr3Vy0Zv2ACnnAJt28Kf/gTz5/84j7y4GN59F4YN83pC\nTZrAHXdAv37Rxi41zwzuugtuvNH/jUh6pTpr6BDgAWBH4EPgYqA2MBZoASwDzgshrEmcPwC4BNgM\n9A0hTK7gddUiyCMffAAdO/qHfNOm/sE/YwaMHQvjxkG9er7StGStQEGB3zp3hhYtoo5eMqlnT9h7\nb//iINtKtkWgWkMSuXPOgfbt4YYbtn2sJCl8+CEcc4w++PPdZ5/5LmazZ2s6aXmUCCQrvf66bzC/\neLF/8xepzJ//7K3HsWOjjiR+lAgk6xQXe5dQv35eLE6kKjZu9K7Cxx6D446LOpp4UfVRyTqjRkHt\n2tC9e9SRSDbZeWcYPBiuvhp++CHqaHKDEoFEYuNGHxO4447Kp4mKlNW9u+9TMHJk5edK5dQ1JJG4\n5RZYsADGjIk6EslW06f72pOPP9aXiRIaI5Cs8fnn8POfw6xZmvkhyQsBWraECRN83YlojECyyE03\nwe9+pyQgqTGDrl19o3tJjRKBZNScOV4Sorw1AyLVdeKJ8NJLUUeR/dQ1JBmzahWcdBJceqnXjBFJ\n1Vdfecvyiy9gp52ijiZ66hqSWHv/fS8advrp8PvfRx2N5Iqf/MTXFEybFnUk2U2JQGrcjBlw7LFw\n3XW+KlQzPCSdNE6QOiUCqVHPPw+nngr336+WgNQMjROkTmMEUmPuuw8GDoTx46GDNiWVGvL997DH\nHl6YcPfdo44mWhojkNgIwfcSGDLEi8opCUhNqlPHux5ffjnqSLKXEoGk1YYNvn/wpEnw1ltwwAFR\nRyT5oGtXdQ+lQolA0mbJEq8mWq8eTJ0KjRtHHZHkixNP9AFj9SgnR4lA0uKpp7x5fvXV8OCD2ltA\nMqt1a69EunRp1JFkpx2iDkCy25Yt0L+/J4IJE+Dww6OOSPJRSbmJl16CVq0qP1+2phaBJG3FCujS\nxauIzpqlJCDRKukekupTIpCkLF3qH/zHH+9rBfJ92p5E74QTfGxq8+aoI8k+SgRSbevXw1lneRXR\nwkLfZUwkanvsAfvtB2+/HXUk2UeJQKqluBguugiOPtqLx4nEiaaRJkeJQKrlttt8bOCee6KORGRb\nGidIjkpMSJW98IKXj545E/bcM+poRLb13Xe+fuWTT6Bhw6ijyTyVmJAatWQJXHwxPPmkkoDEV926\nXu781VejjiS7KBFIpdat88HhW2+FI4+MOhqR7VNZ6upT15BsV3ExnH02NG0Kw4dHHY1I5ebPh1/9\nyjdDyjfJdg1pZbEQAnz7rX/zL7mtXet/TpkCq1fD2LFRRylSNT/7GXzzjZel3m+/qKPJDiknAjOr\nBcwCPgshnGFmjYAngH2Aj4HzQghrE+cOAH4LbAH6hhDUgIuB//1feOwxH1yrX99vDRr4n7vvDuPG\nealfkWxg5rOH7r0XjjsONm3yQeRNm/xWXAw9e+bnYHJF0tEi6AssBOonfu4PTAkhDDGz64EBQH8z\nawucB7QBmgNTzOxA9QFF6623vGT06tWw665RRyOSHhdfDLff7pMcdtpp69u8eb7p/cCBUUcZHymN\nEZhZc2AEcCvQL9EiWAx0DiGsNLOmQFEIobWZ9QdCCGFw4rkTgcIQwjbrADVGkBnFxb5pzDXXQI8e\nUUcjkhnvvAPnnONjCLVybLpMVNNH7wSuBUp/ajcJIawECCGsAEqq0jcDPi113vLEMYnIiBH+DenC\nC6OORCRzDjvMW7+vvx51JPGRdNeQmZ0KrAwhzDWzgu2cmtRX+8LCwv/eLygooKBge28h1bV2rdcK\nev5571MVyRdm0KsXPPwwdO4cdTSpKSoqoqioKOXXSbpryMxuA36ND/zWA3YD/g84HCgo1TX0agih\nTTldQ5OAgeoaisYf/uDJ4IEHoo5EJPNWrvTNbD79NLfGxjLeNRRCuCGEsHcIYT+gO/BKCKEn8BzQ\nK3HaRcD4xP1nge5mVsfMWgIHADOSfX9J3uLF8MgjXjdIJB81aeI76o0bF3Uk8VATQyWDgBPNbAnQ\nJfEzIYSFwFh8htEEoI++9mdeCD44PGCA9hSW/FbSPSRaWZx3nn8e/t//8yl0Whsg+ez776FZM9+/\nIFcWnqnonFRq0yZvDdx5p5KASJ06PmPukUeijiR6SgR5ZOhQOOgg6NYt6khE4qFXLxg50tfU5DMl\ngjyxYAEMGeKtARFxhx7qpVReey3qSKKlRJDjVq2CPn2goAD++lc48MCoIxKJj5I1BSNGRB1JtJQI\nctTGjT49tG1bXz28eLHXXxGRrfXoAePHw/r1UUcSHSWCHFNc7H2eBx0Ec+bA9OneHbT77lFHJhJP\njRt7izmf1xRo+mgWC8E3kl+40McAFi70+in168Pf/+5b9olI5Z55xr8wTZ0adSSpSXb6qBJBFnrg\nAf/Wv2CBV09s1867gNq1g0MOgWOOUf0gker4/nto3hymTYP99486muQpEeSJzz7zD/vRo33Gwx57\n6ENfJB2uvto3ZLr55qgjSZ4SQZ647DLv+hk8OOpIRHLL3Llwxhne0t5tt6ijSY4SQR746CM4/HB4\n7z0N/orUhD59YNYsmDgxO/+PqcREHvjzn+GKK7LzH6hINvjHP+D4430W0X/+E3U0mZOOPYslA5Ys\n8YJxS5dGHYlI7jKDQYN8rODYY2HKFNh336ijqnlKBFmisBD69YOGDaOORCS3mcENN/g4wXHHweTJ\nvolNLtMYQRaYPx9OPNE3286l3ZRE4m7kSOjfHyZM8L2O4y7ZMQK1CLLAwIFw/fVKAiKZdtFF/v/u\npJPg6ad9jU4u0mBxzM2eDTNmwKWXRh2JSH76n/+BRx+Fs8+Gf/3LV/TnGnUNxdwpp8Bpp/m0NhGJ\nznvveVJo3x7++U/YeeeoI9qWpo/moDff9PpBl1wSdSQi0qqVF3HcsgWOPNLH7HKFEkGM/fGP8Kc/\neRlpEYneLrvAY4/B73/vRR2ffTbqiNJDXUMR+Pxzr3JYVOSlomvX9g/7klvdul5O+t//hkWLYAcN\n6YvEzrRpcP750LOnL/asXTvqiFRiItZWrYKXX/YP/qIi+OIL6NzZVy8efrjPW960Cb77zv8suXXs\n6PsKiEg8rVoF553nhSCHDo06GiWC2Jo3z9cAHHWUf/AXFMDBB3v5aBHJfl9/7YngwQf9/3qUlAhi\n6IMPfGXiHXd4E1JEctOUKb4V7L//DT/5SXRxKBHEzOefe62S666D3r2jjkZEalrfvt5VNHp0dDFo\n+miMfP21r0S85BIlAZF8MWiQ72kQZSJIlloEafbNN95PeOSR8Le/afcwkXwya5YvAn3nHd/6MtPU\nNRQD33/vOxztuSc89JCSgEg+uuUWeO01ePHFzE8KyXjXkJk1N7NXzGyBmc03s6sSxxuZ2WQzW2Jm\nL5pZg1LPGWBmS81skZl1Tfa94+iHH+A3v/E1APffryQgkq8GDID1632Tm2yRdIvAzJoCTUMIc81s\nV2A2cCZwMfBlCGGImV0PNAoh9DeztsAo4AigOTAFOLC8r/7Z2CK49toft7irWzfqaEQkSkuXevfw\n669DmzaZe9+MtwhCCCtCCHMT9zcAi/AP+DOBkYnTRgJnJe6fAYwJIWwJIXwMLAU6JPv+cfLCCzB2\nLDz1lJKAiMCBB8Jf/uKrjjdvjjqayqWlB8vM9gUOBaYDTUIIK8GTBdA4cVoz4NNST1ueOJbVli/3\n2UGjRkU7f1hE4qV3b9/y8v77o46kcikngkS30Digb6JlULZPJ7v6eKrhhx+gRw/fUD5XN6wQkeSY\n+czBW27xMYM4S6mcmZntgCeBR0MI4xOHV5pZkxDCysQ4wqrE8eVAi1JPb544Vq7CwsL/3i8oKKCg\noCCVUGvEX/7ihaYGDIg6EhGJo8MOgxNO8IRw883pf/2ioiKKiopSfp2Upo+a2SPA6hBCv1LHBgNf\nhRAGVzBY3BHvEnqJLB4snjoVunf3+cJ77hl1NCISV8uW+WY2775b858VGV9HYGZHA68B8/HunwDc\nAMwAxuLf/pcB54UQ1iSeMwC4BNiMdyVNruC1Y50IVq/2TH///XDyyVFHIyJxd+21sG6db3VZk7Sg\nLENC8EVjbdrAkCFRRyMi2eDrr32Hs9deq9nppKo1lCFDh3phqb/8JepIRCRbNGoE/fv7LY7UIqiG\nkjoib78NLVtGHY2IZJPvvoPWreGRR7w8fU1Qi6CGffUVnHsuDB+uJCAi1Ve3Ltx6q5emj9v3XCWC\nKigu9hWCZ5/tNxGRZFxwgRenfOqpqCPZmrqGquC222DCBHj1Vdhxx6ijEZFs9vLLcOmlsGAB1KmT\n3tdW11ANefllGDYMnnhCSUBEUtelC+y/P9x3X9SR/CillcW5bvly7xJ69FFolvVVkUQkLoYMgV/+\n0gvSXXWVVyiIkrqGKrB5Mxx/PHTrBjfeGGkoIpKD3n8ffvtbH4N86CFfZ5AqdQ2l2YABXjlQdYRE\npCYccAAUFcH558NRR8Edd3ghyyioRVCOp5+Gfv28jpBKS4tITfvgA28dbNkCI0Yk3zpQiyANZs6E\niy+G3/0OnnxSSUBEMmP//X1WYvfu3jp4+OHMvn/etwi+/dZnBN17L3zxBVx2mWfmn/40I28vIrKV\nJUt8fHLYsOqvW1LRuWrYssXn8D72mGfeDh2gTx+vJBr16L2IyDvvwEknwbhx0Llz1Z+nRFCB777z\nOuBz5vjFfecd/7lZMzjzTF/Ysf/+aXkrEZG0mTIFLrzQ1zIdfHDVnqNEUI6hQ33Wz4EH+sYQhx3m\nfx5yCOy2WxoCFRGpQWPG+F4Gb7wB++xT+fnJJoKcXVD23HO+aGPxYth776ijERGpvu7dYcUK77Z+\n4w3YffeaupswAAAFgklEQVSaeZ+cbBHMm+fLuJ9/Hjp2TGNgIiIRuP5639RmyhTYZZeKz1PXUMLK\nlf7hP2iQZ1MRkWwXAvTqBV9+Cc88AztU0JejdQT4wPBZZ/kFUxIQkVxhBg884J9xf/tbDbx+rrQI\nQoBf/9qnho4eDbVyKsWJiMCHH/p09zlzoEWLbR/P+xbBrbfC0qW+LkBJQERy0X77wRVXeAmcdMqJ\nFsG4cX5h3n4b9tyzBgMTEYnYt99Cu3a+bW7Xrls/lrctgmee8bIQ48crCYhI7qtXD+6+G668EjZt\nSs9rZm0iCAFuv90vxsSJvlhMRCQfnHYaHHSQl65Oh6zsGvruO68QumSJtwi0e5iI5JuPPoIjjvCy\nOSWLZvOma2jFCigo8B3Epk5VEhCR/NSypW9zmY6B46xKBHPm+GKxU07xGhw77xx1RCIi0bnuOpg7\nF158MbXXyZquoaefht69fd+Ac8+NKDARkZh54QW45hqYPx/q1s2SriEzO9nMFpvZe2Z2fVWfN306\nTJqkJCAiUtqpp0LbtvD3vyf/GhlNBGZWCxgGnAS0Ay4ws9ZVee6QIfCLX9RkdOlTVFQUdQhJy+bY\nQfFHTfFH4667fK/jZGW6RdABWBpCWBZC2AyMAc7McAw1Llv/MUF2xw6KP2qKPxr77uu7LiYr04mg\nGfBpqZ8/SxwTEZEU1KmT/HOzataQiIikX0ZnDZlZJ6AwhHBy4uf+QAghDC5zXvymMomIZIHYb0xj\nZrWBJUAX4D/ADOCCEMKijAUhIiJbyeiexSGEH8zsCmAy3i31oJKAiEi0YrmgTEREMieyweKqLCwz\ns7vNbKmZzTWzQzMd4/ZUFr+ZdTazNWb2TuJ2UxRxlsfMHjSzlWY2bzvnxPnabzf+OF97ADNrbmav\nmNkCM5tvZldVcF4sfwdViT+uvwMz28nM3jazOYnYB1ZwXlyvfaXxJ3XtQwgZv+EJ6H1gH2BHYC7Q\nusw53YAXEvc7AtOjiDWF+DsDz0YdawXxHwMcCsyr4PHYXvsqxh/ba5+IrylwaOL+rvi4WTb9+69K\n/LH9HQA7J/6sDUwHOmTLta9i/NW+9lG1CKqysOxM4BGAEMLbQAMza5LZMCtU1YVx1R69z4QQwhvA\n19s5Jc7XvirxQ0yvPUAIYUUIYW7i/gZgEduup4nt76CK8UNMfwchhI2Juzvh46Rl+8dje+2hSvFD\nNa99VImgKgvLyp6zvJxzolLVhXFHJpqWL5hZ28yElhZxvvZVlRXX3sz2xVs3b5d5KCt+B9uJH2L6\nOzCzWmY2B1gBvBRCmFnmlFhf+yrED9W89hmdNZRnZgN7hxA2mlk34BmgVcQx5YusuPZmtiswDuib\n+GadVSqJP7a/gxBCMXCYmdUHnjGztiGEhVHHVVVViL/a1z6qFsFyYO9SPzdPHCt7TotKzolKpfGH\nEDaUNOFCCBOBHc3sJ5kLMSVxvvaVyoZrb2Y74B+ij4YQxpdzSqx/B5XFnw2/gxDCOuBV4OQyD8X6\n2peoKP5krn1UiWAmcICZ7WNmdYDuwLNlznkW+A38d0XymhDCysyGWaFK4y/dp2hmHfCpul9lNszt\nMiruR4zztS9RYfxZcO0BHgIWhhCGVvB43H8H240/rr8DM/upmTVI3K8HnAgsLnNabK99VeJP5tpH\n0jUUKlhYZma9/eFwXwhhgpmdYmbvA98AF0cRa3mqEj9wjpldBmwGvgXOjy7irZnZ40ABsLuZfQIM\nBOqQBdceKo+fGF97ADM7GugBzE/09QbgBnwWWux/B1WJn/j+DvYERpqXxK8FPJG41lnx2UMV4ieJ\na68FZSIieU7VR0VE8pwSgYhInlMiEBHJc0oEIiJ5TolARCTPKRGIiOQ5JQIRkTynRCAikuf+PwAl\nDW7IkOHfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff0fdf31e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "val,bi = np.histogram(dis,bins=50)\n",
    "plt.plot((bi[:-1]+bi[1:])/2.,val)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:mylab3]",
   "language": "python",
   "name": "conda-env-mylab3-py"
  },
  "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.5.2"
  },
  "widgets": {
   "state": {
    "935f20d6242f42228ed24ab01e39a96a": {
     "views": [
      {
       "cell_index": 62
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}