{ "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 1D example\n", "and we show how to define local potentials attached to atoms, which allows for\n", "instance to compute forces.\n", "The custom potential is actually already defined as `ElementGaussian` in DFTK, and could\n", "be used as is." ], "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\n", "a Gaussian potential." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "struct CustomPotential <: DFTK.Element\n", " α # Prefactor\n", " L # Width of the Gaussian nucleus\n", "end" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Some default values" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Main.var\"##419\".CustomPotential" }, "metadata": {}, "execution_count": 3 } ], "cell_type": "code", "source": [ "CustomPotential() = CustomPotential(1.0, 0.5)" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "We extend the two methods providing access to the real and Fourier\n", "representation of the potential to DFTK." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function DFTK.local_potential_real(el::CustomPotential, r::Real)\n", " -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)\n", "end\n", "function DFTK.local_potential_fourier(el::CustomPotential, q::Real)\n", " # = ∫ V(r) exp(-ix⋅q) dx\n", " -el.α * exp(- (q * el.L)^2 / 2)\n", "end" ], "metadata": {}, "execution_count": 4 }, { "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": 5 }, { "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": [ { "output_type": "execute_result", "data": { "text/plain": "2-element Vector{Main.var\"##419\".CustomPotential}:\n Main.var\"##419\".CustomPotential(X)\n Main.var\"##419\".CustomPotential(X)" }, "metadata": {}, "execution_count": 6 } ], "cell_type": "code", "source": [ "x1 = 0.2\n", "x2 = 0.8\n", "positions = [[x1, 0, 0], [x2, 0, 0]]\n", "gauss = CustomPotential()\n", "atoms = [gauss, gauss]" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "We setup a 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", " LocalNonlinearity(ρ -> C * ρ^α)]\n", "model = Model(lattice, atoms, positions; n_electrons, terms,\n", " spin_polarization=:spinless); # use \"spinless electrons\"" ], "metadata": {}, "execution_count": 7 }, { "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 `gauss` 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 log10(ΔE) log10(Δρ) Diag\n", "--- --------------- --------- --------- ----\n", " 1 -0.143581349776 -0.42 8.0\n", " 2 -0.156036092192 -1.90 -1.10 1.0\n", " 3 -0.156768710397 -3.14 -1.56 2.0\n", " 4 -0.157045890563 -3.56 -2.32 2.0\n", " 5 -0.157055546819 -5.02 -2.99 2.0\n", " 6 -0.157056387622 -6.08 -3.71 1.0\n", " 7 -0.157056406098 -7.73 -4.37 1.0\n", " 8 -0.157056406912 -9.09 -5.40 2.0\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380295 \n AtomicLocal -0.3163467\n LocalNonlinearity 0.1212608 \n\n total -0.157056406912" }, "metadata": {}, "execution_count": 8 } ], "cell_type": "code", "source": [ "basis = PlaneWaveBasis(model; Ecut=500, 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": 8 }, { "cell_type": "markdown", "source": [ "Computing the forces can then be done as usual:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "2-element Vector{StaticArraysCore.SVector{3, Float64}}:\n [-0.05568392329096156, 0.0, 0.0]\n [0.05568743546713474, 0.0, 0.0]" }, "metadata": {}, "execution_count": 9 } ], "cell_type": "code", "source": [ "compute_forces(scfres)" ], "metadata": {}, "execution_count": 9 }, { "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": 10 }, { "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+gvaeTAAAgAElEQVR4nOzdd3hTVf8A8HNHdtKM7kX3bilt2VCg7GIZMkUBReWFF1FBVFAcKG7llVdcCAguUBFEoWwQCjLbQksL3XumTZo97/j9EX6Vl1E60iRNzufx8UnDyc3Jyc353nPuGQhN0wCCIAiCXBVq7wxAEARBkD3BQAhBEAS5NBgIIQiCIJcGAyEEQRDk0mAghCAIglwaDIQQBEGQS4OBEIIgCHJpMBBCEARBLg0GQgiCIMilwUAIQRAEuTSHC4T5+flms7mTiSmKgkvE2RFJkvbOgkuD5W9fsPzty4rl73CBcMqUKVKptJOJjUYjRVG9mh+oAzqdzt5ZcGmw/O0Llr99WbH8HS4QQhAEQZAtwUAIQRAEuTQYCCEIgiCXBgMhBEEQ5NJgIIQgCIJcGgyEEARBkEuDgRCCIAhyaTAQQhAEQS4Nt3cGoN5lpkCxki5R0pVqQFDARAEhE3ixQYQQ6S9BGPBCCILsh6RBkYIuVtINOtBmBEwUIAjoxwORQiRahHBh9WwrsKSdU5Me/FxOHaun/m6ifblItAgJFQAmBpgoKFeB882gSEGVqegB7siMYHR2MBIsQOydZQhyFS0G8HsVtbeSuiilfblInBjx4wIxC2gJQFIguwWUKKlyNZ3kjkzwRx8JRSKE8OfZu2AgdDbH6+nPCsm/m+mHg9Cno9CfxqBi1r1T6ghwvpneW0UN+ZPsL0FWxWPpgQj8wUHQHT7++OMrV67c/TxJkhiGdelQShOo0dAyI/BgAy82GMdCcBQAAJoBaL4tWSQAYTRQGMFeI/2FHnBxEMAD3hxX/3mOGTNm+fLlvXFkxNEWrQ4ICLh06ZK/v39nEuv1eiaT2dVz0VkdrqXX55IaM3ipPzonFOV1+iLHTIFfK6j/FFBmCnw0GJsc0Nmfm1qtFggE3cwu1GOw/G0jNTV14sSJ0dHR9s6IS8vNzS0oKDhw4ED7M1Y8/2GL0BkUKehVF8kqDXh3IDojCEW7eN3IQMFj4ehj4egf1dSqi+SnPPDFCCzczcWvPiHoH2lpaSNHjrR3Llwah8MpKCjopYPDwRJ9m5kCb1+lRh0kJgWg+TPxmcFdjoK3mx6EXp+Jpweiw/4kNl6nSMfqLIAgCOoVMBD2YQVt9OA/iEtS6urD+Mp41CpDQHEUrIxHL03HM2uosZlEvRYGQwiCnJx1AuGPP/4YERHh5eX11FNP6fX6+yU7dOjQwIEDDx48aJU3dXHbiqmxmcRzcWjmJNyfZ+VuzFABcmIKPjEAHbifOFoHYyEEQc7MCoHw5s2bzzzzzM6dO4uLi6urqzds2HDPZEqlcs2aNXK5XCaT9fxNXZmOAAtOk5sLqawMfHFkb7XpUQSsG4D+Mg5/6iy58Trc/RiCIKdlhWr022+/nTFjxogRI8Ri8WuvvbZ9+/Z7jkRdvXr1ypUrPT09e/6OrqxGQ6ceJDAEXJyGR4t6fTzLKB/kwjTspzLqySzSBKMhBEHOyAqBsKioKDEx0fI4MTFRKpW2tbXdkebUqVPl5eVPPvlkz9/OlZ1vpof+SSwIR78bjXFsNeA3kIecnYq3GUHGUUJtttGbQhAE2YwValO5XN4+mcPNzQ0A0NraKpFI2hNotdrnnntu//79SCfmg2q12oCAgPY/f/jhhxkzZtwvsUvNI9xfi67Kxr8ZSkzwNWo0tn73nUPBC9n46APU3lFmT/atFr9Wq+3Mdwr1Elj+tkFRTtIZIpVKk5OT6+rqLH9u3bq1sbHxjTfeuF/6I0eO7Nu375tvvrFVBh+AJEnNbXVfJ89/LpeLog9o8lkhEEokEpVKZXlseeDu7n57gvXr1w8aNEipVObk5Gi12qqqqurq6qCgoHsejcfjFRQUdHJCPYZhLhIIPy2gPr1OnXgIS5Qw7ZWH7WlgfS455TR6Ih2zDM+haZrP59srPxAsf9t4YDXaV1AU1draanms1+vffvvtnJycDtJPmjTppZdeKiwsjIuLs0kGHwDDsNtPeCue/1YIhBEREYWFhZbHBQUF7u7utzcHAQAURV2/fn3p0qUAgKqqqp07d2IY9tprr/X8rV0BDcCrV8gDNfT5aViAtUeHdtX6ZIzPoEZnkienYEF82BaBIDuorq7Oy8sLCwvbs2fPqFGjxo4dW1RUdPDgQbPZPG3aNEvQ0mq1mZmZhYWFHA5n6tSpd0eyvXv3pqSkeHl5AQBOnTrl5eUVHx//3XffzZw5U6vVZmVlzZ07F0GQxx577Kuvvvr888/t8DltyApXOosXL963b9+1a9d0Ot0HH3ywePFiS3N1w4YNhw8fBgBs3Lgx+//FxcWtX78eRsFOImmw5Cx5ppHOysDtHgUtXkxAV8ajYzLJKjWcVgFBdpCfn798+fIVK1aIxWIURTMzM9PT0xEE4XK56enpWVlZAIDc3NysrKx+/frRND1u3Ljs7Ow7DpKZmZmWlmZ5vH379lOnTgEAVq1aJZPJSktL169fb/mntLS0zMxM2302O7FCizAhIeHjjz/OyMjQaDTp6elvvvmm5fni4uLAwMA7Evv7+8PVETvJRIEFf5FKE318Ct75hUNtYEUsigIw9hCZOQaJgV8m5HpyW+k1V0jbvNdQL2RDyp13f1Qq1f79+4VCIQAgMjJy+/btY8eOBQB4e3u///77o0aNSk1NTU1NtSRmsVjffvvtwIEDbz/C9evXFy1a9MB3j4yMrKqqUqlUlvEfzso69euSJUuWLFlyx5M//vjj3Sn37dtnlXd0enoCzD5JMFHkz4k4y/HugS6PRUkaPPQXM2sq7SBNVQiymWABsqa/jX6WHux7PBkVFWWJgiqVqrS0dMOGDe+//z4AQK1WS6VSAEBjY+OyZctu3rzJYDD0ev3dXaNqtZrL5T7w3S034dRqNQyEkK1pCTDtGOHLRXaOwnBHvU//bByq1pMTDqNnHsK9OPbODQTZkIQFxvvb8/qPxbq1uZpltOCWLVvap2hbhva89tpr0dHRlrH6X3zxxaFDh+44gqenZ/s8NzabbTAY2v/JYDCw2bfCr0wmQ1H0jvGPzsdRa1kXpiXA1KNEiAD5frTjRkGLZ6OIeaHIhMOE3GjvrECQS2Kz2Wlpabt27RL/P5PJBABobm4ODQ1FEMRkMu3evfvuFw4aNCg/P9/yODEx8dixYwRBWP48cODAgAEDLI/z8vISExPb46KzcuyK1vWozGDCISJahGxNxXqyj4TNrE/Gxvkh044ROsLeWYEgl7R169bDhw8nJyc//PDD8fHxH3zwAQBg2bJlr7766sMPPzxw4MB7zkabM2eOZTAjAGDp0qVCoTAyMlKj0YwfPz4nJ6d9pcwjR47Mnj3bZp/FXmDXqANRmsDkI0SiBPliBNYXguAtG4diS86S048TmZNwJryygqBeNnny5DFjxrT/GRwcfPHixaqqqubm5pCQEG9vbwBARkbGzZs3y8vLIyMjBQKBpZno5eVVW1treVVaWprRaMzPz+/fvz+Lxdq7d69UKo2IiPjxxx+HDh1qSaPT6Q4cOHD+/Hlbf0Kbg/WWo2gzggmHiaFeyFcj+1IUBAAgAHw9AuNgyFNZcAdDCOp1DAbjjrH3CIKEhIQMHTrUEgUtvLy8hg0b5u7uzmQyLWNeUBRtv5WIIMjmzZsvXbp0e3oMw3x8fNqfuXLlyquvvmqZa+jcYCB0CJYomOqDfDq0j0VBCxwFv4zFqjT0S5dsNKYcgqAeGjFixB2j/Z9//nnLYFSL0aNHu8gC0TAQ2p/cCMYfJtJ8kY1DHG+eRKdxcPDnBPxIHQ33bIKgPurNN98Ui8X2zoUdwEBoZzIjGH+IGO+HfNyXo6CFmAWOTMY+K6R2lcNYCDkPRxsIdvjw4QkTJvTSwQ8ePLh8+fIOEhQVFU2cOPGee+31XTAQ2lObEUw+TIz3Rz4c3OejoEUADzk8GXvhInms3ql+J5DLKlPRNRrHOpnj4+NXrlzZG0emaXrNmjUrVqzoIE10dDSbzT5w4EBvZMBeYCC0G6kejMkkJgUgHzlLFLSIFSG/jsMXniby5I5VfUBQVzXqwMTDpA/HsW7cEwSh0+kAACaT6cMPP6yvr1+/fv3atWurq6t1Ot2nn3760ksv5ebmWhK3trZu3bp11apVGzZsKCkpaT9ITU3N22+/vW7durKysk2bNimVSgDAmTNnOBxObGwsACArK8uyAOnOnTtramrUavXGjRstr124cOFXX31l40/dq2AgtI8mPUjLJB4ORt4Z6FRR0GKUD/LFcCzjKOlol9IQ1HlqM3joKPFUFCpi2Tsr/6ugoOCzzz4DABiNxrVr1z711FMBAQFKpXLChAmLFy8GAEgkkrS0tJaWFgDAhQsX6urqBg8ezGQyR44cWVZWBgBoaWkZMmSIwWCIjo5etmzZK6+8IpfLAQCHDh1qn5hx9OhRy3o0n332WUVFhUKheP311y3/NGbMmL/++kuv19vhw/cOOI/QDuq19LhD5KII9NUBTnshMjsErdeC9CPkuam42MHqEQh6IDMF5pwkBnsi6wagR+76V1N1keK3L2yTE2ZovOjhpR0keOeddwYOHPjkk096eHiMHz/eMhD02LFjp0+fnjNnztSpU6dOnUrTtEKhqKur27179+uvv75t27bRo0e/9957AIDk5OT4+HjLoQoKCqZNm/bALHl6enI4nLKysoSEBGt8RPuDgdDWqtT0+MPkshj0xQSnjYIWz8ejtVp6xnHiaDrOdsJ2L+S0aACWnCWZKPLFiHufuLh3P9Hc52yTGZTD6ziBJYyhKOrh4dEe0ry8vGQyGQCgsLDw6aeflsvlAoGgqakpIyMDAFBSUpKUlGRJGRsb275yqU6n43A6tXAwl8tVq9Xd+kCOCAZCmypV0hMOky/1R5+JdfIoaPHRYOyx0+TC0+QvY/vGinEQBAB4PZssUtKnpuD3m9WLsrnMwAjbZuq+MOyfaG1ZcRsAgCCIZWDn888//8QTT1j2RX/llVcs0VEkErWvuK3Vai3rzgAAPD09LX2kAAAOh6NQKNqPrNfr22MkRVFyufz2yft9nUtUxw4iX06nHSLfTHaVKAgAQBGwcxTWaqBfgBPtoT7i65vUr5X0gYk41ymaCe0RS6FQ/Prrr5YnH3rooV27djU0NAAA2ofAAACGDBmSl5dneZycnHzixIn2Zt9vv/2WkpJieVxSUiIQCEJDQ232KXqbq9TIdndJSk88THw6FF0c6VplzsLA7xPwUw30J3CiPeTw/qymNlyljkzGPJ1lu4XVq1c//fTT6enpqampycnJlifHjx+/dOnSpKSk4OBgiqI4HI5lb8K5c+e2b0MxZcqUadOmRUVFFRcXL1iwYP/+/Zs3b7a8PDMzc86cOQjiRJ08tIPx9/evq6vrZGKdTkcQRK/mxyqO1VFeP5oO1VD2zoiVqVSqTqas01BBu80/lZG9mh9X0/nyhzrj7ybK60dTdsudv9ORI0eePXvWLlm6J4Ig9Ho9TdOWLsr25xUKhdlstjzWaDQGg8HyWCaT5eXl6fV6vV6v1Wrb01MURVFUYWGhUCikqFufes6cOXv37m1Po9PpYmJifv/99/ZnSJJMSEgoKirqtc93bwcOHMjIyLj9GSue/07R+Hdsv1dRy/8m947DR/o40QVUF/nzkEOTsbGZhDsLmRTguuUAOawyFT3nJPn9aDzFw9HPTwzDLPcFEQS5fUW025cJ5fH+GWIjkUgkEskdB1myZEl8fLzZbP7666/XrVvX3rz74IMPbp8sz+Fw2Gy2SCRqf6asrGzZsmVRUVFW/Ux2BgNh7/pvAbXxOnV8Ch4vdvRfV2+LFSF7x+MzTxCZk/CBDl/XQC6lTkuPP0R+OBh1nau0xx57LCcnB0GQ77//fvjw4e3Ph4aGPv/887enfOKJJ/r169f+Z2RkZGRkpO0yahMwEPYWGoC3csk9FfTZqVgQ31V+XR0b4Y1sT8WnHiVOZ+BRQlgmkENQmkDGUfLZOHRBuAvdvx8zZsztmxp24LnnbDRRxI5c6Iu3JT0B5p0k/2qgz03FYRS8XUY/5N1BWPoRskEHF52B7E9HgIeOEhP8kdXOPq8X6gD87q2vSQ/GZBIsDBxLh4uq3MOTkejSaHTSYbLNaO+sQK7NsnxMuBvyUd/f+wXqCdg1amVXZfTDx8kno9DXk+AM8vtak4i2GuiHjhLHp+A8eA5C9kDR4IkzJI4i21L73m7YFEV98803V69eDQwMXLFixe2DWaBugC1Ca/q1gpp8hPhkCPoGjIIP8tEQLEaEzDxOGOFUe8genr1A1uvoX8ZieB+sBZcuXXrw4MH09PTq6uq0tDSShL+iHoFX49ZBUGBdNrmnkj6ejveXwCD4YAgA36Rij5wiH/2L7KOVEdR3vXKFvNJCn5jSzVVwb7SWbLrytbUzdW/9vWJXpDx9+zN1dXV79uypra0VCATTp0+PjY09fPiwZRFRqHtgILSCZj145BTBwcGVGbg7vCnYaRgCfkrDZhwnnswid46Gi5FCNvLeNepgDX06A3djdPMIQcKA1UM62sbdiviMOxfdLi4uDgsLEwgEAAAEQZKTk4uLi2Eg7AkYCHvqRD39RBa5JAp9PQmFVXlXMVGwdxw+5Six7G9yy8i+d6sG6nM+LaC+K6VOP9Sja1YegxslCbdeprpGp9OpVKr2P1UqlZubm70y4xxgh1T3mSmwLptcnEV+Pxp7MxlGwW7i4ODARPxGG/38BXifA+pdX96gPi+kTk7BfLn2zkrPlJWVnT59GgBQWVl55syZtLQ0e+eob4OBsJuKFPTwP4l8OZ37MD7WD8bAHuEzQOYk/JKUfuEijIVQb/mmiPoonzo5BQvg9fkf7IABAzZs2DB8+PBBgwa9++674eF2a546B9g12mUkDT4rpN6/Rr4zEPtXNLySsA4hExxNxyccIlZfIjfCSV2QtW0tot69Rp2aggUL+nwUBADweLyTJ09WVVV5enrevqwo1D2wHu+awjZ65AHiQDV1YRoOo6B1iZjg+BQ8q5FefYmEq85AVrS1iHrnGnVqChbm5gxRsF1wcDCMglYBq/LO0hHg1StkWiaxOBI9+RDuZL8oByFigmPp+Lkm+rnzMBZC1vH5Deq9PKeKgjExMcuWLbN3LpwKDISd8lslFbeXqNGC/FmMf0XDYTG9SMwCx6fgV2X00nMkBYMh1DOfXKc2FVCnH3KeKAgACA8PX7BgQbdfTtO0wWCwSk4sWx5a5VD2BQPhA1xpoUcfJN69Ru0Yhf04BvPh2DtDLsCNAY5MxsuU9MLTpBluaw9115s55PZi6sxDTrj9S0FBQV5eXsdprly5UlJScvfz165di46Otko2YmNjr1692u2XEwSxZ88eirL/jxwGwvu6oaBnnSBnniAXRaA5M/Axvs72W3JkfAY4NBlXm8GsE6QBjiSFuogGYOVF8kANnZWB+/f9MaJ3+/nnn3/44YeO03z55Ze3b7HrgPR6/dy5cx1hfTg4avQe8uX0e9eo043US/2xH8dgHFhI9sDGwN7x2BNnyMlHiD8m4EKmvTME9RFmCjyVRVZq6FMP4SJnPG1u3rx56tQpgiDWrl0bFha2ZMkSgiC2bNmSnZ3t7++/YsUKHx+fS5cu5eTk1NTUtLS0JCUlzZs3756HMplMX3311dWrV4OCglasWOHp6Wl5/vTp0/v375fL5cnJyc8//zxJkjt37rx8+TJBEKmpqY8//jiK3rcRtXHjxrS0tD179kil0tmzZ6enp1ueLykp2bZtm0wmGzNmzIIFCxAE+fLLLwEA69atQ1H0mWeeCQwMtHZRdRZsEf6PE/X0lKNE+hFykCdSPo+xOgGFUdCOGCj4YQyWKEFGHyQadfbODdQXaAkw/TihMIFjk50zCgIAhEKhr6+vt7d3SkqKZbP4BQsWHDhwYM6cOSRJpqSkKBQKDw8Pd3f3gICAlJSUkJCQ+x1q1qxZp06dmjt3rlqtHjRokFarBQB8++23ixYtSkpKevTRR1tbWymKMhgMJSUl06ZNmzlz5vbt2995550Osrd169ZHH300Li5u0qRJTz/99J9//gkAKC8vHzZsmJeX1/Tp0zdt2vTKK68AACw9tMnJySkpKfYd/gqreQAAUJjAD6XUliIKAWBVAvr7eJQFZ7I5BhQB/x2GvZ9HjThAZE7CYkRO2M0FWYtUD6YdJ+JEyJaRvbuMu6pSV7Krrhff4DaiCF74XP/bn/Hz84uJiTEYDHPmzAEAlJeX//nnn/X19WKxeMqUKdnZ2Tt37ly5cmVwcHB8fLwlzT1dv349Kyurvr6ez+dPmTLl4sWLu3btWrJkyWuvvfbdd99NmDABADB58mQAAJ/P/+ijj/R6fXNz85o1a9auXfvGG290kOdly5ZZhvNotdpPPvlk2rRpmzdvnjVr1osvvggACA8PT0pKevPNN8eOHQsAmDVrFoPR3VVfrcSlA6GZAsfr6R/LqMO11ORA9PPh2GhfBFa0DuiVRDSQB8ZmEr+Mw0f5wK8IuodiJf3QUXJhOPpGcq+P6+YHsOOXBffym9yCMh8Q0svKyoKDg8ViseXP5OTke46RuVtpaWlkZCSfz7/9hQqForGxcejQoben1Gq18+fPLyoqCgoK0uv1jY2NHR85MTHR8mDAgAGWzJSWlrYvCx4bG4thWHV1tb+//30PYVuuGAi1BDhRT+2vpg9UU1EiZEE4+vlwhgTuGuHYFoSj/jxk7kni4yHYwnDYpQ/9j9ON9PxTxAeDsccjbHFuoAyU7W7njtf2eQsikUihULQ/39bW1n6rr2Mikaitre32F/r5+fH5fCaTKZPJLLtbWOzYsQPDMEtIu3LlysSJEzs+cnt+FAqFJUKLxeL2J3U6ncFgaI/cjsBVKhQjCc410e9eoyYcJnx/Mn9+g0pyR67NxP+eiv87BoVRsE9I80X+egh/K5d69QqcYgj9Y3sxNf8UsWssbpso6Ag8PT1ramosj/v374+i6J49ewAAtbW1+/bts4xPuT3NPQ0aNEilUmVmZgIAKioqMjMzJ0+ejOP4tGnT3nnnHctgzpaWFpqm9Xo9giAAAJIkN23a9MDsffPNNwRBUBT19ddfWzpX09PTv//+e0ss3Lx5c3Jysre3N5/P53A4HWfSNpyzRUjSoF5Ll6lAYRtd0EbntNI3FXScGBnlgzwfh6VNQHjO+bmdX4wIuTQdn3mcmHWC/H4MJrDznQXIzggKvHSZPFRLZ2XgEUIX6jOfN2/e3r17g4ODR44c+eOPP+7evfuJJ5546623pFLp2rVrR44cCQB48sknFy1aFBwcPHPmzP/85z/tr0VRlM1mAwAEAsGuXbueeuopgUDQ0tLyzjvvJCcnAwC++OKLxYsXW7pbDQbDzZs3H3/88e+++y4mJoaiqDlz5ohEIsuh2Gz2PYePRkRExMfHm0ym4ODgr7/+GgAwf/787OzsqKgosVjMYrF27doFAEAQZMOGDePHj0dRdN++fe0dqraHONq6AAEBAZcuXepk33GD0vDfYhxBEI0ZqExAZQZSPd2kBw062pONhLuBGBGSIEGS3JEB7kj3tqKGOqBWq2/vP7EZEwWeO0+ea6b3T8DCnWjFkK6yV/k7CJkRzDtJMFCwe2zvDhBNTU19//33LdHFkbW0tEgkEgzrck0nlUrd3d3veKFl10Nvb2/k/wdONDQ0SCQSSxDtQHR09I4dO5KSkjQajYeHx+3/ZDabVSqVu7t7V3MIADh48OBzH3w9bMMfPlxgWZffiud/324ZYQgtYgAURQJ4wI0BhEzgyUZ9ucCPi8Bhn06MiYKvR2JbiqiRB4htqXhGP9eNhS4rt5Wec5KcHYK8NwiDGzpbdPLW4N28vLzufpLL5XK5/7Nto5+fX+ePyWaz7w6ZDAaje1HQwpcLpgQivbEMcd8OhG4M8HJ/BMNc5cYAdLul0WiiBJl3irzUgqxPhrWhC9leTL2aTX45HJsVAn/7Duftt98ODQ3tjSNLWMhjvTNQDp5GUB821AvJnoFflNITDsEZ9y5BYwaLTpObCqisDBxGQcc0d+5cb29ve+eia+CZBPVtnmxwdDKe5oem7DcfrXOsG96QdeXJ6UF/EEwMXJqOR7nS0Biot8FACPV5KAJeT0J/Hov/6xy56iJcpNsJ0QB8WkBNPEy8kYRuS8W4ffuWDuRwYCCEnMQoH+Taw3iDDgz+g8iTw6ah86jV0pMOE3srqUvT8PlhsMqCrA+eVZDzELPAL2Oxl/qjEw8T716jCPtvcwb11Hel1MD9xBhf9EwGHiyA3aFQr4BdDJCzWRiOpvkiT58l91dR20dh/SWw9uyT6rT0snNkvQ4cT8ft+yWmpKTMmzfv7skANE0jcHFiW9HpdCNGjOilg/ftCfV6vZ7JZHZjAilkFY48oZsGYGcJtfYK+a9odN0AzCmXU3Dk8u8JigZbiqg3c8jn4rA1iSjD3v1W91tmWqvV2nfzIFfj4+Nz+9RGOKEegh4AAWBxJDo5AH3uApmwl/hqBDbeH1689wH5cnrpORJHwekMPNYxdt3icDj3nBjnrBciLggGQsiZ+XLBnnFYZi39r3PkIE9k4xA0gOcQdSt0N6UJrM8ld5VT76RgT0f3+lZKENTO3p0OENT7HgpECmfh0UKQ9Dvx3jVKT9g7Q9D/omiwvZiK+c2sI8CN2YwlMApCtgUDIeQSODh4KwW7PB2/JqNj9xK7y+E+To7iVAM9cD/xfSl1YCK+ZSTmDvdEg2wOdo1CLiREgPw6DjvXRK++RG68Tn0wCN44tKdrMvqVK2SZCrw/CJ0N10uD7AcGQsjljPRBLk7Hf6uknjlP+nPBOwOx4d4wHNpUkYJen0udaaReS8L+FW3/caGQi4MnIOSKEADmhKCFs/AFEehjp8n0I8T5ZthXagtFCnrRaXJ0JjHAHSmbx3gmFkZByP7gOQi5LhwFT0aixXPwmcHogtPkuEPE8XoYDieE10cAACAASURBVHtLbiv9yClyTCYRLUJK5zLWJqI82CEFOQZ4JkKujomCJdHo4kh0Vzn1wkWSiYEXE9DZIbClYh00AMfq6I3XySIFWJWAbktl8Bn2zhME/S8YCCEIAABwFCyKQBdGoJk19KcF5MuXqWdi0aeiUM8719WCOktLgB/LqM8KKCYGVsWj88PgtQXkoGAghKB/IABk9EMy+uF5cnpzIRW1x5wRiC6NQUfA0TRdUdhGf1NE/VRGjfZFPx+BpfnC0oMcGgyEEHQPiRJkWyr28WBsRwm15CyJAPBkFPpYOOrDsXfOHJjKDH6toHaUUNUa8GQkcnUmHgjX8YH6AhgIIei+xCzwQgL6QgJ6toneUULF/mYe7oU8Go5OD4IDPf5hpsCxenpXGXWolhrrh65NRKcEohiMgFDfAX/NEPRgqT5Iqg+2mcB+r6J+KqOe+ZucFIDOCUHSA1GX3S2doMBfjfSeSmp/FRUpRB4NQz8bzoDrwkB9kav+iCGo63g4WBCOLghHZUbwexX1TRH11FkyzRedFoRMCUS9XaPXVGkCx+qpP6vpw7VUhBCZE4Jmz8D78WEDEOrDYCCEoC5zZ4Gno9Cno9A2IzhUS/1ZQ6++ZA4TIJMDkYn+6DAvBHeu4ZEUDa7J6OP19NE6KruVHumDTOuHfjgY9+PC+Ac5A2sGQqPRyGJ11DNiNpsZDDiHCHIeYhZ4LBx9LBwQFPZ3M32snnrhIlmipId7I6N90VQfZKAHwuqbewKTNMiX02eb6DON9JlGyouDTAxAXuyPjfFFXLY3GHJW1jmj6+rq5s+fn5+fz2QyP/300wULFtyR4M0339yxY0dLSwuPx1uxYsX69eut8r4Q5CBwFIz2RUb7Yu8OBG1GkNVEnW6kV16gipR0ogQZ7IkM9kRSPJBwIeLIbahaLZ3TSme30Bel9JUWOoCHjPRBZoUgnw9n+HIf/HII6qOsEwhXrlwZHx9/5syZnJycsWPHpqWl+fv7357A09Pz+PHjUVFRRUVFo0aNiouLmzNnjlXeGoIcjZgFpgeh04MAAEBjBtmt9CUpvaeSfiWbkhvo/u5IghhJkCDRIiRGhNhxPobcCIqV9I02+oaCzpfTeTIaRUCKBzLIE1mdgA3xQiRw5AvkGhCa7uniigqFwtPTs7S0NDg4GACQkZExatSol19++X7pZ8+eHRcX99Zbb93zXwMCAi5dunRHHL0fvV7PZDIxrG/2PfV9arVaIBDYOxd9idwI8uT0dTl9o42+qaBvKmgDCcLdkFA3JIQP+vGRQD7w4yJ+XODJfnCfamfKn6CA1EA360G9FtRp6RoNXaUBlWq6VEkTNIgSInFiJFqEJEqQBAmA9/y6BJ7/9mXF8rdCi7C6uprJZFqiIAAgNja2vLz8fonlcvm5c+eWLVt2vwQ0TSuVSi73VkcMj8djMpk9zyQEOQIJC6T5IrevtKIwgTIVXammK9WgVEWfaABNOqpBB6R6moUBDzYiYQEREwiZCAcHPBxwcMD+/wBpMuFMJml5bKaAxgwMJNAToM1EK01AYQKtBlptBp5s4M1BAnjAn4sE8JCHAkGIAA13Q7xcY5grBD2QFQKhQqHg8XjtfwoEgsrKynumNJvNixYtmjRp0vjx4+93NJ1ON2zYMBS9Nepuy5YtU6ZMuV9i2CK0L61Wizj0Pa8+AAcgmg2i2QB43vlPKjMiMwKFCShMiMoM9CTQk4iOACbqVpnjtIkFbj3GMDqYDZgozcWBkAGETCBk0BIWkDDv0+VDAo2mtz6Ui4Dnv311svy5XG57QLkfKwRCDw8PpVLZ/mdbW5u3t/fdyUiSXLhwIQBg69atHRyNx+N1vmsUwzAYCO2Ipmk+n2/vXDgtPgB+HSZQq80CAVwU3G7g+W9fVix/K0x3Cg4OxnG8sLDQ8ufVq1djYmLuSENR1OLFi+Vy+W+//Qa7OiEIgiDHYYVAyOPxHn300XXr1jU2Nu7evfvq1avz588HAOTn52dkZFjSLF++/OTJk88+++y5c+dOnDhx8+bNnr8vBEEQBPWcdaZPbNy4cdWqVcOHD/f19f3jjz8kEgkAgKZpgiAsCTQaTWxs7GeffWb5c8qUKXe3GiEIgiDI9qwwfcK64PSJPgQOH7cvWP72BcvfvqxY/s61JCIEQRAEdRFcNNBptRkUBS03C1qKqlW1jZpmmb7NTBEGwuDGFAjZbl5cjzBxSIQ4JMVngDtHbO/MQpBr0Zi0uc35xbLSckVVg6ZZaVCpTCocZbAwppDl5sv3DnTzi/WISvCM9eF52Tuzzg8GQmfTqpMdrzpztvZijaou3jM6ziN6WsRkX563O1fCRBlsnK00qpRGdZO2ubyt6u+6y5tztnlxPUYFDp8cOtabd9dcNgiCrEdlVB+vOvNX9dkKRXV/r9gY96iM8EkBAj8hSyBkuZkpwkSa2gzKJk1zlar2XO2lL3N3CFluowKHjgsaFSQMtHf2nRa8R+g8cpryfi/JzGsuHNVv2Jh+I5K8++Pog0uGoqmbspLjlWdOVZ+Ncg+fHzMz2ad/J98R3iOxL1j+9tWl8i9rq9x9Y9/FhuxhfoMmho5J8kpgYA/eioei6SJZyZma88erTgcK/KdHpo/pNwJF4C0tAKx6/sNA6AzO11/+oeBXnVk/J3rauODRHLw7k6zNpPlEddbPN/axcfaSxIUDfQc88CWwIrYvWP721cnyL22r2Hrth3JF1ZzoaVPDJ/EY3dnIg6DIc3UX9xYflOvbFsTNnhiahiGuXu/BQHgLDIQ3Wou/zP1WTxgeT3hkZMBQtMcLPtGAzqq5sC3vBy+e5zPJT4WKgjpIDCti+4Llb18PLP8WXeuWq9/nNuc9nvDIlLAJDNQKt6LypIXfXf+5RSdbmvT4yIAhPT9g3wUD4S2uHAjl+rYvc3fkSQueSlwwMSSt5yHwdiRNHiw7tiN/18SQtMX9H71fExNWxPYFy9++Oih/kiZ/KzrwU+FvMyKnzI+d2b1Omg5cbsj96uoOEVu4atCyfm4B1j14XwED4S2uGQgpmt5fcui76z9nRExcGDeXjffWrnEKo/Lr3J25zdfXDH02xSfx7gSwIrYvWP72db/yL1dUvX9+k5gtWjloqb/At5fenaKp30sOfX/9l4zwCY8nPMLEXG7pShgIb3HBQFijqv/o4mcogr44ZEU/t06VUg/lNOV9eHHzUL+U5cmL2f97YQsrYvuC5W9fd5c/RVM/Ff62t/jAsqTFk0PH2iAPcn3b5pxtZW0VLw99LsHTtZbrgoHwFmsFQtJAmVRms5YkjRRpJAEAKAPFmChLzGCJGQjqEDut0ID+rejAjwV7nuj/yPSIKdbtC+2Y1qzbnL21sLXojZEvRYhD25+HFbF9wfK3rzvKv1nb8s75/zAxxitDn/fgutsyJ2drL27K3jIuKPXpxAWO0zQ0KswGmYky05ZKFWOiKBNl8HEGH2fwrNB6gYHwlm4EQrOG0DUb9S0mQ6tR32IyyExGuZmmaaYbg8HHUCaKszEAAGmiKBNlkJvMaoLrxxZH8sUxAreQ7gz3sgqprvX9C5vMJPHq8JV+fB+75OFk9dnN2d8sjJ87K2qq5RlYEdsXLH/7ur38z9Ze3Hj5i3kxD8+LediWF6ntVEb1xstfVqvqXhu+Kvy2q1Ub0zYa5IVqRbFGXavHWSjbg4mxUIz1T6Vq1hAmFUGRNFvMYLsz2R5MjieL48nkeLJYIgboSsnBQHhLB4GQImiT0mxsMxvaTAaZ2dBq0rcaDVITwADHk8X1YnE8mWwPFtuDyZYwcM59QylF0OpqnaJEI8tXIRjil+ruNVCEYDY90U/X/L3pytdzoqfPj51p3ylEjZrmN89+6Mv3fnnoszwGF1bE9gXL374s5U9Q5JZr32XVnF+f+nKMe6R9s3Ss8vSXudsfjZs9J3oa0qWo0kM0kBWoGrJkBpnJvb+bOIrvFsrDWPetrEgjZWwz61uNhlaTvsWkbzHqpUbCQHE8mBxPJtudxfZgssUMlpjBFDEw5r2PAwPhLWqZtumkAkEQmqJJA0WaKEJPEjrSrCYoM8V0Y1j6Nm9dd3gwOR4svNtNchooSjX1p1uNbebwuf62aR0aCOPmnK3XmgteH7E62j3CBu/4QGbSvDlnW25z/oZRr3igYlgR2xEMhPalVqvNOPHG2Q/5TO6rw1e5MR3iu2jUNG/4eyOPwX11+EoxW2SDd9RLjWV7GkgTFTDWwz3Brdv3kkgjpW8xGlpN+laTQWYytpmNCrNJYQYAMAQ4g4fhXAxjYQwBHjbTF8BA2E6j0KoK9SiKAgTBOSjKQHEuxuBiDAHeQSOvh2T5qor9je793UKm+fTq7cNKRfX6cx9FuUesGrTM6sOve8hy4bk0/vH0qPH2zovrgoHQvi5X536U+/m0iEkL4+fatPn1ICRN7sz/+VDFiVeHrbzneG8rasiS1Z5o6TfR03eEey+VAWmkTCqC0JGEniSNJE0BzyQhgIGwnb1GjRIGsvj7WoAg0YsCO2j+98SBsqPbrv24PHnxJJuMPeuGsrbKdaffHRucumTAQrjmk13AQGhHll/oq8NXDvFLsXde7u1q8/V3z386KTTtyf6P9sYyNDRFV/zeqKrUxT4VxBI/eLk4q4OB8BY7Tp+gKbp8X6O6Wpfw7xCca80MaM26jZe+qFbVvTnyZdtMkOi2BlnjJ9e+xFH8jREv8pk8e2fH5cBAaBcERX6W/U2etGDdwFWRPuH2zk5HFEblu+c/NRDGN0as9uR6WPHINEnf3FlDk3T04/16qTHwQHA/QvtDUCR8tp8okn9jezVFWO1iokRevuTwKj6T/9Wkjx08CgIABEz+x2PXB7r5/fvoizWqentnB4J6ncKofOHkazK9/KtJn/jyvO2dnQcQsYQfpb051C/lX0dWX6jPttpxaVD6az0AIPbpIHtFQetyhs9gRyEZPmx3ZtF3NTRlhVi4rzjz5b/WL0lc+MLgfzvOZKCOYQj2bMqS+bGznjv+yuWGXHtnB4J6UbmiatmRF/t7xW0Y9SqXwbF3djoFAchjcbPfTl376ZWvvsrdQVBkz49ZfbhZ12yMWhjoIHOsew4Gwp5BQPhcf9JI1R5v6clh1CbNa1nvH6k8+eWkj9OCRlordzYzJWz8O6Ne/fDiZ7/e3G/vvEBQrzhbe2H1ydf/NWDR04kL7DJTsCcSPGO2pW+qUdU/e3xNo6a5J4dqzVO2XFPGPR10v1kNfZHzfBJ7QXEkelFg0wW5ulrfvSNcb7n51KGVvjyvLyZ+ZK/J8j0X7xn91eRPjledef/CJjNptnd2IMhqaEB/d/2XzTnbPkpbPzYo1d7Z6SY3luC9MevGBY3+99EX/6o+172DmFRE+b7GqAWBDL5TbeoOA6EVMPh42Cy/4h9rSSPVpRdSNPXd9Z/fOPvBqkHLnkl5yiq7tNiRF9dj84QPjKTp+RPrZPo2e2cHgqzAQBjWn/3ocmPO15M+iZSE2Ts7PYIAZHb01I/S1m/P/+nDi58ZCEPXXk+D0l/q/UZIBP36Rrdw58FAaB3uCW5uIdyqzC70OTRppc8dfzVPWrgtfdMw/4G9lzdbYuOsN0e+NMx/4LIjq4tkpb33RhRB65qNihKN5T9dk7GrVyFQX0QRtF5qVJZpLd+7tt5Amnrxe2/SSp85tobL4Gwa966EI+69N7KlSEnY1vRPaZpecnhVsays8y9svCAndGTAeM/ey5u99O0miEMJfdg354NSn2Finu+DJ78frTj1Ze6O+XEz50bP6HP3GzqGAGRh/NxQUfDa02//O8nK8yD1UmNrnkqWr9Q1G1li5q3FCWlgUpkNbWacgwmCOG7BXFEUvzPfAtQnGOQmRZFGValT1+iMbWaWmMEU3loK36wm9C1GlpjpniDw6C/kW7WlcrX5+tt/f7Igbnb74rpOg4Oz1w57/q/qc2tOvzU7atqjcbMeOBWY0JE1R6QJy0OcZoDM7eA8QmtqOCuT31DHLw3uII3CqNx46ctadcNrw18IF4fYKmu9ouN5PNXK2nVZ7w32TX4m5cmez+c1yE3Vh6XKEo1HktDS/r77B2mQm9RVelWVtu2mhiZp9wQ3jyShWxDXkRb9sCbnnkeoazK2XFXI8lVmHSmO5gvDeIIgLteLdee3SQNNg0GWr2zJVbLdmcEZ3vwAK4TDPUV/7rqx940RLyZ5J9wvjROUf4uu9f0L/zUQxleHrwwQ+HWQsmJ/I0XQ4bM7SmNjcEL9LY4WCGmSzv24LHS6jzjm3l/PmZrz/83eMilk7JP9H2VgdliLwboeeCJqTNp3zm/UmfVvpa7pybKH9Wdaa0+0+KW6+4/26OS8Jb3U2HpN2XJNSZlor0Eir4EitnvfmJHSeU5QEd/NrCGkOQrpFQWhJz0GCD0ThfxATmcuZWiSbr7UVnNM6pEoDJnm0+2V8Q2E8eNLn9eo6t4Z9Yo3z6uDlM5R/jSgfy8+tPP67scT5j0cmXHPDip9iyn/s/LkNREONUYGBsJbHC0QAgDkheqqzKakF8PvaK+0GRSbrmypVNasGfpsnEe0vbJnXZ05ESma/u76z5nlx9ePfDnes8sfnDRRZb/U62WmmCf6sUTduXTQ1huar7S15Cp5fmyfYWL3eDcbbx7Se5yjIr6FBopSTdMFuaJEK4kXeA8WC0N53WjKkwaqZFedWUdGPx7IFHS51q5XN76e9X6EJPSFwctZD5rL60zlX6du+PDiZpqmXh763N1LedzcUSMI4gSMday7gzAQ3uKAgRAAkL+5wi/V3WOA0PInDehD5Se2Xvv+obAJjyc80ldmyndG50/EC/XZH1787LG42bOjp3Z+eWLSSBV8Xcn1ZofN9kPxHkUviqDlBaqmC3Jds9F7sNhnmMQuqyNal3NUxISWbL7S1nhejrNQn2ESz2QRxu7ZID4a1J5oab7UlrAipEsXT2drL2y8/OXi/o9Oj0jvTHrnKP92FE3/UXpoR/7uWVFTH42b1T6IXdtgKNxaPXBdZA9/g1YHA+EtjhkI5TfUNUekA14IAwBUKmv+c/krgjKvHvxMX78jeLcunYhNWukbZz/w4XmtGfocj/HgTawogr6xrZotYYTP8bfiTT691Nh4Xi7NUbgFc31HuIuj+H33DmJfr4jVNfrGv2XyArUkXuA73F0QZM2hLvVnWpsvtiWsCO3MZuhmivjm6ndZtRfeSl3T+f3O+nr531OztuW/2Vvq1I0vDFo2wDsBAFDyUx3Xlx0w1ppLlVoFDIS3OGYgBDTI/bjUL0OyV7//eNXpxf0fnRo+2cmGhlp09UQ0k+Yvr+64WJ/9xsgXH7CFKQ2Kvq8BAIlaGNAbo9QoE9VyVdn4t5zQkz7DJd6DxZ2pLh1NH62ISRPVkqtsOi8nDKTvMIn3YHH3dwntUPWh5rZiTcIzIR2vgdKoaX7r3MfuHPGaYc91aU/BPlr+nXG29uLnOdviPKKWhC6q/ko+aF1UT5vpvQAGwlscMxBSNHU88+/WXGXjpOolAxYKWW72zlFv6d6JeLb2wn+ufD03evq8mIfvd33QkCVrvaZMeCakt+/nqWv0TeflsusqcQzfZ5ike/el7KXPVcTaBkPTBXnLVaUwlOczXGKD5njJrjqUgYbPue9Yx5PVZzdnb10QN2dWdEZX9xTsc+XfJUbS9FPhb+qjpjCPoPELh7EdbEtUAANhOwcMhOfqLm279oMHy33W33MTloTx/R3u7LGibp+IUl3ru3//B0XRV4etvHt3GG2DoeDrqsSVoWyJje6nEnpSmq1ouiCnSeA9ROw1SNSNcRa211cqYsJAtl5VNl9qM6kI7yFi7yHi7o176gbSSF3dWBYy1cc94c7rUa1Zt+nKlmJ52esjVkeIQ7tx8L5S/t1GaMnL7xWfG/vXZVXOEwnz0sPG98a+ht0GA+EtDhUIc5rytuf9ZCSNSwYsHOo3sO6vVr3UGDHP0bdS6omenIgUTe+6sfe3oj+WJz81MWTMP88T9LX/lAeM8/BK6f50i25TV+uaLrbJ8lVuIVyvQSJJnJujDRC4nYNXxDRFK0u1zdmKthtqUSTPe4gtmoB3U1frb35bPWB1ONPtn4ubq83XP7jw3yF+KcuTn2TjrG4e2bHLv+fqT7fqGo0R8/2LZKXfXPteqmtdnDA/LSjVQW70wEB4i4MEwiuNV7+7/rPSqH4i4ZH2s8SsIXLeLx30RpRz7Nd1Tz0/EUvbKt49/2mgwO+Fwf+2TDSsPiLVS43RiwKtlMfuIE2ULF8lzVZo6vTuCW6eSUJhOM8BF9Rw0IqYBuoafctVRes1JUvE8EwReSWLeukuYCfVHJXqGg3RT/QDABgI49a878/UXHh5yIrBfsk9OayDlr/15H5YGj7P3y341tC2nKa8b/N/Upu0i+LnpgWNtHvrEAbCW+wbCEmaPFNzfveNfWaKWBA3e2xQ6h3LFBXtrBFFC3yGOskShXezyoloJs07r+/OLD+xPPnJ0eIR1zaWDVgdbrOus46ZVISlQjfIzO4Jbh6JbsIwnuNMQ3SsipgG6hpda75KlqdCGYjHAKFnspDj2c3GlnVRBJ37QWnEfP8KfsVHFzfHe8Y8O/DpLo2LuSfHKn9rU1Xqyn6tT15z5xja7MZrPxT8KtW1zouZMTl0XLfb0z0HA+Et9gqEKqM6s/z47yWZvjzvebEPD/MfeM/b7G031TXHpInP9+0V6ztgxROxRF7+4cXPJhRNSgiPjpvWnRs2vcogN8nyVK35Kn2rURwtcI8TiKL4OMd5roi7jTJRynKtrEAtv6HGOahHf6F7fzeen8PdGq+90lR4pHJHwvZVg5cN9bPOGveOUP69p3R3PdeP5T/63rMmbrQW77qx77r0xpSw8TMi0ztegqeXwEB4i40DIQ3oPGnhwdJjFxqupAYOmxWV0fE9dpqis98piVsSxHXSBaCtWxEoqrR520s/S/jvw7FT5sXMcMwl6EwqQn5DLS9UKcu1PF+2OIoviuTz+3Hs0nFqt4qYBtomg6JEoyjWqCp1/ECOJE7gHufG9nDExSJoQB+pOLX16g9Li5fGjQsLGGK1KtuJAyFpoK5sKE555QFrqjVqmveVZB6pOBnrHjUtYvIQvxQctd3VIQyEt9gsENapG05UnTla8RcLZ2WETZwUmiZg8jvzwuojUtJAhs7w7e0c2oV1K4LrX1Z6DRLRsebPc7ZVKmqeHfi0ta7cewNF0KpybVuxRlGiMbaZ3UK4wjCeWwiXH8ixWd+pLStimqJ1TUZVhU5ZoVWWazEWKorkiyP5oki+A84wa1csK9uUvQUAsHLQUj+1f9F3NQPXRVrrC3LiQNh0Qa4o0UQ/3q8ziY2k6a/qcwfLjtWpGyYEjx4fMjpKEt7bOQQwELbr7UBYo6o/V3fxr+pzMr18TL+Rk0LTuvoFG+SmvP9WDH4zygGHWvScFU9EdbW++IfalFcjLAV1uTF3c/ZWH5738uTFIaIgq7xF7zFrCGWFTlWuVVZoDa0mnh9bEMQV9OPwAzlsCbP3xkn2dkVsUpo1dQZ1jU5do9dU6xluuCXeC8N5DnITtwMtutateT9mN15dMmDR5NCxlpsXhd9UeQwQeg+2zm17Jw6E+ZsrAsd73m/zgPupVzcerTx1oioLAciYoBGjAoZFuod1dXZm58FAeEtvBEITabrecvNSQ86F+mw9YRgZMGR0v+GJXnEP3K/rfvI2lQc95C2K6FQLsm+x4ol4c2eNKJznO9K9/RmCIv8oPfRDwa8jAoYsTpjvwXXv4OWOgzRS6hq9xhI8avWkkeL5sXl+bK4vm+vN4vqwrHhn0boVMWmk9FKjrtmoazRoGgzaBgOgAT+Aze/HFfTjCIK4fWXxHa1Zt/vGvj9KD0+LmPxY7Gwu45+V2xSl2op9DckvR1ilcnbWQGhSmq9+UjZ4fXS3m87F8rIzNefP1l40EIYhfilD/FJSfBJv/yKsAgbCW6wVCA2EsUhWkie9kSctuCkrCRUFWb68SIkVLmfq/2rVy0wOtY+XtVjrRNS3GPM3Vw56LRK9ayksjUm768beA2VHHwqbMD92Zp9bpsesJbUNem2DQddk1DUZ9M0mBAMcTxbbk8l2Z7IlTJaYwRIzmG6MbkxY7F750xRtVhMGudmoMBvlJoPMpG816VtMhI7kerM43iyeD8sSuR2/2XcHA2HcX3ro5xv7hvkPWtz/Ua+71moAAFz7tLzfRC9JnBXOW2cNhA1ZMm2jwSpzoGtV9Rcasi815NxoLQ4W9kv2Tkjwio33iOEzeT0/OAyEt3Q7EBoIQ4WiplxRWSIvL5KV1qjqwkQhCZ4xid7xA7zirXvl4sS9o9Y6Ecv2NDDd8H6T7juKoVUv/6Hg11PVZ6eGT5oTPV3MFvb8Te3FpCb0UqNBZjLITEa52SA3GRVmk4rAORhTgDMs/3ExnIfhHAznYBgLxVgozsYACiytSQRDLJNTNRoNn88HAFBmmjJTAADSSNEUTehJykSTBpIwUISeNGsIQkeataRZZTapCbOWZPAwloTJEjHYEgbbncl2Z3I8WSwRow8tL3cHA2H4s/TI7pu/J3rFLU6YHyS87zzU1mvKhrOy/s9aYWSyswbC/M0VgRO8xNHW7MQykabC1uKrzdevt9wokpV6cj1i3CMiJWFh4pBwUUj34qIVy78PrCPVQ0qjqkXXKtW1Nqib6zWNdeqGGmWdwqgMcgsMEwdHiMPSQ8dHiEN6b4wiW8JkixnKcp0owgoXQc7HrCFa85Qpr3S0BrcHR7Jq0LLH4mb/VPjbogPLJ4SMmRczw5vnWLujdRJTgDMFuDDszpPBrCZMGsKsJkxqgtCSZh2pbzGRBpI0UKSJIvQkoAGhJwEANEmTRgoAQNM0giAAAJSBoAwUAICxUARFcA6KMlGMheFsFOdiLDGD789h8DGGG4MpwBl8zJmuyVQm9e/Fh34vOTjAO+E/Y99+4B1lFchImAAAIABJREFU9/5uVZnN6mq9dTe7cBompVkvNVq9smJizCTvhCTvBAAARVOVypoiWWmxrOxU9dlKRQ0HZwcJA/0Fvn58Hz++jzfP04vrIWKLbLaETd9uEcrU8sPVJ1EU1Zn1JEUaSIOBMOrMerVJozSqlUaVwqDgMDieHHdvnqcP38uP7xsg8AsSBvjwvG25SlDdX60GZ+wdtcoVWf2ZVl2jMeKRzvbDyPVtvxb9kVl+fJBv0ryYGbYZn+aYnLVF0kkNmqY9RX+cqMwaGThkfuysu7eTvZ+6v1oNLcbwuT3t+nPK8rdiv2jnSXWtNcq6OnVDvaapUdMs1bZIda1qk1rIchOy3NxYbm5MPpfBYWEsHoPLxlmPJzwCYIuwHUVTGpMWQRAug4MxMA/MnY2zeAwun8ETst1ELDcRW9S+vaQdeSS65W2qCJvp60xX4tbSfKktfE4XfnUSjnhZ0hML4+ceLDv2etYHHhzJrKiMUf2GO8IXDdkARdPZjVf3lRy82VqaET5hZ8bn7pyujQL1HijK+aA0ZLqvEy9/2G2tecrACbaeHe/F9fDiegz0HXD7kyRNKgwqhVGpMqqVRpWeMBgJo47QszHrT8vu23UHn8FbkrjQ7muNPpBlTISqUnd3h5iLU1fraAq0L2bYeTwGd17MjDnR0/6uu/x7SebmnG2TQ8dmhE8MEDhbsxtqJ9O3Ha44mVl2jM/kPRz50Fupa1lYd6bwMwS4MJzXek3pPcRplz/sHpOa0DVbv1+0ezAEc+eIu3qV0z19OxD2IZJYQdtNNQyEd2i+1OY9RNztMRoogqYGDk0NHFqnbjhYduzZ46/4830mh44b02+EVYalQY7ARJrO1185UnGyoKVodL/h60e+HOXe0/5w7yHiupMtMBDeoe2GWhTJd5zVdG0GBkIbEUcLSn+pD86wdz4cCWmiWvNVyS/fuapvNwQI/JYlPbFkwMJLDTlHKk59mfttsk/i+OBRQ/0G2nFRYKgnCIrMbc47VXX2XN2lSEnY5NCx60eusda3KY7ml/1ar2s2cr3h6fGPtiK1JLaPzVCyChgIbUTQj2PWEMY2M0vcx+Zm9Z7WPKUwlHf7LnE9hCHYcP/Bw/0Ha826rJrzmeXHP7q4eZBv0qjAYUP9B/IYXe6AhWzPRJqym66drb34d93lQDe/Mf1GLhmwyOr9YwiKeA0SSy+3BU/1se6R+y6apBUl2rCZrnhzAQZCW0GAOJrfdlPtM1xi76w4ipZcZS/tUcVjcNPDxqeHjVcZ1X/XXTpRlbXx8pfR7hFD/QcO9RvY+eGFkM206Fotyzldbb4eKQlLDRx6vxnx1uKVLCzcWh2c4dN3Z09al6pSx/ZkMgSuGBRc8TPbizha0HJVAQOhBaEnNdV68eJOrerbbW4sgSUiGghDTlP+hforvxX9iQBkkG9Sik9isk//PrdUjTPRE4ZrzQW5TXlXmq7J9W2DfJPSgkauGfZcz3cK7AyuLxtjoeoaOKHwlrYitaSLi4s6DRgIbUccwy//rYEi6G4spuV8ZPkqURQPu2tNtV7CxtkjAgaPCBgMAKhW1l5pvHas8vTHlz735nkO8I7v7xmX4BnTV5Yz7dNURnVB68186Y08aWGlojrGIzLFO3Ht0OciJeG2nNpr4Z4obM1XwkBoIb+h6fx0XicDA6Ht4ByM68tSlWtFUU64AHdXtearvAeJ7PLWQcLAIGHg7OipFE2VyMvzpIXHq05vurKFhbPiPaJiPKJi3CPDxSFwlI1VEBRZoagqkpXeaC2+IStp1cliPCL7e8YuS3oixj2C2a35D9bi0d/t5rc1IbB3FACjwmzWEIJAF70mgIHQpsQxAvlNNQyEhJ5UV+qiF913QUjbQBE02j0i2j1iXswMAECtqv6mrORGa8nJqqxKZY0v3ztCHBIhDg0Th4SKgvv0Aqe2pDXrKhRV5W1VpW0VpfKKalWdP98nyj08zjN6dvS0UFFQtzdysTqeHxvBEE2dnu+qAaBd2021OJrvshcEMBDalCiCX7an3t65sD9ZgUoYwXO0dT0C3fwD3fwnhqQBAAiKrFRWl8orStsq/q6/UtFWhSJoqDion1tAkFtgPzf/ADc/L66n7XvzHE2rXl6nbqhV1deo6qsUNVXKGo1ZGyzsFyYKjpSETQkbHyZy6La1R6Jba54KBkJFqdZlbxACGAhtjB/INraZzRqCwXfpkpflqzwHOHQDC0exCHFohPifPQpk+rYqZU2VsrZKWXO27kKdqkFhVPnxvf34vn4Cbx+etw/Py5vn6cXzELEc+qN1j9ask2pbGrXSZq20SSNt0DQ1aJrr1Q1snBUg8A8SBgS6+Sd79w8R9fPmefbeXqxW554oLNpZE5zhbe+M2BUNlGXaEBeeSeLS1bHtISjiFsJVlms9Ep2wruwkykwpy7WRjwXYOyNdY1ntKcUnsf0ZA2Fs0DRZ/mvUNF9tzm/Wtkp1LQbC6MX1cOeIPbkeErZIwhFL2CIRWyhiCUVsoZAlYOPWXyyxh0ykSWVUK4yqNoNCYVS2GZRyfZtM3ybTy1v1cqm2BUVQL66HN8/Lm+fpy/eO8Yj04/v4C3z7+uxMvj8b0MDFZ9brmgwYG3XlKc4wENqaMIKnLHPpQKgs0/IDODjb0VeIfSA2zgoVBYXete+PkTRJtS0yQ1uLrrXNoJTp5JXKmjaDok2vUBpVKpOaomk3Jl/A5POZfB6Dy2NwuAwun8nj4GwmxuQzeUyUwcJZTIzJwpgYgnL+f4NMAfN/7i5rdVo1or39GZ1ZR9IUAMBAGAmKIChCTxjMlNlAGPVmvZE06cx6rVn3///XqU0alUmjNqlJihKyBEKWm4gtlLDFYrZQwhGHioIkbLEH192L62H17cUdhyia33ZT7cqBUFGmFUW49MAFGAhtTRTOL7pQa+9c2JP8pkYc7cx3I1gY03K78X4JjKRJbVSrTRq1SaM163VmnY7Qq00aA2FQGJT16kYjaTKRJiNpNJFmkqb0Zj0AgAa0xvQ/YY+iKBT9n/usXAYXQ1AAAAtnMVAcQzEuzmGgDDbO4uBsJs4UsPg+fC8OzuExuDwml8/gCZh8N5aA43iNVJuRxAgasmT+Y3px8r6DU5ZpPRz7VkVvg4HQ1nh+bEJLmJRmptBFOyLaitQxvTyP3sGxMCaL697zaYtOuR+e7QkjeMU/1ZJGytFGb9kIDVQV2jCn2y21S1zyi7cvBLiF8ZTl2gendEZ6qZEmaJ6P67Y/IEeDMf+vvTuPkqK8+wVevdfS3bMw09M9PRvLAIIKCKKIVxRE8SgR405cIiaiSI5rztHgi74mOWrwzVGTq0YCekxuomLkRPDeaDyKkgAiqCwRGJZh1u7Zp7u6qnqrrvtHk8k4a890ddf2/fzV01PT9dB09beqnt/zPGZXFd17PKJ0Q5QRaRFsTqvdkDOr9UEQKqBwChM6YdAg7DkaKTrLpZ2iQjCEorOcPUcMGoSh41yBsTsICQShIgpqnb3HDRqE3elxuwBqUnyWq+coq3QrlBE6yRVMMfrinQhCBdAehxhLxXoTSjck38R4im3gC6ciCEFdKI/DZDbxwZjSDck7iQif5gsmaXsMTPYQhEowEa5qim3glW5HvoVPcs4KyqAlCaBuRdOcPccMd1HIB6M2xmLw+T0IBKFS3BOZcL3hgjB0gis0/E0YUKf0AF+lW5Fv4XrePRGHJIJQIe6JtAGDsPcEuuVBpQomM+FTvJSSlG5IXoXreVeN0e+LEghCpTgrKaE9JsZSSjckf5JRUWiPYXZjUCeb02ovsHEtUaUbklfh07x7IoIQQagQs9XE+MhIk6B0Q/InfJJ3VdNYlBhUq7CWMVQ5d4JNJnmR9hh3brk+sgUhz/ONjY2p1LCXOKIoNjQ0RKPGOuEagctgd0d7j0cKa9EbAepVMIUJnTDQaMJQPe+eSGNQLyFXEL700kt+v3/x4sXTp08/cuTI4A327ds3adKkpUuXlpeXv/nmm7LsVOvcNXS43kCnn6ETGK4EqlYwmQmf5iXRKN2EbD2H+6JpMgRhQ0PDz372s127dp04cWLlypUPPPDA4G1Wr179yCOP1NXVffTRR2vWrOnq6sp+v1rnnsSwpwWDdM4neTHWnXBWoIMQ1MtKW8gJduN0WKBSpo8MQfjWW29ddtllZ511FkEQa9as+eSTT9ra2vpvcOTIkW+//fZHP/oRQRDz5s2bNWvW1q1bs9+v1tkYi81l4dsMMYY3dIJzTaRNFtyFAVUrnML0GmMQRSqR4oMxF4rXCIKQZfWJ+vr6qVOnph97PB6Xy9XQ0FBWVtZ/A7/fT9NnTj1qa2tPnz493KulUqmDBw8Gg8H0j5MmTSoqKhpuY0lMJpobRbNWS36Y0mTP1w22pFbbL/J8nM7ojLL7m6SzxBRvOp7rJhlK5u8/ZIgpTAW/Fsum9Waysabff7ZVIoulZNtJpRsyZiazxeafJO9ryhCELMuWlPxnKS+GYUKh0IANKIoaYYP+YrHYT3/6U5vtzBJFTz/99KWXXjrcxnx7a+ovL4675crjz+6uL7Ec3qF0O8YplUoJmZ2F9HbdVOr6vOt4MNdNMpTM33/IUEpyRNpv7/zz/zYRo/dZaPr97+XOtYgFXX/eqXRDxsxEOZkf/hdBEJFIRpVNNE1bLKMsAy5DEHo8np6enr4fe3p6+l8Opjfo7e3tv8GMGTOGezWKoj788EO/f9hFTfuzWqvsj/xm1H+karkaheNbWnyP3KR0Q8Ypw/XwxHjq9PqjVY/9F8ZOyAvrEeZC23PHC37wP4x/9JXCNP3+h//Y5J3u8sy7VemGZEWu91+G05lZs2bt3bs3/fjw4cNms3ny5Mn9N5gxY0Z7e3tLS0v6x717986aNSv7/eoA4yejHXExrvNh9WwDz/hJpCBogqvGEOOa2AbBVYUOwjNkCMIbb7yxsbHx+eefP3DgwIMPPnjXXXcxDEMQxGOPPfbss88SBFFWVnbDDTesXbv24MGDTz75pCRJV111Vfb71QGTxUR7HbqfzIKtx+wVoBnuiXT4tM6DMMmLSV6kSjGU/gwZgpBhmI8//njXrl2rV6+eN2/er371q/TzFRUVPp8v/fiVV16pqqr68Y9/XFdX9+GHH1qtRp/svI+zimYbdX7UhU+jShs0w11Ds3oPQrZRcFZQGErfR55AOuecc957770BT65du7bvsdvtfvFFLVe15IyritL5iqASwTYKU1ciCEEbqFKHGE/FehOOQpvSbcmVSCPvxH3RfrRa8qQbzkqKbdTzAF4uGLUxVix4BpphgOVC2SYBIwj7QxAqjPY4kpyY5ESlG5Ir7GnejfuioCm6XyUt0iTgirA/BKHSTISzgmT1O6tTuJ53oVIGNEXf3YSxngQhETq+8TsOCELlOStpHU9vGMYVIWiNs5Lig7GUTsc1sbgcHARBqDxnFaXXwtEkJyYjIl2GKm3QErPNTHsdEZ2Oa4qgg3AQBKHynBWkXocSsk2CsxJV2qA9zipKr/dpIumjEvpBECqPLLKL8VSCTSrdEPmxqNIGbXJV0not5+ZaBMaPo/I7EIQqYCKcflKX92EijZjGCTRJrx0WsZ4EYTbZ3RjO9B0IQlVg/BTXosPTz0gzbsKAJul1XFOkWcD62IMhCFXB6ScjzXq7IkSVNmiYiWD8ZKRZb6enXEvUWTH6whpGgyBUBaaCiujuipBtFFzVGDgBWuWq0mE3YaRFcKKDcBAEoSrQHkeCTSajuroPE2lCpQxomC67CSPNUQZXhIMgCNXBRDA+vQ2iYFEpA1rmqqIi+roiTESSqXiKLLIr3RDVQRCqBaOz0YQSwTVHcRMGtMtRaCNMRKw3oXRDZBNpjjorSIzrHQxBqBZOP6Wnnnm+PWZzWayMRemGAIyfs1JXF4Vci8CgZHQoCEK1YCooPRWORppQpQ2a56yk9DQhfqQl6ixHB+EQEIRqwXgd0e54KqGTeX4jzTj3BM1zVuhqgC+uCIeDIFQLk8VEldj5YEzphsgDw5VAB5x+UjczjorxVDyUpEpRKTMEBKGKMOUk16qLu6MSwbWgUgY0z15gM5lN+qiX4VujlNdhMqNUZggIQhXRTRAKnTErbbHSqJQBzWMqKH2Uc3Ot6CAcFoJQRZhykgvo4pBrwaBd0AmnXiZa41qjtA9H5dAQhCqimyvCCEYQgl4wehnXxLVGGVwRDgNBqCI2p9Vs0UOHRKRZQKUM6INO1s2WCC6IIBwWglBdGL8eLgqx8ifoBlmsh3Wzo91xK2WxUui2HxqCUF0Yn+aDECt/gq7oYt1srjXKoINweAhCddFBNyFW/gSd0cG62eggHBmCUF1o7QchhtKDzuhg3WwE4cgQhOpCexyx3oSmJ1qLoIMQ9EUH62ZzAQThSBCE6mKymKhSbU+0xrVGnX4ccqAfVKk9EU6KUa2enoqxVCKcJEswudqwEISqo+l6maQgJnmRLMYhB/phMpsor4MPavWo5ANRqgyTq40EQag6jI/U7iF3pjgNRxzoi6ZPT3FfdFQIQtWhfSQX0OqtUfTJgy4x5WREs0HIB2O016F0K1QNQag6tNfBa3bGUa41yqCDEHRH0+OauAAGEY4CQag6jkJbSpQSEU3OZMG14IoQdIgpJ/lAlJCUbse48MEYptseGYJQjWgvqcXCUSklCe0x2otDDvTGSlmstCXaFVe6IWMWZ5OEJNldmOlpJAhCNWK8Di2uxyR0xO0FVosDHyrQIY3eHeUDWH1pdPjOUiPap8krQlTKgI4xfkqTQRiMoYNwVAhCNWJ8mqyX4VqjTDnmlAF90ugVIReIomR0VAhCNaJ9JKfBnnlcEYKOObUZhLg1mgkEoRpZKYuFtGhuhV6uRUAQgl6RE+wJLpkURKUbMhYSwbfF6DJcEY4CQahSjE9j9TJJTkwlJEehTemGAOSGiaB9pLb6LKLdcSuN9XhHhyBUKc0dclz6DgwmVwP9Ynwkp6kqNgylzxCCUKUYrxYPOdyBAT1jfA5eU92EGEqfIQShStFaKxzlgxhKDzpHe0ltdVjwQZSMZgRBqFK0xyF0xKWUZipHuVZUyoDOMeUkH4hpqJybD8QYnJ5mAEGoUma72e62amZKp3RxGs49QdestMVsN8VC2ijnllKS0BWnPFgcdHQIQvWivQ6tzC8T7Y5bKRSngf6dmX1bC6KdcUeB1WzDl/zo8B6pl4am3kZxGhgErZ0VerkAuu0zhSBUL8br0MpS9XwAxWlgCIx21s1GpUzmEITqpbUrQhxyoH8aGuCLQu7MIQjViypzCF3aKBzFfIZgEHT6qBS1cFTiijBjCEL1MltNjgKb0KH2wtFUUor2JGgPDjnQP7PV5Ci0Ce1qv1UjiVK0O0GV4qjMCIJQ1TRROCq0xchiu8mC2dXAEDTRTSh0xB1FNrMVR2VGEISqRmuhXoYLooMQDIT2aeCo5INRBvdFM4YgVDVN1MugZBQMhdHCRGscZhkdCwShqmniihB98mAotE8DHRY4KscEQahqtMcR7U6ovEQNo+nBUMgJ9gSbFGMppRsyEoydGBMEoaqZLCayWNUlamIsleRFshjzGYJRmMwmyuPg29R7VKaSUqwnQZXgqMwUglDt6DJVH3J8IEqVObAeLxgK7VX1KmlCe4ycYEMhd+YQhGpHeUk1ByEXxDovYDgqr2Lj22J0GY7KMUAQqh1dpuqeefTJgwExPgffpt4rQqyJNlYIQrVTeeEo+uTBgGivqsfU4/R0rBCEaqfywlEe022D8TgKbal4ShRUWjjKB3FrdGwQhGpnspgcRTahU40zjib5VCop2QtsSjcEIL9MBF3miLar8aiURCnWk6BKUTI6BghCDVDt3dFoWxyzV4Ax0T5SaE8o3Yoh8O0xshglo2ODINQA1dbLxDoSmM8QjIn2OmIdagxCAZUyY4cg1AC6TKUjKIT2BCplwJhoLxltU2MQon5tHBCEGqDqW6M49wRDYnwOQZV9hHxblC7DUTk2sgWhIAj79u1rbGwcboPOzs79+/c3NzfLtUfjoDyOaJcaC0ejHQkEIRiTzWk1mUxxNql0QwbigzEE4VjJE4T79++fPHnygw8+eP755z/++OODN7jlllumTp26Zs2aOXPmXHfddfG4Gs+kVOvMotgqKxyNhxJmi8nmtCrdEABlkB6b2iZaSy9MT2Jh+jGSJwgfffTRBx544B//+MdXX331yiuvfPvttwM2WLVqVTAY/OKLL06dOnX48OHXX39dlv0ahwqXqucCMbIMJdpgXI5Sm9qOSqE9RhZjYfoxkyEI29raPvvss7vvvpsgCL/fv2zZsnfeeWfANldccYXdbicIwuVyzZw5s62tLfv9GooKuwn5tqijBJeDYFxUmU1tK/TybbgvOh4yfJE1NTW5XK6SkpL0jxMnThyhp/DEiROfffbZU089NdwGoiju3Lmz79XOPfdcj8eTfSO1ji4juw6HlW7Fd/CBGOnFFSEYF+mxhw6p7KhEyei4ZBqEzzzzzKFDhwY8OW/evIcffpjn+fTVXhpJkhzHDfki3d3d119//aOPPjp79uzhdhSPx3/zm984HGdOah577LGLLrpouI0FQbDb7RaLJcN/hYa5xUirEIlElG7Hf7AtXHEtpaomGQ3HcSYTboIpJuVMcK1ChI2oZxmycEukcCZjkKMyw88/TdNm8yj3PjMNwoULF06ePHnAk36/nyCIsrKy3t7eVCqV3llXV5fX6x38CqFQaNmyZZdffvm6detG2BFFUe+88076lUdlsVgMEoR0jXS8J8hQjFomjJCIWGeysNrtdDqVbopxSZKE919BkiRZqbAt6XAUqWWWwXhnoLi6gHYa4qJQxs9/pkF4ySWXDPermpqawsLCPXv2pC/d/vnPfz700EMDtuE47nvf+97s2bOff/75cbfVyPoKR1XSARDtjltpi8WBcahgaIzPwQejKglClIyOmwxfZA6HY82aNWvXrv3b3/72+OOPd3R03HDDDQRB7Nq1q6qqKr3NDTfc0NjYOHfu3I0bN7722muff/559vs1GlUVjqIrAoAgCLpMResxoWR03OSp+lu/fv2ECRNefvllv9//2WefkSRJEERZWdnNN9+c3mDRokWzZs2qr69P/9hXCwOZo8vSa4G6lW4IQRAEH8TqSwAE7XWETg5dEpF/KBkdN5MkqWu+koqKii+++CLDPkIDFcsQRMdXoa7D4el3VCrdEIIgiGP/p7lompOaZnG5XEq3xbhYlsX7ryCWZYke68l3W2c/PLB+QhGNf2snCKJqmVHK7GX8/KOPRzPUdWs0gCWwAQja6+DbY1JKFZcTfBuOynFCEGoG5XFEu+JqmHFUSklCR5zy4JADo7PYzXanNdatimUoMMvouCEINcNsVctS9dHOuL3AarHjwwNA0D6HGuaXSZeM4vR0fPBdpiUqWaGXC8QYLEwPQBAEQdBeUg1HpYCF6bOAINSSfxeOKowPRmkEIQBBEATBeFVxRYiS0WwgCLVEJeeefCDKoE8egCAIgqB9pBomxMfQ3mwgCLVEJWtQcDjkAP5NJetmo2Q0GwhCLVHDIZdKSrGeBFWKdScACIIgzFYTWWzj2xW+VcMFcGt0/BCEWnKmcLRDycJRPhilSu3okwfoQ/tIZZeqP3N6ipLR8UIQakx6kl8FG4CuCIABFJ/sQmiPkRNwejp+CEKNob0kp+ghh1lGAQZgfKSyhaOY6SlLCEKNob0OZW/CcAFcEQJ8h+Ll3FwwhkLubCAINYb2KlyrzQeiNK4IAfqhSuyJSFKMpZRqAIb2ZglBqDFUqT0WSqYSyhxyyaiYFESyCCWjAP2YCKpUyW5C9NxnCUGoMSaziSqx823KHHJ8IEZ7HQS65AG+i1auik2Mp+JskpxgU2Tv+oAg1B4FOyT4YBSzjAIMxijXZ8EHY7THYTLj/HT8EITaw/gUq5dBpQzAkGgfybUqdnqKktEsIQi1h/aSSt0a5VqjTDmCEGAgppzkWgVFdo0OwuwhCLWHVm62ewxXAhiS3W0lTKY4m8z/rjG0N3sIQu0hi+1JXkxGxTzvN9aTMNtMNqc1z/sF0ATG5+BbFThDRYdF9hCEGmQiqDKHkPd6GS6A+6IAw6J9JJf3IEwKohgTHYUoGc0KglCTFJnSiWtFySjAsJQ8KlExmh0EoSYpcsjxAcxeATAsnJ5qF4JQk5hyBWq1cWsUYAS0zyF0xPO8XCgfiNI4KrOGINQkppzkW6NEHo+4VFKKdmHBM4BhmW1mR2G+lwvlArgilAGCUJOstMVsN8V6E3nbI98Wo0rsZiv6IgCGxeS5XkZKDyLE6Wm2EIRaxZTn9ZDjW9FBCDAKxpfXMb7RrriNsVopS972qFcIQq1iyvPaM88FMGgXYBS0j8zn9IdcK9ZEkweCUKvyPGgJlTIAo1Lg9BRHpRwQhFrF5PfcE7dGAUZFFtuTgpgU8jTrE49KGZkgCLWKLnNEuxP5WaE3Hk6mUhJmrwAYhYlgfCTXkqcz1AhOT2WCINQqk8VEleZphV6uRXD6qTzsCEDrnBVUpCUfy1CI8VQinKRK7XnYl+4hCDUsbzNZRFqiTj9OPAFGx5Tn6YqQD8QorMcrEwShhtE+Mj+z3XMt6JMHyAjjz1MQYii9jBCEGub0k5G8HHKRFoGpwK1RgNExPlLoiqeSOZ/2iWsRGNynkQmCUMOYCopryflEa2IUXREAmTJZTFSJPQ8V3ZGWKIJQLghCDbMxFovDHO3J7dyGXKtA+0h0RQBkyOmncn6rRiL4AIJQNghCbWP8JNec20MOJ54AY8L4SS7HhaN8e8zmslpJTK4mDwShtjn9Oa/V5lpRMgowBnmol8GIJnkhCLUtH4dcs8DgkAPImNNPcYHcdt5HWqLOCpyeygZBqG3OCirSnMMrQkmU+I44g3VeADJmIc02p1XoyOFkFzg9lReCUNscRbaUKMXZZI5en2+LkUW4Lpv0AAASoElEQVQ2sx2fE4AxyPXQpkgreu7lhC84zXOWU1zOLgojTYKzEieeAGPjrMzhrZpod9xsMdld1hy9vgEhCDXPWZHDc0+2UXBVIQgBxsZZSUUacxWEXEvUiQkuZIUg1DzGT+WuVjvSJDgr6Ry9OIBeOauoSIsgpXJSMBNpiWKmJ3khCDXPWUFGcjOUMJWUhPYYU45KGYCxsZIWu8sqdORksguuWcCIJnkhCDWPKnUkuGQu1gLlWqNUqd1sw4cEYMyclXSO7o5GmgXcGpUXvuO0z0Q4K3LSIRFp5J1VuC8KMB6uKopt4mV/2VhvgpAIRxFWyZYTglAPXNU02yD/IYdKGYBxy1G9DNsguKpxeiozBKEeuKooNhdXhBg7ATBeTj/JBWOyr8cUaeSdOD2VG4JQD1zVtOxBKMZSsd4EjTllAMbFbDdTJXY+KHMhG9souKoRhDJDEOqB3W01W03RbjlL1CJNAlOO1ZcAxk/2u6NSSkKlTC4gCHXCWUWxDXIecizuiwJkx1VFsU1yHpV8MOYotFkprL4kMwShTriq6EijnPUykUYelTIA2XBWyVzFxjYK6CDMBQShTrjkviIM1/PuiYyMLwhgNIzPEetNJDnZxviyDTxKRnMBQagTziqKC0QlUZ4SNaEjTphMGKsEkA2T2eSqosOnZbsojGBEU24gCHXCYjeTxXYuIE+JWrieK5iME0+AbLkn0eF6eYJQjKWiXXHGh8nV5Icg1A9XNcXKdO7JnubdNbgvCpAt90QmfIqT5aXYRoHxkyYLCrnlhyDUD/dkJnRSniAMneLdk3BFCJAtdzXFBaKyDKsPn+QKJuP0NCcQhPpRMJkJneSIrI+4BCcmwkkMpQfIntlupsscslR0hxCEOYMg1A9Hoc1iN/PtsSxfJ3yKc02kMZQeQBbuiUz23YSppBRpFlw1uE+TEwhCXSmYzIRPZtshEa7n3RNxvAHIwz2RDp3KNgjZBp72OiwOfGPnBN5WXXFPoUMyBCGHIASQi3sSzdbzWa5Wj/uiOYUg1JUz3YRZEGMpPhjDWCUAudicVpvbygey6rMIn+QRhLmDINQVsthuspiEjvHPvh06wbmqaaxKDyCjomnOnmORcf+5JEpsI+/CfZqcwfed3mR5UdhzjC2a7pSxPQBQNN3Ze4wd95+zjQLlcVhJzLWdKwhCvck2CI9GCqchCAHk5J7MsE2CGEuN789DJ7mCSbgvmkMIQr0pnObsOcqOr2c+2h1PxVOMF3M4AcjJYje7KsdfyNbzLYvT05ySJwhTqdTu3bs/+OCD3t7eETaLxWL79+/v6uqSZacwJEehzVFgG99KFD1HIkXTXQQGEALIrXC6s+foeLoJk5zIB2MFU3BFmEMyBKEoitdcc80999zz2muvTZ8+/dChQ8NtuX79+vnz52/fvj37ncIIime4ur8dT4dEz1F0EALkRNG0cXYTdh9hC6YyZivOT3NIhiDcvn37iRMn9u7d+9e//vWee+558sknh9xs7969O3funDNnTvZ7hJEVz3R3/ys81r+SRCl8iseJJ0AuMD4yvXzEWP+w+1/h4pmuXDQJ+sgQhFu3br3uuusoiiIIYuXKldu2bUsmkwO2icfj995778svv2yxoPAp51xVVIITo91jO+TC9TxVarc5rTlqFYChmYiiaa6x3h2VRKm3jiuejiDMLRm+9Zqbm+fNm5d+XFVVlUwmg8FgRUVF/21+8YtfLF++fPbs2aO+WjKZfP/994uLi9M/XnjhhQNeqj9RFEVRttWf9aRourPzUMh3cXHmf9JxoLdwhnNM7yfef2Xh/VfWWN//whlM4B/dngsLMv+T0HGOLLWbaRP+owfL8P03m80m0yg3ljMKwgMHDjz88MODn9+8eXN1dXU8HrfZzixlnn4Qi31nDoWDBw9u3br1yy+/zGRfoihu27YtfX1JEERZWVlpaelwG8diMUmScJU5mLPW0fFFuPj8TO9zSimp80B42uryAf93I4vH42PaHuSF919ZY33/6Uk2fkuUbefsBZlegXQeCrmnUfhfHlKG7z9JkvIEYU1Nzfr16wc/n44or9fb2dmZfqajo4MgCJ/P13+zX/7yl36//6mnniIIoqmpacuWLS6X6/vf//6Q+3I4HBs3bvT7/Zk0zGQy2e12BOFgjnPIxvc67SaHlcrozek5wlIljiL/GM5VCYIQRZGmMduFYvD+K2sc7/+Ecwq4o/HCy9wZbS0R4WPNZ91dRdMY0TQEGT//GQVhQUHBokWLhvvtwoUL33///XXr1hEE8emnn86ZM2dA41atWtXQ0JB+7HA4SktLS0pKsmgzjM7iMBdOd3Z8FfItzOjuaMdXIc95Y0tBABir0vMKT28L+i/L6AswdJKzOMyMDymYczL0Ed5xxx3PPffcI488MmPGjCeeeOKFF15IP79w4cLrrrvu0UcfvfLKK/s23rRp06WXXnrJJZdkv18YmffCovr3g5kEYSqe6v6Wnfg9bx5aBWBkhVOYOJvk22J02egLXwf39JRdWJSHVoEMVaNFRUV79uyx2+379u3bvHnzzTffnH7+vvvuGxx4999//9y5c7PfKYyqsNYpxlORptFH1ncdZl3VlM2FelGAHDMRpbPdnV+HRt0wyYk9R1jPvMI8NArk+e6rqqp65plnBjx52223Dd7yjjvukGWPMDoT4b2gKLinZ0rlKGsqte/rKT0PxxtAPpTOLTz6RlPlFaUm80gVHO37e4tnujLs44csYa5RPfPML+r8JjTyVL+RZoFvi5XMRgchQD44KyhHsa3jq1EuCoN7ur24L5ovCEI9s7usBVOY9n0jTQDb+Lf2iiWlmMAJIG+qrvQ0/b19hJnxQyc4KUW4J2KapzxBEOpc1TJP40ftCW7oYaeRZoELRMvm48QTIH8KJjP2gmEvCqWUdOqvgeplHkx/nzcIQp1jfGTJrIKG/9c25G9xOQigiOqryho/GvqiMPCPbhtjRW9FPiEI9a/6Kk/34XCkcWD5aNveHr49hstBgPxzT6SpEnvDBwPPUONssunjjsnf9w35V5AjCEL9s1KWmmu8dX9u7j/zffe3bMP/bZv542pcDgIoYtptld3fsi2fdfY9k+TEuj82lV1QRHlGH2UIMsLQMUPwzCtMCuKBF09NWuF1FNkjTULTxx0zflRNleJ4A1CGlbbMXF1z8DenkkKqYDJNSMTxt1tKZhVUL/Mo3TTDQRAaRfn/muCeSB9/q8VsNTF+6qxVVa6qUcYXAkBOOQptZ6+uCfyzu+mjjlgoMfn68uIZWHFJAQhCA3FWUHMenaJ0KwDgPyiPY9J16BFUGPoIAQDA0BCEAABgaAhCAAAwNG0HYWdnJ8uySrfCuAKBQDQaVboVxtXU1CSKQ88ZBLmWSqUaGxuVboVxxWKx1tZWuV5N20G4bt269957T+lWGNedd965f/9+pVthXIsWLeru7la6FQYVCoUuvvhipVthXN98882QCxyNj7aDEAAAIEsIQgAAMDQEIQAAGJpJkoZdE0sR11577cGDB83mjBKaZVmbzUaSZK5bBUPq7e1lGMZmsyndEIPq6uoqKirK8GABeUmS1NXVVVJSonRDDCqZTLIsW1Q0+poB27dvP+uss0beRnVByLJsR0eH0q0AAAA9qKiosNvtI2+juiAEAADIJ9xUAQAAQ0MQAgCAoSEIAQDA0BCEAABgaJpZj7Cjo2PTpk0dHR1XX3314sWLB28giuKbb7556NCh6dOn33XXXajpl1FHR8e2bduOHDlSWFh400031dbWDtgglUr9/ve/7/vx7LPPvuiii/LbRj3bsWNHXV1d+rHVal21atXgbVpbWzdv3tzb27tixQpM/SWvjRs39i8qHPzxrqur27FjR9+PK1as8HiwynxW2tvb9+/f39TUdOmll06dOrXv+ebm5tdffz0cDl9//fUXXnjh4D+Mx+ObNm06fvz47Nmzb7vttgwHF2njipDn+QULFhw7dqy6unrlypVvv/324G3uv//+l19+uba29o9//OMPf/jDvLdRz9auXfvhhx/6fL729vbZs2fv2rVrwAbJZHL16tVHjx49derUqVOnurq6FGmnXr3xxhtvv/12+r2tr68fvEFvb+/8+fNbWloqKiquvfbaDz74IP+N1LH6+vpT//bQQw8dOnRowAa7d+/esGFD3zaxWEyRdurJkiVLnn766ccff3z37t19T3Z2dp5//vkdHR0+n++qq676+9//PvgPb7311rfffru2tvbFF198+OGHM92fpAWbN28+//zzU6mUJEl/+tOfzj333AEbBAIBh8PR1NQkSVJ3dzdJkidPnlSgoTolCELf49WrV69atWrABukjn2XZ/LbLKO68885f//rXI2zwwgsvLF68OP341VdfXbhwYV7aZTj79u2jKKqnp2fA82+88cby5csVaZJeiaIoSdIFF1zwxhtv9D353HPPXXXVVenHL774Yt9nvs+RI0domg6FQpIk1dfXUxTV2dmZye60cUX4+eefL1261GQyEQSxdOnSgwcP9vT09N9g9+7dU6ZMqaioIAiiqKjovPPO27lzpzJt1aP+c/dEo1Gn0znkZr/73e9eeumlr776Kl/tMpA9e/Zs2LBhy5YtiURi8G8///zzK664Iv146dKlu3fvHnIzyNKmTZtuvPHGwsLCwb9qbm7esGHD5s2bMR+ILIa8pTngc75z585UKjVgg/nz57vdboIgampqKisrv/jii4x2l3WD8yEQCJSWlqYfT5gwwWq1BgKB4TYgCKKsrEzGpaqgz+7du7du3fqTn/xk8K+WLFnS1dV15MiRxYsXb9iwIf9t07GqqqrS0tLu7u5nnnlm/vz5PM8P2KD/59/j8aRSqWAwmPdm6pwgCG+99daQHbQFBQVnn312KBTaunXr9OnTB987BVkM+JwnEonOzs7+GwSDwf5B4PF4MgwCbRTLWK3WZDKZfiyKoiiKA6bMsVqt/VcoTSQSo86pA2N17NixG2+8cePGjVOmTBnwK7vd/vHHH6cf33rrrZdffvmaNWsYhsl7G/Xp6aef7ntw3nnnbd68ee3atf036H+ApB/g8y+7v/zlL0VFRZdccsngX61YsWLFihXpx/fdd99///d/v/vuu/ltnSGM+jkfdxBo44rQ7/f3BXv6gc/nG7BBS0tL348tLS3l5eX5bKHu1dXVLVmy5Nlnn73ppptG3nLBggXJZLK5uTk/DTMUm802f/78wfUy/Q+QlpYWm82GyaBlt3nz5rvvvjvdQTOCiy666NSpU/lpktEM+JzTND3gNvW4g0AbQbh8+fJt27ZFo1GCIN59993Fixenrza++eabhoYGgiAWLVrU2dmZXi39+PHjR48e7buVDNlraGi48sor169fP2BJ6C+//DL9sRMEoe/J7du30zRdU1OT50bqWN/by7Lsjh07Zs6cSRBEPB7/5JNP0mVKy5cv37p1a7pfcMuWLVdffbXFYlGwwfpTX1+/c+fO22+/ve8Znuc/+eST9HVJ33+QJEkffPDB2WefrUwr9W758uXvvfde+ppvy5Yty5cvTz//5ZdfpgPyyiuvPHDgQPpMcc+ePTzPL1y4MKOXlqXCJ9eSyeQVV1wxd+7c22+/fcKECbt27Uo/v3jx4p///Ofpxy+88ILP51u1alVlZWXfkyCLZcuWMQwz99/uvffe9PNz5sz57W9/K0nSxo0bZ86c+YMf/GDZsmVut/sPf/iDou3Vm+Li4uXLl69cudLn811zzTWJREKSpPSRf/r0aUmSotHoxRdfvGDBgpUrV5aUlHz99ddKN1lv1q1bd/XVV/d/5ujRowRBdHV1SZK0bNmyJUuW3Hbbbeecc05tbW1DQ4NCzdSPhx56aO7cuQzD1NTUzJ07N/2dz/P8BRdccPHFF99yyy0ej+fw4cPpjc8999xXX301/fiJJ56orq5etWqV1+t95ZVXMtydZlafEEVxx44dnZ2dixYt8nq96SePHTvmcrn6Ln4PHz6cHlA/Z84c5VqqQ3V1dSzL9v3ocrnSQ1z/9a9/lZaWpnut9+/fX19fX1BQMG/ePIwmltfp06e/+eabeDw+derU2bNnp59MJBJff/317Nmz070giUTi008/7e3tveyyy/rXC4As6urq3G533zcPQRDRaPTgwYNz5861WCzd3d179+7t6enx+/0LFizAbB7ZO3nyZG9vb9+PtbW16VrQ9I2QcDi8ZMmSCRMmpH97+PDhsrKyvo/9/v376+rqZs2aNWPGjAx3p5kgBAAAyAVt9BECAADkCIIQAAAMDUEIAACGhiAEAABDQxACAIChIQgBAMDQEIQAAGBoCEIAADA0BCEAABgaghAAAAwNQQgAAIb2/wGJDUkqnfj94AAAAABJRU5ErkJggg==", "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": 11 } ], "cell_type": "code", "source": [ "ρ = scfres.ρ[:, 1, 1, 1] # converged density, first spin component\n", "ψ_fourier = scfres.ψ[1][:, 1] # first k-point, all G components, first eigenvector\n", "ψ = ifft(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": 11 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.1" }, "kernelspec": { "name": "julia-1.8", "display_name": "Julia 1.8.1", "language": "julia" } }, "nbformat": 4 }