{ "cells": [ { "cell_type": "markdown", "source": [ "# Custom potential" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "We solve the 1D Gross-Pitaevskii equation with a custom potential.\n", "This is similar to Gross-Pitaevskii equation in one dimension and we\n", "show how to define local potentials attached to atoms, which allows for\n", "instance to compute forces." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using DFTK\n", "using LinearAlgebra" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "First, we define a new element which represents a nucleus generating a custom\n", "potential" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "struct ElementCustomPotential <: DFTK.Element\n", " pot_real::Function # Real potential\n", " pot_fourier::Function # Fourier potential\n", "end" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "We need to extend two methods to access the real and Fourier forms of\n", "the potential during the computations performed by DFTK" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function DFTK.local_potential_fourier(el::ElementCustomPotential, q::Real)\n", " return el.pot_fourier(q)\n", "end\n", "function DFTK.local_potential_real(el::ElementCustomPotential, r::Real)\n", " return el.pot_real(r)\n", "end" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "We set up the lattice. For a 1D case we supply two zero lattice vectors" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "a = 10\n", "lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "In this example, we want to generate two Gaussian potentials generated by\n", "two nuclei localized at positions ``x_1`` and ``x_2``, that are expressed in\n", "``[0,1)`` in fractional coordinates. ``|x_1 - x_2|`` should be different from\n", "``0.5`` to break symmetry and get nonzero forces." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "x1 = 0.2\n", "x2 = 0.8;" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "We define the width of the Gaussian potential generated by one nucleus" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "L = 0.5;" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "We set the potential in its real and Fourier forms" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "pot_real(x) = exp(-(x/L)^2)\n", "pot_fourier(q::T) where {T <: Real} = exp(- (q*L)^2 / 4);" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "And finally we build the elements and set their positions in the `atoms`\n", "array. Note that in this example `pot_real` is not required as all applications\n", "of local potentials are done in the Fourier space." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "nucleus = ElementCustomPotential(pot_real, pot_fourier)\n", "atoms = [nucleus => [x1*[1,0,0], x2*[1,0,0]]];" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "Setup the Gross-Pitaevskii model" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "C = 1.0\n", "α = 2;\n", "n_electrons = 1 # Increase this for fun\n", "terms = [Kinetic(),\n", " AtomicLocal(),\n", " PowerNonlinearity(C, α),\n", "]\n", "model = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,\n", " spin_polarization=:spinless); # use \"spinless electrons\"" ], "metadata": {}, "execution_count": 9 }, { "cell_type": "markdown", "source": [ "We discretize using a moderate Ecut and run a SCF algorithm to compute forces\n", "afterwards. As there is no ionic charge associated to `nucleus` we have to specify\n", "a starting density and we choose to start from a zero density." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n Energy Eₙ-Eₙ₋₁ ρout-ρin Diag\n", "--- --------------- --------- -------- ----\n", " 1 +0.293977096351 NaN 4.29e-01 8.0 \n", " 2 +0.442071143377 1.48e-01 5.35e-01 1.0 \n", " 3 +0.260880214052 -1.81e-01 1.72e-01 2.0 \n", " 4 +0.243058294552 -1.78e-02 2.70e-02 2.0 \n", " 5 +0.242660926714 -3.97e-04 5.98e-03 2.0 \n", " 6 +0.242635893634 -2.50e-05 1.85e-03 1.0 \n", " 7 +0.242633412289 -2.48e-06 3.10e-04 2.0 \n", " 8 +0.242633389813 -2.25e-08 1.45e-04 2.0 \n", " 9 +0.242633376491 -1.33e-08 1.54e-05 2.0 \n", " 10 +0.242633376330 -1.61e-10 6.30e-07 1.0 \n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown:\n Kinetic 0.0304520 \n AtomicLocal 0.0972467 \n PowerNonlinearity 0.1149346 \n\n total 0.242633376330 \n" }, "metadata": {}, "execution_count": 10 } ], "cell_type": "code", "source": [ "Ecut = 500\n", "basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\n", "ρ = zeros(eltype(basis), basis.fft_size..., 1)\n", "scfres = self_consistent_field(basis, tol=1e-8, ρ=ρ)\n", "scfres.energies" ], "metadata": {}, "execution_count": 10 }, { "cell_type": "markdown", "source": [ "Computing the forces can then be done as usual:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "2×1 Matrix{StaticArrays.SVector{3, Float64}}:\n [-0.038689357456715534, 0.0, 0.0]\n [0.03868888948444092, 0.0, 0.0]" }, "metadata": {}, "execution_count": 11 } ], "cell_type": "code", "source": [ "hcat(compute_forces(scfres)...)" ], "metadata": {}, "execution_count": 11 }, { "cell_type": "markdown", "source": [ "Extract the converged total local potential" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1" ], "metadata": {}, "execution_count": 12 }, { "cell_type": "markdown", "source": [ "Extract other quantities before plotting them" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=4}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ2BUZfYw8HP79Jn0TgIEpEuVIhCk2BCkWBBXXewKllWx7rruuq6rvur6t66dBRUVBQErKIqAKKAUCb2Elp5Mn7n9/XDZGEPJJJly78z5fZqZ3Jl5MjPPPfdp5yFUVQWEEEIoVZGJLgBCCCGUSBgIEUIIpTQMhAghhFIaBkKEEEIpDQMhQgihlIaBECGEUErDQIgQQiilYSBECCGU0jAQIoQQSmkYCBFCCKU0HQXCLVu2SJIU4cGKomByuMRSFCXRRUhpsiwnuggpDU9BCRfFU5COAuHYsWO9Xm+EB/M8jyeCxAoGg4kuQkrDzz+x8BSUcFGsAjoKhAghhFD8YSBECCGU0jAQIoQQSmkYCBFCCKU0DIQIIYRSGgZChBBCKQ0DIUIIoZSGgRAhhFBKw0CYKkI1fP2vkeYrQAjFk2dfwFeBGRIShk50AVA8qIq6c/5hKSBXr2/sOj2fS2MSXSKEEACAFJAPLK1y7/GrKgy8r5Q2UYkuUSrCQJgSjn5bz9jo/nd2Pfpt3eZn9w3+c3eKxc4AhFonSdK1114bDodbPK4oCkEQBEF08PV9B4O0lTJlsKFaAb4BSw7XwRdMYhMnTrzmmmti8coYCJNfuEE4+m3dmXd0ISiicFyWZ3+wYZs3a5Ar0eVCyAACgcBHH3309ttvJ7ogqW716tVff/01BkLUTgeWVBWek2nKYLW72YNcNRsbMRAiFCGGYS699NJElyLVhcPhFStWxOjFsX8sySmC4t7jzx2R3vRIRh+7ryIk+iLd8QohhJIbBsIk59kXsBWaKe63L5pkyfTe9trNngSWCiGE9AMDYZJz7/a7zrC1eDBroKt2kzsh5UEIIb3BQJjkGncF0rq3DISu7lbeLYZq+IQUCSGEdAUDYTITvJLgFa2FphaPEySR3tvRuMufkFIhhJCuYCBMZu5dflc3G0GeZKmTvZPZfzgU/yIhhJDeYCBMZo0nGyDU2DuZfYcwECKEEAbCJKaCZ4/f1d160j+aczjBI0ohOc6FQgghvcFAmLSCNTzJkKZ09qR/JUjCWmD2H2mZOAohlKx27tw5c+bMprv//e9/P/roo9Mcv2jRovnz58e+XImHgTBpBY6GbEXm0xxgLzL7D2PCe4RShc/n27Bhg3bb6/U+/PDDY8eOPc3x48aNe/jhh30+X1xKl0gYCJNWoJK35rWcL9qcrZPZj8OECBnNkSNHtm3bVllZuWDBgv379wNATU3NkiVLPv30U7//+FRwRVE2bdo0f/78zz77LBAInPgiCxYsKCsrS0tLA4CNGzfW1NTIsrxy5UpVVauqqn755RcASEtLGzVq1DvvvBPHfy4xMBAmrWBl2Jp/ukBo72T24cRRhIxm6dKlV1111eTJk9esWXPw4MEvv/xyyJAhK1aseP/99wcNGlRVVQUA8+bNe+SRR7Zs2fL222+feeaZDQ0NLV5kyZIl5513nnb7rrvuWrt2bSgUmjBhgqIo33zzzf3336/9acKECUuWLInnf5cQmHQ7aQWOhS15p9vSxZTOyoIi+CTWjj8DhNrgowPKV0fV+LzX2Hzi8i4tWyzHjh3bt2+f3W6XJKlz587vv//+iBEjAODBBx986qmnnn766VmzZs2aNUs7+Prrr583b96f/vSn5q+wadOmJ598stV379Onz6ZNm6L0r+gXngGTkxSSpZBsSjv5TJnjCG2YMJTeyx6vciGUDAqtxKDMOL1XkfUk64CHDRtmt9sB4MCBA5WVlfPmzZs3b552V5ZlAKitrf3b3/62YcOGuro6t9ttNv9uuoCqqm632+l0tvruTqezoaFBVdWO77yoZxgIk1OwMmzJM0FrP10bBkKE2m5oNjE0O5GBwWo9vixKlmWGYS655JKmQKWFt9mzZxcVFX3++efp6ekPP/xwTU1N86cTBOF0Or1er3aXpmlJ+m07GlEUGYbRbnu9XpfLldxREHCMMFm1OlNGowXCOJQHIRQLXbt2TUtLU1V1/P8MHDgQAHbs2DFp0qT09HRFUb788ssTn9i/f/8dO3Zot7t37/7jjz82/Wn9+vXdu3fXbm/fvl17weSGLcLkFDgWjiQQWnK4YDWm3kbIqBiGeeutt2bNmnXBBRekp6fv2LFj0KBBjzzyyNSpU2+++eZLLrlkzZo1FEWd+MSLL7545cqVM2bMAID7779/7NixWlycPHny3r17v/nmG+2wb775ZvLkyfH8jxICA2FyClaGswa2PgBgymAFj6hIKkknedcHQklj8uTJo0aNarp73nnnbdmyZe3atW63e8qUKUOHDgWAv//972PHjt2/f/+VV15pNptDoRAA9OjR491339WedfXVV/fv39/n89nt9pKSkvLy8pUrV3722Wd33333qFGjtK5Rr9e7cuXKp59+OhH/ZVxhIExGKgSqImoREiTBpbPhesGSc7r5pQgh/SgsLCwsLGz+SGZm5sUXX9zisDFjxowZM6b5I3a7fciQIdrttLS0+++//9NPP9UahSaTSTu4rKysqQX56aef3n///S6XKyb/hp5gIExC4UaBNlG05ST9IScyZ7GhGh4DIUKp5pZbbml+l6Ko8ePHN58Xc8UVV8S9UImBgTAJBSt5SwTNQY05m8MdehFCZrN5xYoViS5FYuCs0SQUqIyoX1RjyeZCNUJMy4MQQnqGgTAJhWp4c85pl9I3Y85ig7XYIkQIpS4MhEkoXC+YMyIOhNg1ipChfP311y+99FKMXnzDhg2nz7L9008/JV8abgyESShUJ5gyI538wthoABADuEMvQsbA8/xJN5SIitmzZ/fq1es0B/Tp0+evf/1ri1Q1RoeBMNnIvCLzSpvyaGsTR2NXJIRQFI0ePfraa68FAEEQKioqFEXZuHHjli1btL/u3Llzw4YNzVOmHTlyZPXq1bt27Wr+IpIk/fLLL9u2bVMUZf/+/aqqAsDq1asBYMCAAQBQU1Pj8XgAoKKiQhTFcDh8+PBhALBYLJMnT37ttdfi9N/GBc4aTTbhesGUzraaZbQ5cxYXquEdnS0xKxRCKGr++9//rlmz5t133927d29ZWdnw4cNVVd28efPMmTPD4fCOHTsqKyuLi4uXL18OAP/v//2/999/v0uXLjt27OjUqdPSpUtJkvR6veedd56qqtnZ2RzHLVq0SBRFmqYXLVp0wQUXaO9yzz33DBw48M477xw8ePD69eurq6tnzZqlRdOJEyfefffdDz30UCI/hajCQJhswvWCOTPSAUKNOZsL1eLEUYQixe/ZIhza1fpx0cAWdee69z/VX+vq6u65557Ro0fv3r27V69eTz/99HPPPcfzfGFhYXl5ea9evebMmXPPPfcAgKqqI0eOXL58+eTJk1944YWsrKxPPvmEIIh///vfixYt0l7txx9/vPfee1stUr9+/bZt2xYKhVpsamFcGAiTTahOMEU8U0ZjzuZqN7ljVB6Eko8qCkrQH6/3Cp/mr06nc/To0QDQvXt3kiQvuugiAOA4rmvXrocOHerVq1cwGHziiSc2b97s9/sPHTq0c+fOyZMnr1u3rmnDimnTpjVtVVhXV6ftWX96Wi7vurq6oqKiKPyHOoCBMNmE64TTb0x/IhwjRKhNTL2GmHoNSXQpAACat8koimq6S1GUNkw4c+bMHj16/OMf/zCbzQ888EA4HAYAWZabMsg0T8lts9mCwaB2myRJRVGa/iTLctOR2lQdbUPE5ICTZZJNuL7tLcIsLlwvqEqcdtxGCMXN+vXr77jjjt69excXF//666/ag4MHD27am+mLL75oOrh379579uzRbpeUlDQdDwBbt24tKSnRbu/evbuwsDCZcpBiizDZtGOMkKQJ1sHwjWJbIyhCSOfKyspmz5598cUXL1u2TNu8HgDuuOOOUaNGTZo0KTc399ChQyRJag3EiRMnfvDBB1pP6U033TR06NDrrrsuFAo99thjn3/++cKFC7Wnf/vttxMnTkzUfxQLGAiTiiqrgkfk0pi2PpFLZ8INGAgRMoDx48f369cPAAoLC1999dWmx995552mEb5HH320Z8+eAPDee+/Nnz+/srLyz3/+M8dxHMcBQGZm5saNG9esWUOSZFZW1kUXXaR1e06fPv2BBx6oqanJzs7Oy8vbvn370qVLFy5c2L1798ceeywvLw8AVFVdsGDB22+/Hff/O4YwECaVcIPIOhmCavPmglwawzcKANZYlAohFEXdu3fXdpB3OByTJk1qenzatGlNt8eOHavdsFgsN91004kvsmzZsrKyMrfbfeedd1566aXagyaT6eGHH3755Zf/+te/AoDdbr/yyiu1A7QoCABff/31wIED+/c/5URWI8JAmFTaMUCoMaWzfKMY9fIghPRpxYoVzz33nMlkOuecc5ovmbj++utbHHnmmWdq7UjN+PHjx48fH6dSxgsGwqQSrhNMbRwg1HBpjGdfrJI2IYT05o033ojwyJUrV8a0JHqAs0aTSriejzzddnOmdJZvwBYhQigVYSBMKqF2twjTmXADJpdBCKUiDIRJJVzf3kDoYkSfpMq4lBAhlHIwECYVvkE0pbcnEBIkwdhp3oO9owihlIOBMHmIAZlgCIpr53eKw4QIGUgwGDxw4IAoYp2NAgyEyYNvFDhXm5fSN+HSmHAjDhMiZABLly7t1q3b5ZdfXlBQsHjx4kQXx/Bw+UTy4BvFjgRCbBEiFKF97oOHvUfj816F9rzStC7NH2loaLj66qvXrl3bu3fvdevWTZw48ZxzzkmmzJ/xh4EwefBukUtrf440Lo3x7g9GsTwIJat9jQfXHvkxPu81rGBwi0C4atWqgQMH9u7dGwBGjBjRqVOn9evXn3/++fEpT1LCQJg8eHeHWoRcOhvGXQkRisC5ncec23lMot69trbW4XA03XU6nR6PJ1GFSQ44Rpg8+Mb2pNtuYkpnsGsUIf1TFGXr1q3aZoHhcLi8vPyMM85IdKGMDQNh8uDdHQqEnIsRvCLuSoiQ/gWDweuvv3758uVXX3314MGDkywFdvxhIEweHZwsQ1AEY6MFjxTFIiGEYuGcc86ZOHHi4sWLBw4ciLNGOw7HCJOEKqtSQGIdHfpCuXQ23CB0pFmJEIqP6dOnT58+PdGlSBLYIkwSvEdk7DRBtnknwuZwmBAh/eM4zmrFrUOjKdIGRHV19caNG3fs2DFs2LCRI0ee9Jjt27e//fbbkiRdeeWVgwcP1h6UZXnevHkbNmzo1q3bzTffbLFYolNw9HtCx6aManBNPUL6d91111133XWJLkVSibRFeNNNN/3zn/988cUXT7U31e7du0eMGOFwOAoLC8eNG7dx40bt8bvvvvull14666yzVq5c2XwDZRRd4UaRa1eW0eY4FyO4sUWIEEotkbYIlyxZAgAzZ8481QHPP//8jBkz/vKXvwBAY2PjM8888+677zY0NLz66qtbt24tLS2dMWNGXl7ezz//PHDgwKgUHTXXwZkyGs7F1P/qi0p5EELIKKI2RrhmzZpx48Zpt8eNG/f9998DwKZNm7Kzs0tLSwHAbDaPGDFizZo10XpH1FwHV9NrWCe2CBFKchUVFcuWLYvKS73++uvhcDgqL5VYUZs1WllZmZmZqd3Ozs6urq5WFKWqqiorK6vpmKysrMrKylO9giiK11xzDcv+rn/vlVdeOemwcCgUYhiGpnHW63HBurC1KxsMdihHmsIp4UYhwhcJBoMkiZOtEiYUClEUlehSJL9QKJToIpzEk08+WVpaevrBpnvvvXfixIllZWUtHt+yZctTTz01adKkjhfjlltuufjii00mU/ueXlFRce+9977//vsRHi/LcvOzU4SnIJPJ1OphUQskLMtK0vElaIIgMAxDkmTzB7XHOY471SuQJDllyhSbzdb8QZvN1iI0ahRFwUDYnOSVbVmW03y8EeEAFKCBiWQvJ1EUO/p2qANOX5tQtJz0/JNw27dvb7VgP//8c9OkRX0KBAKrV6+O/HiSJJv/5iM8BRFE63PpoxZICgoKjhw5ot0+cuRIQUGB9uDRo0dVVdWKcuTIkbFjx57qFSiKmjp1anp6eiRvR/1PNMqeDAS3ZM7gOv6BsE5G9iuspfVeVvz8Ews///jQ4Ye8aNGiTz/9dNWqVQsWLLjoooseeeSR77//fvbs2bW1tU6n89lnn73gggteeumln3766cCBA08++eQf//jHOXPmnPSlvv7669tvv72hoSE9Pf3555/Xzs/V1dW33XabNrw1ePDgZcuWHT16dPr06UeOHJEkacSIEa+//vppTtR9+vS55pprXnrppUAgMGPGjGeeeYamaUmSHnroofnz5wPAyJEjX3755YyMjFmzZtXX12vR+osvvmjqUzwVgiCafx1RrAIdCoRut3vdunUXXnghAEyZMmXhwoWzZs0iCGLhwoVTpkwBgKFDh1IUtXLlygkTJhw8eHDTpk2Rt4JR5OSwoqoqbY7Cb4Jz0bxHNGdjUwOhU+LdouiPUw4mxka3GP6/5JJLli1bNmDAgDvvvBMAGhsbp06d+uabb06ePPm7776bNGlSeXn5rbfe+vHHH994442XXXbZqV65pqZm+vTpH3744YQJE7744otp06bt2bMnKyvrD3/4Q7du3Q4fPkzTdHl5OQBYrdYFCxaUlpZKkjR79uzHHnvs6aefPtXLHjhwYP369Xv27AkGg2PGjHnzzTdvvPHGV1555auvviovL7fb7VddddWdd945f/78t956q/kSgwSKNBA+/fTT77333oEDB7799tvly5fff//9l1xySXl5+cSJE1VVBYCbbrppwYIFY8aMMZlMe/fu1SbFMAzz5JNPzpw5c8KECd9///3cuXPz8/Nj+N+kqnBj1NLBsE5GcGOWNYROp36rt2ZjnLZqyRrgLDjndE2lb7/9tqioaPLkyQBQVlZ21llnff755zfccEOrr7xixYo+ffpMmDABAM4///yePXuuWLFi3Lhxq1atWrJkiTbw1KtXLwBwuVz19fWvvPKKx+NRFOXHH1vZguqee+6hadrhcMyePXvRokU33njjxx9/fNttt2mbJj700ENnnXXWf//734g/g5iLNBBedtllY8aMabrbqVMnAOjfv//WrVu1R1wu18aNG1evXi1JUllZWdPC+auuumrkyJG//PLLQw89pG2ghaIuKqvpNZyL4T04cRSh08kfnZE/OiPRpTiurq4uJyen6W5ubm5tbW0kT6ytrT3xiVVVVQ6Ho8UUxW+++eaaa66ZPXt2YWFhMBhsNRA2zZHMzs7WCtP8vXJzc4PBYCAQiKSQ8RFpICwqKioqKmrxoMVi6du3b9NdjuO0i4sWOnfu3Llz53YXEbWK90jRCoSsiwkcS4b50AglMZqmtW2YAKCkpGTnzp2KomhzI7dv3z5x4kTtGFmWT/MiJSUlr776qjaHQ1XV8vLyWbNmde7c2efzVVRUFBcXNx354Ycfzpkz57777gOAF198sdXi7dq1S1s1t3PnTu3k37lz5x07dmgF2759e1ZWls1ma7WEcYOzLpOB4BFZZ5RahE6msRzX1COkaz169Fi2bNkZZ5xRUlJyzjnn2O32e+65549//ONHH33U0NCgTdHo0aPHokWLXC5Xt27dtLDUwoUXXjh37tyHHnpoxowZ77zzjqqq559/Psuys2fPnjlz5r/+9S+73f7TTz/deOONJSUly5Ytu+CCCyoqKp5//vlWJ6w+9thj6enpHo/nmWeeeeeddwBgzpw5V1999RlnnJGZmXnnnXfedtttAFBYWMjz/Ouvv15SUjJq1KgEzoLGdWDJIIqBkMWuUYR0b86cOZMnT/700083bNhA0/SqVasAYO7cufX19evWrTObzQDwyCOPDBkyZOnSpdu2bWv+3OLiYm0RIcuya9as8fv9c+fODYfD33//vRbhnnnmmauuuurpp59+8MEHPR4PANxxxx2jR4++6667li5d+tprr11yySXaS11//fXae7Xwl7/85dlnn3355ZdffvllbSbqueeeO2/evAULFjz++OM33XTTn//8ZwCwWCyfffbZjh07Pvzww8QuzCe0qS56kJGRsWfPngiXT+CC+ua2v1aRPzI9rae94y8l+qWfn9gz9NGerR7p9/tbLPpE8eTz+ez2KHzj6PQ8Hk9xcbHbHaepMUnAarXu2rWrsLAwui87f/78FStWNJ9iE8VTELYIk0EUW4SMlZYFVRGVqLwaQgjpH7aokgEfvUAIBLAOWvBIpkw9JtRACOnc0aNHHQ5HokvRNtgiNDxFUhVeYSxRy3/BuRgeU28jhNrF5XIZLguxwYqLTiR4RNbBQIe2pv8dnC+DEEopGAgNT/CIrDOaXdwcbsaEEEolGAgNj/dIURsgBID/pRuN4gsihJCeYSA0PMEjclENhJhuFCGUUjAQGp7glaLcNYqTZRBCqQSXTxge7xbtnU6S3KHdWCcjYNcoQv+jqmpjY2OiS5Hqmu9NH3UYCA0viqvpNaydFoOyqqgEGb2pqAgZk9lsLisr69q1azueq4gKSZNRnNENKiiSQjIp2pN3qu2FOw4DoeEJnih3jQIBjI0WvFHb0QIh42JZdunSpSc+3nqWRxXW3bd92D97kXTUIqGqqOvuKx/xRC+8SI2uFL2ySB4qCD6JdUQ5YnFOWvDifBmE2k8MSJSJimIUBACCJBgrJfqwbkYZBkJjE4MyxRLRrWyAw4QIdZjgkVhH9LvcWBfDezAQRhkGQmOL+gChRks3GvWXRSh1CN6Y1E0OL1JjAAOhscUqEDoZwYuVDaH242PUInTSGAijDgOhsfEeiYvuTBkAwBYhQh0meKOc6ULDObFrNPowEBpbzLpGGR5bhAh1QKzGCLFrNAYwEBpbzLpGsUWIUIfEaIwQu0ZjAQOhsQnemFx1cjhGiFDHxKhusg7sGo0+DITGFs296ZuhLZQiqoqoRP2VEUoRgjtmLUK8SI02DITGJnglLgZXnaDNl8E19Qi1i6qoYlBmbFTUX5k2UaCAzONFajRhIDQwVVGloExbo1/Z4PjEUbzwRKg9BK/E2OgYJUJjHDQml4kuDIQGJngl1h6rysY6GWwRItQ+Ud8crTnWweDW2dGFgdDABK/ExKZfFAA4J+5Tj1A7RX277OZw2CLqMBAamOARuWin227COLBFiFA7xWgRoYbD+TLRhoHQwARfLCsbJpdBqL0Erxj1PWGaMA5GxIvUqMJAaGAxWqikwXSjCLVbjMcIaR4DYVRhIDQwwSMyMbvqxAQWCLUb74lhi5B1MCJepEYVBkIDE7wxybitYTG3L0LtJXhi3CLEuhlVGAgNLKbjEBRLEiRIYTlGr49QEovR1hMaTC4TdRgIDSymM9Pg+BageOGJUNsokqoICm2OSaYLwOQyMYCB0KhURZVCMmOLYSBkHXjhiVCbHe+qiUmii+MYXEoYVRgIjUpLKxPTysY6sEWIUJvFdIBQwzpwUnc0YSA0qpiundCwTrzqRKjNYjp4r8HkMtGFgdCohFjOz9bgVSdC7RDrwXvA5DLRhoHQqOLRIsSrToTaLg51E5PLRBcGQqMSvCIGQoR0SPDGZEve5jC5THRhIDQqwSvFvrIxmFwGobaKS28NJpeJJgyERhWnrlHc/xOhNorPZBlMLhNFGAiNKg6VjWRJkiQwuQxCbRKHyTKYXCa6MBAaVRzWKoG2bhcvPBGKmCKpihjDtDIaTC4TXRgIDUmVVSkkM9aYB0IOlxIi1BZxSCujYRy0iCMXUYKB0JAEX8zTymgYXEqIUFvEoV9UwzoYHueyRQkGQkMSvBITl8rG4QoKhNoiplvyNsdiizB6MBAakuAVuRjPlNEwDhpnaSMUuTjMYtPgMt8owkBoSHFrEbIO3J4XoTaIw7omDWZAjCIMhIYUx8qG3S8ItUFcW4R4kRolGAgNSYxfZcMBeYTaIH6TZXBGd/RgIDSkuLUIcfkEQm2CXaNGhIHQkARPzDNua0iWJEjA5DIIRQgnyxgRBkJDErxSfCobHE/vi/UNodYpkqoICm2JbVoZDW2mtLeLw3slPQyExqMqqhiUGVs8Khtgel+EIha3tDIa1o5p8aMDA6HxiH6ZsVIEGafaxjpxwxeEIhK3AUINzpeJFgyExhO3QQgNbgGKUITiNmVUwzpwD4rowEBoPPG+6sTkMghFRvSJ8cl0oWEdDLYIowIDofHEbcqoBisbQhHiPRLnjGtvDdbNqGjD+bSysnLx4sUkSU6bNi07O7vFX7dv315eXt78kalTp9I0vX79+sOHD2uPsCx78cUXd7DEKJ5TRgErG0IRE7yiOdsat7djHUywxh+3t0tikQbC/fv3Dx06dOrUqbIsP/LIIxs3biwsLGx+wPbt2xctWqTd3rt377Fjx6ZNmwYAzz333O7du7t27QoAVqsVA2HHCV7JVmCK29thchmEIiR4pPhkw9dglrVoiTQQPvvss1OnTn311VcB4KqrrnrhhRf+9a9/NT/gsssuu+yyy7Tb06dPHz9+PEUdn99/3XXX3XrrrdErc6oTvCLbyx63t8PkMghFSPDGddiCwd6aKIl0jPDLL7+cNGmSdnvSpElffPHFqY6sq6v79NNPZ82a1fTI9u3bFyxY8NNPP3WkoKiJ4JVYe/wqG8mSJElgchmEWiV4JTaOY4ScgxF82FsTBZGeT48dO5aXl6fdzs/PP3bs2KmOnDdv3uDBg3v27KnddTqdlZWVX3311b333jt48ODFixc3tRRbkCTpmWeeMZvNzR+cM2eOyXSSbkCe5xVFkeVUPDvzHlE1KTzPx+0dGQflrw2as9nfFYPnGSZ+dR61wPM8y7KtH4di48RTkJbnRSYlmY9XK40ChVdCgTBJx2sNv55EeApiWZYgWvl8Ig2EBEGoqqrdVlX1NK/71ltv3XXXXU13X3nlFe2Gx+Pp16/fe++994c//OFUz3W73S3O701vio5TQdZ3A3EAACAASURBVI5jWhkNY6dFn9QiECKEmhO9EmOn45ZWBgCAANpGSX6ZdcWviygpRfrx5ebmVldXa7erqqqaWoctrF+/vqKi4tJLLz3xT06nc+TIkS1mlv6uKDT997//PT09PZLyKIrCMAxNp9zXL/okykSZLPGbLAMAJherBgmO435XElFs8QiKJ0EQ8PNPoBNPQXxY5lxMnL8UzslCiOByUvGXEMVTUKRjhOedd97y5cu128uXLz/33HO123v27AmFQk2Hvfnmm5dffrnd/ttUjqYmHc/zGzZs0KaPonYTvBLrjHf4Z52YwAKhVsQ55ZMGVzdFRaSn1DvvvHPYsGEEQciy/MUXX2zcuFF7vHfv3itWrCgrKwOAQCDw/vvvf/75503PEgShX79+48ePN5lMX3zxRWZm5mn6RVEkElTZmHCDEOc3RchYBE+iLlIxEHZUpC3C0tLSrVu39u7du3///lu3bi0qKtIe/+CDD3r37q3dDgaDb7311ogRI5qexTDMa6+91qNHj4KCgieffHLNmjXYmdNBcZ4yqsHlSgi1KkF1E7fnjYI2fG35+flz5sxp8eCUKVOabmdlZWmL6JsQBDFq1KhRo0Z1pIioOcEjJuKqEysbQq0QvKIl1xbnN2XttPdAMM5vmnww16jBxDm/mgbHIRBqVZy3ntDgRWpUYCA0mDhvPaFhnYzgEQFXsiB0aoJXjOdqeg0OW0QFBkKD4RNR2UiaIBlSCqVi+gKEIpSoFiGPLcIOw0BoMAmpbKBNTsPU2widgiIoiqzS5rhmugAAxkIpvKJI2F3TIRgIDUUF0Z+gQIi7EiJ0agkZswAAIHDiaBRgIDQSwScxFoogE5BXkHXS2AOD0KnwngSMWWhwLlvHYSA0koSMxmtYByNiZUPoFBLWIjw+bIF1s0MwEBpJogYIAa86ETotwStycV/gq8Gu0Y7DQGgkiWwRaisoEEInI3olJu4LfDV4kdpxGAiNJJHdLw6ax8qG0CnwXolLWN3Ei9SOwkBoJIkNhFjZEDoVwSMmrEWIebc7DAOhkQieBGw9oWEdjOiTMLkMQieV0ItUbBF2FAZCI0nIZoQakiYojhQDeOGJ0EkIHpFzYYvQqDAQGomQuLVKcHy+DNY3hFqSwwoQQHGJOZ3SZkqRVVlQEvLuyQEDoWGoiiqFZMYa7xxOTVgnw2MPDEIn4D0il7grVABg7TQu8+0IDISGIXglxkYnJK2MhsN0owidTGK7agA3Y+owDISGkcDReA12jSJ0UoInYYP3GlxK2EEYCA1DD1ed2DWK0IkS3jXKOWmsmx2BgdAwBK/E2hN51YldowidVMIvUhlMBdwxGAgNQ/CK2DWKkA4lvGuUw8RPHYOB0DAEr5TYq04Ou0YROpmEd41iKuAOwkBoGIInwS1C2kopoqqIuFwJod8R3Ikev8fJMh2DgdAwEt4iBKxvCJ1AVVQxKDO2hC3wBWwRdhgGQsPg3QnufgFt4qgb6xtCv9FmsSVwgS/A8aQ2chh7a9oJA6ExKKKiiCptSeRVJxyfOIotQoR+I3jExM6U0XAuHMJvPwyExnB8b/pEXnQCYA8MQifgPYkfswCsmx2DgdAYeHfCcts3x+K6XYR+T0j0lFENh8MWHYCB0BgSvmJXw+FSQoR+TyddoywOW3QABkJj4D0Sp4vKht0vCP2OXrpGcYywAzAQGoNOWoTYNYpQC/rpGsWL1HbDQGgMvEdkdTBGyDkZ0SeBmuhyIKQbOuoaxTHC9sJAaAyCR9RD1yhBEZSJEgM4FIHQcYJHYh26uEjlcYywvTAQGgPv1sU4BBzvHcX6hhAAgBSSCZLQ1rMnFmOj5bCsSNhd0x6J//5Qq1RFFf0J3pW3CQ5FINQk4ftO/IYAxk7jPvXtg4HQAES/zFipxOZwasI6GRyKQEijk1lsGs6Fq5vaCQOhAfBuXcyU0XBpuG4XoePCjSKXppe6iaub2g0DoQEIOki33YRzYSBE6DhBHymfNJhcpt0wEBoAr7PuF6xsCGl0kvtQw7pobBG2DwZCAxD0kVZGw7kYvhErG0IAWiDUTdcorqBoNwyEBiDoYzW9hnUxgkfENfUIgd5ahDhG2F4YCA2A10cOJw1JE5SJEv144YkQ8G4dDVuwThqHLdoHA6EB6CSHUxMcJkQIAKSATNK6WE2vwQyI7aaXrxCdhqCP9PZNcAUFQqAN3uumXxSaMiBib03bYSDUOykoExRBsTr6prBFiBDobIBQw+FmTO2io9MrOim99YuCtvMZBkKU8gSvpJ9ZbBrcnrd9MBDqXbhR5NLYRJfid3AFBUIAwDeKJt2sndBwaQzfKCS6FMaDgVDv9Nn9gulGERI8umsR4rBF+2Ag1DteT8kMNVjZEAKd5VfTcC4We2vaAQOh3umwRcg6adEvqQpO00YpjW/UXd3k0nDYoj0wEOqdDluEBEnQVlr0yYkuCEKJo4Lo09e6JsDemvbCQKh3fKOgt0AIAJyLkbwYCFHqEgMyZaZIWhe7hDY53lsjY29N22Ag1DVVUUWfpJ/8ak1wC1CU4gS3bvamb4YgCcZO41LCtsJAqGuCV6KtNEHp66oTjk8cxUCIUpfe0so0wUnd7YCBUNd0OFNGw6UxInaNohQmeiTWpbsWIQBwaThxtM0wEOqaDudna1jsGkWpTfDoccwCMBVwu2Ag1LVwg8il67GymVyM6MEWIUpdglvUaYsQEz+1HQZCXdNv12gGdr+glMY3SKZ0feU+1GCLsB0wEOqabrtGGSulyqoUxkYhSlG8W2TTsEWYJDAQ6po+FxFqWBfDN2B9Q6lI5hWQgLHpMhCmMWHMu91Gbfgia2trf/jhh/z8/EGDBhFEywn91dXVR44cabrbp08fjuO02xUVFZs3b+7WrVuvXr06XuKUEnbrbuuJJqyL4htFa74p0QVBKN7CDQLjohJdipOjzRQASGGZNum0hDoUaYtwzZo1vXr1mjdv3pVXXnnVVVedeMA777xz3nnn3fQ/NTU12uPvvvvu4MGD33333XPPPfdvf/tb1AqeAhRBUQSVsej018ylM+F6vPBEqYhv0F3iw+awd7StIm0RPvDAA3/5y19uv/12j8fTo0ePdevWjRgxosUxkyZNeuutt5o/Iori3Llz58+ff/755+/fv79Pnz433nhjXl5edMqe7MJallHdLaY/jnFSuPMZSk3hBoHVcyBMY/lG0ZqHvTWRiqhFWFtbu2bNmpkzZwKA0+m88MILFy9efOJhfr9//fr1FRUVTY/8+OOPgiCce+65ANClS5dBgwYtX748SiVPfjrMbd8cl8aEcYwQpSS+QeB0OVNGg8ll2iqi7/LIkSMmkykzM1O726lTp127dp142M8//zx37tzy8vKhQ4d++OGHVqv16NGjhYWFJHk83BYVFR09evRU7yLL8uLFi202W/MHJ0+ezLInGSSTZZkkyROHKpNJqIFnXbQs63RmJuUgQvW8bouX9GRZxg8/UUL1giPPLMuyPk9BjJMKNwpJ//OIsApEEikiCoSCIDDMb00ThmHC4XCLY26++ea77roLAHw+37hx4/71r389+uijgiDQNH36JzZRFGXJkiUtwt7YsWOtVuuJB/M8ryhKcn/TofowZSN4nk90QU5Otch8vajb4iU9QRDww0+UcD3vtHI8r9MLQcpGePeEkv7nwfN888B0KiaTKTqBMDc31+/3h8Nhk8kEALW1tfn5+S2OsVgs2g273X7FFVd8/vnn2hPr6uqajqmtrR04cOCp3oVhmHnz5qWnp0dSJIIgGIZpHmWTj+xrcHazNn2weqOkKwTRyBKcNksNxZksy7r9bSQ9wS078+wmC6fPU5CYozZs8if9z0NRlGj9jxGNERYVFXXq1GnVqlUAoKrqqlWrzj77bO1PqnqSja927NihzYgZOHBgVVXV/v37ASAcDq9bt27kyJFRKXcqCNcL+kxd0YRLx2FClHKksKwqKmXW7yJsUwaLM7rbJKLLGZIk77nnnltvvfWvf/3r2rVrg8HgtGnTAGDdunVnn322Fguvv/760tLSrKys9evXv/fee2vXrgWAjIyMG264YcaMGXPmzFm0aNGQIUMGDRoU0/8nmYQbRFOGrgOhKZ3lGwRbAU5OQymEr9d7xeScjBSUFVEhGf1Ga12JtF0/Z86cvLy8r7/+Oj8/f82aNdpi+ZKSkieeeEI7YNq0ad9+++3BgweLi4vLy8s7deqkPf7ss8+++eaba9euHTVq1K233hqL/yEpKZIqBSTWoceOlyZcOhNuwAtPlFrCDfrN93QccXwFhTmbS3RRjIE4ad9mQmRkZOzZsyfCMcJQKJTcY4ShGr78jYpBD3RPdEFOye/3e3/mw/VCl6m4MDQBfD6f3W5PdClS0dHv6vhGMe88l55PQdtfPZg/KiOtZzL/Qvx+f4tVBu2GDWedCtcLOu9+geNLCbFFiFIL36jrtDIaUzobrsfx+0hhINSpcL3I6XumDACYMljMu41Sjf5nsQEAl8HiRWrkMBDqVLjBAJXNlI6VDaUcvlE0Rt3EiaMRw0CoU+EGA3SNUiaSoAgxoMc1xQjFSLhB4NJ13zWagcMWbYCBUKfC9YIpQ++VDbQVS3VY31CqEP0SQRH6TyKBSwnbBAOhTvG6X0SoMWexobokz+SEUJNQrWDONMCaBNpMEQQhBbG3JiIYCPVICsqggv6vOgHAnMmFa/HCE6WKUB1vzjLAFSoAmHDH0IhhINQjo/SLAoApE1uEKIWEawVTpkECIU4cjRgGQj0yxEwZjTmLDeEYIUoZoTpjdI2CtoICW4SRwUCoR+F6gTNKIMzkwjVY2VCqCNcaqGsU19RHCgOhHoXrDbBQSUNbKSABV1CgFBGqw67RJISBUI8M1DUKAOZMNlyLw4Qo+Qk+iWRIQ8xiAwBTOsNjIIwMBkI9CtcLJt2v2G1iyuJwmBClgnAtbzZIcxAAuHSWbxRVRS/bKugZBkLdURVVcBsg0WgTcyYbwhYhSgEh40wZBQCSJhgbzbtxmLB1GAh1J1wvsE6GpIlEFyRS5iw2hEsJUQoI1QnG2uHPnMWFcC5bBDAQ6k6o1mCVzZTJYZY1lAqM1TUKAOZs7K2JCAZC3QnVGGZ+tsachZUNpQQDTRnVmLM4rJuRwECoO6FawZxlpBYhbaYImhB9UqILglAsqRA2XCDMxq7RiGAg1J1QLW/ONlJlAwBzJk4cRUlO8EkkR9ImY6yd0GBvTYQwEOpOqIY3VosQcA8KlAJCRhsgBABTOiv6ZUVQEl0QvcNAqC8yr8hhhXMaZhGhBntgUNIz3Cw2AAACTBkM9ta0CgOhvoRqeVMWC4ZZOnGcJYcLVoUTXQqEYihYGbbkGi0Q4nyZyNCJLgD6HSP2iwKAJdcUrMLK1mZhGRp41S2AmwevCD5RbeQhKEFIBo+gCjIEJJAU8IkAAAqAR/gtS4gssxT12wQlG00wJACAlQGWBBMFZhrsDGGiwM6AgwE7SzgYcLHg4iCDIyxY9dsoWMWn97InuhRtZs7mcJlvq7A26Ishu18AzJms6JNkQaFY7GP4jaBAZVA9EoCjAbUyCJUhtSoItWG1KgS1IajnVVmFdA5cLOFiwcGCgyFcHFhpMFOQxhIsBVYaKAIcLAAACeBs9vEGg7zF8tuQlU9UJRUAICCCoEBIgrAMPlFt4MEvgkcAn6h4BPAI0ChAA6+qKmSYiEwOcsyQYyayzFBoJXLMUGQl8i1QYCUMNSkkHoJVBm0Rsp69gUSXQu8wEOpLqJZP62m8q04gwJzFhqp5W5E50UVJjKMBda8X9vnU/V71gA8O+tWDPqjn1RwzUWSFfAuRZ4F8C9HTBVkmMssE2WbINBHWDtQ/n0+x25v3obetPz0kQT2v1oWhOgQ1YbUmBIf96sZaOBxQjgbgaFBNY6GTjSixE53t0NlOlDqIrnboZCNIo/XbR4UUkBVRZY02eA8A5iyu6oeGRJdC7zAQ6kuoRsgfbbyrTgCw5JkCVeFUCISKCvt96rYGdYcbtjeqOz3qbo9qpaHUQZQ6iK4O4oIi6Gwni22QZyEovYYNMw2FNFFo1e6dpJSVQajwqwd96gE/bKxVF+5T9nmhNqx2cxI9nESvNOjlIvqkE90dBJ0CvQCBqrAlz5AVEyeyRQIDob6E6gyWVqaJJZdL1mHCsAxbG9Sf69Rf6tUtDeqvDWq2meiTRvROg/MKiTv7kN2dhNOQX9rp5Fkgz0IMy/5djAxIsNuj7nKr5W71gwPqXzYphwNqDydxZgbRP50YlEn0zyBsxms1tS5YxVtyTYkuRXsw2o6hfomx4dn+lPCj0RHBa6Tdzlqw5Joq1yZPD8xer/pDjbq+Rv2pRi13q92dxMAMYmAmcVUp2S+DcCTjuT4SVhoGZBADMn6LjiEJfm1UN9erv9Sr7+1Xfm1QuziIIZnE8BxieDbR05UkXakGHSDUmLO4UK2AgfA08KPREYNOGdVYcrlgpYFXUKgA2xrUbyvV1VXqmiqFIYmzc4hh2cSVXcmBmThz5JTMNAzJIoZkHQ93ogLbGtQfa9Xvq9Qntyq1IfXsHGJULjkmjxiUaeBO1GAVn9HPmehStJM5iw3V8I7OlkQXRL8wEOpIsCpsNexVpymNlcKyFJKN1aI94FNXHFVXHlVXVSoZHFGWR0wpJp4ZSneyJUVDJu4YEgZmEgMziVt6AgBUh2BNlfJ9tXrTGuWgXx2VS4zPJycUEr1cBvt4DV03LbkmXOZ7ehgIdSRQyVsLDDkOAQBAgCWbC1Xz9hK9X3gKCqyuVD87rHx2WPUI6vgC8qJOxL+H0/kWg52d9S/HDNM7k9M7AwDU8/DNMWXlUfXf2xVFhQsKiQuLiPEFpP5XNGoJ5Rm77gt6CtY809Fd/kSXQteM+tUmpcCxcPZgV6JL0X6WXFOgSr+B0CvC54eVxQfVL48oPVzERZ3Id88hB2QSGP3iI4ODSzuTl3YGANjpVj8/oj63XbnqW3l0HjGlmJxcTGbp9SIwUBW25Om1cBGw5psCx7BFeDoYCHVDhWC1gQfk4fjEUd3VN58In1QoHx5Qv6tURuUSU4rJ/xvOZCf/Kg9d6+EieriIP/UhPQJ8flhZXKHe86M4MJO4rAs5rUR3ETFYxRu6YrIOGgBEn2TcRm2s4eeiF+EGgTZTxhpga8GSa2rcqZceGFGBzw8rC/aqXx1VyvLIyzoT88cwKTvbU7ecLMzoSs7oCmGZ+uKI8uF+9f6fxOE5xMyu5LQSvfSaBquMPGYBAACWXC5QGXbZbYkuiE7p44eGAALHwlYjd7+AbnpgtjSob+5SFu5XejiJP5SS/xnJpBn4aj5VmCiYUkxOKYagRH1SobyzV7n9B3lKMXltd3JkboJ7r40+ZgEA1jxToDLs6o6B8OQwEOpFoDJszTd2IGQdNEEC3yhyaQloeflFeHef8tpOpTYMs7qT6yfTne04/Gc8Fhqu6Epe0ZWsDsE7e5Wb1sgKwHVnkH/sRmYmon6oshqsChu+RZhn8h0IJroU+oWBUC+ClWHjLlRqYisy+w6H4hwId7rV58uVhfuUsjzyH4OpCQVJsog7xeWY4a6+5F19yR9q1P/sULp9IF7UibytN3lWVly/3UAVz6WxRs8mb803Va1LnnwXUYeBUC8Cx/hO5xv7qhMA7EVm/+FQZj9HHN5LBfjyiPrsNnlrg3pjD3LbdFz/kJyGZxPDs6lGnnpzt3LFN3KWGf7Uh5xeQsZneb7/UNDeyfBzq6y5XLCGVxUVLxJPCgOhLiiiwntEg2YZbc7WyXL027pYv4ugwDt7lae3KQwJf+pDLu1CcgaeY4QiksbB3X3JP/Uhlx1Snv1Vue8n5Y4+5A1nkLFObeo/HEqCVPIkS3JOOlQrWHJwwPwkjBoIvSI8vVWlSNnFEgDgYMFCg4OBdI7IMEGmicgxG2mb92Alb8lik+BizVZk9h8OgdrWTYEiFZLg1Z3K09uUHi54bjg1Lt/wnxhqE5KAi4vJi4vJTXXqE1uUxzeLs3tRt/cmYzcfyn84lDssPVavHkeWPFOwMmysQFgXhtqwWheGBl71ihCUwM1DOgc39Ihyb4BRAyEJkMYSBEk0CioAVPghIIFXgAZeqQtDPa828JBtJoptUGwjOtuh1EGc4SR6pREuXTa6ApVhi8FnymgYK0WbqVCdEPXWbVCCl3coT2+Th2eTiydQgzIxBKa0QZnEB+OoPR7yX1uUbh+IN/ck7+pLpUf7JK9IaqhWSI66qU0czeyvx4kIAQnKG9XdHnWPV93vhYN+9XAAKoOqjYYsM5FpgnSOcDBgpcHFARGDHBhGDYQ2Bub2BfrUowSiAtUhtcIPB33qfh98c0x9ZYeyw606WaJvOgzIIAZlEkOziAKrLs6ngUrDr51oojUKoxgIwzL8Z4fyxFb57Bzyy/Ppvum6+MqQHnRzEm+Mpv7iI/+5Ren+gTinN/mnPlQUt8QKHA2ZszmSToafnDXPVL2hMdGlOK4uDD/Wqpvq1F/q1C0NanVI7eEiujmIbk4YV0CU2MhONsi3EHEb8jBqIGwVQ0KhlSi0wtk5v/sRV/jVLfXqL/Xw9m71lrUySxIjc4lROURZHtE7LWE/98DRcHovA25MfzK2IrP/cDBrYBQuPGUV5u9R/vqz0j+D+Px8+kwMgehkSuzEqyOpB84k//az0v1D8b5+1Oxe0Rk29h9Onr2mrQWmwJJELvPd71O/O767i1obVodmEYMyiZmlxFMZZGd7grewTtpAeCrFNqLYRkwu1u5R+33qmip1dZX69DYlKKnjC8gLiojzCsmMOHakq4rqP5IMA/IaW5H50Je+jr/Ol0fUe36UXRwsHEsNz8YQiFrR2U68XUaVu8kHNyj/t116bDA5s7Sjo+6+w0FnF2t0ypdopgxWkVTBI7LO+K1u8ouw8pjy+WF1xVE1LKvj8smRucTcfmQPp75mRKRcIGyhi53oYieu7gbwvx15Ptyv3rpW7JtOTCkmp5YQXWK/KDtYyXNpjKGTqzVnKzIFjoY6MlG73K3evV7e74MnzyIvLjb2+i0UZ71cxJIJ1Joq9a4f5f/brjw7jBqR0/4q7D8UKijLjGLxEsveyew7FMroG/NAWBmETyqUJRXKD9XqsGzigiLy9t5kArvcWpXqgbC5znbixh7EjT2Al6lVleqSg8qIpXKRlbi8K3lFlxiOJvoOBe2ddLpjQzvQJop1MqFqvh0J+90CPPKz/O5e5aH+1K29SCY1gqAsKIJbFHySGJBFvyQFZDEgyWFFCstyWJF5RZVVKSQDgBSSQT3+LFVVm88aoM0UEECZSJIiKBNFcSRlImkTRVspxkLRVpq104ydZh100lxyncbIXOLHi+l39yozvpHL8ognziLbscZU5hXeLRo63XYL9mKLryKY0TdWy3wbePhgv/L+fmVrgzqxiLzhDPKj8aTVCEHGCGWMO46C8wuJ8wupF8+mvqtUF+5XzvxY7p9B/LE7Ob2ENEf7M/NVhOzFSdIvqnGUWDwHgm0KhCrA/D3K/RvkSZ3I7Zcwett/ICqksByuE8J1QrheCDeIfKPAN4q8WwQVWBfDWCnGRjN2mrFQXBpDmyjKTNEmkmJJgia06EWbKPjfxYHf77fZfksdqcVIOawokirzsswrUliWQ4oYlIM1vOgLin5J9EuCR1JklUtjOCfDpTGmdIbLYM2ZrCmTY6xJFSAJgCtLySkl5OOb5TM/lh44k7q9d9vW4HsPBG2FZn114XWMvdh8eEVt1F9WUuDTw8pbu9XvKpULi8g/9SHPKzTY0l4MhKdDETA2nxibTz0/nFp+SHlrt3LnD/LlXclbepJ9otfM91UE80dnROvV9MBZam3Y7s0bEenqqx1u9eY1ckiGpefSg5NlXYQcVgLHQsEqPlAZDtXwwWpe5hVTJmvOZE0ZrK3AlNHHzqWxnIuhTO1p+VIS2bxtF3k7TxEU3i3ybpFvFMP1QuMO37E6IVwrAAmWbM6Sw1nyTJZczpZvpo0fGq00/GMwdU038rYf5Lf3KK+c3YaeUs9ev7M0SQYINbZOZv+RDg1btHDQp/5npzJvj1LqIK7tTs4fw9iNucELBsKIcJS20TZ5NKC+sVs9/wu5qx1u601OLSE7ONlJCstJ1v0CAM5S64GlVZEsq+dleHyL/FK58shA6uaepKEvvmVe8R8O+Q6F/EdC/sMh0SdZ8kzWPM6SZ8ro6zBnc5xLFycJkiXN2Zw5u+VPTvRLwWo+VM0HKsP1W72Bo2HKRNoKzbZCk63IbC+20BajxsVuTuKL8+kPDyiXfSNfXEz8c3BESyzcewNdJufGvnTxQ5sozsUEq/gO5vdXAVYeVf9vu7y+Rr2mG7lqIn2G08hVFwNhWxVYiYcHEA+eSX5SoTy3XZn7k3J77w7lefIfCiVZ9wsAcC6GNlOt7qexvka9drXc00VsnmbUNKGCR/TsC3gPBL0HguF6wZpvsneyZPSxF5+fbco0WKogxkY7bbSz629toHCD4D8S9h8OHf223nf4MOtgHCUWRxeLs4vVlKnLzBSndWlnckIB+cAGuc9H0stnUxd1Ot23I4XlUDVvL06ewXuNNkzY7kCoZTd8ZptCEnBHH/KDsdEfKkqIpPgn4o4mjzcQN9WpT21VHt8s3tSTvLMP1Y5FF8k3QKhxllo9ewOnqm8hCR7aKL+/X31uOHlJZ4NNiZECsnuv37074Nnrl0KKo4vF0dmSPSTNVmBK8GKoaDOls6Z09ngKdRUCVWHvgaB7T6Di8xqCAGep1dXd5upu0zZANwQXCy+fTc3ool7/vbxwP/F/w0+ZjMa7L2gvtiTZFwoA9mKzryKUO7zNTwxI8J8dyrO/Kr3T4N9Jl93QML9gfRqUSSwcS+3zkk9tVc74QLy+B3l3X6pNEz18FcGcoWkxK2DCuLpZaza5Tzr2ub5G/eN38qBMYut0Op7rNTtEBf/RUEO5XymAeAAAIABJREFUr3GHL1QtOLpYXN1t+SPTLbkmI+W07QgCrHkma55JG/oN1wnuPf6G7b79n1RyTiathz2tp83R2WKIRnBZHrFlGv3njXLfj6RXRpKTOp3kUsy9x+/qllQDhBp7seXY6vo2PcUvwgvlyr9/lcvyyGXnUv0zDPAVtxUGwijo6iBeGUk9NIB8YovS80Px5p7k3X2piLIAq+A7FCq9rCDmRYw7Z6l174fHWgzLCwo8skl+e4/ywghqWokBGoKqrHr2Buq2eRu2+2gTmd7LXnxhrrNLEjYU2sqUyeZmpucOT1cV1X843LjDd+CTKr5RSOtpz+jrSOth0/naFwsNzwyjppaos1bLSw6q/x5OtZjl4dkbKL00CSumJZfj3aIUlmlT6yO+IQle3qE8uVUel0+umkj3dCXtzx4DYdQUWYkXRlBz+5GP/qKc8aE4tx91W2/y9D+2YA1PsaSBepYix9ho1skEjv6Woarcrf5hldzJRmyeymTruzNYVVT37kDdFk/9r15zFpfR19FvdmcjjorFAUES9mKzvdjc6fxswSPWb/dVrm3Y/d6RtB72rP7OtF52PSfqHJVLbJ5K3/2j3P9jaV4ZNTL3eFHFgMw3iraiJFzEQ5CEvZPFuz94+pyOsgpv71Ye+VkZkkWsvJCO4iR5fUrCU3BiFduI10dRO/uSD25UXiiX/nnaPE/uXX5X9yTsftG4ulndewK2IrMK8FK58ref5ceHUNedoeuGgu9QqHaTu3azx5TBZvV3Fp+fHc98VEbHOpm8Eel5I9KlgFy3zVu5rmHPB0cz+jiyBrlcpVZ99iHbGPjPSGr5IfXyb+RrzyD+OoCiSfDsDRilm7cdXGdY3bv9pwmEXxxR7/lRzjTBonHU0NTIboiBMCZ6uIiPx1PrqlvJ89S4y59zliv+xYuPtB62I9/UscMzr10t1YZh3WS61KHTSiX6pJqN7uqfGlUFsgc5z7y9iykD23/tR1up3GFpucPSBJ9U94vn4LIqMSDlDEnLPstlStfjB3tRJ+KXqfS1q6WRy6V3xlDKDl9aD1vrTzMmV3fb7neOnPRPWnbDg3546qxWZtUmGQyEMTQih/hhMv3ePuVyLc/TELJ5njZVVr0HAmdcWZjAEsaUs5vt1/lHZnwYntKLeWQgpccxIxXce/xVPzS4dwcy+jlKLytwlFj02XAxKNZO54/OyB+dETgWrv6pccu/99sKTTnD0jP62PXW3so2w7Lz6Be2K6M+ERfu8A67IDvRJYoVW4FZDMi8W2y+sNUtwN9+lt/dp/y5P3VzT52P8EYfBsLYIgBmdiUvLib/uVnuv1i6/0zq9t7Hf2TeA0FLNmfcRcqnJynw518Uh9X2YmZgzGDd5c2Rw0r1hsbKtQ0kQ+QNT+92eWH70rugCFnzTV2m5JVclFu/1VP5ff3+xZV5I9Jzh6cxNh2dggiA23qTI4L8tr3sG5vJF0dAu9cH6xoBru5W9y6/Nl9dBfjvHuWBDfLkTuT26UxmEg6Mtk5Hv8IkZqXhscHUrO7kbevkt3crL51NjcolGnf5Xd2Ts/ulwq9e8Y2czsELF7p8P9QB6CgQhhuEyu/rqze607rbul1e4OicbCum9YykiayBrqyBrsCxcOXahk2P78no5ygYndGO/Oyx49jvHXOOaw0Jg5dI74+jknIXTNcZNvdOf87QtF8b1VvWynxyZTdsh7YFwgMHDpAkWVxcfNK/iqJYUVFhNpsLCn6bdhwIBARB0G4TBOFyJe2QWKtKHcTn59MfHVBmrpLPLSBu2eXvdnFSJXDSLKlQbl4jz+1H3dWXVCXbTx8dEX0SY0/8JVfgaPjIqjr3Ll/O0PQBd5fqJOFZarLmm0ovzS+5MKdqfcOv/zloLTAXjs1sntQmUVRZbfjV2/+87Ndd1Lv7lAmfSX8bRN3SM9l6C9K62w4sq37wJ/n13crfB1E39jB2dsOOi/T0FAgEJk+efODAAVmWe/fu/fHHH5tMv7uIe+ONN+677z6XyxUKhfLy8j766CMtXl5//fVLly7lOA4AXC7X/v37o/4/GMv0zuSEQvLRtUJDpfCpbJqZ6PJEES/DvT/Jyw6pS8+lz8oiAICgifSe9rptbUjAHQveA8HDK2uDleH8sozSS/MpLtnOawZFW6nCcVn5ZZm1m9x7PzzGWKjC8VnpPe0JHKZ17wmYs45nhZ3ZlRySScxYJa86pr42KqL0pEbxnZ9ukMnAsfDW6dZcfa9lio9IzwgvvviiLMu7d+/eu3dvY2PjG2+80eKAbt26/fLLL3v37j106FDXrl3vu+++pj899dRTDQ0NDQ0NGAU1Dgbuc4VsnS1P/qpe8IV00Ke2/hzd2+tVz14mHQ3Cz1OPR0FNZn9n3WZPokrlPRD89eWDu989ktHXMfih7gVlmRgF9YakiZyhaYPu65ZfllnxefXmZ/c1bPdBgupE3RZPZn9n091uTmLdJDrHDAMXSxtqk6Ge1oXh6m/lm9bIWT3t99iCGAU1kZ4UFi5ceP3119M0zTDMtddeu3DhwhYHjB49uqioCAAoihozZszhw4eb/iTLcmNjY7RKnBzqtni6DnJsnEKX5ZFDPpGe+1VRjFzLFu5Tzl4mzepOLhpHuX5/4ew6wxasDIcbhDgXyVcR+vWVg7vfPZI12DnogW65w9IwHYyuEZB5pmPAXaVF52Yd+qJ687/3Ne70x7kIMq80/OrN6Pe7fWs5Cp4fQf2/oeRFX0nP/mroagrv7VP6fiRmmmDbdHrQMEf9loRdoepNpF2jBw8e7Nq1q3a7a9euFRUVpzpSkqT58+dPmjSp6ZEHH3zw73//u6Iojz766K233nqqJyqKsmXLFofjd7/C/v37U1SyzauUw4p7V6D0sgKahPvPJKeXEDeskRfuV14fRfU2WgaHoAR3/CB/X6V+eT590iSEJE3kDEs/trq+y5S8OBWpij/4aVXgaLjo3Oycs1x6m6aPToeAjD6OjN6O+m3eA0srD6+kSybmxG1CU/VPjc5utpMOHk8tIQdkEFeskr85prw1mjbc1MojAfWWtfIhPyw9lx6iddh0tQp+OVjNW3KMkvA3hghVjegSh+O49evXDxgwAADWr19/wQUXnKqRd/vtt2/atGnVqlUsywLAwYMHO3XqRJLk6tWrL7zwws8++2z06NEnfaLdbu/cuTNN/y42f/rpp3b7STIghEIhhmFaHGwUDb/4PTuDna/4baGSCvD2PurvW6kbuyl395JYI/TeBQKBA6Jt1jpmYIb6zCDJSp/yhyR6pZ0vHuv1p5gvURA8UtU3bu/uUPYoR+ZZDj0n9+q4FjvUJyEVGrb4q75xm3PZvAlppqzYTm5SFXXHc0dLLs22FJ5yMFBU4NFt1IcV9H+GiUPsAUOcglSAN/dS/9hG3dJd+VNPqfkCwWNfNhIMkTfWqBMYA4GA1dr6BCuLxUKSrZx5Iv0Wc3JymiJffX19Tk7OSQ974IEH1q5d+/XXX2tREABKSkq0G6NHj77wwgu//vrrUwVClmW//fbb9PSIZlVQFGWIX+FJHdxelzc8o8VZ7LYzYVqpeus6pWwF/dooapi+MxupAK/spp4sp58ZRv2hlAQ47UWlDTJ6+33bhMJzMmNUHplXjnxdW/lDQ97w9O4PFUWSUNjoVFVN8kAIYBtpKxyWU7m2ft9b1Rn9HMXnZcdu+nHdVq/JxWb3aOX888zZcEGJ+sfviD90of42kDKxuj4F7fKoN3wvyyp8N4nqdULK7Pyh1M55h7tNKjRuEoloVYFIr9AHDBjwww8/aLd/+OGHgQMHnnjMww8//Nlnn3311VenWiNRXV3douczBQleyX84dNJEfwVW4pMJ1J8HkNNWSnf8IPvE+JcuIlUhuOhLaeFB8ofJ9B9KI/oJFZRlVn5fr8rRH2FRFbXqh4ZNj+8WPNKAu0uLL8xJhSiYOkiaKCjLHPRAN4olf35yz5GvaxUpJuN0x76tKyiL6EJtQgHxy1R6SwOM/kzd7dHpoKGowD9+UUYtky7tTH5/EX1iFAQAW6GZoAjfoWD8i6c3kV7O3H777Zdffnnfvn1lWX7ppZeWL1+uPd6zZ8958+adddZZL7744uOPP/7II4989NFHAGCz2WbOnCmK4ty5c8ePH2+1Wj/++OOtW7e+8847sfpXDKJusye9j/00KYwu70JOKCDv+VHu85H0woiTb5aWQEsqlFvXytefQd41XHQ5Ih1dsBaYzFlszUZ3dDdfdO/xH/ikirZQvW4osRUYbdzm/7d339FxVfe+wM/03rvajIrVLEuWLFu2bMvdhmAblhcmBJzFC4/Y1NxQcsN9QBZZcENyCVxISACTC8YvIXlcDAFsXADjbtlWL1YbtZE0vfd2zn5/zI1iZCOrjGakmd9n8cdotKXZyHv2d845+/w2mDIqi5K7U6laLR46bGz6dZ9mu0JaIYjjcYxL6wt7ouKyyTZkuJachX28gXi3n7Lm8+gL1ZS9xfPrQvRFM9p7FtfwsIY7qDncybomqxJYmlw8dbqXlZhqEG7atOntt9/ev38/iUQ6ePBgbW1t7Pn169fHjv/4fP7999+v0+l0Oh2GYRKJBMMwCoXCZrPfeuutSCRSWFjY0NCQlZWypTWnyHzFodlxk/voxQzs3TrKNwb00Dn8vV70+ipyNif5bzRHCPvJRfyyBR3aTF0lJ3mnuaZPs1PVuX9IvJgXl6paQWt48DOjzxDM3amULEn30wxpgimhF9+X4+r3DX5q1J+15d2hGt/kazaIKNJ+pM+9XTmtdVUkDHuwGNucRb3vNP7xIPHOWsrkkZMYjhD29BX8yAh6tYZ8V97NP0PLq4Stvx/Q7FCm9gX1m5rqYpkEkEgkfX19U7xGuEAXyzh7vQOfGKr+ddEUP8yGcOw/2ojfdeJPLaE8voScxEU0hwaJf6kndueSflVNYVExbEaLNYYOG0POSNGe7Nn0JBrAR760mBucmeulGXWStH0DezyeG64jSwsIM112DB8ziwq56tsUs9zRc/ioKWAKFf+vnGn91PgUhCPs5TbiPzvw56so+5JXooVA2Hu9xLMN+O488gvLpnH7/9U/DYsX85SrklnyYmbiuF5soQah1+Eb+dRGIpGobAqJglHoFCqHQmNTaDwqQ0CjC6jzqpjvuPY/DipWiOTV01umNeBBj9cTPU7025ok7I2i86LHLhL9brR/zbc2k5rBKCTCRNPL2vxdKlHJTGZwhCPjRbvuS4ukjK++VT4//4kTJq2DEMOwf6yQMl60q9ZIMjdIKTP6nOg3BNvfHKp8qmC6aTrhs3iXE+09hxMIe3sNJfHb2F4wocfrcRoZe6OWcsO7mCbhHvT3fjC67N8WzcMbjaIBPOyMhJyRsDca9eERH06ECQqLor5FjsU1CBfqVEKiYOJyLoVMiQZwhCM8RES80YApFPFGQ85IyBVBUcSU0FlyBlvJYCuZ3Exm0neY8wz7Q/aIrEpw86bflscjfbqFcnQEPXUJf60De7mGUjnNsT4zgSj223bid534T8so/71pZvPMt5Dp5IK7Mnv/MlL2UO707l5CmLXVNXzUzJTQyh7UcOZTjWaQLBQGWf09hXKVeOiIqfGlvpytMsWK6ZVNCHuiXe+PaLbP9pgSw7ASIenMduo73cSmL6J355Gfr6KIEnJ7Xr8b/Z8G4qIJ/fty8p7v3gN8EvxcNkNIsza7ZMuSfB9F2BXxjgb9xqDfGPKbQ0FrGEMYXUhjCKk0LpXGoVA5VDqfSuPEP7YW6hHhTU+NRoN40BYJmEJ+Y9BnCPrGgniI4KlZPDWbn8fha1iJ33Hr6n8Ni0p4s6m6GSWwd3qIF5uJ1QrSL5eRS260EiwucIT93z7i+SZihYz02xryDS9+zPjjmKXJOfi5aclDGpZ8ClMFwmydbt1xM4lC0nxPkar7dcwAHBFeyzsSGPrCFLSGc7bKZVWCqcRhxBNt/+OgvFqYtUk2g1f8rinIFsKea8A/HiL+tZzycCl57pYwj/jQi83Ex0PEvyymPLGEzJ5FOji6vYOfGat+VpDg+ygQjjwjAfeAzzMUcA/7MYRxs5icDCZbwWQpGEwJncaZ7M8Hp0Znco0w4o16hgPuIb97wOfTB7lZLFExV1jE5WayEvDP79EFut4drn62aPbXtPxR7I2rxCvt+Fol+efl5OWyePY+SmB/GyB+1ULIWdhLyymrvvt2xtmMQnODc/gLU8HujEnOkeJhwtrsGjtlJdPJ2VtkksX8hXu301yAILyeq9838qUlYA5lrJMolosm2ezTM+zv+9uYrFKQvXWGG/BOPgVddaJnG4gGC/p5Bfn+QjIrrscwPS70chvxyRCxt5j8s3KKOB6Hni2v9mdtlF5bZ3XuBMwhR7fX0e1xD/lZUoYgn83TsHlq9nQ3hIEgnO1iGTxMuAf8jm6Po9uLhwjxYp6kjC9cxJmjcpREFLW8qs3ZKo/jOPNFsf/qIV5tJ7K52IPF5DtzyYzZffa0BrGDfcTvrxK5XOzfllK2ZN7kTzHLUejs9fZ/bGCKaZkbZDw1a/zEa9gddQ/67Z1ue6eHn8fOXCcVFCR/d555CILwu3hHA/rTNvtVj6CAIynn83PZTPH/XBYhosg7GjCct7v7fZrtClnVzE8GTmUKumJBL7USF0zEgyXk/10027XfOMKOjqA3u/BGK3q4lPJoKTkuERjjGQ50vTdc+VTBXF16R5h70Gfr8Ng63AhHomKeqJgrKOBQWTOftiAI47lqNGgN2zrctja33xKSlPKllQJhISe+142HDpuCtnDxfbNaLXlDOMI+1xFvdRENFnSHhnxXHnmdkjStRPRGsKOjxEeD6MQosVNNfqSUvGJqh5izH4UIR4YLdmuLyzcWpAtoiEBEmCAIxNewhYVcWaUgzZfDTA6CcHJ4kLC2uhzdXvegj8ARlUkhUUhhZ4QlZ4jLeJnrZ7i4ZtzUp6BuJ/rDVeKv/USNnLQ7l7xTPb0AwxFWb0YfDhD/PUhouKSHSsm7c+fkjOucTFMIcw/5Lc0uW5uLxqVKy/nixXxOnG75hSCck9snwq6ItdVtaXYFbWFZpUC+TMjNicNdSp4hf9cB3Rx+1MIwDMPGfOjDQfTRINFuR6sVpDVKcqWEtESMZbAnHuXiCBvyoB4XdslMnDWiRiuqVZB2aci788jC6SwniuMoRDgK2sIkColMI9N5VDgFOhUQhFMX8UbxEIFwRBfS4rDoC8Ow6U9BgSj2yTDx8RD6cpQoE5PWKEgr5aRSESmPR5qwXAFhmCmAddhRix1dMKFTBkLNJe3SkO/OIy0SzOF7I3biKnuLXFYZhxNXflPI0ug0NzopDLKsSiitELBkcV6uCEE4t/cRBm1hc6PT0ujCSJhiuVC2TDjj3cx9+mDn/qGCuzJvWFNtLjjD2Dd6ot6Mmm2o04FZg0jGIjHIGJeGhXDMGUbOMKZkkQoFWLWUtEZJXqMk8Wf0P5f6RZ/nNwjC5JrxFBSIYvUWdMaArliIHhc24kM8GiZikDhUzB3GIgRmCiAxAysSkpZKSMulpM2ZZEWidg30jgY69w8X7cma8aq0qA+3NDtNDc6wKyqvEsiqhXO3xhuCMEE31HuG/KYGp7XVxc1kyauFkiX8aW3r6h0LXt0/lLcrQ1qRtNInUQIzBVCIwLwRjEXBeHSSiI7N8mpiDARhckEQJle8pqAIgTnDmDOEAjjGpWFUEqZkk5JYOsM96O96T1d4T5aoeBrvbiKKHF0ec4PTpfWJSnmKaqGwkDvXp3YgCBNaWYaIInunx9zgcPf7hcVc6VKBqIg7eSIiHBnO20e+shTszkjVAmAQhMkFQZhcC7S41VR4hv1d7+oUNaKszbLJzyQTUeTS+qwtLluHm5PBlFcLpRWCaR0tzAYEYXJGYdSPW9vc1laXZ9jP17AFBRx+LoeTwbz2Hz5oDTt6vIZzNoaIlneHakq3yi1MEITJBUGYXCkchBiGhV2RwcMmd79PtVYiKuZee3qTiCK/Iege9Lv6fa4+H1vFkFYIpEsFsy9KMF0QhEkehXiIcPZ4XQM+95DfbwiR6SQ6l0pEUTSAk2lkUTFXWiGY1omFhQiCMLkgCJMrtYMwxjPsNzc6Hd3eiCdK41LJVFLEj+MBnCVj8HLZgly2sJg3+T3vcwpKrCUZhUGWlPMl5f9zzjPiw6O+KIlCojDIsOIfAJAaeGp2bIemaACP+nGEIwqTkvgjvwRIwf+lxKNxKEn8WAQAAHOKyqLM5s73+W9+bfoKAAAAJBgEIQAAgLQGQQgAACCtQRACAABIaxCEAAAA0hqsGgUgmRCGvGGfPxII4eFANBCIBnEC94Z9CEP+SABHOIEIX8Q/3j4YDUXwCIZh4XCYTqdjGEYhU9i0fxajZFNZFDKFhJG4dA6GYVwah0wic+hsJoXBoDI4NDaHxiaT4BMwAP8EQQjAnPBHAtaAzRl0OYIuZ8jlCrldIbc75PWEPZ6w1xP2+cI+X8QfiAZ5dC6LymRQGWwqi0VlUsgULp1DwkgsGotKopBIJC7tnzsyMqkMJoOLYVgICzEYDAzDokTUE/KONzB6zQQiYvmKYZgn7CUQ4Y8EgngoGA35I35fxEen0Nk0NpfG5tI5PDqPT+fxGFwBgydg8AUMvpglEjIEEpaIR4eCCSAtQBACMHOOoMvkM5v9VrPPavJbrH6bxW+zBez2oIOEkWRsiYDBFzIFIqZQwOArOPJCUT6PwePRuVw6h0fjsGnsaw/mpmU2lWWC0aAvEvBFfJ6w1xP2xuLZFXIPuUacQZcj6HSG3PaAIxgNiVkiGVssYYllbImcLZOzpXKOVMGRS1giEmyXBVIFBCEAN0cgwugzj3r0Yx7DmMcw5jEavEaDz8SgMJQcuZwjlbNlCo6sWFwgY0skLLGEJWZS52+ZWSaVyaQyJSzR5M0ieMQedFj8dmvAZvHbzH7rVWuP2W81+syekFfJlau4ikyuKoOnzOJlZPEyMrhKKjmVb7sGqQqCEICJogQ+4hkbcuoGXbph14jOPTrqMYiYwiyeKouXkcVTVcjLMnkqFVfBos7VXmvzAY1CU3DkCo78+m+F8bDBa9J7TXqvcdSjbzC0jLjHrAG7giPTCLJz+FkaQU6uMEfDz6ZRZriXJwAJA0EIAOYN+3rt/VrHQJ9jcMA5POIeVXBkuUK1RpCzLqc2h5+Vzc9kUOK8v/aCRqfQ1YJstSD72icjRFTvMQy5RnTu0fqxK3+9emjMY1ByFXlC9SJR3iJR3iJxnogpTFafAfguEIQgHQWjoV67ttvW12Xr67b1OUOufGHuInFepWLJncU7NIIciL0ZoJGpE9IxSuDD7pEB55DWPvi3rk967f1MKrNInF8iKSyRFhZLFnFo7CR2GIAY2IYJzNCC24bJ7Le2m692WLs6LT3D7pE8oaZYsqhEUlgkzs/mZ5FJC2zpxwLdhknvNfba+7usvV223j7HgIItWywrLpOVlEmLs/mZye7dNMAUlHSwDRMAU6L3GltMHS3mjlZTRwgPl8lKymUlmzXrCkX5cO0qKTK4ygyucn3OagzDcIQPOIfbzV0Nhpb32j6IENEK+eIK+eJKRbnm2yddAZhTcEQIZmjeHhE6Q64GQ0ujsbXR2IYT0UpF+VLF4nJ5Wc6COuC4qQV6RDgJk8/Sau5oMXU0m9oD0WCVsnyZsmK5qlLOlia7azcAU1DSwQ71MAqTb14FIY7wTkv3JX3jZUOzwWtaqliyTFm+TFmRw89KdtfmSuoF4bVMPnOjsa3R2NpobOEz+CtUlTUZyyrki+nz5totTEFJB0EIozD55kMQusOe+rHGi2NXGgwtKq5iRUZVTcayUmkhhZT6d7OldhCOQxjqsw9cNjRf0jf2OwaXKspWZS5flblcyprSRDF3YApKOghCGIXJl8QgNPnMZ0cunRut77X3VykrVmVWr8yovunt4SkmTYLwWp6w97K+6cLYlcuGpkyuak1WzdrsleokXU2EKSjpIAgxg9P0+KlnSSQShmGxwox0Cp1OoXFobBaNxaGxeHQul84V/qN2opgpEjEFUGs4jhIfhGMewynd+dO6CyafZXXWijXZK5cpK9L2Poc0DMJxOMJbTZ1nR+vPjdQzqcz1OavX5dQWiHIT2QcIwrhzBJ2OoMsasLlCblcwVpLX648GApGAJ+yNEtFgNIRhGIfO/s9NL2IQhBiGef1eR8RFoVAwDIuV6g9FQxEi6g37AtGgL+L3hr2esNcVcjtDbnvA6Qg63CGPgClQcmQKtkzJVcSqQ2XylHKODKomzkDCgtDkM58cPvfN8DlLwLY2a+UG9ZoK+WL4TJPOQTgOYajb1ndad+GU7jyNTN2gXrtRvTYxK04hCGfMFnDEqhUafSaD12Twmi1+qy1g59I5IqYw9p+AwefRuTw6l0NjsWgsHp1LIVNihZyEDH6s2hEE4UxGIY5wW8Bh8llMPrPRazb4TKMew5jH4Al7snmZakF2rlCdK8gpEGluWFMKTDDXQegKuU8On/166MyIW78upxbybwIIwgm6bX0nh89+M3yOz+Bt0tRt1qyb0+WmEIRTZA84tM7BAefwkFMXqzpEp9AyeaosXoaKq1BxFEquQs6WStkSGnl6f0wIwniOwkA0qHOPxgpLDjiH+x2DITxcKM4vFOfHbrhWcGSzf5XUM0dBGMbDF8auHB/4pt1ydVVm9SZNXbWyEko5Xw+C8IYIhNotV78eOnNKdz5XqN6Wu2FdTu1c1K+BIPwu9oCjy9bXY+/rsfX32LU4wheJ8vKEao0gJ0+ozuZnxmt7LwjCuR2FzpCr197fa+vvsvV22XoRhi2WFpXJSpbISorEi2BSjol7EF619hwd+PqU7nyhOH9b7sa67JXMlC5pPUsQhJOLENFL+sZjAyebTW2rMpbfmr+pUlEex/pBEITjCIQGnIOt5qud1u4OS3cgEojVzysSFxSK82Rum+tJAAATDklEQVRzdlwOQZjQUWj2W9stXZ2W7nbL1VGPvkRSuFRRVqkoL5EUpnMoxmsUOoKuE4PffNH/JY6IW/I2bctdP3fvnFQCQThF7pDnq6EzRwe+coU8t+Zt+l7+prhc+0jzICQQ0joGmk3tLab2NstVKUu8RFZaJitZLC1KWKk8CMKkjUJfxN9mvtpiam82tY969OXy0uWqyuWqqhSrWjIVsxyFBEJXDE2HtSeaTe1rslfelr95iaw0jt1LeRCE09XvHDqi/fLroTMF4tzt+VvXZK+c7kWpa6VnEJr91sv6xiuGliZTm5gprFKWL5UvqVAsFjIEie8MBOG8GIXukKfJ1HbF0HzF0EzCSDUZy1ZmVlcpyufzjqxxNONRaPFbv+j/6kj/l2KmaHvBlo3quhnv0p7OIAhnJoJHzozWH9GeGHAObc3dsL1g28w+xc6HKSgxogTeZu68pG+s1ze4Qu5qZeVy1dJlqqVJr2kAQTjvRuGQa6Re31A/1tBr76+QL16dtaI2c4U4pW/xnu4oJBBxSd/0Wd+xDkvXJk3d9oKtCb7xK8VAEM6S3ms8rD1xbODrbF7mzkW31GWvmlYd9vk2BcWdJ+y9pG88N3qpwdCSxcuIla1YJM6fP/u0QBDO31HoDfsuGZrOj166om/O4mfUZa+qy16VyVMlu1/xN/VRaA3Yj2i/PNJ/QsIS7Si4ZaN6bZocNM8pCMK4iBL4udH6z/qODTiHbsnbtL1gaxYvYyo/OG+noFmy+m1nR+vPjtR32/oqFUtqs1bUZi6fn9spQxAugFEYJfBWc8dp3YVzo/UCBn9dTm1ddm2eUJ3sfsXNTUdh7Crg59rjrabODeo1OxdtKxDlJax7KQ+CML7GPIbPtcePDZzME6p3FGy76RXE+T8FTYveazwzcvG07sKoR78qc3ld9qrlqsp5XrYJgnAhjUICoavWntMjF87oLtAptHU5q9fnrE6Bs4KTjEJrwP5F/5dHtF8KmPydBds2adax4EaIeIMgnAsRInp25OLn2uODTt0teRt3FGz7rtM5C2gKmsSYx/CN7vxp3XmL37o2e9W67NqliiULZTE8BOGCHIUIQz02bawcFJlEXq9evSFn9cI9SLp+FBKIqNc3HtGeaDNfXa9evaNgW6E4P1ndS3kQhHNq1KM/rD1xbOCkRpC9vWBrXfaqCTtALcQpaNyoR39Kd+HU8Dl70FGXXbs+p7Z8AZZtgiBc2KMQw7Bee/8p3flTuvMYhq3Lrl2fs7pIUpDsTk3PtaNQ5x472v/V8cGTSo5ie8HWjeo1cC/8XIMgTIAIET0/eumI9stue98mdd2t+ZuKxP/zPl2IU5DOPRorW+8MuupyajfkrC6Tlc6fxS/TBUG4IEfhDfXa+0/rzp/SXYgS0bqc2nXZtaXSogUxNL1eL0FDJ4fPnhg8ZfSZtmjWfy9/c7L2xElDEISJZPJZjg18fWzgJINC35a3cbOmjkviLJQpSOsYODNSf0Z3wRvx12WvWp9Tu6DzbxwEYeoE4bgB5/Bp3YWzIxcdIdearJq1WSsrFUumtZ47YfyRwIWxy8e1JzvtPSszq7flbliuqlxw51UWOgjCxEMY6rB0HRs4eWbkooafs1lTt169WsDgJ7tfN0Agot3SdW700tmRiySMVJe9al1ObYm0MJV22oEgTMEgHBdbvnVh9LLWMVilLF+ZsawmY9l8qDrmDnkujF05N1rfZGyrkJetUlRvWbQeVsEkCwRhEkXwyJmhi+f09ZeNzcWSRXXZq1ZnrpCyJcnuF+YKuS8bmurHGq8YmpVc+erMFWuyV+YLNcnu15yAIEzlIBznDnkuGRpjY1rMFFarllYpy8tli7l0TsL6QCCizz5wxdBcr28YcA4vV1WuzqqpzVzOpXOSuEM9wCAIky02BeEkInbX+cWxKwq2rCZzWY2qqlRanMiFl8FoqNPa3WRsvWJsGfMYKhXlKzOqVmZUz4dgnlMQhGkRhOMIhPoc/VcMLS2m9qvWHhVXuURWslhWVCxelMXPiPu5jjAe7nMMdFi628ydbearEpaoWrV0ZUZ1hXzxtadqIQiTC4IwuSZMQQQiOq3d9WMNV4wto279YlnxEllJuay0UFwwFxUETT5zt03bYe3utHQPOIcKRHlVyvJqZcViWTGFtDBufpg9CML0CsJrRQm8167tsHZftfb02LSukDtflJsryNEIcrJ4qkyeSs6WTuvKYhgP670mvdcw7Bodcum0jsER95hakF0mKymTlVTKy76rUBwEYXJBECbXJFOQJ+xtNXe2m6+2W7r6nUNytnSRKDdPqMkRZGXyVBlc5bQuKOAIt/ptsV3Eh1wjsTcpjUIrEueXSorK5CXF4kXpWaoJgjB9g3ACd9jT7xgadOqGXLoxr2HMY7D6bTwGT8QQ8Bk8AYPPprEYFEbsjUcgwhfxE4jwRnzesM8edDqDTm/Yr+LKVVylRpCtFmQXCHPzhOqpRCkEYXJBECbXFKcgHOFDrpF+x9CQSzfsGhnzGg1eI41Mk7BEIqaQS+dy6RwqicKisagkCoZhQTwUxsPesM8X8TuCLmfQ6Qq5RUxhBk+VxVOp+Vm5QnW+UJPadYynKI5TULoHyULHp/MqFUsqFUvGn0EYsgecjqDTFXJ7wl5fxB/Gw4FokISRKGRKFjWDhJG4dA6PzhUxhSKmUMRMwv4pAKQJComSL9RMWK7iDntib1JP2OsN+3CEB6JBAhEEIqQUMYNC59K5bBpLwOBLmCIxSwRLsucaBGGqIWEkCUskgQ+MAMxXfDqPT+dp4KbbeQM+aAAAAEhrEIQAAADSGgQhAACAtAZBCAAAIK1NLwgjkcjkDQiCwHF8Bj8IAAAAJMVUg9BoNG7cuFEikchksvfff//6BgihJ554QigUikSiBx54IBqNxp7v7OysqKiQSCRqtfqrr76KW8cBAACAeJhqED755JMajcbhcBw/fvzRRx8dHh6e0ODDDz/8/PPPBwYG9Hp9c3Pz/v37Y8/fd999d999t8vlev311+++++5AIBDP7gMAAACzM6Ug9Hq9hw4devrppykUSlVV1ZYtW/785z9PaHPgwIG9e/dKpVIul/vYY48dOHAAw7COjo6urq6f/vSnJBLpjjvuUCgUhw8fjvv/AwAAADBjU7qhXqfT4Ti+aNGi2JelpaVarXZCG61W++ijj05ooNVqc3NzWSzW+PN9fX2TvJDT6SR9e7tIoVBIuuEGkggRfi+R9iXWkggFfASstUoeFPARlNTZW27BQcEgEaHCFJQEZBKZGecdeKb0r+h0Ojkczngg8Xg8u91+fZvxsm88Hs/pdBIEEfvB8TY8Hs/hcHzXq4RCoaqqqgmx19HRIRDcoAaY32KM7P/5VDoP5ghCyL3wN7leuBBCHvj7g/RDFil4D/8HhmE+n28q7dlsNpl8k8/sUwpCqVTq8XgIgoj9OofDoVAorm/jcrlij51Op0wmI5PJUqnU7XaPt3E4HCUlJd/1KgwGY+pFtymUTNoL/w+KbicRFN1OLii6nVxQ938+iNcUNKVzW9nZ2Ww2u62tLfZlS0tLcXHxhDalpaXNzc0TGhQXFw8MDIwHZEtLyyRBCAAAACTelIKQxWL98Ic/fO6554xG48cff3zx4sU9e/ZgGNbV1XXLLbfE2uzdu/ett95qaWnp6+t75ZVX9u7di2FYQUHB2rVrn3nmGavV+tprrxEEMd4eAAAAmA+melz/m9/85vHHH1+xYoVSqTx06JBcLscwjCCIUCgUa7Bt27ZnnnnmBz/4QTQavf/++++5557Y8wcPHnzkkUeqqqoKCgoOHz4MZxIAAADMK7AxL5ghuEaYXHCNMLlgCkq6OE5BC3X9+xtvvPHFF18kuxfpq6mp6bnnnkt2L9LaPffc4/V6k92L9AVTUHI1Nzc/++yz8fptC/XjTHd3N3wcTiKbzTa+eAokxblz56CEbxL19PTAGZEkstvtra2t8fptC/WIEAAAAIgLCEIAAABpbR4tlrnlllt6e3tvXFDtOh6Ph0ajMZnMue4VuKFIJOL3+29Y9AckhtVqlUgkU3y/gLiDKSi5IpGIz+cTCoU3bXn48OGb3r8+j4LQ5XLZbLZk9wIAAEDqyMrKotPpk7eZR0EIAAAAJB5cIwQAAJDWIAgBAACkNQhCAAAAaQ2CEAAAQFpbAEH4/vvvZ2Rk8Hi8nTt3Xr8hMIZhnZ2dy5cvZ7FYZWVlFy9eTHwPU5jZbL7rrrvkcjmDwaipqblw4cL1bV555ZX8a4zvugXi4oknnhj/21ZUVNywzWuvvSaXy/l8/j333OP3+xPcw9Q2YXgXFhaGw+EJbVatWjXe4KGHHkpKP1PMCy+8sHXr1vz8/OPHj1/7/IsvviiVSgUCwQMPPHD9PwSGYSdPniwsLGSz2bW1tVqtdoovN9+DsL+//9FHH/373/9utVpZLNbTTz99fZt77713165dPp/vySef3L17N9SdiiOfz7d69erW1lav17tr164dO3YEg8EJbRwOx4YNG778Byh9F18Wi2XPnj2xv+2nn356fYPLly+/+OKLZ86cMRqNFovlpZdeSnwnU9iPfvSj8bG9efPmwsLC69fiDw8P//GPf4y1ef7555PRzVRDIpEeeOABhNC129CfOHHizTffbGhoGBkZ6ezs/N3vfjfhp/x+/1133fWrX/3K4/Fs3Ljx/vvvn+rrofntF7/4xZ133hl73NLSwuFwgsHgtQ0aGxt5PN74kxqN5vDhw4nuZXrw+/0kEqmnp2fC888888yTTz6ZlC6lgz179rz++uuTNNi3b98jjzwSe3zixImMjIyE9Cvt4Diek5Pz8ccfX/8tlUrV0dGR+C6lvNLS0kOHDo1/+f3vf/+ZZ56JPf7www9LS0sntP/ggw/Kyspij91uN4PB6O3tncoLzfcjwr6+vrKystjjxYsX+/3+sbGxCQ0KCwsZDMZ4m97e3kT3Mj0cO3ZMoVDk5uZe/62DBw+qVKqampq//vWvie9Yyvv1r3+tUqnq6upOnDhx/XevfY+UlZXp9XrYlWIuHD9+3O/333bbbTf87m233aZWq++8886pn44D0zVhqGu1WvTt++D7+vqWLFkSe8zj8XJycvr6+qbym+d7EDocjvES71Qqlc1mT7hM6HA4OBzO+Jd8Pv+G1xHBLA0MDDz00ENvvfUWjUab8K077rjjm2++aW1tfeqpp3784x8fPXo0KT1MVfv27Tt9+nRzc/O99957++23X19x/9r3SOy8tMPhSHQv08C7775733333bBGyTvvvHP+/PlTp06JxeJt27YFAoHEdy8dTBjq4XB4wme+GcfBfN+GSSKRuN3u2ONYfUuZTPZdDTAMczqdExqA2dPpdJs2bXr++edvv/32679bXV0de7B79+5z58599NFHt956a2I7mMrWrFkTe7Bv376jR49++umnE5bMSCSS8QVKTqcTwzCpVJrgTqY8m832+eefNzU13fC744eJb775plQqbWpqWr16dQJ7ly4mDHUmkzlhJyyJRKLX68e/nHoczPcjwuLi4vGPwG1tbTweT6VSXdugqKiot7c39hEMIdTe3l5UVJSEjqausbGxjRs3Pvzwww8++OBNG0MN6Dl1wz9vcXHx+N6QbW1tarWaxWIltl+p7+DBg5WVlaWlpVNpjKBu5dyYMNSLioomvCOubeB0OkdGRgoLC6f0q2d9OXNu6XQ6Lpd7/Phxu92+ffv2xx57LPb8L3/5yw8++CD2uKam5umnn/Z4PK+++mpubm40Gk1ef1ONxWIpLCy89957G/7B7XYjhC5fvrxnz55Ym/fee6+vr89sNh86dIjL5R45ciSpXU4pBEG8/fbbQ0NDRqPxT3/6E5PJbG5uRghZLJadO3darVaEUHNzs1AovHjxotlsrq2tfeGFF5Ld6xRUXl6+f//+a5956aWXDhw4gBDq6en57LPP9Hr98PDwgw8+mJub6/f7k9TN1NHT09PQ0JCbm/vyyy83NDR4PB6E0KlTp2QyWUtLy9jYWEVFxe9///tY4717954+fRohFAwGFQrF/v373W73T37yky1btkzx5eZ7ECKEPvnkk7KyMoVCcd9998X+HAihxx9//J133ok97u/v37p1q0QiWbt2bWtra/J6moKampqWfVtjYyNC6OzZs9u2bYu1eeCBB/Ly8mQyWU1NzV/+8pek9jfVEASxY8eOnJwchUKxbt26Y8eOxZ43mUw1NTVmszn25YEDBwoLC1Uq1SOPPBIKhZLX39Sk1WprampcLte1T/785z9/4403EEKdnZ1r165VKpWxxTLXL6sGM7Bv375rp52WlpbY83/4wx/y8/MzMzN/9rOfjR/z7Nq1a/ytcenSpZUrV0ql0ttuu210dHSKLwe7TwAAAEhr8/0aIQAAADCnIAgBAACkNQhCAAAAaQ2CEAAAQFqDIAQAAJDWIAgBAACkNQhCAAAAaQ2CEAAAQFqDIAQAAJDWIAgBAACkNQhCAAAAaQ2CEAAAQFr7/3OXZqU8yFs2AAAAAElFTkSuQmCC", "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "execution_count": 13 } ], "cell_type": "code", "source": [ "ρ = scfres.ρ[:, 1, 1, 1] # converged density, first spin component\n", "ψ_fourier = scfres.ψ[1][:, 1]; # first kpoint, all G components, first eigenvector\n", "ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n", "ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\n", "\n", "using Plots\n", "x = a * vec(first.(DFTK.r_vectors(basis)))\n", "p = plot(x, real.(ψ), label=\"real(ψ)\")\n", "plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n", "plot!(p, x, ρ, label=\"ρ\")\n", "plot!(p, x, tot_local_pot, label=\"tot local pot\")" ], "metadata": {}, "execution_count": 13 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.1" }, "kernelspec": { "name": "julia-1.6", "display_name": "Julia 1.6.1", "language": "julia" } }, "nbformat": 4 }