{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Wang Landau algorithm\n", "\n", "\n", "\n", "\n", "\n", "The algorithm is very useful for studing the phase transition phenomena because it does not suffer from the critical slowing down and gives direct access to thermodynamic potential and its derivatives\n", "\n", "\n", "\n", "Consider the 4-site Ising model.\n", "\n", "\n", "\n", "\n", "Many energy levels are highly degenerate.\n", "\n", "In case of periodic boundary conditions, the allowed energy levels of\n", "the ising model are: $-2J N$, $-2J (N-4)$, $-2J(N-6)$,$-2J(N-8)$\n", "... $2J N$. The number of\n", "energy levels is of the order of $N$ while the number of all possible\n", "states is $2^N$. Therefore many states have to be heavily\n", "degenerate. In Markov-chain simulation at certain low temperature $T$, most of\n", "states visited have the same energy close to $\\langle E \\rangle$ ( the\n", "probability for a state is $e^{-E/T}$).\n", "\n", "\n", "The classical MC simulation generates canonical distribution at a\n", "given temperature $P(E) \\propto g(E)e^{-E/T}$. The idea of Wang-Landau is\n", "to estimate the many-body density of states $g(E)$ directly. The\n", "temperature is not required for the simulation and by analyzing $g(E)$\n", "one can study thermodynamics at *any* temperature.\n", "\n", "\n", "If the density of state $g(E)$ is known, the free energy is simply given by \n", "\\begin{equation}\n", " F(T) = -k_B T \\log (Z) = -k_B T\n", " \\log{\\sum_E g(E) e^{-E/T}}\n", "\\end{equation}\n", "\n", "In the Wang-Landau algorithm the random walk in energy space is\n", "performed and the probability to visit a state with energy $E$ is\n", "taken to be $P\\propto \\frac{1}{g(E)}$. The resulting histogram \n", "$H(E) \\propto g(E)P(E) \\propto const$ is therefore flat in the\n", "simulation.\n", "\n", "In our example above, a simple random walk would visit a state with $E=0$\n", "six-times more often than the other two states since there are\n", "six-times more states at $E=0$. If we assign probability for the state\n", "at $E=-4J$ to be proportional to $1/2$, at $E=0$ to $1/12$ and at\n", "$E=4J$ to be $1/2$, we will visit all three states on average the same\n", "number of times. Hence the histogram of visits is flat.\n", "\n", "\n", "Let us sketch of the algorithm which generates the density of states\n", "$g(E)$.\n", "\n", "The transition probability is taken to be\n", "\\begin{equation}\n", " P(E_1\\rightarrow E_2) = min\\left[\\frac{g(E_1)}{g(E_2)},1\\right].\n", "\\label{Ptr}\n", "\\end{equation}\n", "We always accept moves with smaller density of states but we often\n", "reject moves that would bring us to the state with high density of states.\n", "\n", "\n", "\n", "\n", "- Start with $g(E)=1$ for all energies of the model.\n", "- Try to flip a random spin and accept the move according to\n", " probability $P(E_1\\rightarrow E_2)$.\n", "- Every time we accept a move, the density of states for the\n", " current energy is updated by multiplying the existing value $g(E)$\n", " by a modification factor $f>1$, i.e., $g(E)\\rightarrow g(E) f$. \n", " \n", " This modification factor serves like a temperature in the simulated\n", " annealing algorithm. At the beginning it is taken to be relatively\n", " large $f\\sim e^1 =2.71828$ and later it is reduced to unity. Because\n", " of this large factor, the density of states will grow very fast for\n", " highly degenerate states and we will quickly visit all possible\n", " energies. However, the density of states obtained by large $f$ is\n", " not very accurate. The error is equal to $\\log(f)$.\n", "\n", "- The simulation at certain modification factor $f$ is running for\n", " so long that the histogram becomes flat, i.e., each energy level was\n", " visited equal number of times. Here equal means that the deviation\n", " is within certain bound, like 80-90\\%.\n", " \n", "- At that point, the histogram is reset to zero $H(E)=0$ and a new\n", " random walk is started with finer modification factor (like $f_1\n", " =\\sqrt{f_0}$).\n", " \n", "- The simulation is again running so long that the histogram\n", " becomes flat. After that, the histogram is again reset and\n", " modification factor is reduced.\n", "\n", "- When the modification factor $f$ becomes sufficiently close to\n", " unity, the density of states does not change anymore. \n", " **If we are able to obtain a flat histogram using converged $g(E)$ in\n", " $P(E_1\\rightarrow E_2)$, the density of states $g(E)$ is obviously the true\n", " density of states of the model.**\n", " \n", "\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wang-Landau simulation" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from numpy import *\n", "from numpy import random\n", "from numba import jit\n", "\n", "@jit(nopython=True)\n", "def CEnergy(latt):\n", " \"Energy of configuration for the 2D Ising model\"\n", " N = shape(latt)[0]\n", " Ene = 0\n", " for i in range(len(latt)):\n", " for j in range(len(latt)):\n", " S = latt[i,j] # Spin, can be either +1 or -1\n", " WF = latt[(i+1)%N,j]+latt[i,(j+1)%N]+latt[(i-1)%N,j]+latt[i,(j-1)%N]\n", " Ene += -S * WF\n", " return Ene/2.\n", "\n", "def RandomL(N):\n", " \"Random lattice corresponding to infinite temperature\"\n", " return array(sign(2*random.random((N,N))-1),dtype=int) " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def PrepareEnergies(N):\n", " Energies = (array(4*arange(-int(N*N/2),int(N*N/2)+1),dtype=int)).tolist() # -2 N^2...2N^2 in steps of 4\n", " #Energies = range(-2*N*N,2*N*N,4)\n", " Energies.pop(1) # take out E[1]\n", " Energies.pop(-2) # take out E[-2]\n", " Energies = array(Energies) # make array of energies again\n", " Emin, Emax = Energies[0],Energies[-1]\n", " #index array to energies\n", " indE = -ones(Emax+1-Emin, dtype=int) # index table to get index to particular energy g(E)~g[indE[E]]\n", " for i,E in enumerate(Energies):\n", " indE[E-Emin]=i\n", " # Ising lattice at infinite T\n", " ## g(E) we know E = -2*N^2,2*N^2 in steps of 4...\n", " return (Energies, indE, Emin)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Emin= -200\n", "[-200 -192 -188 -184 -180 -176 -172 -168 -164 -160 -156 -152 -148 -144\n", " -140 -136 -132 -128 -124 -120 -116 -112 -108 -104 -100 -96 -92 -88\n", " -84 -80 -76 -72 -68 -64 -60 -56 -52 -48 -44 -40 -36 -32\n", " -28 -24 -20 -16 -12 -8 -4 0 4 8 12 16 20 24\n", " 28 32 36 40 44 48 52 56 60 64 68 72 76 80\n", " 84 88 92 96 100 104 108 112 116 120 124 128 132 136\n", " 140 144 148 152 156 160 164 168 172 176 180 184 188 192\n", " 200]\n", "[ 0 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 2 -1 -1 -1 3 -1 -1 -1 4 -1 -1 -1\n", " 5 -1 -1 -1 6 -1 -1 -1 7 -1 -1 -1 8 -1 -1 -1 9 -1 -1 -1 10 -1 -1 -1\n", " 11 -1 -1 -1 12 -1 -1 -1 13 -1 -1 -1 14 -1 -1 -1 15 -1 -1 -1 16 -1 -1 -1\n", " 17 -1 -1 -1 18 -1 -1 -1 19 -1 -1 -1 20 -1 -1 -1 21 -1 -1 -1 22 -1 -1 -1\n", " 23 -1 -1 -1 24 -1 -1 -1 25 -1 -1 -1 26 -1 -1 -1 27 -1 -1 -1 28 -1 -1 -1\n", " 29 -1 -1 -1 30 -1 -1 -1 31 -1 -1 -1 32 -1 -1 -1 33 -1 -1 -1 34 -1 -1 -1\n", " 35 -1 -1 -1 36 -1 -1 -1 37 -1 -1 -1 38 -1 -1 -1 39 -1 -1 -1 40 -1 -1 -1\n", " 41 -1 -1 -1 42 -1 -1 -1 43 -1 -1 -1 44 -1 -1 -1 45 -1 -1 -1 46 -1 -1 -1\n", " 47 -1 -1 -1 48 -1 -1 -1 49 -1 -1 -1 50 -1 -1 -1 51 -1 -1 -1 52 -1 -1 -1\n", " 53 -1 -1 -1 54 -1 -1 -1 55 -1 -1 -1 56 -1 -1 -1 57 -1 -1 -1 58 -1 -1 -1\n", " 59 -1 -1 -1 60 -1 -1 -1 61 -1 -1 -1 62 -1 -1 -1 63 -1 -1 -1 64 -1 -1 -1\n", " 65 -1 -1 -1 66 -1 -1 -1 67 -1 -1 -1 68 -1 -1 -1 69 -1 -1 -1 70 -1 -1 -1\n", " 71 -1 -1 -1 72 -1 -1 -1 73 -1 -1 -1 74 -1 -1 -1 75 -1 -1 -1 76 -1 -1 -1\n", " 77 -1 -1 -1 78 -1 -1 -1 79 -1 -1 -1 80 -1 -1 -1 81 -1 -1 -1 82 -1 -1 -1\n", " 83 -1 -1 -1 84 -1 -1 -1 85 -1 -1 -1 86 -1 -1 -1 87 -1 -1 -1 88 -1 -1 -1\n", " 89 -1 -1 -1 90 -1 -1 -1 91 -1 -1 -1 92 -1 -1 -1 93 -1 -1 -1 94 -1 -1 -1\n", " 95 -1 -1 -1 96 -1 -1 -1 97 -1 -1 -1 -1 -1 -1 -1 98]\n" ] } ], "source": [ "Ene,indE,Emin = PrepareEnergies(10)\n", "print('Emin=',Emin)\n", "print(Ene)\n", "print(indE)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "from numpy import random\n", "def WangLandau(Nitt, N, flatness):\n", " \"Wang Landau in Python\"\n", " (Energies, indE, Emin) = PrepareEnergies(N)\n", " latt = RandomL(N)\n", " return RunWangLandau(Nitt,Energies,latt,indE)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "@jit(nopython=True)\n", "def RunWangLandau(Nitt,Energies,latt,indE):\n", " N = len(latt)\n", " Ene = int(CEnergy(latt))\n", " #min,maximum energy\n", " Emin, Emax = Energies[0],Energies[-1]\n", " # Logarithm of the density of states\n", " lngE = zeros(len(Energies))\n", " # Histogram\n", " Hist = zeros(len(Energies))\n", " # modification factor\n", " lnf = 1.0 # f = exp(lnf)=e\n", " N2 = N*N\n", " for itt in range(Nitt):\n", " t = int(random.rand()*N2)\n", " (i, j) = (int(t/N), t%N)\n", " S = latt[i,j]\n", " WF = latt[(i+1)%N,j]+latt[i,(j+1)%N]+latt[(i-1)%N,j]+latt[i,(j-1)%N]\n", " Enew = Ene + int(2*S*WF) # the energy change if we flip the spin\n", " # P = g(E)/g(Enew) = exp(log(g(E))-log(g(Enew)))\n", " # P = exp(lngE(E)-lngE(Enew))\n", " lgnew = lngE[indE[Enew-Emin]] # log(g(Enew))\n", " lgold = lngE[indE[Ene-Emin]] # log(g(Eold))\n", " P = 1.0\n", " if lgold-lgnew < 0 : P=exp(lgold-lgnew) # P = g_old/g_new = exp(log(g_old)-log(g_new))\n", " if P > random.rand():\n", " # accept the step\n", " latt[i,j] = -S\n", " Ene = Enew\n", " Hist[indE[Ene-Emin]] += 1\n", " lngE[indE[Ene-Emin]] += lnf # g(E) -> g(E)*f hence log(g(E)) -> log(g(E))+log(f)\n", " \n", " if (itt+1) % 1000 == 0: # checking for flatness of the histogram\n", " aH = sum(Hist)/N2 # average\n", " mH = min(Hist)\n", " if mH > aH*flatness: # histogram is flat\n", " Hist[:] = 0\n", " lnf /= 2.\n", " print(itt, 'histogram is flat', mH, aH, 'f=', exp(lnf))\n", " return (lngE, Hist)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "68173999 histogram is flat 60078.0 66576.171875 f= 1.6487212707001282\n", "79598999 histogram is flat 10059.0 11157.2265625 f= 1.2840254166877414\n", "99577999 histogram is flat 17561.0 19510.7421875 f= 1.1331484530668263\n", "114655999 histogram is flat 13267.0 14724.609375 f= 1.0644944589178593\n", "130326999 histogram is flat 13789.0 15303.7109375 f= 1.0317434074991028\n", "143994999 histogram is flat 12019.0 13347.65625 f= 1.0157477085866857\n", "160480999 histogram is flat 14495.0 16099.609375 f= 1.007843097206448\n", "181231999 histogram is flat 18260.0 20264.6484375 f= 1.0039138893383475\n", "198768999 histogram is flat 15433.0 17125.9765625 f= 1.0019550335910028\n", "229439999 histogram is flat 26985.0 29952.1484375 f= 1.0009770394924165\n", "254823999 histogram is flat 22318.0 24789.0625 f= 1.0004884004786945\n", "295870999 histogram is flat 36091.0 40084.9609375 f= 1.0002441704297478\n", "349000999 histogram is flat 46699.0 51884.765625 f= 1.0001220777633837\n", "489810999 histogram is flat 123760.0 137509.765625 f= 1.0000610370189331\n", "712426999 histogram is flat 195665.0 217398.4375 f= 1.000030518043791\n", "910861999 histogram is flat 174411.0 193784.1796875 f= 1.0000152589054785\n", "Done\n" ] } ], "source": [ "from numpy import random\n", "\n", "flatness = 0.9\n", "N = 32\n", "Nitt = int(10e8)\n", "\n", "#N=10\n", "#Nitt = int(10e6)\n", "\n", "(lngE, Hist) = WangLandau(Nitt, N, flatness)\n", "print('Done')" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(80508.46934509277, 80507.8903503418)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lngE[0],lngE[-1]" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "# Proper normalization of the Density of states\n", "# g *= 4/(g[0]+g[-1]) = 4/(g[-1]*(1+g[0]/g[-1]))\n", "# g = g * 4/(g[-1]*(1+g[0]/g[-1]))\n", "# log(g) = log(g) + log(4)-log(g[-1]) - log(1+g[0]/g[-1])\n", "# log(g) += log(4)-log(g[-1])-log(1+exp(log(g[0])-log(g[-1])))\n", "# \n", "if lngE[-1]>lngE[0]:\n", " lgC = log(4)-lngE[-1]-log(1+exp(lngE[0]-lngE[-1]))\n", "else:\n", " lgC = log(4)-lngE[0]-log(1+exp(lngE[-1]-lngE[0]))\n", "\n", "lngE += lgC" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3.999999999984037, 2.5633442630349297, 1.4366557369491075)" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(lngE[0])+exp(lngE[-1]), exp(lngE[0]), exp(lngE[-1]) # g[0]+g[-1]" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAta0lEQVR4nO3deVxVdf7H8dfnsioiiCAioKKCCiou5NZqZmrLqJVlq02LNdU01fymdGpqamqq6dfUVGOTbWObji2mWdqYWtmUGu4bCooKioAiCCL79/cHx/mRS6yXc5fP8/Hgce/93nO47y/qh+P5nvP9ijEGpZRSnsVhdwCllFItT4u7Ukp5IC3uSinlgbS4K6WUB9LirpRSHsjX7gAA4eHhpnv37nbHUEopt7J27dpDxpiI073nEsW9e/fupKam2h1DKaXciojsPdN7elpGKaU8kBZ3pZTyQFrclVLKA2lxV0opD6TFXSmlPJAWd6WU8kBa3JVSygPVe527iPQG/lWnqQfwKPCO1d4d2ANcbYw5Yu0zA7gVqAbuNcZ82aKplWphFVU1ZB8pZV9BKUdKKygpq+JoWRUVVTX4+Qh+Pg7aBvgS0S6AiOAAokPbENk+ABGxO7pSp1VvcTfG7AAGAoiID7AfmA9MB5YZY54RkenW64dEJBGYAiQBXYCvRCTBGFPtnC4o1ThV1TVs3l/Eqt0FbNlfxJYDRewrKKWxSxsEB/oS36kdybGhDO0eRkr3MCKCA5wTWqlGauwdqqOBXcaYvSIyAbjAap8NfA08BEwA5hpjyoFMEckAhgI/tEhipZqgqLSSL7ceZOn2XFbtOkxxeRUAXcPa0i+6PRMHRtM1rC2xYW2JCA6gXYAvwYG++Ps4qKoxVNXUUFJWRX5JOXnF5WQVlLIzt5gdB4v5YPU+3v7PHgB6RwYztl9nxvfrTJ/OwXpkr2zT2OI+BZhjPY80xuQAGGNyRKST1R4NrKqzT7bV9hMiMg2YBtC1a9dGxlCqfpXVNfx7ay7z12fzzc58KqsNMR3acFlyFGf3CmdEj450bFf/kba/Q/DHQVt/Xzq1DyTppPcrqmrYcqCINZkFrEjL45Xl6by0LJ1uHdty5eAYrk6JpXNIoHM6qdQZSEOX2RMRf+AAkGSMyRWRQmNMaJ33jxhjOojI34EfjDHvWe1vAl8YYz4+0/dOSUkxOreMaim5R8t4f/U+5qzZR35xOZ3bB3J5chSXJ3ehf3SI04+mD5WUs3RbLos2HeA/GYfxcQgX9unEzSO7M7JnRz2aVy1GRNYaY1JO915jjtzHA+uMMbnW61wRibKO2qOAPKs9G4its18Mtb8UlHKqrIJS/vHNLj5MzaaypoYLEiK4aUR3zk+IwOFovYIa3i6Aa4d25dqhXdl7+Bhz1mTxYWoWS7flkhwTwq8u6MnFiZ1bNZPyPo05cp8LfGmMedt6/RxwuM6Aapgx5kERSQI+oPY8exdgGRD/cwOqeuSumiPvaBkvfLWTD1OzEYGrhsRy5/k96NYxyO5o/1VWWc3H67J57Zvd7CsopU/nYGZc0pfz4sP1SF412c8duTeouItIWyAL6GGMKbLaOgLzgK7APmCyMabAeu9h4BagCrjPGLP4576/FnfVFKUVVcz6djezvt1NZXUN1w3typ0X9CQqpI3d0c6oqrqGRZtyeH7pDrIKjnNOr3Cmj+9Dv+gQu6MpN9Ts4u5sWtxVYxhjWLLlIH/8bCu5R8u5pH9nHhrXx6WO1OtTXlXN+6v28fLydAqPV3Ld0K48OLYPIW397I6m3IgWd+UxDhQe59EFW/lqey6JUe15YkISKd3D7I7VZEfLKnlxaTr//D6TsCB/Hrk0kQkDu+ipGtUgWtyV2zPG8N6qvTyzOI0aA/ePieeWs+Pw9fGMGTS27C/i4U+3sDGrkDGJkTx9RX/CG3CZpvJuWtyVW8svLufBjzayYkc+58aH8+dJ/YkNa2t3rBZXXWN4+z+Z/OXLHQQH+PL0Ff25OKmz3bGUC/u54u4Zhz3KYy3bnsu4F7/l+12HeWJCEu/cMtQjCzuAj0O47dweLPr1OXQOCWTau2uZ/vEmyip15g7VeC6xQLZSJ6uqruEvX+5g1re76RvVnrlTBhIfGWx3rFaREBnM/LvO5sWvdjLz611szC5i5vWDiQt3nwFjZT89clcu51BJOTe+uYZZ3+7mphHd+PTukV5T2E/w93Xw4Lg+vP3Ls8gpOs4vXv6OJVty7I6l3IgWd+VS1u87wuUvf8e6fUd4fnIyT0zoR4Cvj92xbDOqdyc+v/dcenRqx53vrePpxduprrF/nEy5Pi3uymV8un4/17y2Ch+H8PGvRnLlkBi7I7mE6NA2fHjHCK4f1pXXvtnNHe+mUmLNaqnUmWhxV7YzxvDiVzu5718bGNQ1lM/uOUfv2DyJv6+Dpyb154kJSazYkc9Vr35P9pFSu2MpF6bFXdmqvKqaB+Zt5MWv0rlicDTv3jqMDkH+dsdyWTeN6M4/f3kWBwqPM/Hv/2HL/iK7IykXpcVd2eZoWSU3vrmG+ev389sxCTw/ORl/X/0rWZ9z4yP45K6zCfD1YcqsVXyfccjuSMoF6b8kZYvDJeVcO2sV6/Ye4cVrBvLr0fF6y30j9OrUjo9/NZLo0Dbc/PaPfL5Jr6RRP6XFXbW6A4XHmfzaD2TklfD6TSlMHHTKQl2qATqHBDLvjhEkx4Zwz5x1vLtqr92RlAvR4q5aVeahY0z+xw/kHy3n3VuHMapPp/p3UmcU0taPd28dxug+kfzh0y28sXK33ZGUi9DirlpN2sGjTP7HDxyvrGbOtOEMjXPf2RxdSaCfD6/eMJhL+0fx5OfbefXrXXZHUi5Apx9QrWLHwWKue301/j4O3rttOL06tbM7kkfx83HwtykD8XEIzy5Jo7K6hntHx9sdS9lIi7tyuvTcYq57fRV+PsKcacN1jhQn8fVx8MI1A/F1CH9dupPK6hoeGJOgA9VeSou7cqqMvBKufX01Dofwwe1a2J3NxyE8NzkZXx/h5eUZ+Doc/OYiPYL3RlrcldNkHjrGta+vAmDO7cPpGaGnYlqDj0N45ooBVNfAC1/tJCjAh9vO7WF3LNXKGjSgKiKhIvKRiKSJyHYRGSEiYSKyVETSrccOdbafISIZIrJDRMY6L75yVTlFx7nhjdVU1xjm3D5Mz7G3ModDePbK/ozv15knP9/O3DX77I6kWllDr5b5G7DEGNMHSAa2A9OBZcaYeGCZ9RoRSQSmAEnAOGCmiHjvtH5eqOBYBTe8sZqi45W8c8tQr5uu11X4+jj425RBnJ8QwYz5m1m48YDdkVQrqre4i0h74DzgTQBjTIUxphCYAMy2NpsNTLSeTwDmGmPKjTGZQAYwtGVjK1dVXFbJ1LfWkH3kOG9OTdEJwGzm7+vgHzcM4azuYTzwrw0s255rdyTVShpy5N4DyAfeFpH1IvKGiAQBkcaYHADr8cTdKNFAVp39s622nxCRaSKSKiKp+fn5zeqEcg1lldXc/k4q23OO8uoNgxnWo6PdkRTQxt+HN6emkNilPXd/sI61e4/YHUm1goYUd19gMPCqMWYQcAzrFMwZnO66q1NWFzDGzDLGpBhjUiIiIhoUVrmumhrD/3y4kVW7C3j+6mQu7BNpdyRVR3CgH2/dfBad2wdy2+wf2ZVfYnck5WQNKe7ZQLYxZrX1+iNqi32uiEQBWI95dbaPrbN/DKAn+zzcs1+msWhTDjPG92HCQJ0rxhWFtwtg9i1DcYgw9a015BWX2R1JOVG9xd0YcxDIEpHeVtNoYBuwEJhqtU0FFljPFwJTRCRAROKAeGBNi6ZWLuXdH/bw2je7uXF4N6adp5fcubJuHYN46+azOFxSwS/f/lFXdPJgDb1a5tfA+yKyCRgI/Bl4BhgjIunAGOs1xpitwDxqfwEsAe42xlS3cG7lIr7alstjC7dyUd9OPHZ5ot4N6QaSY0OZef1g0g4W86v31lJZXWN3JOUEYoz9i+2mpKSY1NRUu2OoRtqUXcg1r60iPrIdc6cNp62/3hPnTub9mMWDH2/i6pQYnr1ygP5idkMistYYk3K69/Rfo2qS3KNl3DY7lbAgf96cepYWdjd09VmxZB0p5eXlGSREButdrB5G/0WqRiurrGbaO6kcK6/i47tGEhEcYHck1UT3X5RARl4JT32xnbjwIEb31aucPIXO564axRjD9I83sTG7iBeuGUifzu3tjqSaweEQnr86maQu7bl3znp2HCy2O5JqIVrcVaO89u1uPt1wgP+5OIGLkzrbHUe1gLb+vrx+UwpBAb7cOvtHDpeU2x1JtQAt7qrBlqfl8uySNC4bEMXdo3rZHUe1oKiQNrx+Uwr5xeXc+d5aKqr0Chp3p8VdNUhGXjH3ztlAYlR7nrsqWa+s8EDJsaE8NzmZH/cc4anPt9kdRzWTDqiqehWWVnDr7FQC/Xx4/aYU2vjrJJ+e6hfJXdicXcjrKzMZEBPKlUNi7I6kmkiP3NXPqqqu4Z4P1pNTWMZrNw6hS2gbuyMpJ3toXB9G9OjI7+dvZsv+IrvjqCbS4q5+1vNLd/JdxiGenNSPId061L+Dcnu+Pg5euW4QHYP8uePdtRQcq7A7kmoCLe7qjJZtz+XVr3dx7dBYrk6JrX8H5TE6tgvgHzcOIb+knHvnrKe6xv472VXjaHFXp5VVUMoD8zaSGNWexy5PsjuOssGAmFCenNCP7zIO8dyXO+yOoxpJi7s6RXlVNfd8sI6aGsOrNwwm0E8HUL3V1WfFcv2wrvzjm10s2XLQ7jiqEbS4q1M8/UUaG7OLeG7yALp1DLI7jrLZo5cnkhwTwu8+2si+w6V2x1ENpMVd/cTnm3L45/d7uPWcOMb1i7I7jnIBAb4+vHLdYAS4Z846yqt0Bm93oMVd/dfu/BIe+ngTg7uGMn18H7vjKBcSG9aW5yYnsym7iKe/SLM7jmoALe4KgOMV1dz1/jr8fIRXrhuMn4/+1VA/NTapM7ecHcc/v9/D4s05dsdR9dB/wQqARxdsYUduMS9cM1BvVFJnNH18H5JjQ3nwo01kFej5d1emxV3xYWoWH67N5p5Rvbigdye74ygX5u/r4JVrB4HAb+aup0qX6HNZWty93K78Eh5dsJXhPcK476IEu+MoNxAb1panJvVn3b5CXlqeYXccdQYNKu4iskdENovIBhFJtdrCRGSpiKRbjx3qbD9DRDJEZIeIjHVWeNU8FVU13Dd3AwF+Dl68ZhA+Dp3pUTXML5K7cOXgGF5Zns6azAK746jTaMyR+yhjzMA6i7FOB5YZY+KBZdZrRCQRmAIkAeOAmSKid8G4oOeX7mDz/iKevXIAnUMC7Y6j3MzjE5KIDWvL/f/aQNHxSrvjqJM057TMBGC29Xw2MLFO+1xjTLkxJhPIAIY243OUE3yXfojXvtnNdcO6MlZXVFJN0C7Al79NGUTu0TJ+P38zxuj8M66kocXdAP8WkbUiMs1qizTG5ABYjydG4qKBrDr7ZlttPyEi00QkVURS8/Pzm5ZeNUnBsQoemLeBnhFB/OHSRLvjKDc2MDaU+8ck8PmmHD5am213HFVHQ4v72caYwcB44G4ROe9ntj3didtTfqUbY2YZY1KMMSkRERENjKGayxjDgx9torC0kpeuHaQLb6hmu/P8ngzvEcZjC7ey59Axu+MoS4OKuzHmgPWYB8yn9jRLrohEAViPedbm2UDd+WFjgAMtFVg1z3ur9/HV9lweGt+HpC4hdsdRHsDHIbxwzUD8fBzcO3e9rr/qIuot7iISJCLBJ54DFwNbgIXAVGuzqcAC6/lCYIqIBIhIHBAPrGnp4KrxMg8d48+fb+fc+HB+ObK73XGUB4kKacOzV/ZnU3YRL3y10+44ioatoRoJzLcWRPYFPjDGLBGRH4F5InIrsA+YDGCM2Soi84BtQBVwtzFGZxqyWXWN4X8+3Iifj/DcVck49LJH1cLG9Yvi2qGx/OObXYzq3YmhcWF2R/Jq4goj3CkpKSY1NdXuGB7ttW928fTiNF64JplJg3TRY+Ucx8qruOSlldQYw+LfnEe7gIYcP6qmEpG1dS5P/wm9Q9UL7Mwt5vl/72RcUmcmDjzlwiWlWkxQgC/PT04m+8hxnvp8m91xvJoWdw9XWV3DA/M2EBzoy5OT+mGdXlPKaVK6h3HHeT2ZsyaL5Wm5dsfxWlrcPdzfV2SwZf9RnprUj/B2AXbHUV7i/jHx9OkczEMfb+bIsQq743glLe4ebHN2Ea8sz2DSoGhdVUm1qgBfH/569UAKSyt45NMteveqDbS4e6iyymoemLeBju38+ePlSXbHUV4osUt77rsogc8357Bwo97q0tq0uHuoF5buJD2vhGevHEBIWz+74ygvdcd5PRjcNZQ/fLqFg0VldsfxKlrcPVDqngJmrdzNtUO76uIbyla+Pg6ev3ogFdU1OrlYK9Pi7mFKK6r47YcbiQ5tw8OX9rU7jlLEhQfxu7F9WJ6Wx/z1++2O4zW0uHuYZxansfdwKf87OVlvIFEu4+aR3RnSrQOPf7aNvKN6eqY1aHH3IN+lH+KdH/Zyy9lxDO/R0e44Sv2Xj0P4y1UDOF5ZrVfPtBIt7h7iaFklD360kR4RQTw4rrfdcZQ6Rc+IdjwwJoF/b8tl0aYcu+N4PC3uHuJPn23j4NEynp+cTKCfztGuXNNt58SRHBPCYwu3crik3O44Hk2LuwdYmZ7Ph2uzueP8ngzq2qH+HZSyia+Pg+cmJ1NSVsWjC7faHcejaXF3c6UVVcz4ZDM9woP4zeh4u+MoVa+EyGDuHd2LzzflsGSLnp5xFi3ubu6v/95J9pHjPH1Ffz0do9zGHef3JKlLex75dKvOPeMkWtzd2MasQt76TybXDevKML06RrkRPx8Hz12VTGFpBU8s0qmBnUGLu5uqqKrhoY83EREcwPTxfeyOo1SjJXZpz12jejF//X6WbdepgVuaFnc3NevbXaQdLObJif1pH6hzxyj3dM+oXvSODOaRT7dQUl5ldxyPosXdDWXklfDSsgwuHRDFmMRIu+Mo1WT+vg7+fEV/Dh4t43+/3GF3HI/S4OIuIj4isl5EFlmvw0RkqYikW48d6mw7Q0QyRGSHiIx1RnBvVVNjmPHJJtr4++hUvsojDOnWgRuHd2P2D3vYkFVodxyP0Zgj998A2+u8ng4sM8bEA8us14hIIjAFSALGATNFRC/jaCHvr9nHj3uO8MilfYkI1pWVlGf43djeRAYHMv3jTVRW19gdxyM0qLiLSAxwKfBGneYJwGzr+WxgYp32ucaYcmNMJpABDG2RtF4up+g4zy5O49z4cK4aEmN3HKVaTHCgH49PSCLtYDFvrMy0O45HaOiR+4vAg0DdX6mRxpgcAOvxxMTh0UBWne2yrbafEJFpIpIqIqn5+fmNze11jDE8Mn8L1TWGP0/qrwtdK48zNqkzY5MiefGrnew9fMzuOG6v3uIuIpcBecaYtQ38nqerOqdMAWeMmWWMSTHGpERERDTwW3uvRZtyWJaWx28vTiA2rK3dcZRyisd/0Q8/HwcPz9eZI5urIUfuZwO/EJE9wFzgQhF5D8gVkSgA6zHP2j4biK2zfwygCyg2Q2FpBX9cuJXkmBB+eXac3XGUcprOIYE8NK4332Uc4pN1urBHc9Rb3I0xM4wxMcaY7tQOlC43xtwALASmWptNBRZYzxcCU0QkQETigHhgTYsn9yJPf5FG4fFKnrlyAD4OPR2jPNv1w7oxpFsHnvx8m84c2QzNuc79GWCMiKQDY6zXGGO2AvOAbcAS4G5jTHVzg3qrNZkF/Cs1i9vOjaNvVHu74yjldA6H8PQV/Skpr+Kpz7fXv4M6rUYVd2PM18aYy6znh40xo40x8dZjQZ3tnjLG9DTG9DbGLG7p0N6ioqp2UeHo0DY646PyKgmRwdx5fk8+Wb+flel6wUVT6B2qLuz1lbvJyCvhyYn9aOuv66Eq73L3qF7EhQfx8PwtHK/Q//w3lhZ3F7Xn0DFeWpbOpf2jGNWnU/07KOVhAv18+POk/uwrKOVvy9LtjuN2tLi7IGMMf1iwBX8fB49enmh3HKVsM6JnR64aEsMbK3ezM7fY7jhuRYu7C1q48QAr0w/xu3G9iWwfaHccpWz1+0v60i7Ql4fnb6amRq99bygt7i6mqLSSPy3aRnJMCNcP62Z3HKVsFxbkz+/H9+XHPUf4aF223XHchhZ3F/Psl2kcKa3kz1f012valbJcNSSGs7p34OkvtlOgy/I1iBZ3F7J2bwEfrN7HL0d2J6lLiN1xlHIZDofw5MT+FJdV8cxivfa9IbS4u4iq6hoenr+FqJBA7h+TYHccpVxO787B3HZuD+alZrMms6D+HbycFncX8c4Pe0k7WMyjlyUSFKDXtCt1OveO7kV0aBse+XQzFVU67/vP0eLuAvKOlvHXpTs5LyGCcf062x1HKZfV1t+XJyYksTO3hDe/03nff44Wdxfw1Bfbqaiq4fFfJOk87UrVY3TfSMYmRfK3ZTvJKii1O47L0uJus+93HWLBhgPceX4P4sKD7I6jlFt47PIkHCI8/tlWu6O4LC3uNqqsruHRBVuJDWvDXaN62R1HKbfRJbQN946O56vteSxPy7U7jkvS4m6jt77LJCOvhD9enkSgn64hrlRj3HJ2HD0jgvjjwm2UVerEYifT4m6TnKLj/G1ZOhf1jWR030i74yjldvx9HTwxoR/7Ckp57ZvddsdxOVrcbfKnRduoMYbHdGIwpZrs7F7hXDogiplfZ+jg6km0uNvg2535fLH5IPeM6qWLXSvVTI9c2hcfh/D4Z9vsjuJStLi3svKqah5buJUe4UHcfl4Pu+Mo5faiQk4Mrubq4GodWtxb2axvdpN56BiPT0giwFcHUZVqCTq4eqp6i7uIBIrIGhHZKCJbReRxqz1MRJaKSLr12KHOPjNEJENEdojIWGd2wJ1kFZTyyooMLu0fxbnxEXbHUcpj6ODqqRpy5F4OXGiMSQYGAuNEZDgwHVhmjIkHllmvEZFEYAqQBIwDZoqIHqICj3+2FR+H8Mhlfe2OopTH0cHVn6q3uJtaJdZLP+vLABOA2Vb7bGCi9XwCMNcYU26MyQQygKEtGdodfbUtl6+253HfRfFEhbSxO45SHkkHV/9fg865i4iPiGwA8oClxpjVQKQxJgfAejyxinM0kFVn92yr7eTvOU1EUkUkNT8/vxldcH3HK6r542dbSYhsxy/PjrM7jlIeSwdX/1+DirsxptoYMxCIAYaKSL+f2fx0M1+dsvChMWaWMSbFGJMSEeHZ559nfp1B9pHjPDGhH34+OoatlDPp4GqtRlUaY0wh8DW159JzRSQKwHrMszbLBmLr7BYDHGhuUHeVeegYr32zm0mDohneo6PdcZTyeHUHV2d9672Dqw25WiZCREKt522Ai4A0YCEw1dpsKrDAer4QmCIiASISB8QDa1o4t1swxvDogi0E+DqYcUkfu+Mo5TVODK7+fYX3Dq425Mg9ClghIpuAH6k9574IeAYYIyLpwBjrNcaYrcA8YBuwBLjbGOOV/zdavOUgK9MP8cDFCXQKDrQ7jlJexdsHV+tdz80YswkYdJr2w8DoM+zzFPBUs9O5sWPlVTzx2TYSo9pz4/BudsdRyuucGFx9ZnEay9NyubCPd03Qp6N7TvLSsnQOHi3jTxP74auDqErZwpsHV7XqOEF6bjFvfpfJ1SkxDOnWof4dlFJO4c2Dq1rcW5gxhj8s2EJQgC8PjdNBVKXs5q2Dq1rcW9jCjQdYtbuAB8f1pmO7ALvjKKX4/8HVJxZ5z+CqFvcWdLSskic/305yTAhTzupqdxyllOXE4OrSbbmsSMurfwcPoMW9Bf3tq3QOlZTzp4n98HGc7kZdpZRd/ju4+tlWrxhc1eLeQjLySpj9/R6mnNWVATGhdsdRSp3kxODq3sPeMbiqxb0FGGP406JttPH34X8uTrA7jlLqDLxpcFWLewtYnpbHNzvzue+iBB1EVcrFecvgqhb3ZqqoquFPi7bRMyKIm0bonahKuTpvGVzV4t5Mb/8nkz2HS3n08iSdzlcpN+ENg6tajZohr7iMl5dncFHfTpyf4Nlz0ivlSbxhcFWLezM8t2QH5VXVPHxpot1RlFKNVHfN1Zyi43bHaXFa3JtoY1YhH67N5pZz4ogLD7I7jlKqCWaM70ONqT1Q8zRa3JvAGMPjn20lvF0A94zqZXccpVQTxXRoy23nxPHJ+v1syCq0O06L0uLeBAs2HGDdvkIeGteb4EA/u+MopZrhrlG9CG8XwJ8WbcOYU5Z7dlta3BvpWHkVTy+unT/mysExdsdRSjVTuwBffjc2gbV7j7BoU47dcVqMFvdGmvl1BrlHy3n08iQcOn+MUh7hqiGxJEa155nFaR5zaaQW90bIKijl9ZWZTBoUrYtwKOVBfBzCI5f1ZX/hcd78LtPuOC2i3uIuIrEiskJEtovIVhH5jdUeJiJLRSTdeuxQZ58ZIpIhIjtEZKwzO9Canl2ShkPgwXG97Y6ilGphI3uGc3FiJDNXZJBXXGZ3nGZryJF7FfBbY0xfYDhwt4gkAtOBZcaYeGCZ9RrrvSlAEjAOmCkiPs4I35pOnI+bdl5PokLa2B1HKeUEv7+kLxXVNTz/5U67ozRbvcXdGJNjjFlnPS8GtgPRwARgtrXZbGCi9XwCMNcYU26MyQQygKEtnLtVnZj1sVNwAHee38PuOEopJ+keHsTNI7szb20WW/YX2R2nWRp1zl1EugODgNVApDEmB2p/AQCdrM2igaw6u2VbbSd/r2kikioiqfn5+U2I3noWbjzAhqxCfje2N239fe2Oo5RyonsujKdDW3+3vzSywcVdRNoBHwP3GWOO/tymp2k75SdkjJlljEkxxqRERLjuvCxlldX8ZckOkrq010sflfICIW38uP+ieFZnFvDl1ly74zRZg4q7iPhRW9jfN8Z8YjXnikiU9X4UcGLuzGwgts7uMcCBlonb+t78LpP9hcd55NJEvfRRKS9x7dCuxHdqx9OLt1Ne5Z6XRjbkahkB3gS2G2P+WuethcBU6/lUYEGd9ikiEiAicUA8sKblIreevOIyZq7I4OLESEb07Gh3HKVUK/H1cfDIZYnsPVzK7O/32B2nSRpy5H42cCNwoYhssL4uAZ4BxohIOjDGeo0xZiswD9gGLAHuNsa45a++F5bupKK6hhmX9LU7ilKqlZ2fEMEFvSN4eXkGR45V2B2n0Rpytcx3xhgxxgwwxgy0vr4wxhw2xow2xsRbjwV19nnKGNPTGNPbGLPYuV1wjh0Hi/nXj1ncOLy7zvqolJeaMb4vx8qreGVFht1RGk3vUD2DZ5ekERTgy68v1FkflfJWvTsHc9WQGN79Ya/bLaitxf00fth1mOVpedw9qhcdgvztjqOUstH9YxJwOOB//+1ec75rcT+JMYZnFm8nKiSQm0d2tzuOUspmUSFtuOXsOBZsOMDmbPe5sUmL+0k+35zDxuwiHhiTQKCf28+aoJRqAXde0JMObf14evF2t7mxSYt7HRVVNfxlyQ76dA7mCr1hSSllaR/ox72j4/l+12G+2enad9SfoMW9jg9W72VfQSkPje+Dj96wpJSq4/ph3ega1pZnFqdRXeP6R+9a3C3FZZW8tDyDkT07ckGC606HoJSyh7+vg9+N7U3awWLmr99vd5x6aXG3vLEyk4JjFUwf34fam3KVUuqnLu0fRXJMCM//e4fLr9ikxR0oOFbBGyt3M75fZwbEhNodRynlohwOYfr4vuQUlfHeqr12x/lZWtyBV7/O4HhlNQ+MSbA7ilLKxY3o2ZFz48OZ+fUuSsqr7I5zRl5f3HOKjjP7h71MGhRDfGSw3XGUUm7gtxf3puBYBW+58HqrXl/cX16egTGG+y6KtzuKUspNDIwNZUxiJK9/u5vCUtecVMyri/ueQ8eY92MW1w3tSmxYW7vjKKXcyG8vTqCkoorXvt1td5TT8uri/uJXO/H1Ee7WycGUUo3Up3N7Lh/QhX/+Zw95xWV2xzmF1xb39NxiFmw8wM0j4+gUHGh3HKWUG7p/TAIV1TXMXLHL7iin8Nri/vLyDNr4+TDtvB52R1FKuam48CCuGhzDB6v3sb/wuN1xfsIri3tGXgmfbTrATSO6E6ZT+iqlmuFe62KMV5an25zkp7yyuP99RQaBvj7cfm6c3VGUUm4uOrQN15wVy0drs8k+4joLenhdcd+dX8KCDfu5cUQ3OrYLsDuOUsoD3HlBTwD+8Y3rnHuvt7iLyFsikiciW+q0hYnIUhFJtx471HlvhohkiMgOERnrrOBN9cqKDPx9Hdx+rp5rV0q1jOjQNlw1JJZ5P2ZzsMg1rpxpyJH7P4FxJ7VNB5YZY+KBZdZrRCQRmAIkWfvMFBGXWfFiz6FjLNhwgOuHdSMiWI/alVIt564LelJjjMscvddb3I0x3wIFJzVPAGZbz2cDE+u0zzXGlBtjMoEMYGjLRG2+v6/IwNch3HG+HrUrpVpWbFhbJg2KZs6afS5x3XtTz7lHGmNyAKzHTlZ7NJBVZ7tsq+0UIjJNRFJFJDU/3/krm+w7XMon6/dz3bCuel27Usop7h7Vi8rqGl53gbtWW3pA9XQToZ92yRJjzCxjTIoxJiUiwvmLY8z8OgMfh3Dn+T2d/llKKe/UPTyIiQOjeW/VPg6XlNuapanFPVdEogCsxzyrPRuIrbNdDHCg6fFaRlZBKR+tzebas2KJbK9H7Uop57lrVC/Kqqp5w+YZI5ta3BcCU63nU4EFddqniEiAiMQB8cCa5kVsvplf78Ih8t/LlZRSyll6dWrHZQO68M73ezhyzL4ZIxtyKeQc4Aegt4hki8itwDPAGBFJB8ZYrzHGbAXmAduAJcDdxhhb16LKPVrGR2uzmJwSQ1RIGzujKKW8xK8v7MWximr++f0e2zL41reBMebaM7w1+gzbPwU81ZxQLemf3++husboHDJKqVaTEBnM6D6deHfVXn51QU8C/Vr/inCPvkO1pLyK91ftZVy/znTrGGR3HKWUF7n9vB4UHKvgk3X7bfl8jy7u837M4mhZld6NqpRqdcPiwhgQE8IbK3dTU3PaiwadymOLe1V1DW9+l8lZ3TswqGuH+ndQSqkWJCLcfm4Pdh86xrK0vPp3aGEeW9wXbznI/sLjetSulLLN+H6diQ5tw+srW/+mJo8s7sYYZn27m7jwIC7qG2l3HKWUl/L1cXDLOXGsySxgQ1Zhq362Rxb3VbsL2Ly/iNvOjcPhON1Ns0op1TquOSuW4EDfVj9698ji/vrK3YQF+XPl4Bi7oyilvFy7AF+uH9aNxZtzyCpovcU8PK64Z+SVsDwtjxuHd7Pl2lKllDrZzSO74xDh7f/sabXP9Lji/s4Pe/D3cXDjiG52R1FKKQA6hwQyvn8UH63N4nhF69y071HF/WhZJR+tzeay5CjCdQk9pZQLuWFYV46WVfHZptaZS9Gjivun6/dTWlHNzSO72x1FKaV+YmhcGAmR7Xh/1d5W+TyPKu7zUrNI6tKeATGhdkdRSqmfEBGuH9aNjdlFbMoudPrneUxx35RdyJb9R7k6Jbb+jZVSygaTBkfTxs+H91rh6N1jivv7q/bRLsCXSYNPu6qfUkrZrn2gHxMHdWHhxgMUlVY69bM8orgfLavk8805XJwYSftAP7vjKKXUGV0/rBtllTV8sj7bqZ/jEcX96x35lJRXMWVoV7ujKKXUz+oXHUKfzsEs3nzQqZ/jEcX9g9V76RrWlsFdQ+2OopRS9bo8uQtr9hSQnlvstM9w++KeU3SctXuPMCYxEl8ft++OUsoLXD6gCwBLt+c67TPcvhou3ZZLZbVhcorOI6OUcg+xYW0Y3DWUT9c7b5UmpxV3ERknIjtEJENEpjvrcz5et5/wdv70jgx21kcopVSLEhHOiY9gZ24JGXklTvkMpxR3EfEB/g6MBxKBa0UksaU/p6yymk3ZhUwcGI2ITu2rlHIfkwbVXra9wkmrNDnryH0okGGM2W2MqQDmAhNa+kMyDx3DGEiODW3pb62UUk4VFx5EWJA/uw+50ZE7EA1k1XmdbbX9l4hME5FUEUnNz89v0oe0b+PHvRf2IlmnG1BKuaEJA7vQq5NzTin7OuW7wunOkfxk+W9jzCxgFkBKSkqTlgaPDm3DAxf3bsquSillu8cuT3La93bWkXs2UHeSlxigdea5VEop5bTi/iMQLyJxIuIPTAEWOumzlFJKncQpp2WMMVUicg/wJeADvGWM2eqMz1JKKXUqZ51zxxjzBfCFs76/UkqpM3P7O1SVUkqdSou7Ukp5IC3uSinlgbS4K6WUBxJjmnT/UMuGEMkHmrOoYDhwqIXiuDpv6it4V3+9qa/gXf11Vl+7GWMiTveGSxT35hKRVGNMit05WoM39RW8q7/e1Ffwrv7a0Vc9LaOUUh5Ii7tSSnkgTynus+wO0Iq8qa/gXf31pr6Cd/W31fvqEefclVJK/ZSnHLkrpZSqQ4u7Ukp5ILcu7q21CHdrEZFYEVkhIttFZKuI/MZqDxORpSKSbj12qLPPDKv/O0RkrH3pm05EfERkvYgssl57ZH9FJFREPhKRNOvPeISn9hVARO63/h5vEZE5IhLoSf0VkbdEJE9EttRpa3T/RGSIiGy23ntJWmpBaGOMW35RO5XwLqAH4A9sBBLtztXMPkUBg63nwcBOahcY/wsw3WqfDjxrPU+0+h0AxFk/Dx+7+9GEfj8AfAAssl57ZH+B2cBt1nN/INSD+xoNZAJtrNfzgJs9qb/AecBgYEudtkb3D1gDjKB2BbvFwPiWyOfOR+6tsgh3azLG5Bhj1lnPi4Ht1P4jmUBtYcB6nGg9nwDMNcaUG2MygQxqfy5uQ0RigEuBN+o0e1x/RaQ9tcXgTQBjTIUxphAP7GsdvkAbEfEF2lK7GpvH9NcY8y1QcFJzo/onIlFAe2PMD6a20r9TZ59mcefiXu8i3O5MRLoDg4DVQKQxJgdqfwEAnazNPOFn8CLwIFBTp80T+9sDyAfetk5BvSEiQXhmXzHG7Af+F9gH5ABFxph/46H9raOx/Yu2np/c3mzuXNzrXYTbXYlIO+Bj4D5jzNGf2/Q0bW7zMxCRy4A8Y8zahu5ymjZ36a8vtf+Ff9UYMwg4Ru1/28/EnfuKda55ArWnILoAQSJyw8/tcpo2t+lvA5ypf07rtzsXd49chFtE/Kgt7O8bYz6xmnOt/75hPeZZ7e7+Mzgb+IWI7KH2tNqFIvIentnfbCDbGLPaev0RtcXeE/sKcBGQaYzJN8ZUAp8AI/Hc/p7Q2P5lW89Pbm82dy7uHrcItzVK/iaw3Rjz1zpvLQSmWs+nAgvqtE8RkQARiQPiqR2ccQvGmBnGmBhjTHdq//yWG2NuwAP7a4w5CGSJSG+raTSwDQ/sq2UfMFxE2lp/r0dTO4bkqf09oVH9s07dFIvIcOvndFOdfZrH7hHnZo5WX0LtFSW7gIftztMC/TmH2v+SbQI2WF+XAB2BZUC69RhWZ5+Hrf7voIVG2W3q+wX8/9UyHtlfYCCQav35fgp08NS+WvkfB9KALcC71F4p4jH9BeZQO55QSe0R+K1N6R+QYv2MdgGvYM0c0NwvnX5AKaU8kDufllFKKXUGWtyVUsoDaXFXSikPpMVdKaU8kBZ3pZTyQFrclVLKA2lxV0opD/R/vL67LL9+hQoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from pylab import *\n", "%matplotlib inline\n", "\n", "plot(lngE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The thermodynamics can be calculated in the following way:\n", "\n", "\\begin{eqnarray}\n", "Z &=& \\sum_E g[E] e^{-E/T}\\\\\n", " &=& \\frac{1}{Z} \\sum_E g[E] E^n e^{-E/T}\\\\\n", "c_v &=& \\frac{-^2}{T^2}\\\\\n", "F &=& -T \\log(Z)\\\\ \n", "S &=& -\\frac{dF}{dT}=log(Z)+\\frac{}{T}\n", "\\end{eqnarray}\n", "we introduce \n", "\\begin{equation}\n", "Z_0=g[0] e^{-E_0/T} \n", "\\end{equation} \n", "so that\n", "\\begin{eqnarray} \n", "Z_w &\\equiv& \\frac{Z_E}{Z_0} = \\frac{g[E]}{g[0]} e^{-(E-E_0)/T}\\\\\n", "Z &=& Z_0 \\sum_E Z_w(E) \\\\\n", " &=& \\frac{\\sum_E Z_w(E) E^n}{\\sum_E Z_w(E)}\\\\\n", "log(Z) &=& log(\\sum_E Z_w) +log(Z_0) = log(\\sum_E Z_w)+log(g[0])-E_0/T \n", "\\end{eqnarray} " ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [], "source": [ "(Energies, indE, Emin) = PrepareEnergies(N)\n", "\n", "def Thermod(T, lngE, Energies, N):\n", " Z = 0.\n", " Ev = 0. # \n", " E2v = 0. # \n", " # = sum_E e^(-E/T)*g[E] E^n /Z\n", " # Zw = Z/Z_0\n", " # = Z0/Z sum_E e^(-(E-Emin)/T) g[E]/g[0] E^N\n", " for i,E in enumerate(Energies):\n", " # Z += exp(log(g)) * exp(-E/T)\n", " # Z/Z_0 = w where Z_0 = g[0] exp(-Energies[0]/T)\n", " w = exp(lngE[i]-lngE[0]-(E-Energies[0])/T) # g(E)/g0 Exp(-(E-E0)/T)\n", " Z += w\n", " Ev += w*E\n", " E2v += w*E**2\n", " Ev *= 1./Z\n", " E2v *= 1./Z\n", " cv = (E2v-Ev**2)/T**2\n", " # Z_correct = Z * exp(lngE[0]-Energies[0]/T)\n", " # Entropy = log(Z_correct) + /T = log(Z) + lngE[0]-Energies[0]/T + Ev/T\n", " Entropy = log(Z)+lngE[0]+Ev/T-Energies[0]/T\n", " return (Ev/(N**2), cv/(N**2), Entropy/(N**2))" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "Te = linspace(0.5,4.,300)\n", "\n", "Thm=[]\n", "for T in Te:\n", " Thm.append(Thermod(T, lngE, Energies, N))\n", "Thm = array(Thm)" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3uklEQVR4nO3deWCU1b3w8e/JZJvsK1kJARJ2WcMiiuCuiFIVFa0LWrXa9r712vZWa6+tfetb22tb21prrXIVi1sXFC1aF1RwAWQJAmHfQwhk37eZOe8fZxKSOIEkM5OZyfw+9z59tjPP88tj+M2T85znHKW1RgghxOAX4usAhBBCDAxJ+EIIESQk4QshRJCQhC+EEEFCEr4QQgSJUF8HcDopKSk6NzfX12EIIUTA2LRpU7nWOtXVPrcTvlJqKLAMSAccwDNa6991K6OA3wHzgUZgidZ685mOnZuby8aNG90NUQghgoZS6nBP+zxxh28Dvqe13qyUigU2KaXe01oXdSpzOZDvnGYCf3LOhRBCDBC36/C11sfb79a11nXATiCrW7GFwDJtrAMSlFIZ7p5bCCFE73n0oa1SKheYAqzvtisLONppvZivfikIIYTwIo89tFVKxQD/AO7TWtd23+3iIy77dFBK3Q3cDZCTk+Op8IQQPtDW1kZxcTHNzc2+DmXQiYyMJDs7m7CwsF5/xiMJXykVhkn2y7XW/3RRpBgY2mk9GyhxdSyt9TPAMwAFBQXS0Y8QAay4uJjY2Fhyc3MxbTeEJ2itqaiooLi4mOHDh/f6c25X6Thb4DwH7NRa/6aHYiuBW5UxC6jRWh9399xCCP/W3NxMcnKyJHsPU0qRnJzc57+cPHGHfw5wC7BNKVXo3PYjIAdAa/00sArTJHMfplnm7R44rxAiAEiy947+XFe3E77W+hNc19F3LqOBb7t7LiH8Tks9bPkrjJgHQ8b4OhohTsuv37QVwq+1NsCLV0PxBkDBTa/CqEt9HZUQPZK+dITor/V/Nsl+4R8hOQ/e/W+w23wdlejGYrEwefLkjumxxx7r2Ldo0SIOHDjAzJkzmTx5Mjk5OaSmpnaUPXToEBdddBFVVVU+/Ak8R+7wheiP1gb4/EnIuxim3AwRsfDarbDrTRh/ta+jE51YrVYKCwu/sn3Hjh3Y7XZGjBjB+vXm1aHnn3+ejRs38uSTT3aUu+WWW3jqqad46KGHBipkr5GEL0R/bP8HNFbAnO+Z9TELIDoVilZKwu/BI2/uoKik+ys67hmXGcdPrhzfr88uX76chQsXnrHcVVddxZw5cwZFwpcqHSH6Y+urkDQScmaZ9RALjLkC9r4LbfKSkT9pamrqUqXz6quvAvDpp58ybdq0M34+MTGRlpYWKioqvB2q18kdvhB9VX0EDn8C5z8EnZvGjb0SNj0PB9fAqEt8Fp6/6u+duLt6qtI5fvw4qakuexH+iiFDhlBSUkJycrKHoxtYcocvRF8VrTTzsxZ13T7sXLBEwMGPBz4m0WdWq7XXLy41NzdjtVq9HJH3ScIXoq92roT0syBpRNftYZGQXQCHP/VNXKJPxo4dy759+85YTmtNaWkpg2EwJkn4QvRFXSkcXQ9jr3K9f9g5cHwrNHv24aTov+51+A888AAAV1xxBR999NEZP79p0yZmzZpFaGjg14AH/k8gxEDa9ZaZj73S9f5hs2GNw3wp5F88cHGJHtntdpfbFy1axPnnn88jjzyCxWIBYMmSJSxZsqRLuRdffJFvfetb3g5zQMgdvhB9sfNN85JVag/dKAydASoEjm4Y2LhEn1mtVh555BGOHTt22nITJkzgwgsvHKCovEvu8IXorcZKOLgWzvk/XVvndBYeDWnjofiLgY1N9Mull565K4y77rprACIZGHKHL0Rv7foXaHvP1TntsqfDsU3gcAxMXEL0kiR8IXprxz8hMRcyp56+XPZ0aKmF8t0DEpYQvSUJX4jeqC+DAx/D+Gt6rs5plz3DzKVaR/gZSfhC9MbON0x1zoRrz1w2eSREJkjCF35HEr4QvbH9n5Ay2jyQPROlTLVO8UbvxyU87r777mPNmjVcffXVTJ48mby8POLj4zva8X/22WcsXryYvXv3+jrUPvNIwldKLVVKnVRKbe9h/zylVI1SqtA5PeyJ8woxIGpL4PBn5u6+t8PKZU+Hkzuhuca7sQmPqqysZN26dZx33nmsWLGCwsJCnn32WebMmUNhYSGFhYXMnj2be++9l1/96le+DrfPPNUs83ngSWDZacqs1Vov8ND5hBg4O14HNEy4pvefGTrdfObYZhh5vpcCCzBvPwCl2zx7zPSz4PLHzlhs2bJlPP744yilGDFiBFu2bOHAgQOEhITQ2NjI6NGjOXDgAH//+9+57LLLzni8OXPmsGTJEmw2W0C9geuRO3yt9Rqg0hPHEsLvbP8HpE+ElPzefyZrGqCkWscP7Nixg0cffZTVq1ezdetWnnvuOSZNmsTHH5tO7t58800uvfRSwsLCet1lckhICHl5eWzdutXb4XvUQH41na2U2gqUAN/XWu9wVUgpdTdwN0BOTs4AhieEC1WH4NhGuOiRvn0uMh5SR8uD2856cSfuDatXr2bRokWkpKQAkJSUxA033MCrr77K+eefzyuvvNLRdUJ/ukzuzReEvxioh7abgWFa60nAH4DXeyqotX5Ga12gtS7o7YUXwmu2/9PM+zOKVXaBSfhaezYm0Sdaa1S3Zy9XXXUVb7/9NpWVlWzatIkLLrgAGPxdJg9Iwtda12qt653Lq4AwpVTKQJxbCLfs+Kd5AJs4rO+fzZ4OTZVQecDzcYleu/DCC3nttdc6RqyqrKwkJiaGGTNm8N3vfpcFCxZ0dJ7W2y6TAfbs2cP48b4Z1KW/BiThK6XSlfMrVik1w3newB8vTAxuZXvMQ8betL13RV7A8gvjx4/noYceYu7cuUyaNIn7778fgBtuuIG//vWv3HDDDR1le9tl8okTJ7BarWRkZHgrbK/wSB2+UuplYB6QopQqBn4ChAForZ8GFgH3KqVsQBOwWGv5O1f4uR0rAAXjvta/z6eOhvBYk/AnLfZkZKKPbrvtNm677bYu2xYtWkT3NDRnzhwefPBBqqurSUhIAGDevHnMmzevS7mXXnqJb37zm94M2Ss8kvC11jeeYf+TmGabQgSOotch52yI6+ddXIgFsqbKHX6A+fWvf82RI0c6Er4rCQkJ3HLLLQMXlIfIm7ZCuFK2B04WwbiF7h0nezqUbofWRs/EJbxu5syZTJw48bRlbr/99oBqf99OEr4QrhS9YebjehjKsLeGzjB98JRscT8mIdwkCV8IV4peh6GzIC7TveNkFZj50XVuhySEuyThC9Fd+T44sd396hyA6GQYMt50rSyEjwVeJZQQ3lb0upl7IuEDjJgLG5dCWzOERXrmmKLXLBYLZ511Vsf64sWLeeCBB3os/9FHHxEeHs7s2bO9HtsTTzxBUlIS69ev59NPP6W1tZWDBw8yevRoAH784x+zbt065s+f3/FymDsk4QvRXdHrpg19fJZnjjd8Lqx7Coo3wPDzPHNM0WtWq5XCwsJel//oo4+IiYlxmfA92VmazWZj6dKlbN68mVtvvRWAQ4cOsWDBgi7xTp8+nbvuussjCV+qdITorGK/edlq/Nc8d8xhsyEkFPZ94LljCrfl5ubyk5/8hKlTp3LWWWexa9cuDh06xNNPP81vf/tbJk+ezNq1a1myZAn3338/559/Pj/84Q8pLCxk1qxZTJw4kauvvpqqqirAtNe/7777mD17NhMmTGDDhg04HA7y8/MpKysDwOFwkJeXR3l5OatXr2bq1Kln/AIZNmwYFRUVlJaWuv0zyx2+EJ21t84Z62brnM4i40x7/r3vwsV97IRtEPnlhl+yq3KXR485JmkMP5zxw9OWaWpqYvLkyR3rDz74YMfbtSkpKWzevJmnnnqKxx9/nGeffZZ77rmHmJgYvv/97wPw3HPPsWfPHt5//30sFgsTJ07kD3/4A3PnzuXhhx/mkUce4YknngCgoaGBzz77jDVr1nDHHXewfft2br75ZpYvX859993H+++/z6RJk0hJSel1z5wAU6dO5dNPP+Xaa/v51reT3OEL0VnRG6Zr44Shnj3uqMtMu/6qw549rjij9iqd9qlzVwrXXGPGOJg2bRqHDh3q8RjXXXcdFouFmpoaqqurmTt3LmDe4F2zZk1HuRtvNO+gnnfeedTW1lJdXc0dd9zBsmVmqJClS5dy++23A/3rmdNdcocvRLuqw3C8EC7+meePPeoyePch2PMOzAy8V/I94Ux34r4QEREBmAe7Nputx3LR0dG9Ol73XjmVUgwdOpS0tDRWr17N+vXrWb58OeCbnjnlDl+IdjtXmrknq3PapeRB6ljn6FnCn8XGxlJXV+dyX3x8PImJiaxduxaAF198seNuH+DVV18F4JNPPiE+Pp74+HgA7rzzTm6++Wauv/76fvfMOWHChH7/TO0k4QvRrugNM7JV0nDvHH/CNXDkczNGrhgw7XX47dPpmmQCXHnllaxYsaLjoW13L7zwAj/4wQ+YOHEihYWFPPzwqSG6ExMTmT17Nvfccw/PPfdcx/arrrqK+vr6juocgMsvv7xLdVBP2tra2LdvHwUFBb35cU9LqnSEAKg5Zjo5u+C/vXeO8dfAh4+aQVVmf8d75xFd2O12l9s719kXFBR0dIs8atQovvzyy459c+bM6fK5yZMns26d6zenr732Wn7xi198ZfvWrVuZNGkSY8aM6dg2bNgwkpOT2bt3L/n5ZvjM3Nxctm/f3uWzb731FosWLfJIc1C5wxcCYOebZt7frpB7IyXPdLWw5UUZBSuIPPbYYz1+ETz22GMcP378tJ+32Wx873vf80gsyp+7pS8oKNAbN8og0GIA/O98aKqCb33u3fNsXgYr/wPueBdyZnr3XH5g586djB071tdhDFqurq9SapPW2mX9j9zhC1F3Ag5/5rmuFE5n/DVmUJQNf/b+ufyEP99UBrL+XFePJHyl1FKl1Eml1PYe9iul1O+VUvuUUl8qpaZ64rxCeMTOlYAemIQfEQPTbjOtdaqPeP98PhYZGUlFRYUkfQ/TWlNRUUFkZN/6ZvLUQ9vnMSNaLeth/+VAvnOaCfzJORfC97a+AkPGQeqYM5f1hFn3wvqn4ZMnYMFvBuacPpKdnU1xcXFH1wLCcyIjI8nOzu7TZzw1xOEapVTuaYosBJY5x7Fdp5RKUEplaK1P/7RCCG8r2w3HNsIlj0K3l2a8Jj4bpt1uetCcdS+k5A/MeX0gLCyM4cO91MxV9NlA1eFnAUc7rRc7t32FUupupdRGpdRGuSsQXle4HJQFJl4/sOed9wCERcF7PxnY84qgNlAJ39Wtk8tKPa31M1rrAq11QW/7mRCiX+w22Poq5F8CMUMG9tzRKTDnP2H3v+DgV1/uEcIbBirhFwOde6PKBuR1Q+Fb+1dDfSlM+bpvzj/rWxCfA29+F1obfBODCCoDlfBXArc6W+vMAmqk/l743Bd/gagUyL/UN+cPs8LXnoLKA/Dvh3wTgwgqnmqW+TLwOTBaKVWslPqGUuoepdQ9ziKrgAPAPuAvwLc8cV4h+u1EkemffuY9EBruuziGz4HZ/wGb/hd2rfJdHCIoeKqVzo1n2K+Bb3viXEJ4xGd/MA9Np3/D15HABT+GAx/Bim/CN96DIQPUPFQEHXnTVgSfmmOw7TWYeitEJfk6GgiNgBtfNlU8L10PDeW+jkgMUpLwRfD57Pem87JZflSzGJ8Ni1+G+hMm6TfX+joiMQhJwhfBpfIAfPGcaZmTOMzX0XSVPQ2uex6Ob4WXbpCWO8LjJOGL4KE1vPMgWMJg3o98HY1roy+Ha5+Fo+vglZugrXdD4AnRG5LwRfDY9nczpuwFP4a4DF9H07PxV8PX/gQHPjZJX+70hYdIwhfBof4kvP0DyJ5ummL6u0mLYeGTcOBDePEaaKr2dURiEJCELwY/hx3e+La5U174Rwix+Dqi3plyMyz6Xzi2CV5YYL60hHCDJHwx+L3/U/OS1WWPQepoX0fTN+O/Bje9AhX7YellQdGHvvAeSfhicNv8ommGOf1O/3jJqj/yLoJbXjft85+9GEq2+DoiEaAk4YvBq/BlM37siHnm7j6Q5cyEO94BSzgsvRyKVvo6IhGAJOGLwWnzMnj9XtNXzeKXTVPMQJc2Du76ANLGw2u3wJr/AYfD11GJACIJXwwuDju8+9/mzn7k+XDjqxAe5euoPCdmCCx5CyYsgtU/l64YRJ9IwheDR32ZeUO1vc7+ptcGV7JvF2Y1L2dd8Ws4uAaePhcOfeLrqIQHaK1ptbfS2NboleN7ahBzIXxr/2pYcY9pr37FbwL3AW1vKWW+1LJnwN+WwPNXQME34KKfQmScr6MLWDaHjWZbM832Zprammixt9DqaKXV3kqLvcWsO5db7ae2tzq+uq1L2U7H6P75NkebWXa0YnPYAEixpvDh9R96/OeThC8CW1MVvPtj2PJXSBkFN/8D0s/ydVQDJ2MifHMNfPgorPsT7H4bLv8ljL1y4AZlH2Ct9lYa2hq+MjXaGmmyNZmEbWumyd7Udd253GRv6rLeuWx7wu2v8JBwIiwRhFvCCbecWm6fR4dGkxSR1GVbWEgYYZYwwkPCO+Yx4TEeulpdScIXgatoJaz6vqnDPuc+58DgVl9HNfAiYuCyX5h6/ZX/YR7o5pwNFz1iWvf4Ca019W311LbWUtdaR21LLbWtzqmllvq2+q5J3NZAQ6uZN7Y10tDWQH1bfZ+ScnhIOJGhkUSGRhIVGmWWLZFYQ60kRSYRGWqWraFWIi2RXdYjLBFEhEYQERJx2iTevhwWEkaI8u9ackn4IvDUlZpEv/NNczd/02uQOdnXUfle9jT45sew5UX46DFYegmMng/n3g9Dp3vsNO2Ju7q5msqWSqqaq8zUUkVNS01HAm9P5nWtdR1zh+65VZFCERUWRXRoNFFhUcSExRAdFk1iZCLRYdE9T6HRRIdHExUaZRK3M2lHWiKxBMpb1QNEmcGo3DyIUpcBvwMswLNa68e67Z8HvAEcdG76p9b6Z2c6bkFBgd64caPb8YlBQmtTdfPuQ6YXyXkPmOEBB0OTS09rbYDPn4LPn4Tmahh2DpzzXfMSl4skaHPYqGiqoKypjLLGMsqayqhsruxI5B1J3bne5mhzedrQkFDiwuPMFGHmseGxHdviI+K7butUJjos2u/vkAOBUmqT1rrA5T53E75SygLsAS4GioEvgBu11kWdyswDvq+1XtCXY0vCFx0qD8Jb95mhAHNmw1W/h5R8X0fl/1rqqP3iL5RsepbS5kpOxqZSljmRsoQsyuwNlDWWcbLxJJXNlWi+mgtiwmJIjEwkMTKRpIgkEiMTSYhM6FhOjEwkKdK5HJGINdSKGqTPDgLF6RK+J6p0ZgD7tNYHnCd7BVgIFJ32U0L0hsMO6582bc6VxbTAmXY7hMidYLva1lqO1h2lpL6k69Rg5vVt9ZBgAVIBUNVbSarcwpDQaFLihjIu61xSYzJItaaSak1lSNQQkq3JJEWah4ti8PBEws8CjnZaLwZcPSk6Wym1FSjB3O3vcHUwpdTdwN0AOTk5HghPBKwTRbDyO6a3yPxLYcFvzFCAQajN0cbRuqMcrjnModpDZqox88rmyi5lo8OiyYzJJCs6i2lp08iKySIjOoOM6AyGRA0hqbGasMLlsPUVOLAHwtbBmPkw4VrInG3G2BWDkieqdK4DLtVa3+lcvwWYobX+j05l4gCH1rpeKTUf+J3W+ox/j0uVTpCytcDaX8Pa35g25Zf/yiSjIKkqqGquYnfVbnZX7mZP1R72VO1hf/X+LvXmSZFJ5MblkhufS25cLjmxOWTFmsQeFx7Xu2oVhx0Ofwrb/wFFb5gmrhHxkHeB+YLNvxiiU7z4kwpv8HaVTjEwtNN6NuYuvoPWurbT8iql1FNKqRSttbwTLro69Cn8634o2wVnXW86PYtO9nVUXtNsa6aoooht5dvYWraVL8u+5ETjiY79KdYURieO5uxxZ5OXkNeR5OPCPfByVYgFhp9npvmPm+cjRa/D3vdgxwpAQdY0GHUpDJ8LmVMgVKp4ApknEv4XQL5SajhwDFgM3NS5gFIqHTihtdZKqRmYLh0qPHBuMVg0lJs+cLa+BPE5cNPfYNQlvo7K45psTWw5sYV1pev44vgX7KrchU2bduVZMVlMHTKVccnjGJU0itGJo0m2DtCXnSXM3NHnX2w6ZCvdCnvehb3/hg//n3mxKywKcmZB7rmQO8d8AUgLqYDidsLXWtuUUt8B/o1plrlUa71DKXWPc//TwCLgXqWUDWgCFmtPtAcVgc/hgM0vmEFKWutNm/HzfjCo+sApqS9h9ZHVfHj0Q7ac3EKbo43QkFAmpkzk9gm3MzF1IhNSJpBi9ZPqk5AQk8wzp8C8H0JDhan6OfSJmT5wtqgOizJ/AWRPPzXFpPo2dnFaHmmH7y1Shz/IlW6Dt+6H4g0w7FzTGdiQMb6OyiNqWmpYdXAVb+x7gx0Vpn1CXkIe52ady8yMmUwdMpWosAD9UmsoP/UFUPyF+e/Y/vZrYm6nL4ACSDtLqoEGmFfb4XuTJPxBqr7MVBFsfgGsSXDpozDxhkHxULaooohlRct479B7tDpaGZ04mvkj5nNhzoUMixvm6/C8o60JSgpN8m+f6o6bfZZwGDLOvAmdMdnMh4yXLwEv8vZDWyF6x9ZiOvha+2toa4QZd5u3Za2Jvo7MLVpr1peu57ltz7Hu+Dqiw6K5Ov9qrsm/hnHJ43wdnveFWWHY2WZqV1NsEn/JFvNlsGMFbHre7JMvAZ+RO3zhfVqbZn/vPQzVh2HUZXDJzwP+TVm7w877R95n6falFFUUkWJN4eaxN3P96OuJDY/1dXj+RWuoOmiS//HCU/PmGrPfEg6pYyBtghnRK22cWY4Z4ruYA5Tc4Qvf0Br2vW+qb0q2mLu6W143I1EFsCZbE2/se4NlRcs4WneUnNgcfnL2T7hy5JVEWOSlJZeUgqQRZppwjdnW/UugdLsZ12DrS6c+F51qvgCGjHd+EYyH1NHB2SuqB0jCF56nNRxaa7pDOLoeEnJg4R9h4mKwBO6vXHlTOX/b/Tde3vUyVS1VTEyZyH9O+08uGHqB9MrYH66+BMC0Cjq5A07sgBPbzRvXG5eCran9g5AwFJLzIDnf/KWYnGemuCzpduM0Avdfn/A/dhvsXGl6aDy2CWIzTd83U24J2PpZh3awoXQDf9v9N1YfWY1N25iXPY8lE5YwdchU6SjMG6KTT70Q1s5hNx3ondwBJ3dBxV4o3wtHl5vmvO1Crc7kPwIShkHiMDNPyIH4oYOquW9/SMIX7muqgsKXzQPZmiPmjm3+4ybRh0X6Orp+OVp7lFUHV/HmgTc5XHuY+Ih4vj726ywatYjc+Fxfhxd8QiyQkmemcQtPbdca6k+Y5F+xF8r3QcU+89fB7rfB3tr1ONFDTPJvn+IyzRSbCXEZZn8A/xV6JoP3JxPe5bDD/g+hcDns+hfYW0y3xZc/Zh7KBmAVR0l9CR8c+YC3D77NtvJtAExLm8Y9k+7h4mEXS/28P1IKYtPNNHxO130Oh/kyqD5iGgtUHzbLVYfNM6Wdb0L3fv1VCMSkQWyG84sgw3wRxGaaB8gxQ8yXQlRyQH4xBF7EwnccdlMnv/NN2PE61JWYJpXTlsCUr0PGJF9H2Cdt9jYKywpZW7yWNcVr2F+zH4CxSWP53rTvcdnwy0iPTvdxlKLfQkJMso7LcD3Uo8MBjeVQW2LeG+iYHze/2xX7zbOo9pZEXSiISjLJPzrl1BdBTKp50Ny+HJVs3jWJiPWL90wk4YvTa6k3v/R73jF38g1lpgndyAvNOKqjLw+Y7nQb2xopLCtk84nNbD65mW1l22i2NxMaEsq0tGlcnX81c7PnSpVNsAgJOXXXzuSey7U2mGE1609Cw0nnvMxM7cslW8wLha11PZwr1NwcWRPNF4A10XxhtG/rWG7flwzxWR7/kSXhi660Ni0j9n1gmlQeWWf+7A2PgfxLYOwCM4/w73bmzbZm9lTtoaiiiJ2VOymqKGJv1V7s2k6ICmF04mgWjVpEQVoBMzNmEhMe4+uQhb8Kj4bkkWY6k7Ym55dAuflyaKwwz7gaK828qdIs1xRD6ZdmW1vjV48TlQz/dcDjP4ok/GDX3ha6vWOsAx9DfanZlzYBZt1rxkHNmeW3d/K1rbXsr97PrspdFFUUUVRRxP7q/di1HYCEiATGJY/jjgl3MC1tGpNSJ0mCF94RZjUtgxL70I1GW/OpL4P2L4f2vok8TBJ+sOme4A99ArXHzL7oVNPtbd6FpsomLsO3sXbT2NbI/ur97Kvex77qfeyv3s/e6r2cbDzZUSYpMomxyWOZmz2X8cnjGZc8jvTodGk+KfxXWCSEZQzIvzdJ+INda4OpXyz+Aoo3mqn9Dj461dm3ubN/85RRPn+wZHfYKWko4WjtUQ7XHeZI7REO1x7mQM0BjtUf6ygXYYlgRPwIZqbPJC8xj7yEPEYljiItKk2SuxA9kIQ/mDTXON9O3GG6rC3ZbN5SdFZtkDQCRsyFoTN8luC11lS1VHG8/jjHG8xUUl9ixmutPUxxfTG2Tn/OWkOt5MTmMDFlItfkX8PIhJHkJ+STFZMlb7cK0UeS8AON1qaOr3K/ecGkYj+cLDIPWquPnCoXmWB6IZxzv+mbPKvA60MFNrY1Ut5UTkVzBeVN5V2m9gRf2lBKs725y+esoVayY7PJS8jjgpwLGBY3jJzYHIbFDSPFmiJ37EJ4iEcSvlLqMuB3mBGvntVaP9Ztv3Lunw80Aku01ps9ce5Bp71tcE2xaRdcWwK1xVBzzNS9V+zr2i5YWcyr5FkFpj182lmmg6m4zH7fvWutabI1Ud1STU1LDTWtNWbeeWqt6djfntSbOvo6OSVEhZAUmURGdAb5ifnMzZ5LRkwG6dHpZEZnkhGdQXxEvCR1IQaA2wlfKWUB/ghcjBnQ/Aul1EqtdVGnYpcD+c5pJvAn53zwcjigrQFa6kxb9tY653KdaarVWGHu1BvKncvlptOouuNfffvPEm7e+EsaDhMWQfJIdNJIbIm5tMSm0YKDFnsLzfZmWu2tNLeU0Vp67NS6vZkmWxONbY3Ut9XT0NZw2uW6trou1SrdWUOtxIXHER8RT0JEQsfwfO1TcmSymVuTSYxIlKoXIfyEJ+7wZwD7tNYHAJRSrwALgc4JfyGwzDmO7TqlVIJSKkNrfdwD5/+KJ16/kVZ7KxptqkCgY9n8n7mLpX2tUxmt27fRXtL5/w60w4HWdtBmrh0OUz+uO2132NEOG9qZMDuPNuA8GnalsAH2kFBsoeHYLaHYIsKwW2OxpSViD7FgC7GYOWAHbNqG3VGDrX4DrTWf0LK3BYd29Ov6RIVGER0W3WXKisnqsp4QkUBCRAJxEXHEh8cTH3Fqki4GhAhMnkj4WcDRTuvFfPXu3VWZLOArCV8pdTdwN0BOTk6/AlpR+SXNChRmotO8fVkBSnddRynQX/2cqW1QznXVUf2glEKFmNJKhZoySqFUSMdEx7LF9NOhQggNjSDUEo4lJAxLiIVQFUpoSCgRymLWQ0IJVaFYQixYlHM9JBSLMuvhlnAiLBFEWCKIDI0k3BJOpCWyY1tEaMSp5U7lYsJiiAqLIkRJ97FCBCNPJHxXla/dh9HqTRmzUetngGfAjHjVn4A+vmNHfz4mhBCDmidu9YqBoZ3Ws4GSfpQRQgjhRZ5I+F8A+Uqp4UqpcGAxsLJbmZXArcqYBdR4q/5eCCGEa25X6WitbUqp7wD/xjTLXKq13qGUuse5/2lgFaZJ5j5Ms8zb3T2vEEKIvvFIO3yt9SpMUu+87elOyxr4tifOJYQQon+kuYYQQgQJSfhCCBEkJOELIUSQkIQvhBBBQhK+EEIECUn4QggRJCThCyFEkJCEL4QQQUISvhBCBAlJ+EIIESQk4QshRJCQhC+EEEFCEr4QQgQJSfhCCBEkJOELIUSQkIQvhBBBQhK+EEIECbdGvFJKJQGvArnAIeB6rXWVi3KHgDrADti01gXunFcIIUTfuXuH/wDwgdY6H/jAud6T87XWkyXZCyGEb7ib8BcCLziXXwC+5ubxhBBCeIm7CT9Na30cwDkf0kM5DbyrlNqklLr7dAdUSt2tlNqolNpYVlbmZnhCCCHanbEOXyn1PpDuYtdDfTjPOVrrEqXUEOA9pdQurfUaVwW11s8AzwAUFBToPpxDCCECWn2LjRO1zTS02JiYneDx458x4WutL+ppn1LqhFIqQ2t9XCmVAZzs4RglzvlJpdQKYAbgMuELIcRg43BoKhpaOVHbTGlNM8drmzlR00ypc73UuV7XYgMgJSaCjT/uMfX2m1utdICVwG3AY875G90LKKWigRCtdZ1z+RLgZ26eVwgh/ILDoSmvb6G4uoljVU0dCbzz/GRdM232rhUWIQqGxEaSHh9JXmoM5+alkB4fSVpcBOlxVq/E6m7Cfwx4TSn1DeAIcB2AUioTeFZrPR9IA1YopdrP95LW+h03zyuEEAOize6gtKaZY86E3nleXNVISU0zrTZHl89Ywyykx0eSHhfJjOFJpMVFkhEfSVqcSfAZ8ZGkxERgCVED+rO4lfC11hXAhS62lwDzncsHgEnunEcIIbxFa015fStHKhs4XNHI4YpGjlQ2UlzVaO7Ya5txdHuamBITQVailfGZ8VwyPp2sBKuZEq1kJliJiwzFeZPrV9y9wxdCCL9nd2hKqps4UtnIoYoGjjgT++HKRo5UNNDQau8oqxRkxpvkPWtEMlmJp5J5VoJJ6JFhFh/+NP0nCV8IMWhUNbRyoLye/Scb2F9Wz/6yeg6UNXC0qrFLHXq4JYTsJCu5ydHMHJ7EsOQo5xRNdqKViNDATOhnIglfCBFQ7A5NcVWjSejdEntFQ2tHufDQEIYnRzM6PZZLxqd3SerpcZEDXn/uDyThCyH8ktaasvoWdpfWsbu0jl3O+d6TdTS3nXpImhwdzsjUGC4Zn8bI1JiOKSvRGpRJ/XQk4QshfK65zc7u0jp2Hq/tSOy7T9RR2emOPSUmgjHpsXx95jBGpcWQNySGESkxJEaH+zDywCIJXwgxoJpa7RQdr2VHSQ3bimvYXlLL3hN12JxNYaLCLeSnxXLx2DRGp8cyJj2W0emxJMdE+DjywCcJXwjhNc1tdnaU1PBlcQ3bjtWw/VgN+07WdzRzTIoOZ0JWPBeMSWVCZjzjM+PJTrQSIlUxXiEJXwjhEVpriqua2Hykii1HqtlypIqi47UdrWNSYiI4KyuOy8anMyErnglZ8WTER/ple/XBShK+EKJfmtvsbD1azeYj1R1Jvry+BTBvmk7MjufOOSOYPDSByUMTSIuL9HHEQhK+EKJXGlttbD5czfqDFaw/UEnh0Wpa7aa1zPCUaM4blcKUnESmDE1gTHosoRYZQdXfSMIXQrhU19zGxsNVrD9QyfqDFWwrrsHm0FhCFBMy47ht9jBmDk9m6rBEkqSlTECQhC+EAMBmd7C1uIa1e8tYu7ecwqPV2B2a0BDFxOx47jpvBDOHJ1GQm0RMhKSOQCT/1YQIYkcrG1mzt4y1e8r5dH85dc02lIKJ2QncO3ckZ49MZmpOItbwwdnVQLCRhC9EEGmzO/jiUCXvFZ3gw10nOVTRCEBWgpUFEzOYk5/K7JHJJERJFc1gJAlfiEGutrmNj3eX8f5Ok+Rrm21EhIZwTl4Kt58znDn5KQxPiZbmkUFAEr4Qg1BlQytvbz/OO9tLWXeggja7Jik6nEvHp3PRuDTm5KcQFS7//IONW//FlVLXAT8FxgIztNYbeyh3GfA7wIIZCesxd84rhPiq2uY23ttxgpVbS/hkXzl2h2Z4SjR3nDOci8elMSUnUToTC3LufsVvB64B/txTAaWUBfgjcDFQDHyhlFqptS5y89xCBL1Wm4MPdp7g9cJjfLi7jFabg+xEK3efN4IrJ2YyNiNWqmpEB3eHONwJnOkXagawzznUIUqpV4CFgCR8Ifpp38k6Xv3iKP/cfIyKhlZSYyP4+swcrpyUyZShCZLkhUsDUYmXBRzttF4MzOypsFLqbuBugJycHO9GJkQAabHZeWvrcV7acIRNh6sIDVFcOHYIi6fncN6oVKmuEWd0xoSvlHofSHex6yGt9Ru9OIer30LtYpvZofUzwDMABQUFPZYTIliU17ewfN0RXlx3mPL6FkakRvOj+WO4eko2qbHSZbDovTMmfK31RW6eoxgY2mk9Gyhx85hCDHq7Smt5bu1B3igsodXu4PzRqdxx7nDOzUuRKhvRLwNRpfMFkK+UGg4cAxYDNw3AeYUISNuP1fD7D/bybtEJrGEWbpg+lCXn5DIyNcbXoYkA526zzKuBPwCpwL+UUoVa60uVUpmY5pfztdY2pdR3gH9jmmUu1VrvcDtyIQaZwqPV/P6DvazedZK4yFDuuyifJbNz5a1X4TFKa/+tJi8oKNAbN7ps2i/EoHGwvIH/+fcuVm0rJSEqjLvmjOCWs4cRFxnm69BEAFJKbdJaF7jaJ6/aCeEjFfUt/P6DvSxff4Tw0BC+e2E+d503QnqiFF4jv1lCDDC7Q/PShiP86p1dNLbauWH6UO67MJ8hMiKU8DJJ+EIMoB0lNfxoxXa2Hq3mnLxkHrlqPHlDYn0dlggSkvCFGAB2h+bpj/fz2/f2kBAVxhM3TGbh5ExpXikGlCR8IbyspLqJ/3y1kPUHK1kwMYOff22CtLwRPiEJXwgv+nRfOd9avpk2u4PHr5vEtVOz5K5e+IwkfCG85MV1h/npyh2MTI3mz7cUMDwl2tchiSAnCV8ID9Na85v39vCH1fu4cMwQnlg8mVhpUy/8gCR8ITzI4dD87K0inv/sEIunD+XRq8+SXiyF35CEL4SHaK15dNVOnv/sEHeeO5yHrhgr9fXCr4T4OgAhBosnV+/juU8OsmR2riR74Zck4QvhAcs+P8Sv39vDNVOyeHjBOEn2wi9JwhfCTa9vOcbDb+zgorFp/HLRREKkzl74KUn4Qrjh/aITfO9vWzl7RDJP3jSFMIv8kxL+S347heinz/dX8O2XNjM+M46/3FZAZJjF1yEJcVqS8IXohy1HqrjzhS/ISYri+dtnSJfGIiC4lfCVUtcppXYopRxKKZcd7jvLHVJKbVNKFSqlZEQTEdCKSmq5bekGUmIj+OudM0mKln5xRGBw97ZkO3AN8OdelD1fa13u5vmE8Kk9J+q45bn1REeEsvzOmaRJH/YigLiV8LXWOwFpgiaCwpfF1dy6dAPhlhCW3zmT7MQoX4ckRJ8MVB2+Bt5VSm1SSt19uoJKqbuVUhuVUhvLysoGKDwhTm/DwUpu+st6YiJC+fs9sxmRGuPrkIToszPe4Sul3gfSXex6SGv9Ri/Pc47WukQpNQR4Tym1S2u9xlVBrfUzwDNgBjHv5fGF8Jr3i07wnZc3k5lgZfmdM8mIt/o6JCH65YwJX2t9kbsn0VqXOOcnlVIrgBmAy4QvhL/QWvPcJwd5dNVOJmTG87+3TyclJsLXYQnRb16v0lFKRSulYtuXgUswD3uF8Fttdgc/WrGdn/9rJ5dPSOe1b54tyV4EPHebZV6tlCoGzgb+pZT6t3N7plJqlbNYGvCJUmorsAH4l9b6HXfOK4Q3Ha1s5LqnP+flDUf4zvl5PHnjVKzh8lKVCHzuttJZAaxwsb0EmO9cPgBMcuc8QgyUt7cd57/+8SVo+ONNU7liYoavQxLCY+T1QCGAqoZWfv6vnfxjczGThibw5I1TGJokzS7F4CIJXwQ1rTWvFx7j/761k9qmNr59/ki+e+EowkOl1xEx+EjCF0Fr0+Eqfvn2LjYcqmRKTgK/uOYsxqTH+TosIbxGEr4IOrtKa/nNu3t4t+gEqbER/L+rz+KG6UNl7Fkx6EnCF0FBa836g5U8/fF+PtpdRkxEKN+/ZBR3nDucqHD5ZyCCg/ymi0GtxWbnne2lLP30EFuPVpMcHc73Lh7FLWcPIyFKerkUwUUSvhiU9pfV88qGI/x9UzFVjW0MS47i51+bwKJp2TJQiQhakvDFoFHfYuPf20t5beNR1h+sJDREcfG4NG6ckcO5eSky1qwIepLwRUBrtTn4eE8ZbxQe4/2dJ2huc5CTFMV/XTaaRdOyGRIr/dUL0U4Svgg4zW12PtlbzntFJ3hnRyk1TW0kRYdz3bShLJycydScRLmbF8IFSfgiIFQ2tLJ610ne3VHK2r3lNLXZiY0I5aJxaVw1OZNz81IIs8jLUkKcjiR84Zdsdgdbi2v4dF85a/eWselwFQ4N6XGRXFeQzcXj0pg5PFneiBWiDyThC79gd2h2lday8VAVn+4r5/MDFdQ121AKJmTG8+3z87hkXDoTsuJkSE0h+kkSvvCJuuY2thypZuPhKjYfrmLLkSoaWu0AZCVYueKsDM7NT2H2yBSSoqW9vBCeIAlfeJXWmhO1LRQdr6GopJai47UUldRyqKIRgBAFY9LjuGZqNgW5iUzNSSQ70Sp38UJ4gSR84RHNbXaOVDZyoKyBA+X1HCxr4GB5A/vL6qlqbOsoNyw5inEZcVw7NZspOYlMGhpPbGSYDyMXInhIwhdnpLWmttlGSXUTx2uaKKludi6b+THnpDsNOZ8aG8HwlGguHZ/O2Iw4xmXGMSY9VpK7ED7kVsJXSv0PcCXQCuwHbtdaV7sodxnwO8ACPKu1fsyd8wr3aK1pbLVT29xGbZON6sZWyutbqWhoMfP6FsrrW6iob6WioZWTtc0d9evtLCGK9LhIMhMimTYskWunZjMiNZoRKTHkpkRJYhfCD7l7h/8e8KDW2qaU+iXwIPDDzgWUUhbgj8DFQDHwhVJqpda6yM1zDxpaa2wOjc2uaXM4sDvnNrvZZnM4sDk0bXazrcXmoKnNTlOrneY2MzU5p+bWTsttDuqbbdQ2t1HnnNc2tVHbbMPu0C5jUQoSo8JJjg4nJSaC8ZlxzB2VSmZCJJkJVjLirWQlWEmNjZDuhIUIMO6Oaftup9V1wCIXxWYA+5xj26KUegVYCHgt4S/4w1qanHek2vk/7elNa91pGdrXtKZLlYR2rmhObdfoTsudy3ff7uKznc7bHk+b3YHdYZK9J1nDLFjDLUSGhhAdEUqcNYyUmHBGpEYTFxlGnDXUOQ8jLjKMeGsYKbHhJEdHkBgVRqi8wCTEoOTJOvw7gFddbM8CjnZaLwZm9nQQpdTdwN0AOTk5/Qokf0gsrTYHOG9AlTlu+ypK0Wn51HYUtK91LdNte6cPqC7n6NjabfupO+HOZcIsilCLwhISQliIItQSQphFYWlfds5DQ0y50BCzPyLUgjU8hMgwS0dyt4ZZiAyzEBEaIi1chBAunTHhK6XeB9Jd7HpIa/2Gs8xDgA1Y7uoQLrb1eEurtX4GeAagoKCgX7e+v71hcn8+JoQQg9oZE77W+qLT7VdK3QYsAC7UWrtK0MXA0E7r2UBJX4IUQgjhPrcqa52tb34IXKW1buyh2BdAvlJquFIqHFgMrHTnvEIIIfrO3adzTwKxwHtKqUKl1NMASqlMpdQqAK21DfgO8G9gJ/Ca1nqHm+cVQgjRR+620snrYXsJML/T+ipglTvnEkII4R5pfyeEEEFCEr4QQgQJSfhCCBEkJOELIUSQUK6bzvsHpVQZcLifH08Byj0YjjcFUqwQWPEGUqwQWPEGUqwQWPG6E+swrXWqqx1+nfDdoZTaqLUu8HUcvRFIsUJgxRtIsUJgxRtIsUJgxeutWKVKRwghgoQkfCGECBKDOeE/4+sA+iCQYoXAijeQYoXAijeQYoXAitcrsQ7aOnwhhBBdDeY7fCGEEJ1IwhdCiCAR0AlfKXWZUmq3UmqfUuoBF/vnKaVqnD15FiqlHvZFnM5YliqlTiqltvewXymlfu/8Wb5USk0d6Bi7xXOmeP3p2g5VSn2olNqplNqhlPquizJ+cX17Gas/XdtIpdQGpdRWZ7yPuCjjL9e2N7H6zbXtFJNFKbVFKfWWi32evbZa64CcAAuwHxgBhANbgXHdyswD3vJ1rM5YzgOmAtt72D8feBszQtgsYL2fx+tP1zYDmOpcjgX2uPhd8Ivr28tY/enaKiDGuRwGrAdm+em17U2sfnNtO8V0P/CSq7g8fW0D+Q6/Y3B0rXUr0D44ul/SWq8BKk9TZCGwTBvrgASlVMbARPdVvYjXb2itj2utNzuX6zDjLmR1K+YX17eXsfoN5/Wqd66GOafuLT385dr2Jla/opTKBq4Anu2hiEevbSAnfFeDo7v6h3O280+8t5VS4wcmtH7p7c/jT/zu2iqlcoEpmLu7zvzu+p4mVvCja+uscigETgLvaa399tr2Ilbwo2sLPAH8F+DoYb9Hr20gJ/zeDI6+GdOvxCTgD8Dr3g7KDX0a7N0P+N21VUrFAP8A7tNa13bf7eIjPru+Z4jVr66t1tqutZ6MGY96hlJqQrcifnNtexGr31xbpdQC4KTWetPpirnY1u9rG8gJ/4yDo2uta9v/xNNm1K0wpVTKwIXYJwE12Lu/XVulVBgmgS7XWv/TRRG/ub5nitXfrm07rXU18BFwWbddfnNt2/UUq59d23OAq5RShzBV0hcopf7arYxHr20gJ/wzDo6ulEpXSinn8gzMz1sx4JH2zkrgVudT+VlAjdb6uK+D6ok/XVtnHM8BO7XWv+mhmF9c397E6mfXNlUpleBctgIXAbu6FfOXa3vGWP3p2mqtH9RaZ2utczH5a7XW+uZuxTx6bd0a09aXtNY2pVT74OgWYKnWeodS6h7n/qeBRcC9Sikb0AQs1s5H3wNNKfUypoVAilKqGPgJ5qFSe6yrME/k9wGNwO2+iLNdL+L1m2uLuVO6BdjmrL8F+BGQA353fXsTqz9d2wzgBaWUBZMcX9Nav9Xt35m/XNvexOpP19Ylb15b6VpBCCGCRCBX6QghhOgDSfhCCBEkJOELIUSQkIQvhBBBQhK+EEIEiYBtlinEQFNKJQMfOFfTATtQ5lyf4ezTSQi/Jc0yhegHpdRPgXqt9eO+jkWI3pIqHSGECBKS8IUQIkhIwhdCiCAhCV8IIYKEJHwhhAgSkvCFECJISLNMIYQIEnKHL4QQQUISvhBCBAlJ+EIIESQk4QshRJCQhC+EEEFCEr4QQgQJSfhCCBEk/j9t69wWQvqInwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from pylab import *\n", "%matplotlib inline\n", "\n", "plot(Te, Thm[:,0], label='E(T)')\n", "plot(Te, Thm[:,1], label='cv(T)')\n", "plot(Te, Thm[:,2], label='Entropy(T)')\n", "xlabel('T')\n", "legend(loc='best')\n", "show()" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Entropy(T=0.500)=0.000932 Entropy(T=4.000)=0.894 * log(2)\n" ] } ], "source": [ "print('Entropy(T=%5.3f)=%8.6f Entropy(T=%5.3f)=%5.3f * log(2)' % (Te[0],Thm[0,2],Te[-1],Thm[-1,2]/log(2)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework\n", "\n", "* Estimate the error of the density of states $log(g(E))==S(E)$ for 32x32 Ising model. \n", "\n", "To do that, you can use five independent runs of WL agorithm, which will produce $S_i(E)$ for i=1..5. \n", "In each WL run use $10^9$ Monte Carlo steps. From the series of $S_i(E)$ you can calculate \n", "\\begin{eqnarray}\n", "&& = \\frac{1}{n} \\sum_{i=1}^n S_i(E) \\\\\n", "&& \\sigma_E = \\sqrt{\\frac{-^2}{n}} = \\sqrt{\\frac{1}{n^2} \\sum_{i=1}^n (S_i(E)-)^2}\n", "\\end{eqnarray}\n", "\n", "Plot $S_i(E)-$ and $\\sigma_E$. \n", "\n", "* According to the publication \"Exact Distribution of Energies in the Two-Dimensional Ising Model\", PRL 76, 78 (1995), the first 10 values for the density of states in 32x32 Ising model is:\n", "$g_{Exact}=[2,2048,4096,1057792,4218880,371621888,2191790080,100903637504,768629792768,22748079183872]$\n", "\n", "Plot $$ and $log(g_{exact})$ for the first 10 coefficients, as well as their difference.\n", "\n", "Verify if your estimated $$ is approximately $\\sigma_E$ from exact result. To do that, plot $-S_{exact}$ and $\\sigma_E$ for the first 10 energies.\n", "\n", "* Calculate thermodynamics from your best estimation of $g(E)$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }