From 94b1bcbb4978b3830cd56ca2ac3837bc2a67e180 Mon Sep 17 00:00:00 2001 From: ggandus Date: Tue, 12 Mar 2019 11:40:39 +0100 Subject: [PATCH] add solutions 3-a and 3-b --- Exercises/Exercise-3/Solution3-a.ipynb | 14551 ++++ Exercises/Exercise-3/Solution3-b.ipynb | 84103 +++++++++++++++++++++++ 2 files changed, 98654 insertions(+) create mode 100644 Exercises/Exercise-3/Solution3-a.ipynb create mode 100644 Exercises/Exercise-3/Solution3-b.ipynb diff --git a/Exercises/Exercise-3/Solution3-a.ipynb b/Exercises/Exercise-3/Solution3-a.ipynb new file mode 100644 index 0000000..6c4232c --- /dev/null +++ b/Exercises/Exercise-3/Solution3-a.ipynb @@ -0,0 +1,14551 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import matplotlib.animation\n", + "from IPython.display import HTML\n", + "import numpy as np\n", + "from scipy.constants import physical_constants\n", + "import copy" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def get_ani(traj,dx,dy,frames=None):#(xdata,ydata,xlim,ylim,frames=None):\n", + " \n", + " xdata,ydata,xlim,ylim=traj_coord(traj,dx,dy)\n", + " if frames is None:\n", + " frames=len(ydata)\n", + " \n", + " fig, ax = plt.subplots()\n", + " ax.set_xlim((xlim))\n", + " ax.set_ylim((ylim))\n", + " l, = ax.plot([],[],'o')\n", + " \n", + " def animate(i):\n", + " l.set_data(xdata[i],ydata[i])\n", + " \n", + " ani = matplotlib.animation.FuncAnimation(fig, animate, frames=frames)\n", + " plt.close()\n", + " \n", + " return ani" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def traj_coord(traj,dx,dy):\n", + " X=[]\n", + " Y=[]\n", + " Xmax=-np.inf\n", + " Ymax=-np.inf\n", + " for i in range(len(traj)):\n", + " x,y=coordinates(traj[i],dx,dy)\n", + " xmax=max(x) \n", + " Xmax=xmax if xmax>Xmax else Xmax\n", + " ymax=max(y) \n", + " Ymax=ymax if ymax>Ymax else Ymax\n", + " X.append(x.copy())\n", + " Y.append(y.copy())\n", + " xlim=(-10,Xmax+10)\n", + " ylim=(-20,Ymax+20)\n", + " return X,Y,xlim,ylim" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def traj_stat(traj):\n", + " nouter=len(traj)\n", + " nalone=np.zeros(nouter)\n", + " ndimer=np.zeros(nouter)\n", + " ncluster=np.zeros(nouter)\n", + " #scluster=[list() for i in range(nouter)] \n", + "\n", + " for i in range(nouter):\n", + " reduced=copy.deepcopy(traj[i])\n", + " while len(reduced)>0:\n", + " cluster=allconnected(reduced,0,nx,ny) \n", + " members=len(cluster)\n", + " if members==1:\n", + " nalone[i]+=1\n", + " elif members ==2:\n", + " ndimer[i]+=1\n", + " else:\n", + " ncluster[i]+=1\n", + " #scluster[i].append(members)\n", + "\n", + " #### I will not check molecules already found in a cluster\n", + " for c in cluster:\n", + " reduced.remove(c)\n", + " return nalone, ndimer, ncluster" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import monte_carlo.py functions\n", + "\n", + "> coordinates
\n", + "neighbors
\n", + "allconnected " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "from monte_carlo import *" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Tasks" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "1. Estimate the dimer binding energy as $de=k_BT \\ ln(n_0/n_{exp})$ where $k_B$ is the Boltzmann's constant, T is the simulation (and experiment) temperature, $n_0$ is the concentration of dimers in the case of zero interaction and $n_{exp}$ is the concentration of dimers found in the experiment.\n", + "\n", + "\n", + "2. Repeat the simulation using as $dE$ your estimate. What do you get?\n", + "\n", + "\n", + "3. Repeat the simulation with coverage 0.1 and $de=-0.02$ then $de=-0.1$ Describe what you obtain. Now try $coverage=0.1, T=400, dE=-0.1$ Comment the result" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Monte Carlo parameters\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "de=-0.1\n", + "T=400\n", + "kb=physical_constants['Boltzmann constant in eV/K'][0]\n", + "beta =1.0/(kb*T)\n", + "nouter=100\n", + "ninner=10000" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Honeycomb lattice define" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "sr3=np.sqrt(3.0)\n", + "nx=50\n", + "ny=nx/2\n", + "dx=10.02\n", + "dy=dx*sr3" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "The figure illustrates the hexagonal lattice. The blue dot is positioned at site $si=0$ of cell [$xi$,$yi$]=[$2,2$] and the red dot at site $si=1$ of the same cell." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## #Molecules" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "coverage=0.1\n", + "nmolecules=int(nx*nx*coverage)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run ...\n", + "\n", + "This Monte Carlo implemets an **outer** and an **inner** cicle, i.e. for each outer step, an inner loop is runned.
Statistics for the system of molecules are sampled for each outer step, by performing statistical averages over the configurations obtained from the inner circle." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "molecules=[]\n", + "i=0\n", + "\n", + "#### Create initial geometry\n", + "while i < nmolecules :\n", + " xi=np.random.randint(0,nx)\n", + " yi=np.random.randint(0,ny)\n", + " si=np.random.randint(0,2)\n", + " newmolecule=[xi,yi,si]\n", + " if newmolecule not in molecules:\n", + " molecules.append(newmolecule)\n", + " i=i+1\n", + "#### END Create initial geometry\n", + "\n", + "#### Initial energy\n", + "totene=0.0\n", + "for m in range(nmolecules):\n", + " n,l =neighbors(molecules,m,nx,ny)\n", + " totene+=n*de\n", + "totene/=2\n", + "\n", + "nacc=0 # acceptance rate\n", + "nrej=0 # rejection rate\n", + "avgene=np.zeros(nouter) # average energy\n", + "traj=[] # trajectory of system of molecules per outer step\n", + "\n", + "#### outer loop\n", + "for i in range(nouter):\n", + " \n", + " maxinner=ninner*20\n", + " nsteps=0\n", + " j=0\n", + " while (j0:\n", + " avgene[i]+=totene\n", + " j+=1\n", + " nsteps+=1\n", + " if nsteps==maxinner: \n", + " print('coverage too big')\n", + " exit()\n", + "\n", + " #### STATISTICS and LOG\n", + " avgene[i]/=ninner\n", + " traj.append(copy.deepcopy(molecules))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Count #alone,#dimer,#cluster\n", + "\n", + "> #alone = unpaired molecule && not in any cluster
#dimer= number of dimers
#cluster= number of clusters" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "nalone, ndimer, ncluster = traj_stat(traj)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Check convergence\n", + "\n", + "1. Visually check that the energy **and** the #dimer converge. \n", + "2. If **not** converged goto **Monte Carlo parameters**, incerase *nouter* and **Run ..** again the simulation" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAEyCAYAAAA4HuM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX9//HXmS0bISQkrEEWBWRHRJCqxQUVV0RLtT9rBdtaFfXrt61ft9a1WluprXWr1rVq64IbtVoV6lIVwaCoLIIICCFAQkL2ZDLL+f1xJyFAVmZCmMn7+Xjkkcy9d+499+bOfO7nnHPPNdZaRERE5MDm6uwCiIiISOsUsEVEROKAAraIiEgcUMAWERGJAwrYIiIicUABW0REJA4oYIuIiMQBBWwREZE4oIAtIiISBzydXYDGsrOz7aBBgzq7GCIiIvvNsmXLdlhrc1pb7oAK2IMGDSIvL6+ziyEiIrLfGGO+bctyqhIXERGJAwrYIiIicUABW0REJA4oYIuIiMQBBWwREZE4oIAtIiISBxSwRURE4oACtoiISBxQwBYREYkDCRuwA9sLKXnmGUJlZZ1dFBERkaglbMCu+3Yj22/7DTXLl3d2UURERKKWsAE7ZfRocLmo+fzzzi6KiIhI1BI2YLtSU0kaPlwZtoiIJISEDdgAKePHUfP5F9hQqLOLIiIiEpXEDtjjxhGuqqJu/frOLoqIiEhUEj5gA2rHFhGRuBeTgG2MecwYU2iMWdFoWpYx5m1jzNeR35mx2FZ7+AYNwp2RQbXasUVEJM7FKsN+Api+x7RrgUXW2qHAosjr/coYQ/L4cdQqwxYRkTgXk4BtrX0fKNlj8gzgycjfTwJnxWJb7ZUybhz+dd8QqqjojM2LiIjEREe2Yfe21m4FiPzu1YHbalbq+PFgLbVfftkZmxcREYmJTu90Zoy52BiTZ4zJKyoqivn6k8eMAWPUji0iInGtIwP2dmNMX4DI78KmFrLWPmytnWitnZiTkxPzQrjT00k65GD1FBcRkbjWkQF7AXBh5O8LgVc7cFstSh43jtrln2Ot7awiiIiIRCVWt3X9A1gMDDfG5BtjfgzcCZxojPkaODHyulOkjh9PqKyMwLffdlYRREREouKJxUqstT9oZtYJsVh/tOoHUKlevhzfoEGdWxgREZF90OmdzvYH38EH4+rWTe3YIiISt7pEwDYuFyljxypgi4hI3OoSARucJ3f516wlXF3d2UURERFpt64TsMeNg1CIiv+809lFERERabcuE7BTjziCpKFDKbj6aoruv1/PyBYRkbjSZQK2KzWVQc89S8aZZ7Dj3vvY/NOLCRYXd3axRERE2iRxA3ZlITzzfSjZ0DDJlZpK3zvvpM9tt1Kdl8eGmWdTvWxZJxZSRESkbRI3YBevg80fw0NTYfVrDZONMWTOmsWg557FpCTz7Y8upPiRR7DhcCcWVkREpGWJG7AHfgd+9j70HALPnQ9v3gChQMPs5BEjGPzii6SfeCKF8/5A/tzLCZWWdmKBRUREmmcOpPG1J06caPPy8mK70qAf3vo1LH0Ico+A//c8pGY1zLbWsvOZv7P9d7/Dm5ND+inTMcY0zE875rukTZ4U2zKJiIhEGGOWWWsntrpcwgfseitehPkXwQk3wjG/2Gt2zZdfUnDtdQTy8xum2VAIQiGyL72U7LmXYdzujimbiIh0WQrYTfnrCRAOOFXlbRCuqWHbrbdR9vLLpE45kv533YUnO7vjyiciIl1OWwN2TB7+ETdGnQVv/crpOZ41uNXFXSkp9PvtHaROnMi2W29l/cyZ9JxzEcbrjVmRXKkpdD/tNFzJyTFbp4iIJJ6ulWHv/BbuGQvTboGjr2rXW2vXrGHLVf9L3YYNrS/cTknDh5N7z5/0JDERkS5IVeLNefg4wMLF77b7rTYUIlReHtPi1Hy2nK3XXYcNBul7+2/oPn16TNcvIiIHNlWJN2fkDFh4k5NtZw7cfd6OddC9H/hSm3yrcbvxZGbGtDjpxx9H8ssvseV/f86Wq/6X6v+3lNQpU2K6jf2uNB/SeoI3pdlFUkaOxNu/v/OiZifUVUNG//1UwARVVQzBWh1HkQTV9TLskg3w5/Fw4m1w1JW7pm9bAQ8fC4OOggtegUa3du0Ptq6Owj/cTcmTT+7X7XYWk5pK31tuIWP6CfCXo52A/T+fg8fX2UWLT3XV8JejoLYM5i6FNHWOFIkXqhJvyUPfBZcXfrrIeR0KwCMnQOFqCNXBGX+Gwy/s+HI0oS4/n3BlZadsOyaWPASf/g3cXvjhS7vd817P+v1s//1d1CxbRo+jhtC73we43MDMh2Hcufu/zIng39fDx/eDywMjzoRZj3d2iUSkjVQl3oRtVdvondobM3IGLLoVSjdDjwHw4T2w9XOY9SR88ojTk/yQaa1XLdZVO1WQTQSlJoXDUJ4PPQ5qdhFfbu7uE8oLoFsfcDUzKF11yW7jpQPQa0Sz1fodqq4aXn0ZRo2Bgs+g7D2YcF2Tiw584nGKfnM9xc+9Rk2/4fSeFMD1yp+AQ/d77Ubc274S/vVXGHEOpGbD+4+B734YfMzuy4UCUFMK3XL2X9lqypyLiKS0/bfNWAvUQqAGUmPbHNaimjKo2Lr7tKwh+14DZa3zv29pH6p37t99jGPJI0diPPs/fHaJDLsuVMddn9zFs2ue5eeH/5w5fb8L906Ak++Ag493Mu7hp8L3n3SC34PfgUFHO6OiNRc8wiF4bLqTlc+4F0bNbLkQlUXw0k9h/Ttw5FyYdnPLH75wCN79Lbw/D6bMhZNv33uZqh3wwBSoKtx9eq9RcPE74ElquUyxlvc4vHYVzH4dPvozbFkGV60AbxO3rAX98NB3qVhbTsHi7oRj3JlPRKSjDPtkKe709JitTxl2RH5FPr947xesKl5Fbrdc7vvsPqYOmMqQ3mNgxUvOCGhJ6XDqPOcNWYPhhJvg39fA58/C+B80veKPH4D8pZA5GF6YDd9+BCf9pukg+e1Hzihr1SVw6OlO1WX+Uvje406Gv6eK7fDij2Hjf531L77fqeY8aPLuy71+tdNh6+xHIDnDmVay3in7+3fB8b/a5+PWbtbCxw9Cn7HOOO7hAPxtBqyYD4f9cO/l3/s9FH1F+uXzOfjXh1H7+afw8s+g92j47i/3X7nj3WfPwKpXnP9133HOtJ0b4Y1rnYvO71wOmz52/jfGQM+hsHU59J/oXAgmdYt9mcIhp1xf/dPJCgM1ULkdxv0ARswAVxzUoFgLX78Nyx53PltpvaBotXOBP/HHHdfXwlp4907YvgKmXL7r+6R0Eyx/Bsb/EEbNaN86i9Y6tYZYyMiFU37vNFnVqylzLrTrKsGTAqffrT4QreiscTMSOmC/s+kdbvjwBrDwp+P+xLiccZz16lnc+OGNPDnyDNzv3OEseM6ju1cTTrrY+RL89zVw8HGQ3mf3Fe9YB//5jZOVz3oSFt4cCcJ5MOM+SGlURf7Fc071e+ZAp828zxhY+TK8egU8dAzMuB/6Tdi1fOFKeOUyqC135o2c4WTRr86FS/67q+f16n/CypfguF/B2Fm7l2/rcvjv3TDijF1f4s2x1sl2m8qC2+ObRbBjDcx8yAkMg6c6mf7iB2D8+bvXVBQshw/+COP+Hww9EQ/Q7YSTwF7oNE+MvRcyB0VXnj0Falvex0Ctc/HTWFoOuA+gj8ieZdyxFt59AU7/IZzZqAMlUyFjC7z3Oyjwwjf/gvGHw6wnIGMALH3YeRjOqs1w1l+coBortaWw4EooXwpnXexcxAb9sOAKWPUY9CiAsx5sezNSU1r7X0bLXwH/vAq2zIejpzl9K5IznBqv/86DtZudz2a33m1fpyepbfu8/B/AYrjwTjjy0l3TrYXnCuHrp2HkJZAzrG3bDdTCQ7+EETlw4q0wfw7Yj2DqjbuWee4C6F0J5z4NL1wIJX+HU+a3v2kqWAfVxbtPS8ve/eKgvaL9XwdqnWMfq2Y2a6HoK4jh4Fnt0eFV4saY6cA9gBt4xFp7Z3PLxrJK/J1N73DlO1cysudI5k2dx4B0J5P95zf/5PoPrufqEXP40eu3OBnvuU/v/Q8t/sapGs8eBt//266R0cJheOJUKFzl9MatD+ar/+kEWn8TVbsjZ8CZ90Fy993X//yFsP3LvZfvOdSpnu89ynm9bhE8fTYcdRWceIuTqd8/2dn2T/+z9weiugQeONLJCi5+p/kPTHUJvHIpbHgfTru7+dqEtnjqbCcruGrFruzj06dgweXwowUwZKozbcsyZ79DAZj7MaQ0ajMr2+IMbDPpZzD9jn0vS2PWOv0S3rweDj3N6VDY+P8AzvF96WKo3rH79OxhzgVZ75GxKUs01i2MlHGPL8Tu/eGyxbtqWOoF65y7HgpXwuRLnLsiGmeF+XlOzVDZ5tiX1dcNzvwzjD5n1zRrYelfnf9Deh/n4iG31RrA3QX9zvvzHoep1zg1Ma4Yj++/faVzfpZ8A8fdAEf/fPf+I1+/7fwfakrav+4jL3MGbWouO6/YBvdPgpxDYc4be+9bxXZnfvYwuOjfbdv3hTc7F8c/fAkOOQFevtRJIn76H+g3Hla+4gTpE26CY37udBp94/9gxgNw2Plt37dvP4L5P4aKgt2n9xjofJf1O6zt6wKnVuaNa+Czp+H4G+Co/22+H09zvnrd+X7rNRK+96hzy240akqdxGndQrjs4zaNltlWB0QvcWOMG1gLnAjkA58AP7DWrmpq+VgG7EAowFOrn+L8EeeT5N5VTW2t5Yr/XMGSrUt4cfwvOeiQ6Xt/2dVb+xa89BOwwFn3Oxnrx39xMu+zHoTx/2/35Us3wTf/cb6c6qX3gWHTm77CC9TAqgUQqN41zZPkbCdpj/aRVy93qsR+shCWPOxUNV/8rpOxN2X1a85jRY+7Aab+397z67+wK7Y5ndS2fQGHXQCn3tXi/dNNKvwKHpjsZPtTr260f7Xwx1HOF/MPnt2V2aX3hXP/1vSHeP6PYe2b8PNVewfW9vJXONneypeg73jY9qWTuX//See4hUNO9eP7dznH4IifgIl8KQRrnVoKf4VTRbjn/3p/adyXoddIOOLHu8oITg1Qc7UR5QXOeAMDm7mvv7oE1ry+22Nno2YMDDm2+TJtWeacd+Vbnex78s/alv2UbHDet3U59D/cWc+Q4+Dsv8auE91nz8C/fuGcd+c8unenvXrlBU7gtuG2r3vr5071ev+JzsXKnk1h1sJzP3SCwSUfQPbQptfz+bNO09HJdzhNGi3Z8ik8Ms05d2fc50yr2elc7KflwA9fhAePcqrJf7LIqU1qnJBctgS69215G+Gw01+lvhZxylwwkQuJUMCpMasqdMp7xE/a9r9unMzU/68PORHOfrhttRShACy6BT66F3JGON/L3hQ4569Ok8a+KPjMKVP5Fqem4sjLYto5tq0BG2tth/0AU4A3G72+DriuueUPP/xwG0vhcNiGw+G9pm+r3GanPDPFXvjGhTYUDrW8kpKN1j50rLU3dbf21Sus/U0fa586x9om1tuhakqtnXeotXcNdcryn9tbf88Lc6y9pae121bsmhYOW7v4AWf6H0dbm59nbTBg7cJbnfU+8B1ri75uX9kWXGntbb2srSzae95/bnfW+9Q5zu9nvm9tVXHz69qc5yy3+IH2lWFPW7+09s8TrL25h7Xv/8HaUMjaDR9Ye9cwp6yLH7T28dOcbb18mbX+qr3XUb5t1zKvNLNMR9pt+3P3//Y7SnWJtX8/z9mvZ3/onNstWfVPa+8YYO1vB1i7+jXnHM57wtpbc6ydN9zajR9GVx5/lXMO3NTdOd7l26JbX3NWvmLtHbnW3jnQ2jX/3n3eFy842//gTy2vIxy29ulZ1t7W29od65pfLlBr7f1HOseneufu81a/5mzrrmF7fz9Y66z3tl7WPnNuy99zVcXO5/mm7tY+9yNra8qaXubp7znLPD/b2trylvdvxUvW3t4/cozedLa/9K/W3ppt7R9GWrtpScvvL8239pETne299nNr62qsLfzK2vsmW3tThrX/ucPaULDldTS22/ZHtL79fQTk2TbE1I7OsL8HTLfW/iTy+gJgsrX28qaWj2WG/eG6HZz/yBJeuGQKRwza+6rs5a9f5saPbmzinSIiIs376Acfke5LvF7iTdUZ7HaFYIy5GLgY4KCDmr8/ub08kZ6odcGmq63OOuQsjDFsqdzS9pWWbHDagzNyW1+2oxStcarZm6vG31NlodNLtLGUHk77eFNVOrXlTlt0uB3VfcYFfcc23+O4ZL3T+7S16rV6VTucavZouNxOmXxN3P9rw87Idt37tb03bOkm2LkpujK1lzGQMzyxe+yWb3WqQFviTXHaW5tqsw36nU6M0VbrZx7U4vgIMRUKOlX7gdpd04wL+o7ZuymsOW05bmk9naaepgT9sONrZ35Tx9WGnWp8f1XL28geCult7HxXVuB8F7TEl+p0lG2qTIFap0wt/a+NcfapuarzojXOLbbtkZwBfUbv9n3pc3fOiIwdHbDzgcaNNbnAbr0SrLUPAw+Dk2HHasNej9POVxdqOvAYYzjrkLNitTkRkcTSzn5i0vHa2e2u3T4BhhpjBhtjfMB5wIIO3iYAPnckYDeTYYuIiMSTDs2wrbVBY8zlwJs4t3U9Zq1d2ZHbrOeLZNiBZjJsERGReNLho0JYa18HXu/o7exJGbaIiCSSjq4S7zReZdgiIpJAEjZgK8MWEZFEkvgBO3TgPI1MRERkXyVuwPYowxYRkcSRsAHb63ZuclcbtoiIJIKEDdgetwuXUYYtIiKJIWEDNoDX7VKGLSIiCSGhA7bP48KvDFtERBJAYgdsZdgiIpIgEjtge1xqwxYRkYSQ0AFbbdgiIpIoEjpg+zyuZh+vKSIiEk8SOmB73S7qghrpTERE4l9CB2xl2CIikigSO2C7DQF1OhMRkQSQ2AHbo05nIiKSGBI6YHvdqhIXEZHEkNAB2+fWfdgiIpIYEjpge9XpTEREEkRCB+wkDZwiIiIJIqEDtldV4iIikiASOmA7vcQ1cIqIiMS/hA7YyrBFRCRRRBWwjTGzjDErjTFhY8zEPeZdZ4xZZ4xZY4w5Obpi7huNdCYiIoki2gx7BXA28H7jicaYkcB5wChgOvCAMcYd5bbazec21AXDWKtqcRERiW9RBWxr7Wpr7ZomZs0AnrXW+q21G4B1wKRotrUvfB5n94JhBWwREYlvHdWG3R/Y3Oh1fmTaXowxFxtj8owxeUVFRTEthNft7J7asUVEJN55WlvAGLMQ6NPErBusta8297YmpjWZ5lprHwYeBpg4cWJMU+H6DFv3YouISLxrNWBba6ftw3rzgQGNXucCBfuwnqgowxYRkUTRUVXiC4DzjDFJxpjBwFBgaQdtq1n1GbZ6iouISLyL9raumcaYfGAK8C9jzJsA1tqVwPPAKuDfwFxrbSjawraXTxm2iIgkiFarxFtirX0ZeLmZebcDt0ez/mjtasNWL3EREYlvCT/SGSjDFhGR+JfQAVtt2CIikigSOmB73c7dZcqwRUQk3iV0wE7SfdgiIpIgEjpgqw1bREQSRUIHbI10JiIiiSKhA3ZDhq2ALSIicS6hA7YGThERkUSR2AFbt3WJiEiCSOiAXV8lHlCGLSIicS6hA7YybBERSRQJHbDrB07RWOIiIhLvEjpg13c686tKXERE4lxCB2xjDF630X3YIiIS96J6vGY88Llduq1LRCSGAoEA+fn51NbWdnZR4kpycjK5ubl4vd59en/CB2yvx6UMW0QkhvLz80lPT2fQoEEYYzq7OHHBWktxcTH5+fkMHjx4n9aR0FXioAxbRCTWamtr6dmzp4J1Oxhj6NmzZ1S1EgkfsL1ul27rEhGJMQXr9ov2mCV8wE7yKMMWEZH4l/ht2G61YYuIJBq3282YMWMaXp933nlce+21nViijpfwAdunDFtEJOGkpKSwfPnymK4zGAzi8Ry4YTHhq8Sd+7A10pmISFcwaNAgbrrpJiZMmMCYMWP46quvAKiqquKiiy7iiCOO4LDDDuPVV18F4IknnmDWrFmcccYZnHTSSYTDYS677DJGjRrF6aefzqmnnsr8+fNZtGgRM2fObNjO22+/zdlnn71f9y2qSwljzF3AGUAd8A0wx1pbGpl3HfBjIARcaa19M8qy7hNl2CIiHeeWf65kVUF5TNc5sl93bjpjVIvL1NTUMH78+IbX1113Heeeey4A2dnZfPrppzzwwAPMmzePRx55hNtvv53jjz+exx57jNLSUiZNmsS0adMAWLx4MV988QVZWVnMnz+fjRs38uWXX1JYWMiIESO46KKLOP7445k7dy5FRUXk5OTw+OOPM2fOnJjud2uizf3fBq6z1gaNMb8DrgOuMcaMBM4DRgH9gIXGmGHW2lCU22s3r9tFRSC4vzcrIiIdqKUq8frM9/DDD+ell14C4K233mLBggXMmzcPcG5N27RpEwAnnngiWVlZAHzwwQfMmjULl8tFnz59OO644wCnh/cFF1zA008/zZw5c1i8eDF/+9vfOnQf9xRVwLbWvtXo5cfA9yJ/zwCetdb6gQ3GmHXAJGBxNNvbF0keF8XKsEVEOkRrmXBnSEpKApyOacGgk7BZa3nxxRcZPnz4bssuWbKEtLS0htfWNt+EOmfOHM444wySk5OZNWvWfm/vjmUb9kXAG5G/+wObG83Lj0zbizHmYmNMnjEmr6ioKIbFcaiXuIiInHzyydx7770NAfmzzz5rcrmjjz6aF198kXA4zPbt23n33Xcb5vXr149+/frxm9/8htmzZ++HUu+u1csDY8xCoE8Ts26w1r4aWeYGIAg8U/+2JpZv8rLFWvsw8DDAxIkTY947zOfRwCkiIolmzzbs6dOnc+eddza7/K9//Wuuuuoqxo4di7WWQYMG8dprr+213DnnnMOiRYsYPXo0w4YNY/LkyWRkZDTMP//88ykqKmLkyJGx3aE2aDVgW2untTTfGHMhcDpwgt1Vl5APDGi0WC5QsK+FjIbX7SKgKnERkYQSCjXdJWrjxo0Nf0+cOLEhQ05JSeGhhx7aa/nZs2fvli27XC7mzZtHt27dKC4uZtKkSbvd7/3BBx/w05/+NCb70F7R9hKfDlwDTLXWVjeatQD4uzHmbpxOZ0OBpdFsa18pwxYRkfY4/fTTKS0tpa6ujl//+tf06eNUMh9++OGkpaXxhz/8oVPKFW2L+X1AEvB2ZIzUj621l1hrVxpjngdW4VSVz+2MHuKgh3+IiEj7NG63bmzZsmX7tyB7iLaX+CEtzLsduD2a9ceCMmwREUkEGulMREQkDiR8wPa53YTCllBYQVtEROJXwgdsr8e5w0z3YouISDxL+IDtczu76FfHMxGRhHTdddfx7rvv8sorr7R4LzbAzTff3DA8abxJ/IDtcXZRGbaISGJasmQJkydP5r333uOYY47p7OJ0mMQP2JEMW7d2iYgklquvvpqxY8fyySefMGXKFB555BEuvfRSbr31Vv76179yxBFHMG7cOM455xyqq6v3ev/y5cs58sgjGTt2LDNnzmTnzp0AHHvssVxzzTVMmjSJYcOG8d///hdwBmu5+uqrOeKIIxg7dmyTA7F0pAP3Sd0x4nUrwxYR6TBvXAvbvoztOvuMgVNartoGuOuuu5g1axZPPfUUd999N8ceeywffvghAMXFxQ0jkv3qV7/i0Ucf5Yorrtjt/T/60Y+49957mTp1KjfeeCO33HILf/rTnwAIBoMsXbqU119/nVtuuYWFCxfy6KOPkpGRwSeffILf7+eoo47ipJNOYvDgwbHd/2YkfMCurxJXhi0ikng+++wzxo8fz1dffbXb+N4rVqzgV7/6FaWlpVRWVnLyySfv9r6ysjJKS0uZOnUqABdeeCGzZs1qmN/4EZ31w52+9dZbfPHFF8yfP79hHV9//bUCdqzUZ9gaPEVEpAO0IRPuCMuXL2f27Nnk5+eTnZ1NdXU11lrGjx/P4sWLmT17Nq+88grjxo3jiSeeaHb0suY094jOe++9d6/gv78kfBt2kjJsEZGEM378eJYvX86wYcNYtWoVxx9/PG+++SbLly8nJSWFiooK+vbtSyAQ4Jlnntnr/RkZGWRmZja0Tz/11FMN2XZzTj75ZB588EECgQAAa9eupaqqKvY714wuk2FrtDMRkcRSVFREZmYmLpdrryrx2267jcmTJzNw4EDGjBlDRUXFXu9/8sknueSSS6iurmbIkCE8/vjjLW7vJz/5CRs3bmTChAlYa8nJyeGVV16J+X41x+x6Imbnmzhxos3Ly4vpOpduKOH7Dy3m6R9P5uih2TFdt4hIV7R69WpGjBjR2cWIS00dO2PMMmvtxNbem/BV4l63RjoTEZH4l/ABu76XuEY6ExGReJb4AVv3YYuISAJI/ICtXuIiIpIAEj5ga6QzERFJBAkfsBsybAVsERGJYwkfsL16+IeISEJr7fGaGzduZPTo0QDk5eVx5ZVX7u8ixkTCB+wkZdgiIgmtPY/XnDhxIn/+859jXob64Us7UtcZ6Sx44AwQIyIi0bv66qt588032bBhA1OmTOGbb75h0aJFfO973+O0007joosuIjU1laOPPrrhPe+++y7z5s3jtdde4+abb2bDhg1s3bqVtWvXcvfdd/Pxxx/zxhtv0L9/f/75z3/i9XpZtmwZP//5z6msrCQ7O5snnniCvn37cuyxx/Kd73yHDz/8kDPPPJNf/OIXHbq/CR+w3S6D22WoC4U6uygiIgnnd0t/x1clX8V0nYdmHco1k65pdbmWHq85duzYhkdnXn311c2u45tvvuGdd95h1apVTJkyhRdffJHf//73zJw5k3/961+cdtppXHHFFbz66qvk5OTw3HPPccMNN/DYY48BUFpaynvvvRebHW9FVAHbGHMbMAMIA4XAbGttgTHGAPcApwLVkemfRlvYfeV1G40lLiKSgJp6vOaej8684IILeOONN5p8/ymnnILX62XMmDGEQiGmT58OwJgxY9i4cSNr1qxhxYoVnHjiiQCEQiH69u3b8P5zzz23I3dvN9Fm2HdZa38NYIy5ErgRuAQ4BRga+ZkMPBj53Sl8bpcEpc+lAAAdt0lEQVQ6nYmIdIC2ZMIdoaXHa77xxhs4eWPr6h+j6XK58Hq9De9zuVwEg0GstYwaNYrFixc3+f60tLTY7FAbRNXpzFpb3uhlGlCfxs4A/mYdHwM9jDF991rBfuLzuNTpTEQkgbT0eM2+ffuSkZHBBx98ANDk4zXbavjw4RQVFTUE7EAgwMqVK2OyD+0VdS9xY8ztxpjNwPk4GTZAf2Bzo8XyI9Oaev/Fxpg8Y0xeUVFRtMVpkjJsEZHE09LjNR9//HHmzp3LlClTSElJ2edt+Hw+5s+fzzXXXMO4ceMYP348H330USyK326tPl7TGLMQ6NPErBusta82Wu46INlae5Mx5l/Ab621H0TmLQL+z1q7rKVtdcTjNQGm3vUO4wf04J7zDov5ukVEuho9XnPfRfN4zVbbsK2109pYjr8D/wJuwsmoBzSalwsUtHE9MacMW0RE4l1UVeLGmKGNXp4J1PftXwD8yDiOBMqstVuj2VY0vG6XxhIXEZG4Fm0v8TuNMcNxbuv6FqeHOMDrOLd0rcO5rWtOlNuJis/j0vOwRUQkrkUVsK215zQz3QJzo1l3LPmUYYuISJxL+LHEIXJblzJsERGJY10iYGukMxERiXddImArwxYRSVytPV6zXuPHbLbXE088QUFBp93sBHSRgK1e4iIiias9j9fcV/sSsGP9yM0uEbDVS1xEJPFcffXVjB07lk8++YQpU6bwyCOPcOmll3Lrrbeybt06pk2bxrhx45gwYQLffPPNbu994oknuPzyyxten3766bz77ruEQiFmz57N6NGjGTNmDH/84x+ZP38+eXl5nH/++YwfP56amhqWLVvG1KlTOfzwwzn55JPZutW5c/nYY4/l+uuvZ+rUqdxzzz0x3d+Ef7wmqJe4iEhH2XbHHfhXx/bxmkkjDqXP9de3ulxLj9ecPHky1157LTNnzqS2tpZwOExhYWGr61y+fDlbtmxhxYoVgPP4zB49enDfffcxb948Jk6cSCAQ6JRHbnaNgK2Hf4iIJKSmHq9ZUVHBli1bmDlzJgDJycltXt+QIUNYv349V1xxBaeddhonnXTSXst01iM3u0TA9rpdBFQlLiISc23JhDtCS4/XbEt26/F4CId3xYXa2loAMjMz+fzzz3nzzTe5//77ef755xsy53qd9cjNLtOGrQxbRCRxtPR4zYyMDHJzc3nllVcA8Pv9VFdX7/b+QYMGsXz5csLhMJs3b2bp0qUA7Nixg3A4zDnnnMNtt93Gp59+CkB6ejoVFRVA5z1ys+tk2CGLtbbNDzUXEZEDW0uP13zqqaf42c9+xo033ojX6+WFF17A5dqVox511FEMHjyYMWPGMHr0aCZMmADAli1bmDNnTkP2/dvf/haA2bNnc8kll5CSksLixYuZP38+V155JWVlZQSDQa666ipGjRrVofvb6uM196eOerzm/e+s464317DmN9NJ8rhjvn4Rka5Ej9fcd9E8XrNLVIl73U5WrdHOREQkXnWRgO3spkY7ExGReNUlArbP4+ym7sUWEYmNA6k5NV5Ee8y6RMBWhi0iEjvJyckUFxcraLeDtZbi4uJ23RO+py7RSzwpkmHr1i4Rkejl5uaSn59PUVFRZxclriQnJ5Obm7vP7+8SAVsZtohI7Hi9XgYPHtzZxehyukSVuM+tNmwREYlvXSJgez3KsEVEJL51iYBdn2GrDVtEROJV1wjYHmfgFGXYIiISr7pGwHY7w5FqpDMREYlXMQnYxphfGmOsMSY78toYY/5sjFlnjPnCGDMhFtvZV15l2CIiEueiDtjGmAHAicCmRpNPAYZGfi4GHox2O9FQL3EREYl3sciw/wj8H9C4vnkG8Dfr+BjoYYzpG4Nt7RPdhy0iIvEuqoBtjDkT2GKt/XyPWf2BzY1e50emNbWOi40xecaYvI4aNUcjnYmISLxrdaQzY8xCoE8Ts24ArgdOauptTUxrsseXtfZh4GFwnofdWnn2hVdV4iIiEudaDdjW2mlNTTfGjAEGA58bYwBygU+NMZNwMuoBjRbPBQqiLu0+8mngFBERiXP7XCVurf3SWtvLWjvIWjsIJ0hPsNZuAxYAP4r0Fj8SKLPWbo1NkdtPGbaIiMS7jnr4x+vAqcA6oBqY00HbaROvW7d1iYhIfItZwI5k2fV/W2BurNYdLWMMPreLOg2cIiIicapLjHQGTju2MmwREYlXXSZge91GbdgiIhK3ukzAVoYtIiLxrMsEbK/bpQxbRETiVpcJ2D6PC78CtoiIxKmuE7DdLgKqEhcRkTjVdQK2x6WxxEVEJG51mYCtNmwREYlnXSZg+9zqJS4iIvGrywRsr0cjnYmISPzqMgFbGbaIiMSzrhOwPRrpTERE4lfXCdjKsEVEJI51mYCtXuIiIhLPukzA1ljiIiISz7pMwPa6NXCKiIjEry4TsJOUYYuISBzrMgFbbdgiIhLPukzA9nlchC0EFbRFRCQOdZmA7XU7uxrQaGciIhKHukzA9nmcXVU7toiIxKOuE7DdBkA9xUVEJC5FFbCNMTcbY7YYY5ZHfk5tNO86Y8w6Y8waY8zJ0Rc1Og0ZtgK2iIjEIU8M1vFHa+28xhOMMSOB84BRQD9goTFmmLU2FIPt7ZOGNmxViYuISBzqqCrxGcCz1lq/tXYDsA6Y1EHbahNl2CIiEs9iEbAvN8Z8YYx5zBiTGZnWH9jcaJn8yLS9GGMuNsbkGWPyioqKYlCcptVn2Op0JiIi8ajVgG2MWWiMWdHEzwzgQeBgYDywFfhD/duaWFWT91NZax+21k601k7MycnZx91onTJsERGJZ622YVtrp7VlRcaYvwKvRV7mAwMazc4FCtpduhjyqQ1bRETiWLS9xPs2ejkTWBH5ewFwnjEmyRgzGBgKLI1mW9FShi0iIvEs2l7ivzfGjMep7t4I/AzAWrvSGPM8sAoIAnM7s4c4NB7pTAFbRETiT1QB21p7QQvzbgduj2b9seRTpzMREYljXWaks/Rk59qkuKquk0siIiLSfl0mYOdmppCV5uOzTaWdXRQREZF26zIB2xjDhIMy+fTbnZ1dFBERkXbrMgEb4PCBmazfUUWJqsVFRCTOdLmADbBMWbaIiMSZLhWwx+Zm4HUbBWwREYk7XSpgJ3vdjOqXoXZsERGJO10qYINTLf55fqnuxxYRkbjSJQO2PxhmZUFZZxdFRESkzbpkwAZ1PBMRkfjS5QJ27+7J5Gam8OkmBWwREYkfXS5gg5NlL/t2J9Y2+YhuERGRA06XDdjby/1sKa3p7KKIiIi0SZcM2BMOUju2iIjEly4ZsA/tk06qz62ALSIicaNLBmyP28VhB/VQwBYRkbjRJQM2wOEHZbJ6azlV/mBnF0VERKRVXTZgTxiYSdjC8s16PraIiBz4umzAPuygTHxuFy9/tqWziyIiItKqLhuwM1K8zD5qEC9+ms+qgvLOLo6IiEiLumzABph77CFkpHi54/XVGkRFREQOaF06YGekerny+KF8sG4H760t6uziiIiINCvqgG2MucIYs8YYs9IY8/tG068zxqyLzDs52u10lB8eOZCBPVO54/XVBEN65KaIiByYogrYxpjjgBnAWGvtKGBeZPpI4DxgFDAdeMAY446yrB3C53Fx7fRDWbu9kvnL8ju7OCIiIk2KNsO+FLjTWusHsNYWRqbPAJ611vqttRuAdcCkKLfVYaaP7sPhAzP5w9trdV+2iIgckKIN2MOAY4wxS4wx7xljjohM7w9sbrRcfmTaXowxFxtj8owxeUVFndOObIzhhtNGUFTh57JnPmVzSXWnlENERKQ5rQZsY8xCY8yKJn5mAB4gEzgSuBp43hhjANPEqprshm2tfdhaO9FaOzEnJyeKXYnOhIMyufmMkXyysYRpd7/HH99eS20g1GnlERERaczT2gLW2mnNzTPGXAq8ZJ17opYaY8JANk5GPaDRorlAQZRl7XCzjxrMyaP7cMfrX3HPoq958dN8fnfOWI46JLuziyYiIl1ctFXirwDHAxhjhgE+YAewADjPGJNkjBkMDAWWRrmt/aJvRgr3/uAw/vHTI0nyuLjwsaUs+PyAv9YQEZEEF23AfgwYYoxZATwLXGgdK4HngVXAv4G51tq4ql+ecnBPXp57FBMOyuR/nv2Mpz7+trOLJCIiXZg5kEb4mjhxos3Ly+vsYuymNhBi7jOfsuirQn550jDmHncITjO9iIhI9Iwxy6y1E1tbrkuPdNYWyV43f7ngcM4a3495b63lF89/zqZi9SIXEZH9q9VOZwJet4u7vz+efj1SeOS/G3j18wLOHNePS489mGG909lR6WdlQTkrC8qoDYQ5dngO43N74HIpExcRkdhQlXg7bS+v5a/vr+eZJZuoCYTI7pbEjkp/w3yXgbCF7G4+jj+0FxMHZlFdF6SsJkhZTYCaQIjMVC9ZaT56dvORlZZERoqX7ske53eKF69bFR8iIl1FW6vEFbD30c6qOp5cvJFNJdWM7Nudkf26M6pvBgDvri1k4epC3l1TSEXtrpHT0nxukr1uymoCBMNNH3e3y3DyqN789JghHHZQ5v7YFRER6UQK2AeAQChMQWkN6clOBu2JZM7WWsprghRX+dlZXUdZTYDySAb+bXE1LyzbTEVtkEmDsvjJMYM5YlAWPVK96uwmIpKAFLDjWKU/yPOfbObRDzawpbQGAJ/bRU56Er27J5Gbmcqg7DQGZ6cysGcaPVK8eFwu3G6Dx2XISvOpWl1EJE4oYCeAYCjM+18XsWFHNYUVtRSV+9lWXsumkmoKSmtopladzFQvMw/L5dwjBjC8T/r+LbSIiLRLWwO2eokfwDxuF8cf2rvJef5giM0lNWzcUUWlP0gwbAmFwwRClsXfFPPUxxt57MMNjB/Qg1PH9OGQXt0Ykt2N3MyUhqr5WKsLhqmuC1LpD1JcWUdRhZ+iSj87KvzkZqVw/PDeZKR6O2TbIiKJThl2giqu9PPyZ1t47pPNfF1Y2TDd6zYMyEqlX0YKfTOS6dsjhX4ZyfTrkRL5SSbJ4+bb4iq+2lbBV9sqWF9USZU/SE0gRG0gTG0gRF0wjD8YJhAKUxcKU+0PURcKt1gmj8tw5JCenDSqNxMHZtG7exJZaT61zcdI/WdZx1MkvqhKXBqUVNWxvqiS9UVVrN9RxaaSKgpKa9laVkNhhZ89TwGPyzT0YncZGJCVSvdkLyleN8k+N8keF0leN163Icnjwut2keJz083nIS3JQ1qSm6y0JHLSk+iV7gTl1VvLeWvVdt5cuY31RVUN2/K6Db3Skxk/oAdXnjBUVfhtVBsI8ULeZj7bVMrWslq2lTv/T4MhNzOFAVmpDMhMoX9mCr27J9OnezJ9MpyfJI+7s4svIo0oYEubBEJhtpfXsrWsloLSGgpKaymrCTAkJ40RfboztHc3kr2x/YJfV1jJ19sr2F5ey/YKP9vKalm4ajuVdUHOGt+fq6YNZWDPtJhusyWBUJiiCj/ltQFq6kLU1IWorgvh87gaLjoyU33tGggnHLYEwxafJ7bNDzV1IZ5Z8i1/eW89Oyr9DbUjvTOS6ds9GYDNO6vZVFJDfkk1Ff7gbu/3uV1MGpzF1GE5HDs8h0N6dWtzRh4OWyr8Qay1GAwY54IuzedJ2EGC/MEQ764pYmtpDSP7ZTCqX3fSknZvSbTWErbOLZki+0IBW+JKaXUdf3lvPU98tIFgyDLzsP6cPSGXyYOz2hQMgqEwXxdW8vnmUj7PL2P11nKq/EHqQmH8Aafq3udxagLSfB5SvG7KawMUVfgprqprdf31NQH9eiTTP9J8kJbkobS6jpKqADur65xb9KoDDbfqWSA3M4XhvdMZ2judYb27MbRXOgfndCPF1/aLoGAozBdbynh/bRFPf/wtOyrrOOqQnlxx/FCOHNKz2fdZaymvDVJY7mTg28pqWbOtgvfWFjU0k/TunsShfboztFc3hvbuRv8eqWyPdGzcXFJN/s4aiqv8lEb2q6mOjh6XabiwyUlPpnuKh2Svm5TIT49UL73qs/zuyfTtkdzmuxhCYcu3xVVU14WoDYSoCYSo8gcpqvCzvdzP9nLnAnPCwExOHNmbg3O6tfm4trTNJRuKefWzAl5fsXW3sRSMgSHZaQzISmVndYAdFX6Kq/z4g2GyuyU11GT075HChIGZHDk4i16RC6l9URsIUVJVR0lVXeQukX1fl+w7fzDE2m2VrCgoY822Cm46Y2RMm54UsCUuFZbXct8763ghL5+aQIje3ZM4Y2w/Jg/pyfbyWjaXVLN5ZzVby2qp9oeoDgSpqQtRXhukLui0oXdP9jCqXwY9Ur0keVz4ItX2gVCYqroQ1f4g1XUh0pO99OruBJpe6clkpHhJ9blJ8blJ9bmpDTiZd2FFLYUVfraX1bKltIaCshq2ltYSDFtSvG4yU71kpvnITPXRI9Xr/KT4cLsM3xRV8vX2StbvqCQQqm9jdgL5wTndyO6W1PD+Hik+LJa6YJi6YJiaQIgv8stYuqGEykimfMzQbP7nhKFMHJQV1XHO31nN+2t3sGRDMV9vr+Sbokr8wV19EIyBvt2Tyc1MJTvd2bf6/XO7DNZC2FrC1rKzOkBhuXOciir8VPqD1Eb6O9QEQoT2iPLpSR6mjezNKaP78N1hOXvV4NTUhfjv10W8vWo7//mqsNkLKrfLkNMtiRSfmw07nGaWIdlpHH9oLzLTfARCYUJhSyBksdZicS5irIUUn7thZMHuyV6Kq/ys3lrOqoJyvtpWQXVdiDSfm5NH92HG+P4M692NVQXlfLmljBVbythWXktWWhLZkRELU7xuiir9TvNEmXPBU13nPKBwcHYahw/MpHf3JLone8lI8ZKe7MXjNriNcTJzAwWlNWyINFtt2FFFYXktVXW7P+TwkF7dOPqQbI46JJuDc9IIW0soDMGw87/zuZ1z3etxUVkbZNXWMlZuKWdlQTlby2rISvPRKz2ZnPQk0pM9bCuvJX9nDVt21lBYUUtWmo++GSn07+H0cWloXslKJTczBYOhqi7Y8NkLhS0uYzAGXMa5rTTZ6yYp0mxWXOlnxRZn2OaVBeWU1gToHbnw6JORTM80X8PyyV43XreLYNj5vwXDlkCkf0ylP0iVP0ggFGZkvwwmDc4iK83X6nleUxdqaPorrPBTWF5LdV2IYNjuVgvWI8Xb8NkNh6GgrMb5rJfWsmFHJWu2VTR8ftOTPSz6xVR6pcfu4kkBW+JadV2QhasLWbC8gPfWFjZ8WHxuF7mZKfTtkUy3JCdTTvF56J7sYUTf7ozNzWBQz7QOr6INRb5M2tpcEAiF2bijymkOiPysL6qkuNLJzBsHy8aG5KTxnYN78p2DszlySM82fUnti1DYkr+zmi2lNfTN2NX5MFr1Wf72SIa/rbyWTzaU8Naq7ZTVBEjzuRnVLwN/KIw/kkFvK6vFHwyTnuzhuOG9OHpoNj1SvKT4Ilm7z01OehI905IaqqG3lNawaPV23l61nY/XFzecLy4DHpcLY5yLEIMTXGoCob36bqRHzqGRfbszcVAmJxzau101IY0FQ2FWbS1nyfoSlmwoYfnmUnZW1+118bKnZK+LQT3TODinG30ykslK85EVuRjcVFLFB+uKWbqhmNpAyx08G/N5XIzok05uViql1XUUljt3b5TXBOgTuSjLzUwhJz2Jkqq6hgvSgrKadm2nJR6XYWjvdLK7+ZymsHI/ZTWBdq/HGBr+b8N6d2PioCz6dE8mPdnjfB/43GwqqWZVQTmrtpazYUfVXv/nem6Xc8EUCIebXMbndtG3RzIHZaUyql8Go/t3Z0z/DA7KSo15x04FbEkYpdV1fFNU6bTVpicnZHtpTV2IspoALkNDjUD970QUCIX5eH0xr3+5jW+KKiNV6C5SvG6yuyVx3KG9mDQ4a5/2PxByvoA9LtPsuVLfHl9eE6CsJkBGitfJIDuwh721luq6EOW1zsiG9eUMRWoqend3+iG0dn77gyE+/baUbeU1uF0uZz+NAZzahEAoTDDkZI4j+nbn4Jy0Jm/ltNa2uL/WWnZU1rF5p9M8sqXU6dSYluQm1ech1edu2G44UuMSDFn8wRD+oHM3SfdkL6P6ZTCsT7e9LgBr6kKUVNfhD+xaPhCyuF1Opu5xG7xuF6k+N92SnA6tYWv5Mr+MJRtKWLqhhE837dytyaJebmYKI/t2Z0Tf7gzKTqVXenJDTVq3ZA8us+tuilDYUlEboLQ6QGlNAGst/TNTyE5L2m/fNQrYIiKS8AKhMFX+IBW1zhgQ/TJS4m68Bw2cIiIiCc/rdtEj1UeP1I5pLjqQJGZ9m4iISIJRwBYREYkDCtgiIiJxQAFbREQkDihgi4iIxIGoArYx5jljzPLIz0ZjzPJG864zxqwzxqwxxpwcfVFFRES6rqhu67LWnlv/tzHmD0BZ5O+RwHnAKKAfsNAYM8xaG2pyRSIiItKimFSJG2fImO8D/4hMmgE8a631W2s3AOuASbHYloiISFcUqzbsY4Dt1tqvI6/7A5sbzc+PTNuLMeZiY0yeMSavqKgoRsURERFJLK1WiRtjFgJ9mph1g7X21cjfP2BXdg3Q1ACsTY6Baq19GHgYnKFJWyuPiIhIV9RqwLbWTmtpvjHGA5wNHN5ocj4woNHrXKCgtW0tW7ZshzHm29aWa4dsYEcM19dV6TjGho5jbOg4xoaOY2zE4jgObMtCsRhLfBrwlbU2v9G0BcDfjTF343Q6GwosbW1F1tqcGJSngTEmry0DqkvLdBxjQ8cxNnQcY0PHMTb253GMRcA+j92rw7HWrjTGPA+sAoLAXPUQFxER2XdRB2xr7exmpt8O3B7t+kVERCTxRzp7uLMLkCB0HGNDxzE2dBxjQ8cxNvbbcTTWqmO2iIjIgS7RM2wREZGEoIAtIiISBxI2YBtjpkcePLLOGHNtZ5cnXhhjBhhj3jHGrDbGrDTG/E9kepYx5m1jzNeR35mdXdZ4YIxxG2M+M8a8Fnk92BizJHIcnzPG+Dq7jAc6Y0wPY8x8Y8xXkfNyis7H9jPG/G/kM73CGPMPY0yyzsfWGWMeM8YUGmNWNJrW5PlnHH+OxJ0vjDETYlmWhAzYxhg3cD9wCjAS+EHkgSTSuiDwC2vtCOBIYG7k2F0LLLLWDgUWRV5L6/4HWN3o9e+AP0aO407gx51SqvhyD/Bva+2hwDic46nzsR2MMf2BK4GJ1trRgBvnllydj617Api+x7Tmzr9TcMYdGQpcDDwYy4IkZMDGedDIOmvtemttHfAszgNJpBXW2q3W2k8jf1fgfDn2xzl+T0YWexI4q3NKGD+MMbnAacAjkdcGOB6YH1lEx7EVxpjuwHeBRwGstXXW2lJ0Pu4LD5ASGZ0yFdiKzsdWWWvfB0r2mNzc+TcD+Jt1fAz0MMb0jVVZEjVgt/nhI9I8Y8wg4DBgCdDbWrsVnKAO9Oq8ksWNPwH/B4Qjr3sCpdbaYOS1zsvWDQGKgMcjTQuPGGPS0PnYLtbaLcA8YBNOoC4DlqHzcV81d/51aOxJ1IDd5oePSNOMMd2AF4GrrLXlnV2eeGOMOR0otNYuazy5iUV1XrbMA0wAHrTWHgZUoervdou0sc4ABuMMF52GU327J52P0enQz3iiBux9eviIOIwxXpxg/Yy19qXI5O31VTuR34WdVb44cRRwpjFmI06TzPE4GXePSJUk6Lxsi3wg31q7JPJ6Pk4A1/nYPtOADdbaImttAHgJ+A46H/dVc+dfh8aeRA3YnwBDIz0gfTidKxZ0cpniQqSd9VFgtbX27kazFgAXRv6+EHh1z/fKLtba66y1udbaQTjn33+stecD7wDfiyym49gKa+02YLMxZnhk0gk4zyjQ+dg+m4AjjTGpkc94/XHU+bhvmjv/FgA/ivQWPxIoq686j4WEHenMGHMqTkbjBh6LjG0urTDGHA38F/iSXW2v1+O0Yz8PHITz4Z9lrd2zI4Y0wRhzLPBLa+3pxpghOBl3FvAZ8ENrrb8zy3egM8aMx+m45wPWA3Nwkg2dj+1gjLkFOBfnTpDPgJ/gtK/qfGyBMeYfwLE4j9HcDtwEvEIT51/kYug+nF7l1cAca21ezMqSqAFbREQkkSRqlbiIiEhCUcAWERGJAwrYIiIicUABW0REJA4oYIuIiMQBBWwREZE4oIAtIiISB/4/HMCXQr8jKMkAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(8,5))\n", + "plt.plot(avgene,label='Energy')\n", + "plt.plot(nalone,label='#alone')\n", + "plt.plot(ndimer,label='#dimer')\n", + "plt.plot(ncluster,label='#cluster')\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create animation\n", + "\n", + "Visualize the trajectory." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "
\n", + " \n", + "
\n", + " \n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Once \n", + " Loop \n", + " Reflect \n", + "
\n", + "
\n", + "\n", + "\n", + "