{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ase.io import read\n",
    "from ase.visualize import view\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma   = 3.405\n",
    "epsilon = 119.8*8.616733e-5 # eV\n",
    "unit_cell_side = 5.269\n",
    "cell_side = 6*unit_cell_side\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def g_step(box, n_bin, d):\n",
    "    box    += 1\n",
    "    box    *= np.sqrt(2)/sigma\n",
    "    del_bin = box/n_bin  #width of a bin\n",
    "    \n",
    "    g    = np.zeros([n_bin,n_bin])\n",
    "    g[0] = np.linspace(0,box,n_bin)\n",
    "    g[0]+= del_bin*0.5\n",
    "    \n",
    "    for i in range(len(d)):\n",
    "        for j in range(i+1,len(d)):\n",
    "            g_index = d[i][j]/(del_bin)\n",
    "            #g_index = (d[i][j]-box*round(d[i][j]/box))/(del_bin)\n",
    "            g[1,g_index.astype(int)] += 2\n",
    "            # we count both for i and j.\n",
    "            \n",
    "    \n",
    "    #nomalized by number of particles\n",
    "    g[1] /= len(d[0])\n",
    "    #and by bin volume\n",
    "    \n",
    "    for k in range(len(g[1])):\n",
    "        g[1][k] /= ((k+1)**3 -k**3)*del_bin**3\n",
    "        \n",
    "    #side of the optimized cell: 5,269\n",
    "    vol = (unit_cell_side/sigma)**3\n",
    "    rho = 4./vol\n",
    "    g[1] /= np.pi*(4/3)*rho\n",
    "    return g;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20 K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory = read('production_bulk_20-pos-1.xyz', index='::2')\n",
    "N = len(trajectory[0])    #number of atoms\n",
    "n_step = len(trajectory)  #number of step\n",
    "\n",
    "#create a distance array\n",
    "dist = np.empty([0,N,N])\n",
    "for frame in trajectory:\n",
    "    frame.set_cell([cell_side,cell_side,cell_side])\n",
    "    frame.set_pbc((True,True,True))\n",
    "    dist = np.append(dist, [frame.get_all_distances(mic=True)],axis=0)\n",
    "    \n",
    "\n",
    "#Lennard Jones units\n",
    "\n",
    "dist  /= sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'radial distribution function')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEOCAYAAACKDawAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztvXeYJHd1r/+e7umenGc2ze5s0moVVmhXWoQSAiRsJAECY2zgIrjG+OpHMCAMxhcnwo978TU2NsY2togmmGsMEpggkBDKQoLdlbS70uaond3Zybmn4/f+UVU9PTNd1dU9UxPP+zz9TMeq0z3dp059vieIMQZFURRl6ROabwMURVGUuUEdvqIoyjJBHb6iKMoyQR2+oijKMkEdvqIoyjJBHb6iKMoyQR2+oijKMkEdvqIoyjJBHb6iKMoyQR2+oijKMqFsvg3IpaWlxWzYsGG+zVAURVk07N69u8cY0+rnuQvK4W/YsIFdu3bNtxmKoiiLBhE55fe5gUo6IvJBEXlORPaLyLdFpCLI/SmKoijuBObwRaQNeD+w0xizDQgDbw5qf4qiKIo3QS/algGVIlIGVAFnA96foiiK4kJgDt8Y0wH8DXAaOAcMGmPuC2p/iqIoijdBSjqNwOuAjcAaoFpEbs/zvDtEZJeI7Oru7g7KHEVRlGVPkJLOK4ETxphuY0wSuBu4duqTjDF3GWN2GmN2trb6yixSFEVRSiBIh38auFpEqkREgJuAAwHub1nz/NkhUunMfJuhKMoCJkgN/yngu8AeYJ+9r7uC2t9y5txgjFd//lF+fuD8fJuiKMoCJtDCK2PMx4CPBbkPBbqG4hgDA2PJ+TZFUZQFjPbSWQIMxCxHn8yYebZEUZSFjDr8JcDAWAKAZEo1fEVR3FGHvwQYciJ8XbRVFMUDdfhLAEe7T6mkoyiKB+rwlwCOhp9QSUdRFA/U4S8BBlXSURTFB+rwlwAq6SiK4gd1+EuAwZiVpaOSjqIoXqjDXwKopKMoih/U4S8BspJOWiUdRVHcUYe/BNAIX1EUP6jDX+SMJ9PEbe1eWysoiuKFOvxFTm7DNG2toCiKF+rwFzmOnAMq6SiK4o06/EWO0zgNVNJRFMUbdfiLHKetQlU0rJKOoiieBDnEfKuIPJNzGRKRO4Pa33LFkXRaa8tV0lEUxZPAJl4ZYw4B2wFEJAx0APcEtb/lyqC9aNtSU66SjqIonsyVpHMTcMwYc2qO9rdsGIglCIeEhsqISjqKongyVw7/zcC352hfy4rBWJKGygjRspBKOoqieBK4wxeRKHAb8J8uj98hIrtEZFd3d3fQ5iw5BsaS1FdGiIRD2i1TURRP5iLCvwXYY4w5n+9BY8xdxpidxpidra2tc2DO0mIwlqS+KkJZWLRbpqIonsyFw38LKucExmDMivCj4RCpjDp8RVHcCdThi0gV8BvA3UHuZzkzMGZp+JFwiKR2y1QUxYPA0jIBjDFjQHOQ+1juDMaSNFRFEdFeOoqieKOVtouYdMYwNJ6kzpZ0kirpKIrigTr8RczweBJjoKHSWrRVSUdRFC/U4S9inLYKDVWWhp/OGDKamqkoigvq8BcxTi98Jw8fUFlHURRX1OEvYgYmRfgCoLKOoiiuqMNfxDiSTn1ldCLC10wdRVFcUIe/iBm0h5+opKMoih/U4S9iJiJ8lXQURSlMQYcvIm8QkSMiMmgPMRkWkaG5ME7xZmAsSXU0TLQslI3wU9oxU1EUF/xU2v418FpjzIGgjVGKY8DuowNMSDrq8BVFccGPpHNenf3CxOqUGQXISjqJlEo6iqLkx0+Ev0tE/gP4PhB37jTGaEO0eWZwLEl9pfUvzEo6umirKIoLfhx+HTAG/GbOfQbtgDnvDMQSbGqpAVTSURSlMAUdvjHmHXNhiFI8VqdMS8MvU0lHUZQC+MnSWSsi94hIl4icF5HvicjauTBO8cYZbwgQVUlHUZQC+Fm0/SrwX8AaoA34oX2fMo+MJ9PEUxnqsxG+SjqKonjjx+G3GmO+aoxJ2ZevAb6Gz4pIg4h8V0QOisgBEblmRtYqWbKdMis1S0dRFH/4cfg9InK7iITty+1Ar8/tfw74qTHmIuByQNM7Z4ncTpmgko6iKIXx4/B/H/hdoBM4B7zRvs8TEakDbgC+DGCMSRhjBko3Vckltxc+qKSjKEph/GTpnAZuK2Hbm4Bu4KsicjmwG/iAMWa0hG0pUxjIaZwGaC8dRVEK4urwReQjxpi/FpHPY+XdT8IY834f274CeJ8x5ikR+RzwP4G/mLKfO4A7ANrb24s0f/kyEMsv6WiEryiKG14RvqO37ypx22eAM8aYp+zb38Vy+JMwxtwF3AWwc+dODU99Mp5MA1AVDQM5ko72w1cUxQVXh2+M+aF9dcwY85+5j4nI7xTasDGmU0ReEJGtxphDwE3A8zOyVsmSsB17pMxy9I6kk9KZtoqiuOBn0fajPu/Lx/uAb4nIXmA78L/9GqZ442j1jpTjtFZIqKSjKIoLXhr+LcCtQJuI/EPOQ3VAys/GjTHPADtnZKGSF0erj0xx+EnNw1cUxQUvDf8sln5/G1aGjcMw8MEgjVIKk0xnCAmEQ5aUEw4JIdE8fEVR3PHS8J8FnhWRe4BRY0waQETCQPkc2ae4kEhnslG9QyQcUklHURRX/Gj49wGVObcrgZ8HY47il2TKZPV7h0g4pJKOoiiu+HH4FcaYEeeGfb0qOJMUPyTTmWyGjkMkLCrpKIriih+HPyoiVzg3RORKIBacSYofkulMNhXToSwc0sIrRVFc8TPx6k7gP0XkrH17NfCm4ExS/JBPw4+GQ9paQVEUV/z00vm1iFwEbAUEOGiMSQZumeJJMp1PwxeN8BVFccVPhA/wYmCD/fwdIoIx5uuBWaUUJJmaHuGrpKMoihcFHb6IfAPYDDwDpO27DaAOfx6xFm0na/gRlXQURfHAT4S/E7jEGKOeZAGRSGemSTpRlXQURfHAT5bOfmBV0IYoxZHMs2irko6iKF74ifBbgOdF5FdA3LnTGFPKUBRllkimTbY1soO1aKsnYoqi5MePw/940EYoxWNF+JFJ90XCIUbivvraKYqyDPGTlvnwXBiiFEciNb3wKqKSjqIoHvjJ0hlmYsRhFIhgNVOrC9IwxZt8Gn4kLKRU0lEUxQU/EX5t7m0ReT1wVWAWKb7IX3il3TIVRXHHb+FVFmPM90Vk2mzafIjISaz++WkgZYzRYSizRP4IP6QRvqIorviRdN6QczOElZdfjFd5hTGmp1jDFG/yF15pHr6iKO74ifBfm3M9BZwEXheINYpvEtpaQVGUIvGaaft/jDF/AtxrjPlOids3wH0iYoB/NcbcVeJ2lCnk0/C1W6aiKF54VdreKiIRwJde78J1xpgrgFuA94rIDVOfICJ3iMguEdnV3d09g10tL9yydDTCVxTFDS+H/1OgB3iRiAzlXIZFZMjPxo0xZ+2/XcA95MnuMcbcZYzZaYzZ2draWsJbWH5kMoZUxqikoyhKUbg6fGPMHxtj6oEfG2Pqci61fnLwRaRaRGqd68BvYvXlUWZI0h5j6NYtU/vcKYqSDz95+KUu0K4E7hERZz//boz5aYnbUnJwdPp83TIBO/qXaa9TFGV5U3Qevl+MMceBy4Pa/nImkbIj/DySDuTX9xVFUdQrLEIcnT5f4ZX1uEo6iqJMRx3+ImQiwp8s2ziSji7cKoqSDz+VttdhtUhebz9fAGOM2RSsaYobjkOPluWXdLS9gqIo+fCj4X8Z+CCwm4mZtso84kg27pKORviKokzHj8MfNMbcG7glim/cNXxL0tGOmYqi5MOPw39QRD4D3M3kEYd7ArNK8SSRzq/hR1TSURTFAz8O/yX239zWxga4cfbNUfyQtBdt8/XDB5V0FEXJj5/Cq1fMhSGKf7Ia/rRF26Ur6RhjsIv4FEUpkYJpmSJSLyKfdRqcicjfikj9XBin5MdNw48uUUnnCw8d45bPPUoms7Tel6LMNX7y8L+CNbXqd+3LEPDVII1SvCmk4S81Sedo1wgHO4d58njvfJuiKIsaPw5/szHmY8aY4/blE4Dm4M8j2Tz8aa0Vlmbh1XjSyga+++mOebZEURY3fhx+TESud27YhVix4ExSClFI0llqrRVitsO/d985YgktBVGUUvGTpfNu4N9s3V6APuD3gjRK8SaZyr9ou1QlnbFEispImNFEmvue7+R129vm2yRFWZQUjPCNMc8YYy4HXgRcZozZYYx5NnjTFDfcNPylKunEkhl2bmhkTX0F96isoygl4zXT9nZjzDdF5I+m3A+AMeazAdumuOCm4S9VSWc8kWZVXTmv29HGXY8cp3s4Tmtt+XybpSiLDq8Iv9r+W5vnUhOwXYoHhdsjL7UIP01VtIw37GgjnTH817Nn59skRVmUuEb4xph/ta/+3BjzeO5j9sKtMk9kJ165FF6llqDDr4iE2bKylm1tddzz9Bneef3G+TZLURYdfrJ0Pu/zvryISFhEnhaRH/k3S/HC6YdfFsqfh59YgpJOZSQMwOu3t7G/Y4gX+sbm2SpFWXx4afjXANcCrVN0/DogXMQ+PgAcsF+nzALJdIZoODSt1UBkyS7apqmMWgezTa2W0tg7mmBdU9V8mqUoiw6vCD+KpdWXMVm/HwLe6GfjIrIWeDXwpZmZqeRizayd3ldmolvm0nH4iVSGVMZkI/zKiBWjjCVS82mWoixKvDT8h4GHReRrxphTJW7/74GPYB0o8iIidwB3ALS3t5e4m+VFMm2m5eDDhMSzlCQdp+iqwnb4VVHrr1N9qyiKf/wUXn1NRKZ5EGOMZ3tkEXkN0GWM2S0iL3d7njHmLuAugJ07dy4dTxUgiXRmWoYOWCmzkbAsqQjfcexV0TL7r+Xwx7TiVlGKxo/D/3DO9QrgtwE/59PXAbeJyK326+pE5JvGmNuLN1PJJZnKTMvBd4iEQ0tKw3daKTgavhPpq8NXlOLx0w9/95S7HheRh3287qPARwHsCP/D6uxnBzcNHxyHv3ROlBxJp1IlHUWZMQUdvog05dwMAVcCqwKzSClIMm3ySjpgZeosqQh/mobvLNqqw1eUYvEj6ezGGmkoWFLOCeCdxezEGPMQ8FCRtikuuGn4sIQlHdvhl9uL1erw8xNPpTlyfoRtbTqjSJmOH0lHSxoXGMl0Jm+WDixBSSer4VsOPxQSKiNhYpqWmZcvPHSMf/zFUZ752G9SU+4nnlOWE34knQrgPcD1WJH+Y8AXjDHjAdumuGAVXuXX8MuWqKTjRPhg6fgx1fDz8pN950hlDIOxpDp8ZRp+vhFfxxpx6LRTeAvwDeB3gjJK8SaRcpd0oktN0klOjvDB0vNV0pnOse4RDp8fAWA0rmdAynT8OPytdj98hwdFRPvhzyOJtKEqujwknXG3CF8d/jR+ur8ze31EHb6SBz/N054WkaudGyLyEuBxj+crAZP0iPCXnKSTmB7hq6STn5/u78wuamuEr+TD1eGLyD4R2Qu8BHhCRE6KyAngl8ANc2WgMp1kOkO0zCsPfwk5fCcts0wlHS/O9I+xr2OQm7dZGdPq8JV8eEk6r5kzK5SiSHqmZQrx5BJy+Ik05WUhQjmtoKuiYXpGEvNo1cLDkXPeeOVafvDMWUbiekBUpuPl8PuNMUNTCq+UBYB34VWIkfGlE91ZrZEnd+OuipYRS8bmyaKFyc+e6+SiVbVcstrqQq4RvpIPL4f/71hRfm7hlYMBNgVol+JBocKrJdUtM2f4iUNFRBdtc+kaHmfXqX7uvOlCqu1UTF20VfLh1R75NWJN2HiZMeb0HNqkFMArD3+pdcvMH+GHtR9+Dvc9dx5j4JbLVlFeFqIsJBrhK3nxzNIxxhjgnjmyRfGJV5bOUlu0HU9Oj/A1S2cyjxzupr2pii0rahARqsvL1OErefGTlvmkiLw4cEsU37gNQIGll4cfy+PwKyJhxpMZMpml8z5nQv9YgraGyuzIy5ryMl20VfLix+G/AviliBwTkb056ZrKPGCMKaDhL708/HySDqBRvs1IPJ3V7gGqy8OMxJPzaJGyUPFTaXtL4FYovknZUa27hr+0JJ2xRJrmmvJJ9+U6/GrtF8NoPEVN+cRB0ZJ09GBYiHODMf78nv189k3bqa+MzLc5c4KfCP9TxphTuRfgU0EbpuTHcebeGv7SkTryafhOb3zN1LEYjacmHfgsSUc1/EI8frSXBw528fzZofk2Zc7w4/Avzb0hImGsISjKPJBMWc582bRWSKazEb2DDkGZzEg8NakzZnVUF2390NFv1XL0jy2fIj6v1gofFZFh4EUiMmRfhoEu4AeFNiwiFSLyKxF5VkSeE5FPzKLdy5aEE+G7LNouuW6ZiXQ2ondQDX+CZDpDPJWZHOFXqMP3w9kBdfhZjDGfNsbUAp8xxtTZl1pjTLM9r7YQceBGu9PmduDm3CZsSmk4ztxLw88YSC+RDJbxZGbaoq1zW3PxJypqVdIpng7H4Y+qw8/lRyJSDSAit4vIZ0VkfaEXGYsR+2bEviwNLzSPFNLwy+wDwVKI8lPpDIl0ZpqGX6kafhbHsddOydIZTaSxymgUNxyH3ze6fDKa/Dj8LwBjInI58BHgFNZQlIKISFhEnsGSge43xjxVsqUKUNjhR+37l4LDH09Z7yFf4RWohg9ks3Emp2WWkc4Y4qnF/x0ICmNM1uEPqKQziZRdcfs64HPGmM8BtX42boxJG2O2A2uBq0Rk29TniMgdIrJLRHZ1d3cXY/uyJFFo0dbuKplaApk6jmRT4SLpqIY/EeFX56Rl1mg/nYL0jCRI2AfEPnX4kxgWkY8CtwM/trN0ikpaNcYMAA8BN+d57C5jzE5jzM7W1tZiNrsscSL3crdK27IlFOEn8kf4KulM4Gj4U7N0ch9TpuNE9yGB/jGVdHJ5E9YC7DuNMZ1AG/CZQi8SkVYRabCvVwKvBA7OwFYFf3n4MJHNs5hxInhNy3Qn36Ktc314CbXJnm2cDJ0tK2qX1aJtwTJF28l/Nuf2afxp+KuBf7PPCELAd4wxPyrVUMUim5bp0S0TloakE8szzxagIhKa9PhyZjhPhO9c1wjfHScH/9K2Ou5/7vw8WzN3uDp8EXnMGHO9nXuf6z0EKwmnzmvDxpi9wI7ZMVNxcKpovZqnWc9bAhG+HcFPzcMXESojYWKalukS4Vuf16h+Pq50DMSoKS9jQ3M1w/GU5xS5pYRXP/zr7b++FmiVuSGZcvLwl76kM56cPsDcweqJrxH+qOeirX4+bnQMxGhrqKSxOgpYxVcraivm2arg8aq0bfK6zKWRygSFNfylI+k4Dn2qpAPWQWAhSjr3Pdc5p1LKSDxNNByivGxy8zRQSceLjv4YaxoqaKyy8k/6l0kuvtc5zG5gl/23GzgMHLGv7w7eNCUfhTX8JSTpuGj4zn0LLUtnf8cgd3xjNz/ae3bO9mk1Tpv8+dRUqMMvRMdAjLbGSpqqJiL85YBXa4WNxphNwM+A1xpjWowxzVhzbu+eKwOVyWQ1/AKSzlLomOk4/Iro9Pe6ECWdx472AFaO91wxtVMmTKRlah5+fkbiKQZjSdoaqiYknWWSqeNnleLFxpifODeMMfcCLwvOJMWLbC8d10XbpdNaYTzhpGVOX2paiJLO47bD75tD5zG1UyZAOGQtamuEnx8nJdOSdJwIXyUdhx4R+XMR2SAi60Xkz4DeoA1T8uM3D38pOPxshJ/n4LbQJJ14Ks2vT/YBcysPjOSJ8MHS8XXRNj9OSubaxkoaHA1/uUs6ObwFaMUaZn6Pff0tQRqluOOUgxfW8JeGpBMNhyjLc3CripYtqG6Ze04NMJ60/jcDcxgt5pN0AGrKNcJ3w6mybWuooiISpjoantOzsvnET+FVH/CBObBF8UHCZ5bOkojwE+lskdVUKqPhrINdCDxxrIdwSLisrX7OJZ21jVXT7rfGHKrDz0fHQIxIWFhRa43ObKiKaoSvLEwKTbxaUpJOngHmDpWR8IKK8B872sPla+tZ11Q1p90XR+PpaVk6YDn8YXX4eenoj7GqvoKQ3WiwqTqqi7bKwiSZzhAOCeFQfknHkT+WQh5+LM88W4eFlKUzNJ5k75lBrrughaaqyJxG+O6Sjkb4bpy1i64cGqoi9OmirbIQsUrA8zt7mJB0lkKlbSw5fbyhQ2U0TDyVIbMAJns9dbyPdMZw3QUtNFRFGRpPkZqDz98Yw2hiepYOBC/ppDOGRw53L8ohK1aV7YQM1lQdXTY98b166XwejwlVxpj3B2KR4kmiQM+PJTUAJc8Ac4dsi+RkOm+EO5c8frSHikiIHe0NHDw3BMBALElLTXmg+40l02QMrou2QWbp/OJgF//j67u4+z3XckV7Y2D7mW2S6Qznh8Zpa5hoo9BYFdVFW6wqW2WBkUxnXPvowBKTdDw0/NypVwvB4V+1sZnysnC2kGdgLBG4wx8Zn944zaE6GmyE39E/BsBzZ4cWlcPvHBwnY6CtcULSaayKMjy+PBqoeTVP+7e5NETxRzJlPL+US03ScfKkp1JpF2ONz3PxVdfQOEe6RnjjlWsB5rSQZyTbGnn6QbGmooxYMk06Y1zXe2ZC13AcgEOdQ7O+7SDJTcl0aKq2vmMDY0laa4M9SM83BUMjEWkF/gS4BMieBxljbgzQLsWFZDpDpMxDww8toQjfS8OPLIy5tk+dsIqtrt3cAlh6MMxNta0zz7amfPpBMdsTP5GirqKoAXW+OD9kOfyD54ZnfdtB4hRdrcmVdHI6Zi51h+/n/OVbwAFgI/AJ4CTw6wBtUjwopOGH7AyepaDhxxLeWTrAvKdmdg6OA7ChxYoYnTOSuVgEzDfP1iHojpldw9b7PnR+eFEt3E60VZgs6cDy6Kfjx+E3G2O+DCSNMQ8bY34fuLrQi0RknYg8KCIHROQ5EdHirVmgkIYPlqyzJBx+0iMPP7ow5tr2jSWIhCUbUU9E+MFLOvnm2ToE7vDtCH94PMVZ+6C3GDg3NE5TdXTSmWPjMuqY6cfhO9/ccyLyahHZAaz18boU8CFjzMVYB4j3isglJdqp2CTT3ho+WMVXS6K1go9F2/luoDYwlqChKoqIJbNVRsJEy0JzEuE7E63csnQguCEo54fHuXi1NfTOyUxaDPQMx2mdspg+lwfp+caPw/+UiNQDHwI+DHwJ+GChFxljzhlj9tjXh7FkobYZ2KpQOA8fHIe/uCP8TMYQT2VcJZ2FouH3jSayPdXBGr/YWBWZk2hxxCvCjwYX4cdTaQbGktywxVq3ONi5eHT8npE4LbXRSfctpwZqfnrpOIPHB4FXlLITEdmANd/2qVJer0yQSBVOHVsKks54yn34CSwcSad/NElj9eRFUSuvew6ydLzSMu37hsdn3+E7cs7m1hraGioXmcNPsKO9YdJ9FZEwVdHwstDwvQqvPmKM+Wu3Aiy/hVciUgN8D7jTGDPt3E9E7gDuAGhvb/dr97Ilmc4UzDuPhEOLPi3TceTuko71Gcy3pNM3luDClTWT7musmpvKTSd6r8pzUKwJUMN3UjJb68q5aFXtokrN7BmJ562PaKyKLoue+F6e44D9t+QCLBGJYDn7bxlj8k7JMsbcBdwFsHPnzsUvPAeMHw2/prwsG/0tVrK98AOWdPZ3DJJIZ0ouHuofTWQX/RyaqqMcmAMnOBJPU1Nelm0Clkt1TlrmbNM1ZC3Srqyt4KLVtTx0uJt4Kj1pru5CZCyRYiyRzu/wq+dGhptvvAqvfmj/LakAS6xVrC8DB4wxny3NPGUqfjT8hqoIA7HFHa3EPAaYA1REQohAbIYO7RM/fI7h8RQ/vfOGol+byRgGYslpDr+hKjInPfHzzbN1cCL8IMYcnrcd/oq6crauqiOdMRzrGuWSNXWzvq/ZpGfYcugtNdFpjy2X9gpeks4P8e6lc1uBbV8HvA3YJyLP2Pf9ae64RKV4CuXhAzRURjnRMzpHFgWD1wBzsBZHKyMzH3N4tGuEZNpgjMlm2vhleDxFOmOyhTsOTjOuTMbkjb5ni5FE/k6ZYB0QQxKcpFMWEpqqoly0qhaAQ+eHFrzD7x6xpKiWPMVVjVVRXugbm2uT5hwvSedv7L9vAFYB37RvvwWr+MoTY8xjQHDf9mWKnzz8hjnKEsklkzGMzGJVZyw7z9ZdJrB64pfu8PtGE1nddmAsOc1xF3y9/Rk3TVm0baiKkjFW2+SGquK2WQyjeebZOoiI3SJ59tc4zg/FWVFbTigkbGypJhoOWRW3O2Z9V7NKj+3wp6ZlgnWQXg4RvqvnsIusHgZ2GGPeZIz5oX35b8D1c2eikkuhXjoA9bakM5cVkN986hTX/dUvZq23TVbD93L40ZnNtT3ePZK9fsYuuS8Gx0FMlXQas2l+wco6I+OpbPplPmrKywKRdLqGx2mts1oTRMIhNq+oWRSZOo7Dd1u0nWlb6x8808EXHzle8uvnAj95+K0issm5ISIbsebaKvNAoV46YEk6iVRmTkcA/nR/J8PjqWxzqpkyXkDSASv6n4mkcyzH4b/QX/zp/EA2wp/i8Oeon47bAHOHoHridw3FWZkji1y0qpaDiyBTx9Hwp/6/gGxq7UzWvr715Gk+fe8BTvcuXGnIj8P/IPCQiDwkIg8BDwJ3BmqV4koinSEa9s6GyPZzic3NKWoskWbXyX5gojnVjLfpw+HPVNI53j2a7SR5pgSH7x7hT7RIDhJr+In751MdYIS/om6ywz8/FF/weew9I3HqKyNEy6a7vdnop9MxECNj4K5Hj5W8jaAp6PCNMT8FtmANMv8AsNUY87OgDVPy4yfCb6yaaPc6F/zqZF827//sLEX4YwXy8J3HZiLpHOseZXNrNXUVZSVJOs46ybRF26q5ifBH42lqKrwlndmO8OOpNP1jSVbWTnSb3Gov3C50WcfKwc+/pjLTttapdIbOoXEiYeE7u87QbdcqLDT8dvvfAmwFLgfeJCJvD84kxYtk2hRctK2vdCLMuXH4jx3pJhq2skJmS9JxHLlbHj5YxVczkXSOd4+wqaWGtY1VJWVo9I0miYZDVE85KDVUz80Bt7CkE571CN+psp0c4VvZOQu9AMut6AomJJ1SD9Jdw3HSGcM7rttIMp3hq48k10i5AAAgAElEQVSfKNnOICno8EXkY8Dn7csrgL8GCqVkKgGQzhjSmcKLto6kMzhHks6jR3rYuaGRVXUVc6rhW5JOaQ4tmc5wum+MzSuqWddUWVqEP5qgsToyLZ2ztryMspBks3iCIJnOkEhlqPFYtK0OIEvHqbJdUTcR4a+sK6e8LDRr//ug6BlJ5E3JhAldv9TsNue9X3dBC7dsW8U3fnmKofGZHfCDmNfsJ8J/I3AT0GmMeQdWlL+0pwQsUJz+OH4d/lyUincNj3Owc5jrt7TQ1lg5qxp+WUjy6q0OM5F0TvWOkcqYbIR/pj9WdFZT/9j0KluwUiIbAm6vMBp376PjEESWjlNluyLHcYoIq+srFnyb5HydMh0aZyjDnc1O0qrg3S+7gOF4in9/6nRphgLD40l2/q+f87f3HZrVbDs/Dj9mjMkAKRGpA7qATQVeowTAhMMvnKUDcyPpPH60B4CXXtDKmoZKzg7OlqTj3inTYSZZOk5K5uYVNaxrrCSWTNNb5I/dzeGDtY7SH2ADNacpmlsePkxk6cymw3Ai/JU5ET7AqvqK7DCYhch4Ms1wPOWq4VdEwtRWlJWsvZ/pnxisctnael66pYUvP3ai5CaGR7tG6BtN8PlfHOVv7zs8a/9DPw5/l4g0AF8EdgN7gF/Nyt6VonB63HtFvWBVWUbLQnOSpfPokR4aqyJcuqaOtoZKzg2Mk56FU9FYMu2Zgw8zy9I51m1VIm9qrWZtozWtqlgdv280kTfFD6yF3CAlHa9e+A415WWk7DbTs8X5ofFslW0ua+orF7TD98rBd1hRW56d5FUsZwdiNFZFsk39Xr+9je7hOKdKTNE82Wt9P2+4sJV/fPAon71/dpy+p+ew++F82hgzYIz5F+A3gP9uSzvKHJNI+ZN0RISGygiDAUf4xhgeO9LDtRe0EAoJaxoqSWXMrGQojCfdxxs6VEbDxFOZkg4wx7tHaK0tp64iwromy+EXq+P3j01vjezQWBWZI0nHIy3TPmDOZqZO13CcVrvKNpdV9RV0Dnkf7E/0jLLjk/dxYB4GpvSMOH103B3+yrqK7KzeYjk7EJs0NtEZeVlqu4YTPWOIwF1vu5I3v3gdn//FUb76+MmStpWLp+cw1iHl+zm3Txpj9s54r0pJ+NXwYW4aeB0+P0LXcDw7CKOt0frCdwzMvPBkLJHyJelAaS2Sj/eMsqmlGoC1tt3FFF9lMoaBscS0SNfBKtUP7vN3JlnVeqVl2m0uZnPh9vzQ+KQFW4fVDZWkMyYbSefjyeO99I8lue+587Nmj196ht376DjMJMLvGIjRluPwnSDiVG9pPa1O9oyypr6SikiY//1bl7GjvYF7nu4oaVu5+JF0nhSRF894T8qMSfjU8MHS8YOWdB490g3A9VuswmvnC98xMPNT+1gy40vSgdKGoBzrHmFTq9XHvrq8jKbqaFER/tB4kozBtVeOs2gbVHsLf4u2zpjD2Yvwu4fjkxZsHVbbB4FzHrLO82etyP6JYz1F7XPvmYEZy4QTko57byMnwi/2f2aMoaN/coTfWlNOVTTMqRIj/FO9o2y0A5JQSLhqQxOHOoezZ/ml4sfhvwL4pYgcE5G9IrJPRDTKnwecCL9QHj7MTYT/5PFeNrZUZx2984WfjUyd8UQ672CPXCqdIShFOvy+0QQDY0k2t1Zn71vbWFnU6beTzeGm4TdVRUllTCCVrjDhxL166QTRE//80Dgr66Y7/FX1tsP3SM10pJynTw/4/p89+8IAt/3j4/xk37kSrJ3Aj4bfWltOIpVhKFbc5zU0nmI0kZ4U4YsI7U2l1XcYYzjRM5qVhQC2tdWTSGc4fH5mxW1+HP4twGbgRuC1wGvsv8ock0xZkYdfSWcw4J74x3tG2bqyNnu7pryM+srIrFTbxpLuA8wdHElnLFncD9TpobO5dWJS1brGqqIOVG5Vtg7Z1NiAZJ1Rj3m2DtWz3BPfqbJdUTtd0nEO9m4RfiZjONg5zIbmKhLpDLtP9fva5wMHLPnnmRcGSrTaomckQW15mWchnyNVFSvrON8bR9J0aG+qKmnRtn8sydB4ig3NEwHJZW31gDWwZyb4aa1wKt9lRntVSiIr6RTI0gFLUgiyRbJzGruuafKXvK2hclYKcGI+F22h+Aj/eB6Hv7axkjMDMd/FLo4+76bhT5TqB/M/8Jpn61Bf6VT8zo4N3dmUzHzdJiOUl4XoHMrvLF/oH2MknuJt12ygLCQ87lPWefCQJRs+d3Zmjq57JO6p3wPZhnDFLtw6AU6upAOWwz/dN1Z0AZUzy8KRdJxt1ZaXsX+Gn4Pf1grKAsBvHj5YP/bxZGbW2hVPpXs4TjyVyS5OOaxpqJyVCH9gLEFdpffs3lI1/GPdo0TLQpMisrVNVSRSmeyQjEI4TbZcs3ScjplBOfxEimg45Jmiu6Z+9iQ2mHCE+RZts8VXLv97R87Zub6R7esaeOJYb8H9dQ2Ps69jkGhZiOfPDs1oPaRn2L2PjkOpEb5Te9I2xeGvb64iXsR3yuGk7fA35Dj8UEi4tK2OfR0zy3AKzOGLyFdEpEtE9ge1j+VGsRo+EJis42S0rJ1yGrt2Fqptx5NpekYSWYflRlbSKSHC39hcne2UCRPvw2/XzKyk41F4BcF1zPQab+hQGQ3TXB2dtZYH3cPTq2xz8Sq+ev7sECGxGq1du7mZfWcGCrYeeMiO7t+0cx1D46mS2l84ePXRcVhRYoTf0R8jWhaieYq8125LMsXKOqd6RwmJJTPmcllbPQfODZVczAXBRvhfA24OcPvLjqLSMgOutnV+fFO/lGsaKhiOp2bUR8TRgaeeIk+l1LTM492jbMpZsIWJ9/FCnz+n0jeWIFoWcp3I1ZTtiR+Uhu/dKdOhrbG0PkH5cBzh1Cpbh9X1la4a/vPnhtnUWkNFJMw1m1vIGPjV8T7P/T10qIuVdeX89pVrAW/9uns4zl/de5BDLh07e0YSBR1+dXkZNeVlxWv4AzHW1FdMq01YX2Jq5oneMdoaK6edvW1rqyeRynC0a8TllYUJzOEbYx4BvP+jSlEkily0heAiTCf7YOpCVVuD9SWfSZTvpolOpaIESSeZznCqb2yawy86wh+1cvDd5uDWVUQISXCf/0jce9qVw9rG2VlTAUvqyFdl67C6voLzLsVXB84NcfFqq6vmjvYGystCnrJOMp3h0cM9vGLrCi5aVUs4JDx31l3OuHvPGf7l4WO86u8f4d3f3D2puCuRyjAYSxZ0+GDn4hcb4Q/Epv0OwPr+hqT44quTPaOTFmwdttkLt/tmsHCrGv4iIivpFOiHDzkLdgFJOmf6Y7TURLOl5A5rGqzobyY6/kQjqkIRvrXvYjpmdvTHSGfMtB9URSRMa225/wh/1HsGbigkgS6ce82zzaWtwZLYZqMe4IW+GKvyRLIOq+srSGUMvVM068GxJB0DMS6xHX5FJMyLNzR55uPvPtXPcDzFy7euoCISZsuKGs+F230dg6yur+D9N17AY0d6uOVzj/LT/Z0A9I46RVeF5wuvqCu++OrsQCyv/BgtC7GmobKoXHxjDCd7Rict2DpsbK6mOhqeUabOvDt8EblDRHaJyK7u7u75NmdBU2ylLRBYe4UX+seyPWhymai2nYnDt35wK+u9I7L6ygiRsGQbevnB6VGyIc8PysrU8ffjHBhLTBtePpXm6mjR0aJfRgv0wndoa6gknspkWwvMhOM9E8Vq+Vhdnz8184DdJ//i1RMpvNdsbuZg5/C0g4PDg4e6iISF6+0q7kvW1LHfI8Lf1zHI9nUN/NFvbuWxP7mRTS3V/Osj1uQpZ7Shvwi/oqjvUyKVoWs47no2ur65uNTM3tEEw/FU3gg/FBIuXVO/uB2+MeYuY8xOY8zO1lYdletFMQ4/O2YvoGrbM/2xaQu2AC3V5UTDM+uNfnYgRmttOeVl3ouS4ZCwrrG4H9RpO9pa3zT9YLWusaooDd+tytZhfXN19gAz2wz7jfDtg/JMZR1jDCe6J9pR5CNbfDWlY6pTYetE+ADXbm4G4JfH88s6Dx3s5qqNTdn3eOmaerqH43mj78GxJKd6x7KSR31VhNuvXs/TpwfY3zHoq+jKYUVtOeeHxn2fEXUOjmPMdGnTob2puihJx9H780X4YMk6z58byg5bv+fpM763DQvA4Sv+SaT9a/hV0TCRsATSEz+dMZwdiE1LyQTsJmoVM9PwB2MF9XuH9c1VRTnVkz1jVNryzVTWNloppX7K+B0N34tNrdWc7B2ble6hU/GTpQMTaxMzzZzqGo4zmkhPW/vIZXV9/vYKB84N0VITnfSZX9ZWT015GU/mcfgdAzEOnR/mFVtXZO+7dI11sMin4zu56S9aW5+977evXEtFJMQ3nzyVTYt064Wfy8q6CsaTGYZ9Fqt1FJAf25uq6B1N+C5+O9FjHRzynYECbGurYzyZ4Vj3KJ2D4/zlD57ztV2HINMyvw38EtgqImdE5J1B7Wu5kEz5T8sUEeoro4Fk6XQOjZNMm7wRPsw8F/+snfXgh/XN1ZzqHfMdkZ3qHWV9c1XexdYNLdWkMqZgRJbOGAZi3ho+WFFaIpUp6rNIZ4yv9zIaT/uTdIpcjHbDqU7e1OIu6TRVR4mWhfJKOhevrpv0mZeFQ+xob2DPqekVtA/b6Zgv3zpxxn+J7fCfz+PwnUXMbWsmHH59ZYTXb2/j+890ZAuZ/Gr4MDHopRCFEgzWNxeXqXOyZ5RwSFx/W5flLNx+9O69RadoBpml8xZjzGpjTMQYs9YY8+Wg9rVcyEo6PhZtwWmvMPuSzhnbIU5NyXSYSbWtMYazA+NFRfgj8ZTv4SUne/NnQADZNhGHCvQrGYwlMQaaqrw1fEf+cBxOIYwxvO3LT/GB//uM5/NG4ilG4ilfEkVdRYTairIZSzrHc+YHuOEUX+U6/GQ6w+HOkUlyjsMV7Y0c7ByaFv0+cayH1fUVkyqh6yoirG+uyqtf7zszyLqmymkH4NuvXs94MsO3njxFVTQ8LcEgH07bCL9rL87nutolQGlvKq5N8oneUdY2VrqexW9qraEyEuYfHjjCg4e6+ZObL/K1XQeVdBYRxWj4AA2VwTRQe8HJwc8j6YAV7XQNx0vq7DcwliSWTPt2+BuyxS2FnWo6Y3ihL8b6lvx2b1lpOZjDLrncDoX66DhsbC3O4e861c8Tx3q5//nzxFPuqaZOrvmFOX2MvHAydWbC8e5RKiNhVrnk4Dusrq+gM0fDP949SiKdyaZk5nLl+kYyBp45PRHlG2N48ngfV29qnnYWdumaurySzt6OgWzkm8u2tnp2tDcwNO7v4Ag5Eb7PhduzAzFaaspde/S0ZyN8fw7fLSXTIRwSLllTx+m+Ma7a0MR/v2aDr+06qMNfRDgafplLWtxUguqYeabfGs7gpGBOpa2xEmMoaQKSU6buX9Lx/4PqHBonkc6wvin/D6oqWkZ7U1XBCD/bVqGAht9aU05NeZlvh3/XI8cBq5Bs10n35mKOw79olT+H78zsnQknekbY0FLtmpLpsLq+MptlBfD8OSsiz+fwt7c3IMKkRmrHukfpGYlz9aamac+/dE09p/vGJlWPD4wleKEvxmVtDXntedvV6wHvtsi5TFTb+vvuuuXgO9RVRGisivhKzfRKycxl5/pGKiNh/vqNLyr4/5iKOvxFRDKdIRoOuRb7TKW+MhpIa4UX+mKsrK1wzaJxFrD8pjjm4jgLvxH+2sYqQgInfTj8U06Pkub8ET7AhStrOHLeu5KxUGtkBxFhY0s1x304/GPdI/z8wHneef1GImHhkcPuKcoHO4eoKS9z1Xmn4hRfzSQX/3jP9OrkfKyyi6+chmEPHuymprws72vrKiJsXVnL7tMTDt/J2rl6U/O051+aR8d39Pt8ET7ArZetpqk6ymqf36ea8jKqomHfEb41+MQ7OGlvrua0j+9nz0iC0UTa8/sJcOcrL+QXH36Z68KuF+rwFxHJVMZX4zSHhoDG7J3pH/N0NhessKQRtzJ3L/xW2Tpki1t8SDrOQWG9xw/lwpW1HOse8ZSj/Eo6YC3cOt05vfjSoyeIhEO862WbuXJ9Iw97OvxhLlxZ4/vA39ZQyUg8VXSfd4d4Ks0LfWNs9uFg1tjFVz2jcU71jvKjvWd560vaXWXIK9Y38vSp/uwB4snjvayur8hq37lcai/K5hZgFXL4FZEw/3HH1fzZrRcXtB2sg7STmlmIVDpDR3+sYIHgertrZiG8akRyqYyGszUPxaIOfxGRTGd8tUZ2aKiMMJpIz3hKzlTO9OdPyXRYWVfBqrqKknqYnx3M34jKiw3N1f4i/F6rS+ZqDx1666paUhnjmepZqDVyLptaq+kYiHl2Le0ZifO9PWf47SvW0lpbzg0XtnKwczhvpogxhkOdw1yURyJxI9s2osTRk6d7x8gYPIuuHFY5xVcD4/zLw8cpC4d45/UbXZ9/ZXsjw/EUR7pGMMbw1PHevPo9WANK1jZWcs/THdk89H1nBmlvqqLeYwF9y8pa3wEEWF0z/UT4BzuHiacy2fx/N9qbqugYiBXMqLl3X2fW3qBQh7+ISKSN7wVbgIbq2S++SqYznBuMsa6AnLB9XUNpDn9gnNUe5fv5sKoZ/UT4o7Q3VXlu21kI9To7GRhLUF4WKjigBawI3xg8I7yv//IUiVSGP3ip5RhvsEdGPnJkeuuBzqFxBmNJ3/o95KZmlqbjH/ORoePgZKs8e2aA7+0+w+9cuTZvO2WHK9c3ApaOf6x7hJ6RRF793uHPbr2Y584Ocdej1nrHvo5BLlvr7XCLZUVtebb3vxd7bCnqivZGz+e1N1dla1fc2HdmkK89cYK3vqS94BnDTFCHv4hwNHy/NFTOfnuFcwPjZAx52yrkcvm6Bk71jmUXOP3i1pfEiw3N1QyMJQu+z1O9YwX10U2tVtvkIx4Lt72jiYL6fXZ7dt66k9Y4lXgqzTd+eZJXXrwym4Z4yeo6WmqieXX8g9kFW/8RfnbWcIkOP99ADjcch/+5nx8hlcnw/92w2fP565uraK6OsutUH7+0u2des6nF9fm3XLaaWy9bxd/ff4Rfn+zjTH/MVc4plRW1Fb4knT2n+rNnHV44Vd1uU75S6QwfvWcvzTXlfKTINMtiUYe/iEimi9fwYXYbqGX74DcVjvABnjlTXJR/dsB/la1DNvWtzz3KN8ZwqneMdpcMHYfysjAbmr0zdfaeGZiUI+6FM5fULVPnqeN99I8lefOL12XvC4WEl25p5bGjPdOmJR08Z9m1tYjT/qbqKBWR0ttdHO8eobW2nNoK77oDZ1/RshC9owluu3xN9n/jhohw5fpG9pzq58ljvaypr5g2RW0qn7htG9XlYe74+i7AXb8vlZV15Ywl0gWrY/ecHuCK9oaCaynb2urZ2FLNh/7zWT75w+endXf92hMn2d8xxCduuzTb9DAo1OEvIvpGE76KRxyC6InvVGy6FV05vGhtPSGZnGNdiFQ6w/mhcdd0TzecvGUvHb97OE4smZ40GNqNratqOeySqdM1NM7h8yPZpl6FqK2I0Fpbzome/Nt74MB5KiKhadu74cIW+kYT00baHeocYnV9hadmPRURYa3HzN6DnUP8yXf3umZ0He/x7qEzdV9OlP/ul1/g6zVXrm/kZO8YjxzudtXvc2mtLefjt12abRuSW2E7Gzi5+F5Rfs9InNN9YwXlHLD67P/ofdfz9qvX85XHT3DL5x7hS48e50d7z/LgwS7+9r7D3HTRCm7ZtmrW3oMb/r2HMq+M27nZb8qJBAsRRE/8F/pihEPiWlnoUF1expYVtUXp+OeH42SM/wwdByej45RH+mM2Q8ejqMXhwpW13Lu/k/FkelpBzWNHLV39+gv8OXywKm7zRfjGGB442MV1m1um7eeljo5/uJsXrZ3IMT/YOVyUfu/Q1pC/E2gileED336GQ+eHGU+l+dybd0x7zvHuEW7ettr3vq7d3MyLNzSx1aedjo4/HE/lTcfMx22Xr+HefZ280D9W1MHPDytzqm3dzuT22PLMFesLO3ywfg+feN02bt62mo/evZdP/fhA9rGqaJhPvn6b76yrmaAOf5Gw62Q/sWSaGy7072jqAxhz+EL/GKvrKyjzsZawfV0DP3u+E2OMry9zsSmZDpVRqwLUK8LPprwVkBjAcvjGwNGukWkZGI8d6aG5Opq3VYAbm1qrue+589PuP3x+hDP9Md6TJxJuqSnn0jV1PHK4hz+8cQtgSXrHukd4eU5TMb+0NVayN4+89k8PHuXQ+WFuvGgFP3jmLK/YuoLX72jLPt4/mqB/LOk7wgf49BteVJRt29rqiYSFZNr4dvgiwj+99QpSmdnNQIPcalv3CH/P6QHKQlK0nHTN5mYe/PDLGYwl6Rwa5/xQnBW15YEu1Oaiks4i4eHDXUTDId8/CIDa8jLCIZllSSd/W+R8bG9vYMBuXeuHicEnxUk6YOn4Xpk6p3vHKAuJrx+Wk6lzeIqOb4zhsaM9XHtBS1FZRBtbqukdTUxbVH7goHUQuPGi/A78FVtXsPt0f3ao9fHuUZJpM6mvvF/aGirpH0symqNLH+wc4p8ePMrrt6/hi2/fyc71jfzF9/dP6vviFI35ydAplYpImG1t9b70+1zCISnYQrsUWn3009lzup9L19S5tlTwQsQajnPRqjpedmFr3irkoFCHv0h46JDVH7wYDd/qmBmZ1alLL/SNFdTvHbILtz5lHafKtpSikg3NVZ7l6yd7R2lrrPR1ZrKhuYpoODRt4fbw+RG6huNcf4H/gy7ARjtT58SUA9IDB7rY1laX7SM/lbdfs55IWPjs/YcBy0EDvqWSXNZOGUyTSmf4yHf3Ul8Z4S9feynhkPB3b9qOAf7oO89kWzo7RWN+cvBnwqdev43P/7cdcyJrFKKuooyKSMg1wk+mM+w9M8AOH/r9QkMd/iKgYyDGka4RXnZh8QNiGiojs5alc+DcEF3D8WyTsUJsWWF19vPv8GPUV0Z8tf2dyvrmarqH45Mi2FxO9Y750u/Bat27eUXNtCZqjx6x0iSv31Lc/2FjtmvmxMJt32iCPaf7uemila6vW1FXwTuu28h/PXuW588OcbBzmLKQeLYodiO3L/54Ms2nfnyAvWcG+cTrLs2mmK5rquKTr7uUX5/s521fforj3SMc7xmlLCQF6y5myqVr6rlyvXv+/VxiVdtWZIe2T+XguWHGkxnf+v1CQjX8RYCTj/2yrSU4/KrIrOXh/939h6mtKONNO9t9Pb8sHOKytfVFOfxi9XuHia6ZY9ne6Q7GWJWzO9rzN9jKx9aVNfx6SgOzx472sKmlumi9tb2pinBIOJGTi//gwS6MgZsu9tbj33XDZr715Cn+5r5DGGO4YEUN0SKqrR2c4fI/2nuOT/7oeU70jHL71e28+rLJi7G/taONWDLNX917kJv//lGaa6K0N1f5OjNaSqxrqmT3qf68C/cTBVf+v08LheX1X1ykPHyom9X1FWxZUXxk11AVnZVK231nBrnv+fP8wfWbisqK2LGugefPDnm2+3U4Ozjuu0vmVLwGTfSPJRkeT/mO8MEqb+8YiDE8bh0s46k0Tx3v852OmUu0LMS6xspJTdR+cbCLFbXlBVMK66sivOvlm/nFwS5+eby3JDkHrOrRSFj43p4zZIzhG++8ik+9/rJpEoqI8NaXrOeBD72Mm7et4tzgOBeuCK7Uf6Hy3pdfQMdAjH9+8Oi0x/ac7p/ThdbZJFCHLyI3i8ghETkqIv8zyH0tVZLpDI8f7eFlF7aWpG82VEboHBz3PWLNjb/7+WHqKyO84/oNRb3u8nUNJNIZDpwr3EhtJhG+U+CTL1PncTuVMt8cWzecwqYnjlndG/ecGiCWTBeVjpnLxpZq9nUM8uiRbk7bOec3XrTC1+LvO67dyIracsaTmZIdfigk3HHDJu585RZ+ducN2bRPN1bUVvAPb9nBD957HR+77ZKS9rmYufaCFl6/fQ1fePhYdtqXw57T/VzR3rgg1huKJcgRh2Hgn4BbgEuAt4jI8vvmzJCnTw8wHE+VpN8DvPbyNfSPJbn9S0+VLO3sOd3PLw52cccNm6jzUW2Zi7Nwm292aS7dw3EGY8mSHX5dRYTm6ihHuiYfWB44cJ4PfedZLl5dx7VFLLbu3NDI6voK3vXN3fzpPfv4yb5zhEPC1ZuLW7Cd2F4Tp3rHeNuXf8UNn3mQ4XiKmy521+9zqYyGed9NVmrmpTMoMvrjV13Ena+8sKjMksvXNZTcmXGx82evvoSKSJi/+P7+bGvpg51DvNAX44r1i0/OgWA1/KuAo8aY4wAi8n+B1wHPB7jPJcfDh7sIh4RrS4wsX3HRCr7w1iv4w39/mjd/8Um+8c6rfE//cfi7+w/TVB3l967dUPT+V9dXcMnqOv7q3oPsOzPIH79q66T2r8l0hm8+eYrPPXAEEXjxhtIXwq5c38jdezroHo7zR79xIZ2D47zv209zyZo6vv77VxVXpVwV5b4P3sDf3X+Ef/vlSdIZwxXtDUUf8Bze8/LNvPHKtRzvHuVEzyjD48lJM1sL8dar2tnYXM21JR5wlOJprS3nI6/ayl/84Dn++aFjHOwc5sd7z1IZCXOjx2L7QkZmMhTBc8MibwRuNsb8gX37bcBLjDF/6PaaCy653Nz6F18LxJ5g3qW97YA+Q4CnTvSxobmK/3zXtTPazqNHuvkfX99Fc3U5W1bWEBYhFBJCYuUzh8S6OMSSaToHxzk3OE7PSJw/vfUi7ijQCMuN0XiKux45zhcfPU4ileH6LS1EwyFCIhw6P8yJnlGuu6CZP7v1kmkLrsUQS6T55pOn+MLDx+gbTSBirSF87fevKtlRgzVw47P3H+a27Wu47fI1JW9HWXykM4Y3/PPjPHtmkOpomNuvWc87r9+YnX27EBCR3caYnb6eG6DD/x3gVVMc/lXGmPdNed4dwB0AzWsSfrwAAAdhSURBVG0br7zszq8EYg9AoIpbQBsX4P03beF129sKPrcQvz7Zx9/87BCxZJp0xpAxkMkY0saQyZhJB8VoOMSq+grWNFjDpN9+zYaSskNy6Roe5/MPHOXXJ62uiMZAbUUZ7375Zm68aMWsaaKj8RRfe+IkJ3tG+dhtl1JTQpqnojic7h3jgYPn+a0dbTT4mIEw1ywUh38N8HFjzKvs2x8FMMZ82u01O3fuNLt27QrEHkVRlKVIMQ4/yCydXwNbRGSjiESBNwP/FeD+FEVRFA8CO9c1xqRE5A+BnwFh4CvGmOeC2p+iKIriTaDipjHmJ8BPgtyHoiiK4g+ttFUURVkmqMNXFEVZJqjDVxRFWSaow1cURVkmqMNXFEVZJgRWeFUKIjIMHJpvOxYwLUDPfBuxgNHPxxv9fLxZrJ/PemOMr8ZMC63m/JDfirHliIjs0s/HHf18vNHPx5vl8PmopKMoirJMUIevKIqyTFhoDv+u+TZggaOfjzf6+Xijn483S/7zWVCLtoqiKEpwLLQIX1EURQmIBeHwddi5NyLyFRHpEpH9823LQkNE1onIgyJyQESeE5EPzLdNCw0RqRCRX4nIs/Zn9In5tmmhISJhEXlaRH4037YEybw7fB127ouvATfPtxELlBTwIWPMxcDVwHv1+zONOHCjMeZyYDtws4hcPc82LTQ+AByYbyOCZt4dPjnDzo0xCcAZdq7YGGMeAfrm246FiDHmnDFmj319GOtHO/N5kEsIYzFi34zYF128sxGRtcCrgS/Nty1BsxAcfhvwQs7tM+gPVikBEdkA7ACeml9LFh62ZPEM0AXcb4zRz2iCvwc+AmTm25CgWQgOP9/kao0+lKIQkRrge8Cdxpih+bZnoWGMSRtjtgNrgatEZNt827QQEJHXAF3GmN3zbctcsBAc/hlgXc7ttcDZebJFWYSISATL2X/LGHP3fNuzkDHGDAAPoWtCDtcBt4nISSw5+UYR+eb8mhQcC8Hh67BzpWRERIAvAweMMZ+db3sWIiLSKiIN9vVK4JXAwfm1amFgjPmoMWatMWYDlu/5hTHm9nk2KzDm3eEbY1KAM+z8APAdHXY+GRH5NvBLYKuInBGRd863TQuI64C3YUVmz9iXW+fbqAXGauBBEdmLFWDdb4xZ0umHSn600lZRFGWZMO8RvqIoijI3qMNXFEVZJqjDVxRFWSaow1cURVkmqMNXFEVZJqjDVxRFWSaow1cURVkmqMNX5gwR+biIfNi+/oTH8xpE5D1zZ9mkfW8QkZjdaGy2tvmE/XfW3peIvN+eAfAtn8+vtIvSEiLSMhs2KIsPdfjKvGCMudbj4QZgXhy+zTG70diskPNeZ/N9vQe41RjzVp82xOz3pH2qljHq8JVAEZE/s6eZ/RzYmnP/iP23WkR+bE9j2i8ibwL+CthsR6SfsZ/3fRHZbU9susO+b4Md5X7Rvv8+u1cMIvJ2Edlrb/cbOfu93Z7+9IyI/Ks9gMfve9mQO3VMRD5sn7W42pH7Xqe+L5f3PnWff2Q/tl9E7rTv+xdgE/BfIvLBPK+5XEQeEZHnRSQjIkanXCkAGGP0opdALsCVwD6gCqgDjgIfth8bsf/+NvDFnNfUAxuA/VO21WT/rQT2A83281LAdvux7wC3A5cCh4CWKa+9GPghELFv/zPw9in7mbZvt8eADwMfd7Mj53kjLq+f9t5dPr9qoAZ4DthhP3bSeX9TXlOB1RjtKvv2/w98hok2Knlfp5flcdEIXwmSlwL3GGPGjNWjPl8X1H3AK0Xk/4jIS40xgy7ber+IPAs8idVOe4t9/wljjKO378ZyqjcC3zXG9AAYY5xpYTdhOdFf2xr9TViR8myQz45CFHrv12N9fqPGmlh1N9Zn6sUrgT3GmF/Zt/diHfC0aZaiDl8JHE9HY4w5zEQk+2kR+cupzxGRl2M5smuMNZf1aaxIFqx5rQ5poAxrqE6+/Qrwb8aY7fZlqzHm40W8lxSTfzMVOdfz2eGJj/eebzhQIbbZ23O4AthTwnaUJYg6fCVIHgF+y84QqQVeO/UJIrIGGDPGfBP4GywHNQzU5jytHug3xoyJyEVYw8q9eAD4XRFptvfRlHP/G0VkhXO/iKwv4v2cB1aISLOIlAOvKeK1MOV9ubz3XB4BXi8iVSJSDfwW8GiBffQCL7K3fyHwBqzBHopSOApRlFIxxuwRkf8AngFOkd9ZXQZ8RkQyQBJ4tzGmV0QetxdI7wX+HHiX3c/9EJas47Xf50TkfwEPi0ga64zg94wxz4vInwP3iUjI3t97bdv8vJ+kiHwSa2buCYocIpLnff186nuf8vw9IvI1wJFnvmSMebrAbr6NNcFpP9ADvMUY01uMncrSRfvhK0oOYg1C/5ExZknOfBVrlN9OZ31DWV6opKMok0kD9bNZeLUQcAqvgAiQmW97lPlBI3xFUZRlgkb4iqIoywR1+IqiKMsEdfiKoijLBHX4iqIoywR1+IqiKMsEdfiKoijLBHX4iqIoywR1+IqiKMuE/wfS3Fo/kbUJYQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_bin = 300  #number of bin in the graph\n",
    "box   = cell_side   #side of the bon in the simulation \n",
    "g_20  = np.zeros([n_bin,n_bin])\n",
    "\n",
    "#average during the thermalized part of the simulation\n",
    "for i in range(n_step):\n",
    "    g_20 += g_step(box,n_bin,dist[i])\n",
    "g_20 /= n_step\n",
    "\n",
    "plt.plot(g_20[0],g_20[1])\n",
    "plt.xlim(0,box/(2*sigma))\n",
    "plt.xlabel('distance [ units of $\\sigma$]')\n",
    "plt.ylabel('radial distribution function')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 400 K\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory = read('production_bulk_400-pos-1.xyz', index='::2')\n",
    "N = len(trajectory[0])    #number of atoms\n",
    "\n",
    "n_step = len(trajectory)  #number of step\n",
    "n_step\n",
    "#create a distance array\n",
    "dist = np.empty([0,N,N])\n",
    "for frame in trajectory:\n",
    "    frame.set_cell([cell_side,cell_side,cell_side])\n",
    "    frame.set_pbc((True,True,True))\n",
    "    dist = np.append(dist, [frame.get_all_distances(mic=True)],axis=0)\n",
    "\n",
    "#Lennard Jones units\n",
    "\n",
    "dist  /= sigma\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "n_bin  = 300  #number of bin in the graph\n",
    "box    = cell_side #side of the simulation box\n",
    "g_400  = np.zeros([n_bin,n_bin])\n",
    "\n",
    "#average during the thermalized part of the simulation\n",
    "for i in range(n_step):\n",
    "    g_400 += g_step(box,n_bin,dist[i])\n",
    "g_400 /= n_step\n",
    "\n",
    "plt.plot(g_400[0],g_400[1])\n",
    "plt.xlim(0,box/(2*sigma))\n",
    "plt.xlabel('distance [ units of $\\sigma$]')\n",
    "plt.ylabel('radial distribution function')\n",
    "#plt.savefig(\"gr.png\",  dpi=300,transparent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(g_400[0],g_400[1])\n",
    "plt.plot(g_20[0],g_20[1])\n",
    "plt.xlim(0,box/(2*sigma))\n",
    "plt.xlabel('distance [ units of $\\sigma$]')\n",
    "plt.ylabel('radial distribution function')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data38  = np.loadtxt('production_38-1.ener' )\n",
    "data150 = np.loadtxt('production_150-1.ener')\n",
    "data450 = np.loadtxt('production_450-1.ener')\n",
    "data817 = np.loadtxt('production_817-1.ener')\n",
    "dataBulk = np.loadtxt('production_bulk_20-1.ener')\n",
    "\n",
    "data38   = np.transpose(data38)\n",
    "data150  = np.transpose(data150)\n",
    "data450  = np.transpose(data450)\n",
    "data817  = np.transpose(data817)\n",
    "dataBulk = np.transpose(dataBulk)\n",
    "\n",
    "\n",
    "#energy in simulation time in LJ units\n",
    "E38   = 27.1442*data38[4]/epsilon #eV / epsilon(eV)\n",
    "E150  = 27.1442*data150[4]/epsilon \n",
    "E450  = 27.1442*data450[4]/epsilon \n",
    "E817  = 27.1442*data817[4]/epsilon \n",
    "EBulk = 27.1442*dataBulk[4]/epsilon \n",
    "\n",
    "#compute average energy \n",
    "e_average    = np.zeros([2,4])\n",
    "e_average[0] = [38,150,450,817]\n",
    "e_average[1] = [np.mean(E38),np.mean(E150),np.mean(E450),np.mean(E817)]\n",
    "bulk_average = np.mean(EBulk)/864.\n",
    "\n",
    "#compute average energy per particle\n",
    "e_part    = np.zeros([2,4])\n",
    "e_part[0] = [38,150,450,817]\n",
    "e_part[1] = [np.mean(E38)/38.,np.mean(E150)/150.,np.mean(E450)/450.,np.mean(E817)/817.]\n",
    "\n",
    "#create an array of crystal energy per particle\n",
    "ecr_part  = np.ones(4)*bulk_average\n",
    "\n",
    "#first plot\n",
    "fig, ax = plt.subplots(1,2, figsize=(12,4), sharex='col')\n",
    "ax[0].plot(e_part[0],e_part[1], linestyle='-', marker='o', color='brown',markersize=12, label='Cluster')\n",
    "ax[0].plot(e_part[0],ecr_part*np.ones(4), color='green', label='Crystal')\n",
    "ax[0].set_xlabel('cluster size')\n",
    "ax[0].set_ylabel('<Energy> per particle [ units of ε]')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#compute surface energy\n",
    "e_surface    = np.zeros([2,4])\n",
    "e_surface[0] = e_part[0]\n",
    "#surface energy = (cluster energy) - N x (energy per particle in a crystal ) \n",
    "e_surface[1][0] = e_average[1][0]-bulk_average*38\n",
    "e_surface[1][1] = e_average[1][1]-bulk_average*150\n",
    "e_surface[1][2] = e_average[1][2]-bulk_average*450\n",
    "e_surface[1][3] = e_average[1][3]-bulk_average*817\n",
    "\n",
    "\n",
    "#from bulk knowledge:\n",
    "vol_part = (unit_cell_side)**3/4.\n",
    "\n",
    "#compute the ray approximating the cluster to a spere\n",
    "r38  = np.cbrt(vol_part*38*3/(4*np.pi))\n",
    "r150 = np.cbrt(vol_part*150*3/(4*np.pi))\n",
    "r450 = np.cbrt(vol_part*450*3/(4*np.pi))\n",
    "r817 = np.cbrt(vol_part*817*3/(4*np.pi))\n",
    "#and the surface\n",
    "S38  = 4*np.pi*r38**2\n",
    "S150 = 4*np.pi*r150**2\n",
    "S450 = 4*np.pi*r450**2\n",
    "S817 = 4*np.pi*r817**2\n",
    "S    = np.array([S38, S150, S450, S817])\n",
    "\n",
    "#plot <Surface Energy> per surface unity\n",
    "ax[1].plot(e_surface[0] ,1000*epsilon*e_surface[1]/S, linestyle='-', marker='o', color='brown',markersize=12)\n",
    "ax[1].set_xlabel('cluster size')\n",
    "ax[1].set_ylabel('<Surface Energy>/S [ meV$/\\AA^2$]')\n",
    "\n",
    "plt.subplots_adjust(wspace=0.3)\n",
    "\n",
    "ax[0].legend(loc='upper right')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}