{
"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\n",
"a Gaussian potential."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"struct ElementGaussian <: DFTK.Element\n",
" α # Prefactor\n",
" L # Width of the Gaussian nucleus\n",
"end"
],
"metadata": {},
"execution_count": 2
},
{
"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::ElementGaussian, r::Real)\n",
" -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)\n",
"end\n",
"function DFTK.local_potential_fourier(el::ElementGaussian, q::Real)\n",
" # = ∫ V(r) exp(-ix⋅q) dx\n",
" -el.α * exp(- (q * el.L)^2 / 2)\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": [
{
"output_type": "execute_result",
"data": {
"text/plain": "2-element Vector{Main.##340.ElementGaussian}:\n Main.##340.ElementGaussian(X)\n Main.##340.ElementGaussian(X)"
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"x1 = 0.2\n",
"x2 = 0.8\n",
"positions = [[x1, 0, 0], [x2, 0, 0]]\n",
"gauss = ElementGaussian(1.0, 0.5)\n",
"atoms = [gauss, gauss]"
],
"metadata": {},
"execution_count": 5
},
{
"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": 6
},
{
"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.143571689200 -0.42 8.0\n",
" 2 -0.156034324285 -1.90 -1.10 1.0\n",
" 3 -0.156768366303 -3.13 -1.56 2.0\n",
" 4 -0.157045437468 -3.56 -2.31 2.0\n",
" 5 -0.157052804997 -5.13 -2.68 2.0\n",
" 6 -0.157056391680 -5.45 -3.74 1.0\n",
" 7 -0.157056406247 -7.84 -4.40 2.0\n",
" 8 -0.157056406918 -9.17 -6.14 2.0\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380294 \n AtomicLocal -0.3163465\n LocalNonlinearity 0.1212606 \n\n total -0.157056406918"
},
"metadata": {},
"execution_count": 7
}
],
"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": 7
},
{
"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.05568147166176689, 0.0, 0.0]\n [0.05568027211637458, 0.0, 0.0]"
},
"metadata": {},
"execution_count": 8
}
],
"cell_type": "code",
"source": [
"compute_forces(scfres)"
],
"metadata": {},
"execution_count": 8
},
{
"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": 9
},
{
"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+gvaeTAAAgAElEQVR4nOzdd3hTVf8A8HPvzR7N6N57t5S2bChQdqEMmYKAIvLCD1FBVFAcKPo6eeUVRREQVMSBDIWyQSjIbAstLd17j+yd3PH7I7x9eRmlI03S5HweH580nNycnNyc7z3nnoFQFAUgCIIgyFmhts4ABEEQBNkSDIQQBEGQU4OBEIIgCHJqMBBCEARBTg0GQgiCIMipwUAIQRAEOTUYCCEIgiCnBgMhBEEQ5NRgIIQgCIKcGgyEEARBkFOzu0CYl5dnMpk6mZgkSbhEnA0RBGHrLDg1WP62BcvftixY/nYXCCdPntzS0tLJxAaDgSTJXs0P1AGtVmvrLDg1WP62BcvftixY/nYXCCEIgiDImmAghCAIgpwaDIQQBEGQU4OBEIIgCHJqMBBCEARBTg0GQgiCIMipwUAIQRAEOTUYCCEIgiCnRrN1BqDeZSJBsYIqUVCVKoCTwEgCAQN4sEC4AOknRujwQgiCbIegQJGcKlZQDVogMwAGChAEBHBBhACJEiIcWD1bCyxpx9SkA7+Uk6fqyb+bKG8OEiVEQviAgQEGCsqV4HIzKJKTZUqqvysyIwidHYQE8RFbZxmCnEWrHhyqIg9UkldbKG8OEitCfDhAxAQaHBAkyGoFJQqyXEUluiLjfdEnQ5BwAfx59i4YCB3N6XrqiwLi72bqiUD0uUj0p9GoiPnwlFocXG6mDlSRg/8k+omRNXFYmj8Cf3AQdJ9PP/30xo0bDz5PEASGYV06lMIIatSUxADcWMCDBcYyERoKAADNADTfkywCgFAKyA3ggIH6Sgc4NODHBZ5sZ/95jh49euXKlb1xZMTeFq328/O7du2ar69vZxLrdDoGg9HVc9FRHa+lNuYQahN4tR86JwTldvoix0SC3yrIf+WTJhJ8Mgib5NfZn5tKpeLz+d3MLtRjsPytIyUlZcKECVFRUbbOiFPLycnJz88/cuRI+zMWPP9hi9ARFMmpNVeJKjX4YAA6IxBFu3jdSEfBU2HoU2HoH9XkmqvE51zw1XAszMXJrz4h6L9SU1NHjBhh61w4NTabnZ+f30sHh4Ml+jYTCd67SY48ik/0Q/Nm0mYGdTkK3mt6IHp7Ji3NHx36J775NknYV2cBBEFQr4CBsA/Ll1GD/sCvtZA3n6CtjkMtMgSUhoLVcei16bSMGnJMBl6vgcEQgiAHZ5lAuHfv3vDwcA8Pj6VLl+p0ukclO3bs2IABA44ePWqRN3VyO4vJMRn4i7FoxkSaL9fC3ZghfOTMZNoEP3TAYfxkHYyFEAQ5MgsEwsLCwueff37Pnj3FxcXV1dWbNm16aDKFQrFu3TqpVCqRSHr+ps5Mi4OF54mtBWRmOm1JRG+16VEEbOiP/jqWtvQisfk23P0YgiCHZYFq9LvvvpsxY8bw4cNFItGbb765a9euh45EXbt27erVq93d3Xv+js6sRk2lHMUxBFydRosS9vp4lpFeyJVp2E9l5LOZhBFGQwiCHJEFAmFRUVFCQoL5cUJCQktLi0wmuy/NuXPnysvLn3322Z6/nTO73EwN+RNfGIZ+PwpjW2vArz8XuTiVJjOA9JO4ymSlN4UgCLIaC9SmUqm0fTKHi4sLAKCtrU0sFrcn0Gg0L7744uHDh5FOzAfVaDR+fn7tf/74448zZsx4VGKnmkd4uBZdk0X7dgg+3tugVlv73fcMAS9n0UYdIQ+MNLmz7rb4NRpNZ75TqJfA8rcOknSQzpCWlpakpKS6ujrznzt27GhsbHz77bcflf7EiRMHDx789ttvrZXBxyAIQn1P3dfJ85/D4aDoY5p8FgiEYrFYqVSaH5sfuLq63ptg48aNAwcOVCgU2dnZGo2mqqqquro6MDDwoUfjcrn5+fmdnFCPYZiTBMLP88nPb5NnpmAJYoat8rArFWzMISafR8+kYebhORRF8Xg8W+UHguVvHY+tRvsKkiTb2trMj3U63XvvvZednd1B+okTJ7766qsFBQWxsbFWyeBjYBh27wlvwfPfAoEwPDy8oKDA/Dg/P9/V1fXe5iAAgCTJ27dvL1++HABQVVW1Z88eDMPefPPNnr+1M6AAeOMGcaSGujwN87P06NCu2piE8ejkqAzi7GQskAfbIhBkA9XV1bm5uaGhofv37x85cuSYMWOKioqOHj1qMpmmTZtmDloajSYjI6OgoIDNZk+dOvXBSHbgwIHk5GQPDw8AwLlz5zw8POLi4r7//vuZM2dqNJrMzMy5c+ciCPLUU099/fXXX375pQ0+pxVZ4EpnyZIlBw8evHXrllar/eijj5YsWWJurm7atOn48eMAgM2bN2f9R2xs7MaNG2EU7CSCAssuEhcaqcx0ms2joNkr8ejqOHR0BlGlgtMqIMgG8vLyVq5cuWrVKpFIhKJoRkZGWloagiAcDictLS0zMxMAkJOTk5mZGRAQQFHU2LFjs7Ky7jtIRkZGamqq+fGuXbvOnTsHAFizZo1EIiktLd24caP5n1JTUzMyMqz32WzEAi3C+Pj4Tz/9ND09Xa1Wp6WlvfPOO+bni4uL/f3970vs6+sLV0fsJCMJFv5FKIzU6cm0zi8cagWrYlAUgDHHiIzRSDT8MiHnk9NGrbtBWOe9hnggm5Lvv/ujVCoPHz4sEAgAABEREbt27RozZgwAwNPT88MPPxw5cmRKSkpKSoo5MZPJ/O677wYMGHDvEW7fvr148eLHvntERERVVZVSqTSP/3BUlqlfly1btmzZsvue3Lt374MpDx48aJF3dHg6HMw+izNQ5M8JNKb93QNdGYMSFJjyFyNzKmUnTVUIspogPrKun5V+lm6shzwZGRlpjoJKpbK0tHTTpk0ffvghAEClUrW0tAAAGhsbV6xYUVhYSKfTdTrdg12jKpWKw+E89t3NN+FUKhUMhJC1aXAw7RTuzUH2jMRo9nqf/oVYVKUjxh9HL0yhebBtnRsIsiIxE4zzteX1H5N5d3M182jB7du3t0/RNg/tefPNN6Oiosxj9b/66qtjx47ddwR3d/f2eW4sFkuv17f/k16vZ7Huhl+JRIKi6H3jHx2PvdayTkyDg6kn8WA+8sMo+42CZi9E4vNCkPHHcanB1lmBIKfEYrFSU1P37dsn+g+j0QgAaG5uDgkJQRDEaDT+/PPPD75w4MCBeXl55scJCQmnTp3Ccdz855EjR/r3729+nJubm5CQ0B4XHZV9V7TOR2kC44/hUUJkRwrWk30krGZjEjbWB5l2Ctfits4KBDmlHTt2HD9+PCkp6YknnoiLi/voo48AACtWrHjjjTeeeOKJAQMGPHQ22pw5c8yDGQEAy5cvFwgEERERarV63Lhx2dnZ7StlnjhxYvbs2Vb7LLYCu0btiMIIJp3AE8TIV8OxvhAE79o8BFt2kZh+Gs+YSGPAKysI6mWTJk0aPXp0+59BQUFXr16tqqpqbm4ODg729PQEAKSnpxcWFpaXl0dERPD5fHMz0cPDo7a21vyq1NRUg8GQl5fXr18/JpN54MCBlpaW8PDwvXv3DhkyxJxGq9UeOXLk8uXL1v6EVgfrLXshM4Dxx/EhHsjXI/pSFAQAIAB8MxxjY8jSTLiDIQT1Ojqdft/YewRBgoODhwwZYo6CZh4eHkOHDnV1dWUwGOYxLyiKtt9KRBBk69at165duzc9hmFeXl7tz9y4ceONN94wzzV0bDAQ2gVzFEzxQj4f0seioBkNBb+OwarU1KvXrDSmHIKgHho+fPh9o/1feukl82BUs1GjRjnJAtEwENqe1ADGHcdTvZHNg+1vnkSnsWngz/G0E3UU3LMJgvqod955RyQS2ToXNgADoY1JDGDcMXycD/JpX46CZiImODEJ+6KA3FcOYyHkOOxtINjx48fHjx/fSwc/evToypUrO0hQVFQ0YcKEh+6113fBQGhLMgOYdBwf54t8PKjPR0EzPy5yfBL28lXiVL1D/U4gp1WmpGrU9nUyx8XFrV69ujeOTFHUunXrVq1a1UGaqKgoFot15MiR3siArcBAaDMtOjA6A5/oh3ziKFHQLEaI/DaWtug8niu1r+oDgrqqUQsmHCe82PZ14x7Hca1WCwAwGo0ff/xxfX39xo0b169fX11drdVqP//881dffTUnJ8ecuK2tbceOHWvWrNm0aVNJSUn7QWpqat57770NGzaUlZVt2bJFoVAAAC5cuMBms2NiYgAAmZmZ5gVI9+zZU1NTo1KpNm/ebH7tokWLvv76ayt/6l4FA6FtNOlAagb+RBDy/gCHioJmI72Qr4Zh6ScJe7uUhqDOU5nAlJP40khUyLR1Vv5Xfn7+F198AQAwGAzr169funSpn5+fQqEYP378kiVLAABisTg1NbW1tRUAcOXKlbq6ukGDBjEYjBEjRpSVlQEAWltbBw8erNfro6KiVqxY8frrr0ulUgDAsWPH2idmnDx50rwezRdffFFRUSGXy9966y3zP40ePfqvv/7S6XQ2+PC9A84jtIF6DTX2GLE4HH2jv8NeiMwORus1IO0EcWkqTWRn9QgEPZaJBHPO4oPckQ390RMP/Kuxukj++1fWyQkjJE74xPIOErz//vsDBgx49tln3dzcxo0bZx4IeurUqfPnz8+ZM2fq1KlTp06lKEoul9fV1f38889vvfXWzp07R40a9c9//hMAkJSUFBcXZz5Ufn7+tGnTHpsld3d3NptdVlYWHx9viY9oezAQWluVihp3nFgRjb4S77BR0OylOLRWQ804jZ9Mo7EcsN0LOSwKgGUXCQaKfDX84ScuzTNAOPdF62QGZXM7TmAOYyiKurm5tYc0Dw8PiUQCACgoKHjuueekUimfz29qakpPTwcAlJSUJCYmmlPGxMS0r1yq1WrZ7E4tHMzhcFQqVbc+kD2CgdCqShXU+OPEq/3Q52McPAqafTIIe+o8seg88euYvrFiHAQBAN7KIooU1LnJtEfN6kVZHIZ/uHUz9UgY9t9obV5xGwCAIIh5YOdLL730zDPPmPdFf/31183RUSgUtq+4rdFozOvOAADc3d3NfaQAADabLZfL24+s0+naYyRJklKp9N7J+32dU1THduK2lEo9Rryd5CxREACAImDPSKxNT70MJ9pDfcQ3heRvldSRCTSOQzQT2iOWXC7/7bffzE9OmTJl3759DQ0NAID2ITAAgMGDB+fm5pofJyUlnTlzpr3Z9/vvvycnJ5sfl5SU8Pn8kJAQq32K3uYsNbLNXWuhxh/H/zUYfTbCucqciYFD42nnGqjP4ER7yO79WU1uukmemIS5O8p2C2vXrn3uuefS0tJSUlKSkpLMT44bN2758uWJiYlBQUEkSbLZbPPehHPnzm3fhmLy5MnTpk2LjIwsLi5euHDh4cOHt27dan55RkbGnDlzEMSBOnkoO+Pr61tXV9fJxFqtFsfxXs2PRZyqIz32Go/VkLbOiIUplcpOpqxTk4E/m34qI3o1P86m8+UPdcbfTaTHXmNW6/2/0xEjRly8eNEmWXooHMd1Oh1FUeYuyvbn5XK5yWQyP1ar1Xq93vxYIpHk5ubqdDqdTqfRaNrTkyRJkmRBQYFAICDJu596zpw5Bw4caE+j1Wqjo6MPHTrU/gxBEPHx8UVFRb32+R7uyJEj6enp9z5jwfPfIRr/9u1QFbnyb+LAWNoILwe6gOoiXy5ybBI2JgN3ZSIT/Zy3HCC7Vaak5pwlfhhFS3az9/MTwzDzfUEEQe5dEe3eZUK53P8OsRGLxWKx+L6DLFu2LC4uzmQyffPNNxs2bGhv3n300Uf3TpZns9ksFksoFLY/U1ZWtmLFisjISIt+JhuDgbB3/Tuf3HybPD2ZFiey919Xb4sRIgfG0WaewTMm0gbYfV0DOZU6DTXuGPHxINR5rtKeeuqp7OxsBEF++OGHYcOGtT8fEhLy0ksv3ZvymWeeCQgIaP8zIiIiIiLCehm1ChgIewsFwLs5xP4K6uJULJDnLL+ujg33RHal0KaexM+n0yIFsEwgu6AwgvSTxAux6MIwJ7p/P3r06Hs3NezAiy9aaaKIDTnRF29NOhzMO0v81UBdmkqDUfBe6QHIBwOxtBNEgxYuOgPZnhYHU07i432RtY4+rxfqAPzuLa9JB0Zn4EwMnEqDi6o8xLMR6PIodOJxQmawdVYg52ZePibMBfmk7+/9AvUE7Bq1sJsS6onTxLOR6FuJcAb5I61LQNv01JST+OnJNC48ByFbICnwzAWChiI7U/rebtgkSX777bc3b9709/dftWrVvYNZoG6ALUJL+q2CnHQC/2ww+jaMgo/zyWAsWojMPI0b4FR7yBZeuELUa6lfx2C0PlgLLl++/OjRo2lpadXV1ampqQQBf0U9Aq/GLQMnwYYsYn8ldTqN1k8Mg+DjIQB8m4I9eY5Y8BfRRysjqO96/QZxo5U6M7mbq+DeaSvZcuMbS2fq4fp5xKxKfu7eZ+rq6vbv319bW8vn86dPnx4TE3P8+HHzIqJQ98BAaAHNOjD/HM7EwI0ZNFd4U7DTMAT8lIrNOI0/m0nsGQUXI4Ws5J+3yKM11Pl0mgu9m0cIFPitHdzRNu4WxKPfv+h2cXFxaGgon88HACAIkpSUVFxcDANhT8BA2FNnG6inLxDLItG3ElFYlXcVAwUHxtImn8RX/E1sH9H3btVAfc7n+eT3peT5KT26ZuXSOZHiMMtlqmu0Wq1SqWz/U6lUuri42CozjgF2SHWfiQQbsoinLxA/jMLeSYJRsJvYNHBkAu2OjHrpCrzPAfWubXfILwvIs5Mxb46ts9IzZWVl58+fBwBUVlZeuHAhNTXV1jnq22Ag7KYiOTXsTzxXQt18gjbGB8bAHuHRQcZE2rUW6uWrMBZCveXbIvKTPPLsZMyP2+d/sP3799+0adOwYcMGDhz4wQcfhIXZrHnqGGDXaJcRFPiigPzwFvH+AOwfUfBKwjIEDHAyjTb+GL72GrEZTuqCLG1HEfnBLfLcZCyI3+ejIACAy+WePXu2qqrK3d393mVFoe6B9XjXFMioEUfwI9XklWk0GAUtS8gApyfTMhuptdcIuOoMZEE7isj3b5HnJmOhLo4QBdsFBQXBKGgRsCrvLC0O3rhBpGbgSyLQs1NoDvaLshNCBjiVRrvURL14GcZCyDK+vEP+M9ehomB0dPSKFStsnQuHAgNhp/xeScYewGs0IG8W/R9RcFhMLxIxwenJtJsSavklgoTBEOqZz26TW/LJ81McJwoCAMLCwhYuXNjtl1MUpdfrLZIT85aHFjmUbcFA+Bg3WqlRR/EPbpG7R2J7R2NebFtnyAm40MGJSbQyBbXoPGGC29pD3fVONrGrmLwwxQG3f8nPz8/Nze04zY0bN0pKSh58/tatW1FRURbJRkxMzM2bN7v9chzH9+/fT5K2/5HDQPhId+TUrDPEzDPE4nA0ewZttLej/ZbsGY8Ojk2iqUxg1hlCD0eSQl1EAbD6KnGkhspMp/n2/TGiD/rll19+/PHHjtNs27bt3i127ZBOp5s7d649rA8HR40+RJ6U+uct8nwj+Wo/bO9ojA0LyRZYGDgwDnvmAjHpBP7HeJqAYesMQX2EiQRLM4lKNXVuCk3oiKdNYWHhuXPncBxfv359aGjosmXLcBzfvn17VlaWr6/vqlWrvLy8rl27lp2dXVNT09rampiYOG/evIceymg0fv311zdv3gwMDFy1apW7u7v5+fPnzx8+fFgqlSYlJb300ksEQezZs+f69es4jqekpDz99NMo+shG1ObNm1NTU/fv39/S0jJ79uy0tDTz8yUlJTt37pRIJKNHj164cCGCINu2bQMAbNiwAUXR559/3t/f39JF1VmwRfg/ztRTk0/iaSeIge5I+Tz62ngURkEboqPgx9FYghgZdRRv1No6N1BfoMHB9NO43AhOTXLMKAgAEAgE3t7enp6eycnJ5s3iFy5ceOTIkTlz5hAEkZycLJfL3dzcXF1d/fz8kpOTg4ODH3WoWbNmnTt3bu7cuSqVauDAgRqNBgDw3XffLV68ODExccGCBW1tbSRJ6vX6kpKSadOmzZw5c9euXe+//34H2duxY8eCBQtiY2MnTpz43HPP/fnnnwCA8vLyoUOHenh4TJ8+fcuWLa+//joAwNxDm5SUlJycbNvhr7CaBwAAuRH8WEpuLyIRANbEo4fGoUw4k80+oAj491Dsw1xy+BE8YyIWLXTAbi7IUlp0YOopPE6EbB/Ru8u4Kyu1JfvqevEN7iEM54bN9b33GR8fn+joaL1eP2fOHABAeXn5n3/+WV9fLxKJJk+enJWVtWfPntWrVwcFBcXFxZnTPNTt27czMzPr6+t5PN7kyZOvXr26b9++ZcuWvfnmm99///348eMBAJMmTQIA8Hi8Tz75RKfTNTc3r1u3bv369W+//XYHeV6xYoV5OI9Go/nss8+mTZu2devWWbNmvfLKKwCAsLCwxMTEd955Z8yYMQCAWbNm0endXfXVQpw6EJpIcLqe2ltGHq8lJ/mjXw7DRnkjsKK1Q68noP5cMCYD/3UsbaQX/IqghyhWUFNOEovDrbEVKM+PFbciqJff5C6U8ZiQXlZWFhQUJBKJzH8mJSU9dIzMg0pLSyMiIng83r0vlMvljY2NQ4YMuTelRqOZP39+UVFRYGCgTqdrbGzs+MgJCQnmB/379zdnprS0tH1Z8JiYGAzDqqurfX19H3kI63LGQKjBwZl68nA1daSajBQiC8PQL4fRxXDXCPu2MAz15SJzz+KfDsYWhcEufeh/nG+k5p/DPxqEPR1ujXMDpaMsVxt3vLbPWxAKhXK5vP15mUzWfquvY0KhUCaT3ftCHx8fHo/HYDAkEol5dwuz3bt3YxhmDmk3btyYMGFCx0duz49cLjdHaJFI1P6kVqvV6/XtkdseOEuFYiDApSbqg1vk+OO490+mL++Qia7IrZm0v6fS/i8ahVGwT0j1Rv6aQns3h3zjBpxiCP3XrmJy/jl83xiadaKgPXB3d6+pqTE/7tevH4qi+/fvBwDU1tYePHjQPD7l3jQPNXDgQKVSmZGRAQCoqKjIyMiYNGkSjUabNm3a+++/bx7M2draSlGUTqdDEAQAQBDEli1bHpu9b7/9FsdxkiS/+eYbc+dqWlraDz/8YI6FW7duTUpK8vT05PF4bDa740xah2O2CAkK1GuoMiUokFH5Miq7jSqUU7EiZKQX8lIsljoe4Trm53Z80ULk2nTazNP4rDPED6Mxvo3vLEA2hpPg1evEsVoqM50WLnCiPvN58+YdOHAgKChoxIgRe/fu/fnnn5955pl33323paVl/fr1I0aMAAA8++yzixcvDgoKmjlz5r/+9a/216IoymKxAAB8Pn/fvn1Lly7l8/mtra3vv/9+UlISAOCrr75asmSJubtVr9cXFhY+/fTT33//fXR0NEmSc+bMEQqF5kOxWKyHDh8NDw+Pi4szGo1BQUHffPMNAGD+/PlZWVmRkZEikYjJZO7btw8AgCDIpk2bxo0bh6LowYMH2ztUrQ+xt3UB/Pz8rl271sm+4waF/t/FNARBNCagNAGFEbToqCYdaNBS7iwkzAVEC5F4MZLoivR3Rbq3FTXUAZVKdW//idUYSfDiZeJSM3V4PBbmQCuGdJWtyt9OSAxg3lmcjoKfx/TuANGUlJQPP/zQHF3sWWtrq1gsxrAu13QtLS2urq73vdC866Gnpyfyn4ETDQ0NYrHYHEQ7EBUVtXv37sTERLVa7ebmdu8/mUwmpVLp6ura1RwCAI4ePfriR98Me/8PLzb4bDAGLHr+9+2WEYZQQjpAUcSXA1wYQMgA7izUmwN8OAgc9unAGCj4ZgS2vYgccQTfmUJLD3DeWOi0ctqoOWeJ2cHIPwdicENns07eGnyQh4fHg09yOBwO53+2bfTx8en8MVks1oMhk06ndy8KmnlzwCQ/pDeWIe7bgdCFDl7rh2CYs9wYgO61PApNECPzzhHXWpGNSbA2dCK7isk3sohtw7BZwfC3b3fee++9kJCQ3jiymIks7J2BcvA0gvqwIR5I1gza1RZq/DG8QWtfnfxQb1CbwOLzxJZ8MjOdBqOgfZo7d66np6etc9E18EyC+jZ3Fjg5iTbGBx1wGD9RB2OhI7sloQYcxhkYuD6dFulMQ2Og3gYDIdTnoQh4MxH9dQxtxSVizVW4SLcDogD4121y4gl8YxK6MwUu/wtZGAyEkINI8UJuPkFr0IJBf+C5Utg0dBy1GmricfxQNXltGu3JUFhlQZYHzyrIcYiY4Ncx2Kv90AnH8Q9ukbjttzmDeur7UnLAYXy0N3p+Ci2ID7tDoV4BuxggR7MoDE31Rp67SByuIneNxPqJYe3ZJ9VpqBWXiHotOJ1Gs+2XmJycPG/evAcnA1AUhcDFia1Fq9UOHz68lw7etyfU63Q6BoPRjQmkkEXY84RuCoA9JeT6G8Q/otAN/TGHXE7Bnsu/J0gKbC8i38kmXozF1iWgdFv3Wz1qmWmNRmPbzYOcjZeX171TG+GEegh6DASAJRHoJD/0xStE/AH86+HYOF948d4H5EqpFZcIGgrOp9Ni7GPXLTab/dCJcY56IeKEYCCEHJk3B+wfi2XUUv+4RAx0RzYPRv24dlG3Qg9SGMHb2cQvFeQHA7Clkb2+lRIEtbN1pwME9b4p/kjBLFqUACQewv95i9Thts4Q9L9ICuwqJqN/N+kJcGc2/TkYBSHrgoEQcgpsGng3Gbs+nXZLQkX/ju8rh/s42YuzDVTyYfyHUvLIBNr2EZgr3BMNsjrYNQo5kWA+8ttY7FITtfYasfk2+fFAeOPQlm5JqPU3iHIl+GggCtdLg2wIBkLI6YzwQq5Opx2oJJ+/TPhywPsDsGGeMBxaVZGc2phDXmgk30zE/hFl+3GhkJODJyDkjBAAZgejBbNoC8PRp84TaSfwy82wr9QaiuTU4vPEqAy8vytSNo/+fAyMgpDtwXMQcl40FDwbgRbPoc0MQheeJ8Yew0/Xww1C5VYAACAASURBVHDYW3LaqCfPEaMz8CghUjqXvj4B5cIOKcg+wDMRcnYMFCyLQpdEoPvKyZevEgwMvBKPzg6GLRXLoAA4VUdtvk0UycGaeHRnCp1Ht3WeIOh/wUAIQQAAQEPB4nB0UTiaUUN9nk+8dp18PgZdGom637+uFtRZGhzsLSO/yCcZGFgTh84PhdcWkJ2CgRCC/gsBID0ASQ+g5UqprQVk5H5Tuj+6PBodDkfTdEWBjPq2iPypjBzljX45HEv1hqUH2TUYCCHoIRLEyM4U7NNB2O4SctlFAgHg2Uj0qTDUi23rnNkxpQn8VkHuLiGr1eDZCOTmTJo/XMcH6gtgIISgRxIxwcvx6Mvx6MUmancJGfO7aZgHsiAMnR4IB3r8l4kEp+qpfWXksVpyjA+6PgGd7I9iMAJCfQf8NUPQ46V4ISle2FYcO1RF/lRGPv83MdEPnROMpPmjHGf9DeEk+KuR2l9JHq4iIwTIglD0i2F0uC4M1Bc5648YgrqOSwMLw9CFYajEAA5Vkd8WkUsvEqne6LRAZLI/6ukcvaYKIzhVT/5ZTR2vJcMFyJxgNGsGLYAHG4BQHwYDIQR1mSsTPBeJPheJygzgWC35Zw219poplI9M8kcm+KJDPRCaYw2PJClwS0KdrqdO1pFZbdQIL2RaAPrxIJoPB8Y/yBFYMhAaDAYms6OeEZPJRKfDOUSQ4xAxwVNh6FNhACexv5upU/Xky1eJEgU1zBMZ5Y2O8EQGuiPMvrknMEGBXAl1qZm60EhdaCQ92MgEP+SVfthob8Rpe4MhR2WZM7qurm7+/Pl5eXkMBuPzzz9fuHDhfQneeeed3bt3t7a2crncVatWbdy40SLvC0F2goaCUd7IKG/sgwFAZgCZTeT5RmrNVbJQTiW4IoPdkUHuSJIbEuaC2PMOQzVqKkdCZbVSV1uoG62UHxcZ4YXMCka+HEb35jz+5RDUR1kmEK5evTouLu7ChQvZ2dljxoxJTU319fW9N4G7u/vp06cjIyOLiopGjhwZGxs7Z84ci7w1BNkbERNMD0SnBwIAgAYH5riyv5J6PYuU6ql4MdJPjMSLkSghEi1EbDgfQ2oAxQrqjowqkFF5UipXStEQkOSGDHRH1sZjgz0QMRz5AjkHhKJ6uriiXC53d3cvLS0NCgoCAKSnp48cOfK11157VPrZs2fHxsa+++67D/1XPz+/a9eu3RdHH0Wn0zEYDAzrm31PfZ9KpeLz+bbORV8iM4BcKXVbShXIqEI5VSin9AQIc0GC+UgwHwTwkAAe8OEgPhzgznp8n2pnyh8nQYueataBeg2o1VA1aqpKDSpVVKmCIigQKUBiRUi0EOknRvqJEdjs6xJ4/tuWBcvfAi3C6upqBoNhjoIAgJiYmPLy8kcllkqlly5dWrFixaMSUBSlUCg4nLu/SC6Xy2Awep5JCLIHIiYY7Y2MvmelFbkRlCmpShVVpQJlSupcA2jUkg1a0KKjmBhwYyFiJhAygICBsGmASwNsGmD9J0AajTQGgzA/NpFAbQJ6AuhwIDNSCiOQG0GbnlKZgDsLeLIRPy7w4yJ+XGRqAAjmo2EuCFw9DoLMLBAI5XI5l8tt/5PP51dWVj40pclkWrx48cSJE8eNG/eoo2m12qFDh6Lo3VF327dvnzx58qMSwxahbWk0GgSx41tefQENgCgWiGIB4H7/PylNiMQA5EYgNyJKE9ARQEcgOhwYyLtlTqOMTPCfxxgVxAIMlOLQgIAOBAwgoFNiJhAzHtHlgwO1urc+lJOA579tdbL8ORxOe0B5FAsEQjc3N4VC0f6nTCbz9PR8MBlBEIsWLQIA7Nixo4OjcbnczneNYhgGA6ENURTF4/FsnQuHxQPAp8MEKpWJz4fNOpuB579tWbD8LTDdKSgoiEajFRQUmP+8efNmdHT0fWlIklyyZIlUKv39999hVycEQRBkPywQCLlc7oIFCzZs2NDY2Pjzzz/fvHlz/vz5AIC8vLz09HRzmpUrV549e/aFF164dOnSmTNnCgsLe/6+EARBENRzlpk+sXnz5jVr1gwbNszb2/uPP/4Qi8UAAIqicBw3J1Cr1TExMV988YX5z8mTJz/YaoQgCIIg67PA9AnLgtMn+hA4fNy2YPnbFix/27Jg+TvWkogQBEEQ1EVw0UCHJdPL81sL81uLqpW1jepmiU5mInE9rndh8AUsFw+OW6goOFwUnOzV35UtsnVmIci5qI2anOa8YklpubyqQd2s0CuVRiUNpTMxhoDp4s3z9HfxiXGLjHeP8eJ62Dqzjg8GQkfTppWcrrpwsfZqjbIuzj0q1i1qWvgkb66nK0fMQOksGkthUCoMqiZNc7ms6u+661uzd3pw3Eb6D5sUMsaT+8BcNgiCLEdpUJ2uuvBX9cUKeXU/j5ho18j0sIl+fB8Bky9guphI3EgYZXpFk7q5Sll7qfbatpzdAqbLSP8hYwNHBgr8bZ19hwXvETqO7KbcQyUZuc0FIwOGjg4YnujZj4Y+vmRIiiyUlJyuvHCu+mKka9j86JlJXv06+Y7wHoltwfK3rS6Vf5ms8uc7B682ZA31GTghZHSiRzwde/xWPCRFFUlKLtRcPl113p/vOz0ibXTAcBSBt7QAsOj5DwOhI7hcf/3H/N+0Jt2cqGljg0axad2ZZG0iTGeqM3+5c5BFYy1LWDTAu/9jXwIrYtuC5W9bnSz/UlnFjls/lsur5kRNmxo2kUvvzoquOElcqrt6oPioVCdbGDt7Qkgqhjh7vQcD4V0wEN5pK96W850O1z8d/+QIvyFojxd8ogCVWXNlZ+6PHlz355OWhggDO0gMK2LbguVvW48t/1Zt2/abP+Q05z4d/+Tk0PF01AK3onJbCr6//UurVrI88ekRfoN7fsC+CwbCu5w5EEp1sm05u3Nb8pcmLJwQnNrzEHgvgiKOlp3anbdvQnDqkn4LHtXEhBWxbcHyt60Oyp+giN+LjvxU8PuMiMnzY2Z2r5OmA9cbcr6+uVvIEqwZuCLAxc+yB+8rYCC8yzkDIUlRh0uOfX/7l/TwCYti57JovbVrnNyg+CZnT07z7XVDXkj2SngwAayIbQuWv209qvzL5VUfXt4iYglXD1zuy/fupXcnKfJQybEfbv+aHjb+6fgnGZjTLV0JA+FdThgIa5T1n1z9AkXQVwavCnDpVCn1UHZT7sdXtw7xSV6ZtIT1vxe2sCK2LVj+tvVg+ZMU+VPB7weKj6xIXDIpZIwV8iDVybZm7yyTVbw25MV4d+dargsGwrssFQgJA2lUmEwagjCQhIEAAKB0FGOgTBGdKaIjqF3stEIB6veiI3vz9z/T78np4ZMt2xfaMY1JuzVrR0Fb0VvDX4kQh7Y/Dyti24Llb1v3lX+zpnXT35uZNMbrQ15y47haMycXa69uydo+NjDluYSF9tM0NMhNeomRNFHmShVjoCgDpfNoDB6NxrVA6wUGwru6EQhNalzbbNC1GvVtBl2rUS8xGmQmiqAYAhqdR0MZKI2FAQAII0kaSb3UaFLhHB+WKIIniua7BNtsA+8WbduHV7aYCPyNYat9eF42ycPZ6otbs75dGDt3VlQ6AhAAK2Jbg+VvW/eW/8XaK5uvb5sX/cS86CeseZHaTmlQbb6+rVpZ9+awNWGiEOtnwEzbqJfeUcmK1aoaHY2JstwYGBPFmHcrVcJAmtS4SYWTBMUS0VmuDJYbg+3OZLsz2O5MppAOulJyMBDe1UEgJHHKqDAZZCa9zKiXmPRtRl2bQd9iBBhguzM5Hky2O4PlxmS5MVgiOo3zyFBK4pSqWisvUUvylAiG+KS4egwQIphVT/TzNX9vufHNnKjp82Nm2nYKUaO6+Z1LH3tzPV8b8gKXzoEVsW3B8rctc/njJLH91veZNZffTVkX5Rpu2yydqjy/LWfXgtjZc6KmIV2KKj1EAUm+suGiRN9mdEtwEUbwXEK4GPORlRVhIA1So05i1LcZda1GXatB12LA9STbjcF2Z7Bc79bMTBGdIaRjjIcfBwbCu1QSTdNZOYIgFEkRepIwkriOwLWESYWTJpLhQjf3bd697nBjsN2Y3W+SU0Beqq4/32aQmcLm+lqndajHDVuzd9xqzn9r+Fqb/8bMTIRpa/bOnOa8TSNfd0NFsCK2IRgIbUulUplo+NsXP+YxOG8MW+PCsIvvolHdvOnvzVw6541hq0UsoRXeUddiKNvfQBhJvzFurvEu3b6XRBhIXatB32bUtd3tqzPITUa5CQBA59PoXIzGwTAmRufTQmd6AxgI26nlGmWBDkVRgCA0NorSURoHo3MwOp9GY/fWCBpJnrLicKNrP5fgaV69evuwUl698dInka7hawausPjw6x4yX3guj3s6LXKcrfPivGAgtK3r1Tmf5Hw5LXziori5Vm1+PQ5BEXvyfjlWceaNoasfOt7bghoyJbVnWgMmuHsPd+2lMiAMpFGJ41oC1xGEgaBI4J4oADAQtrPVqFFcTxT/UAsQJGqxfwfN/544UnZy5629K5OWTLTK2LNuKJNVbjj/wZiglGX9F8E1n2wCBkIb+rP0xK7cn94YtnqwT7Kt8/JwN5tvf3D584khqc/2W9Aby9BQJFVxqFFZqY1ZGsgUPX65OIuDgfAuG06foEiq/GCjqlob/3/BHdxi7AaNSbv52lfVyrp3RrxmnQkS3dYgafzs1jYaSnt7+Cs8BtfW2XE6MBDaBE4SX2R9m9dS8MaA1RFeYbbOTkfkBsUHlz/X44a3h69157hZ8MgUQRXuqaEIKurpgF5qDDwW3I/Q9hAUCZvtI4zg3fmuhsQtdjFRIi1fdnwNj8H7euKndh4FAQB8Bu/TMRsDXHxXnHylRllv6+xAUK+TGxQvn31TopNum/ipN9fT1tl5DCFT8EnqO0N8kv9xYu2V+iyLHZcCpb/VAwBingu0VRS0LEf4DDYUnO7FEtOLvq+hSAvEwoPFGa/9tXFZwqKXB/2f/UwG6hiGYKuSn3sqZtaLp1+/3pBj6+xAUC8ql1etOPFKgkfcppFvcOhsW2enUxCAPBU7+72U9Z/f+PrrnN04SfT8mNXHm7XNhshF/nYyx7rnYCDsGQSEzfUl9GTtmdaeHEZlVL+Z+eGJyrPbJn6aGjjCUrmzmrTQce+PfOPjq1/8VnjY1nmBoF5xsfbK2rNv/aP/4qUJT9lkpmBPxLtH70zbUqOsf+H0ukZ1c08O1ZarbL2liH0u8FGzGvoix/kktoLSkKin/ZsuS1XVuu4d4XZr4dJjq725Hl9N+MRWk+V7Ls496utJn52uuvDhlS0mwmTr7ECQxVCA+v72r1uzd36SunFMYIqts9NNLkz+P0dvGBs46v9OvvJX9aXuHcSoxMsPNkQu9KfzHGpTdxgILYDOo4XO8ineW0sYyC69kKTI72//8vbFj9YMXPF88lKL7NJiQx4ct63jPzIQxpfObJDoZLbODgRZgB7Xb7z4yfXG7G8mfnbv+oJ9EQKQ2VFTP0nduCvvp4+vfqHH9V17PQVKf633GS7mB/SNbuHOg4HQMlzjXVyCOVUZXehzaNK0vHj6jdyWgp1pW4b6Dui9vFkTi8Z8Z8SrQ30HrDixtkhS2ntvROKUttkgL1Gb/9M2Gbp6FQL1RSRO6VoMijKN+XvX1OsJYy9+702aludPrePQ2VvGfiBmi3rvjawpQhy6I+1ziqKWHV9TLCnr/Asbr0hxLeE3zr338mYrfbsJYldCnvDO/qjUa6iI6/34ye8nK85ty9k9P3bm3KgZfe5+Q8cQgCyKmxsiDFp//r3/S7TwPEhdi6EtVynJU2ibDUwR4+7ihBQwKk16mYnGxviBbJcgjjCS15lvAeoT9FKjvEitrNSqarQGmYkpojMEd5fCN6lwXauBKWK4xvPd+gl4Fm2p3Gy+/d7fny2MnT0rcqoFD2sP2DTW+qEv/VV9ad35d2dHTlsQO+uxU4FxLVFzoiV+ZbDDDJC5F5xHaEkNFyXSO6q45UEdpJEbFJuvbatVNbw57OUwUbC1stYrOp7HU62o3ZD5z0HeSc8nP9vz+bx6qbH6eIuiRO2WKDC3vx/8QeqlRlWVTlmlkRWqKYJyjXdxSxS4BHLsadEPS3LseYTaJkPrTbkkT2nSEqIoniCUyw/kcDyY93+bFFA36CV5itYcBcuVEZTuyfOzQDjcX/TnvjsH3h7+SqJn/KPSOED5t2rbPrzybz1ueGPYaj++TwcpKw43kjgVNrujNFYGJ9TfZW+BkCKonE/LQqZ7iaIf/vVcqLn876ztE4PHPNtvAR2zwVoMlvXYE1Ft1Lx/ebPWpHs3ZV1Plj2sv9BWe6bVJ8XVd5RbJ+ct6VoMbbcUrbcUpJHyGCj0GCBkufaNGSmd5wAV8YNMarwlW95yQ47rCLf+AvcEAc+f3ZlLGYqgmq/Jak61uCUIgqd5dXtlfD1u+PTalzXKuvdHvu7J9eggpWOUPwWoQ8XH9tz++en4eU9EpD+0g0rXasz7ojxpXbhdjZGBgfAuewuEAABpgaoqoynxlbD72isyvXzLje2Vipp1Q16IdYuyVfYsqzMnIklR39/+JaP89MYRr8W5d/mDE0ay7Nd6ncQY/UwAU9idSwdNvb75hqw1R8H1YXkNFbnGuVh585De4xgV8V0UkJeqm65I5SUacRzfc5BIEMLtRlOe0JMl++pMWiLqaX8Gv8u1dr2q8a3MD8PFIS8PWsl83FxeRyr/OlXDx1e3UhT52pAXH1zKo3B3DT+Q7TfGvu4OwkB4lx0GQgBA3tYKnxRXt/4C858UoI6Vn9lx64cpoeOfjn+yr8yU74zOn4hX6rM+vvrFU7GzZ0dN7fzyxISBzP+mkuPJCp3tg9J6FL1InJLmK5uuSLXNBs9BIq+hYpusjmhZjlER4xqi+Yas8bKUxkS9hordk4QYq2eD+ChQe6a1+ZosflVwly6ezHsKLum3YHp4WmfSO0b5tyMp6o/SY7vzfp4VOXVB7Kz2QeyaBn3BjuoBGyJ6+Bu0OBgI77LPQCi9o6o50dL/5VAAQKWi5l/Xv8ZJ09pBz/f1O4IP6tKJ2KRpefviR15cj3VDXuTSH7+JFYlTd3ZWs8T0sDm+FrzJp2sxNF6WtmTLXYI43sNdRZG8vnsHsa9XxKoaXePfEmm+ShzH9x7myg+05FCX+gttzVdl8atC6J3Yec1E4t/e/D6z9kqX9hTs6+X/UM2a1n9nba9TNb48cEV/z3gAQMlPdRxvlt8YSy5VahEwEN5ln4EQUCDn01KfdPEB3eHTVeeX9FswNWySgw0NNevqiWgiTNtu7r5an/X2iFeiXSM6SkqBoh9qAEAiF/n1xig10ki23lQ0/i3FdYTXMLHnIFFnqkt700crYsJItuYomi5LcT3hPVTsOUjU/V1CO1R9rFlWrI5/PrjjNVAa1c3vXvrUlS1aN/TFLu0p2EfLvzMu1l79MntnrFvkspDF1V9LB26I7GkzvRfAQHiXfQZCkiJPZ/zdlqNonFi9rP8iAdPF1jnqLd07ES/WXvnXjW/mRk2fF/3Eo64PGjIlbbcU8c8H9/b9PFWNrumyVHJbKYrmeQ0Vd+++lK30uYpY06BvuiJtvakQhHC9homt0Bwv2VeH0tGwOY8c63i2+uLWrB0LY+fMikrv6p6Cfa78u8RAGH8q+F110hjqFjhu0VCWnW2JCmAgbGeHgfBS3bWdt350Y7rO+ntu/LJQnq/dnT0W1O0TsUXb9sHf/0JR9I2hqx/cHUbToM//piphdQhLbKX7qbiOaMmSN12RUgTwHCzyGCjsxjgL6+srFTGuJ9puKpqvyYxK3HOwyHOwqHvjnrqBMJA3N5cFT/Vyjb//elRj0m65sb1YWvbW8LXhopBuHLyvlH+34Rri+j+LL43567oy+5n4eWmh43pjX8Nug4HwLrsKhNlNubtyfzIQhmX9Fw3xGVD3V5uuxRA+z963UuqJnpyIJEXtu3Pg96I/ViYtnRA8+r/P49Stf5X7jXXzSO7+dItuU1Vrm67KJHlKl2COx0ChONbF3gYI3MvOK2KKpBSlmuYsueyOShjB9RxsjSbgg1TVusLvqvuvDWO4/Pfi5mbz7Y+u/HuwT/LKpGdZNGY3j2zf5d9z9efbtI2G8Pm+RZLSb2/90KJtW9JvQWrACDu50QMD4V12EgizGm99n/+rwqB8Om5eamCK+SwxqfHsD0sHvh3pGPt1PVTPT8RSWcUHlz/35/u8POj/zBMNq0+06FoMUYv9LZTH7iCMpCRP2ZIlV9fpXONd3BMFgjCuHS6oYacVMQVUNbrWm/K2WwqmkO6eLPRIEvbSXcBOqjnZom3URz0TAADQ44YduT9cqLny2uBVg3ySenJYOy1/y8n5uDRsnq9L0N2hbdlNubvz9qmM6sVx81IDRzx2MZreBgPhXbYNhARFZNZc2XfngInEF8bOHhM48r4LpaI9NcIovtcQB1mi8EEWORFNhGnP7Z8zys+sTHp2lGj4rc1l/deGWa3rrGNGJW6u0PUSk2u8i1uCiyCUaz/TEO2rIqaAqkbblqeU5CpROuLWX+CeJGC7d7OxZVkkTuV8VBo+37eCV/HJ1a1x7tEvDHiuS+NiHsq+yt/SlJXast/qk9bdP4Y2q/HWD/m/tmol86JnTAoZ2+32dM/BQHiXrQKh0qDKKD99qCTDm+s5L+aJob4DHnqbXVaoqjnVkvBS316xvgMWPBFLpOUfX/1ifNHE+LCo2GnduWHTq/RSoyRX2Zan1LUZRFF811i+MJJHY9u4K8IeKmLSSCrKNZJ8lfSOisZG3foJXPu5cH3s7tZ47Y2mghOVu+N3rRm0YoiPZda4t4fy7z2lP9dzfJi+ox4+a6KgrejnO4dut9yZHDpuRkRax0vw9BIYCO+yciCkAJXXcudI6ckrDTdS/IfOikzv+B47RVJZ75fELgvkOOgC0JatCORVmtxdpV/E//uJmMnzomfY5xJ0RiUuvaOSFigV5RquN0sUyRNG8HgBbJt0nNqsIqaApkkvL1HLi9XKSi3Pny2O5bvGurDc7HGxCApQJyrO7bj54/Li5bFjQ/0GW6zKduBASOjJG5uKk19/zJpqjermgyUZJyrOxrpFTg2bNMQ32ZqjaWAgvMtqgbBO1XCm6sLJir+YGGNK2ISJIamd7FepPtFC6ImQGd69nUObsGxFcHtbpcdAIRVj+jJ7Z6W85oUBz1nqyr03kDilLNfIitXyErVBZnIJ5ghCuS7BHJ4/22p9p9asiCmS0jYZlBVaRYVGUa7BmKgwgieK4AkjeHY4w6xdsaRsS9Z2AMDqgct9VL5F39cM2BBhqS/IgQNh0xWpvEQd9XRAZxIbCOO56osZZafqVI3jg0aNCx4VKQ7r7RwCGAjb9XYgrFHWX6q7+lf1JYlOOjpgxMSQ1K5+wXqpMfffFYPeibTDoRY9Z8ETUVWtK/6xNvmNcHNBXW/M2Zq1w4vruTJpSbAw0CJv0XtMalxRoVWWaxQVGn2bkevD4gdy+AFsnj+bJWb03jjJ3q6IjQqTuk6vqtGqanTqah3dhWaO94Iwrp3cxO1Aq7ZtR+7erMaby/ovnhQyxnzzouDbKrf+As9Blrlt78CBMG9rhf8490dtHvAodaqGU5V/nanKRAAyOnD4SL+hEa6hXZ2d2XkwEN7VG4HQSBhvtxZea8i+Up+lw/Uj/AaPChiW4BHb7SFSuVvKA6d4CsN5FsyknbDgiVi4p0YYxvUe4dr+DE4Sf5Qe+zH/t+F+g5fEz3fjuHbwcvtBGEhVjU5tDh61OsJAcn1YXB8Wx5vF8WRyvJgWvLNo2YqYMJC6FoO22aBt1Ksb9JoGPaAAz4/FC+DwA9j8QE5fWXxHY9L+fOfgH6XHp4VPeipmNof+35Xb5KWaioMNSa+FW6RydtRAaFSYbn5WNmhjVLebzsXSsgs1ly/WXtXj+sE+yYN9kpO9Eu79IiwCBsK7LBUI9bihSFKS23IntyW/UFISIgw0f3kRYgtcztT/1aaTGO1qHy9LsdSJqGs15G2tHPhmBPrAUlhqo2bfnQNHyk5OCR3/ZMwTQqag529nTSYNoWnQaRr02iaDtkmvazYiGGC7M1luDJYrg+XKYIroTBGd4ULvxoTF7pU/RVImFa6Xmgwyk0Fm1EuMujajrtWIawmOJ5PtyeR6Mc2R2/6bfffR4/pDJcd+LTw01Hfgkn4LPB5YqwEAcOvz8oAJHuJYC5y3jhoIGzIlmka9ReZA1yrrrzRkXWvIvtNWHCQISPKMj/eIiXOL5jG4PT84DIR3dTsQ6nFDpaK6TFZZIi0vkpTWKOtChcHx7tEJnnH9PeIse+XiwL2jljoRy/Y3MFxoARMfOYqhTSfdm7//bHXm1LCJc6Kmi1h9LBzey6jCdS0GvcSolxgNMpNeajLIjEYlTmNjDD6Nbv6Pg9G4GI2N0dgYxkQxJoqxMAQF5tYkgiHmyalqtZrH4wEASBNFmkgAAKEnKYrCtQRpogg9getJXEeY1DiuJUwawqQ0GVW4SUPQuRhTzGAK6Swxg+VKZ7ky2O5MppDeh5aXu48e1/9ReuKXwkMJHrFL4ucHCh45D7XtlqLhoqTfCxYYmeyogTBva4X/eA9RlCU7sYyEsaCt+Gbz7dutd4okpe4ct2jX8AhxaKgoOEwY3L24aMHy7wPrSPWQwqBs1UpatK2N6uY6VWO9qrFGWSfVyQIF/mGi4HBRaFrIuHBRcO+NUWSJGSwRXVGuFYZb4CLI8ZjUeFuuIvn1jtbgdmOLVw9cviB21k8Fvy8+snJ88Oh50TM8ufa1O1onMfg0Bp8mCL3/ZDCpcKMaN6lwowrHNYRJS+hajYSeIAwkYSBxPQFIgOsIAABFUISBBABQFIUgCAAApSMoHQUAYCwUQRAaB0PpCMbEaCyUxsGYIjrPl03nYXQXOoNPo/MwR7omUxpVh4qPHSo52t8z/l9j3nvsHWXXfi5VGc2qnEYU1gAAIABJREFUap1lN7twGEaFSddisHhlxcAYiZ7xiZ7xAACSIisVNUWS0mJJ2bnqixXyag6NHSjw9+V7+/K9fXheHhw3D46bkCW02hI2fbtFKFFJj1efRVFUa9IRJKEn9HrcoDXpVEa1wqBSGJRyvZxNZ7uzXT257p5cD1++tx/fJ1Dg58X1sOayCHV/tekdsXfUIldk9RfatI2G8Cc72w8j1cl+K/ojo/z0QO/EedEzrDM+zT45aoukkxrUTfuL/jhTmTnCf/D8mFkPbif7KHV/telbDWFze9r155Dlb8F+0c5r0bbVKOrqVA31qsZGTXOLpq1F26YyqgRMFwHTxYXp4sLgcehsJsbk0jksGvPp+CcBbBG2IylSbdQgCMKhszE65oa5smhMLp3Do3MFLBch00XIErZvL2lDbgkuuVsqQmd6O9KVuKU0X5OFzenCr07MFq1IfGZR3NyjZafeyvzIjS2eFZk+MmCYPXzRkBWQFJXVePNgydHCttL0sPF70r90ZXdtFKjnAGH2R6XB070dePnDbmvLVfiPt/bseHMTcIB3/3ufJChCrlfKDQqlQaUwKHW43oAbtLiOhVl+Wnbfrjt4dO6yhEU2X2v0sVhiBlNEV1ZqH+wQc3Kqai1FgvbFDDuPS+fMi54xJ2ra33XXD5VkbM3eOSlkTHrYBD++ozW7oXYSnex4xdmMslM8BveJiCnvpqxnYt2Zwk/n0wRh3LZbCs/BDrv8YfcYVbi22fL9ot2DIZgrW9TVq5zu6duBsA8Rx/BlhSoYCO/TfE3mOVjU7TEaKIKm+A9J8R9Sp2o4WnbqhdOv+/K8JoWMHR0w3CLD0iB7YCSMl+tvnKg4m99aNCpg2MYRr0W69rQ/3HOwqO5sKwyE95HdUQkjePazmq7VwEBoJaIofumv9UHpts6HPSGMZFueMum1+1f17QY/vs+KxGeW9V90rSH7RMW5bTnfJXkljAsaOcRngA0XBYZ6AieJnObcc1UXL9VdixCHTgoZs3HEOkt9m6IoXtlv9dpmA8cTnh7/JStSiWMcdiPxDsBAaCX8ALZJjRtkJqaoj83N6j1tuQpBCPfeXeJ6CEOwYb6DhvkO0pi0mTWXM8pPf3J160DvxJH+Q4f4DuDSu9wBC1mfkTBmNd26WHv177rr/i4+owNGLOu/2OL9YwiKeAwUtVyXBU31suyR+y6KoOQlmtCZznhzAQZCa0GAKIonK1R5DRPbOiv2ojVH0Ut7VHHpnLTQcWmh45QG1d91185UZW6+vi3KNXyI74AhPgM6P7wQsppWbZt5OaebzbcjxKEp/kMeNSPeUjySBAU7qoPSvfru7EnLUlZqWe4MOt8Zg4IzfmZbEUXxW2/KYSA0w3WEulonWtKpVX27zYXJN0dEPa7Pbsq7Un/j96I/EYAM9E5M9kpI8uonYDpjR5Cd0OH6W835OU25N5puSXWygd6JqYEj1g19sec7BXYGx5uFMVFVDZxQeJesSCXu4uKiDgMGQusRRfPKf28gcaobi2k5HkmeUhjJxR5YU62XsGis4X6DhvsNAgBUK2pvNN46VXn+02tfenLd+3vG9XOPjXeP7ivLmfZpSoMqv60wr+VObktBpbw62i0i2TNh/ZAXI8RhVps93c41QdCWp4CB0Ex6R23l6YP2AwZC66GxMY43U1muEUY64ALcXdWWp/QcKLTJWwcK/AMF/rOjppIUWSItz20pOF11fsuN7UwaM84tMtotMto1IkwUDEfZWAROEhXyqiJJ6Z224juSkjatJNotop97zIrEZ6Jdwxndmv9gKW79XAq/qwmGvaMAGOQmkxrnBzjpNQEMhFYliuZLC1UwEOI6QlWpjVr8yAUhrQNF0CjX8CjX8HnRMwAAtcr6QknJnbaSs1WZlYoab55nuCg4XBQSKgoOEQb16QVOrUlj0lbIq8plVaWyilJpRbWyzpfnFekaFuseNTtqWogw0JqLOnWM68NCMERdp+P5O2kAaCcrVIkieU57QQADoVUJw3ll++ttnQvbk+QrBeFce1vXw9/F19/Fd0JwKgAAJ4lKRXWptKJUVvF3/Y0KWRWKoCGiwAAXv0AX/wAXXz8XHw+Ou/V78+xNm05ap2qoVdbXKOur5DVVihq1SRMkCAgVBkWIQyeHjgsV2nXb2i3BpS1XCQOhvFTjtDcIAQyEVsbzZxlkJpMap/OcuuQleUr3/nbdwKKhWLgoJFz03z0KJDpZlaKmSlFbpai5VHe1VlkvNyi9eZ6+PC8fvpcX19OL6+HJdffguvW5vaI6Q2PStmhaGzUtzZqWJnVLg7qpQd1Ur2pk0Zh+fN9AgZ+/i2+SZ79gYYAn17339mK1ONcEQdGemqB0T1tnxKYooCjTBDvxTBKnro6tD0ERl2COolzjluCAdWUnkSZSUa6JeMrP1hnpGvNqT8leCe3PGAhjg6qxXt3UoG5qVDffbM5r1rS1aFv1uMGD4+bKFrlz3MQsoZgtErOEQpZAyBQIWQIBk8+iWX6xxB4yEkalQSU3KGV6udygkOkVUp1MopNJdNI2nbRF04oiqAfHzZPr4cl19+Z5RrtF+PC8fPnefX12Js+XBSjg5DPrtU16jIk68xRnGAitTRDOVZQ5dSBUlGl4fmway95XiH0sJsYIFgY+uO+PgTC2aFolelmrtk2mV0i00kpFjUwvl+nkCoNSaVSRFOXC4PEZPB6Dx6VzuHQ2h87hMbhsGouBMXgMLgOlM2lMBsZgYgwMQdn/2SCTz/ifu8sarUaFaO59RmvSEhQJANDjBpzEcRLX4XoTadLjBp1JZyCMWpNOY9L+5/9alVGtNKpVRhVBkgImX8B0EbIEYpZIxBKI2aIQYaCYJXLjuHpw3Cy+vbj9EEbxZIUqZw6E8jKNnawvaiswEFqbMIxXdKXW1rmwJWmhWhTlyHcjmBjDfLvxUQkMhFFlUKmMapVRrTHptCatFtepjGo9rpfrFfWqRgNhNBJGA2EwEiaCInUmHQCAApTa+D9hjyRJFP2f+6wcOgdDUAAAk8akozQMxTg0Nh2ls2hMNo3FoDH4TJ4Xz4NNY3PpHC6Dw6Nz+QyeC5PPtr9GqtWIo/kNmRLf0b04ed/OOfmlOYCB0Pq4PixcgxsVJobASTsiZEWq6F6eR2/nmBiDyXHt+bRFh9wPz/oE4dzin2oJA2lvo7eshALKCk2ow+2W2iVO+cXbFgJcQrmKcs3jUzoiXYuBwimul/O2PyB7gzFQfgBHXqq2dUZsQ/3/7d15lBTl3S/w6r2W7p6FmZ7u6dlYBhBUBkEU8YqCKB4lYtyJS8REFMlxzTkafdHXJEcNvjlqctVIQI/JTVSMnAjeG41HURJABJUlAsMyzNo9+3RXV1Vv1XX/aDMZZ+2Zru7avp+/enpquh6arv5W1fN7nqdVsDmtdkPOrNYPQaiAwmlM6IRBg7D3aKToDJd2igrBEIrOcPYeMWgQho5zBbVGH9mMIFRAQa2z77hBg7DnCFs00+hHHahN8Rmu3qOs0q1QRugkVzDN0JUyBIJQEbTHIcZSsb6E0g3JNzGeYhv5wukIQlAXyuMwmU18MKZ0Q/JOIsKn+YIp2h4Dkz0EoRJMhKuaYht5pduRb+GTnLOCMmhJAqhb0Qxn7zHDXRTywaiNsRh8fg8CQagU92Qm3GC4IAyd4AoNfxMG1Ck9wFfpVuRbuIF3T8YhiSBUiHsybcAg7DuBbnlQqYKpTPgUL6UkpRuSV+EG3lVj9PuiBIJQKc5KSuiIibGU0g3Jn2RUFDpimN0Y1MnmtNoLbFxLVOmG5FX4NO+ejCBEECrEbDUxPjLSLCjdkPwJn+Rd1TQWJQbVKqxl+ox0dzTBJpO8SHuMO7dcP9mCkOf5pqamVGrESxxRFBsbG6NRY51wjcJlsLujfccjBp/PEFSuYBoTOmGg0YShBt49mcagXkKuIHzxxRf9fv+SJUtmzpx55MiRoRvs27dvypQpy5YtKy8vf+ONN2TZqda5a+hwg4FOP0MnMFwJVK1gKhM+zUuiUboJ2QYO90XTZAjCxsbGn/3sZ7t27Tpx4sSqVavuu+++odusWbPmoYceqq+v//DDD9euXdvd3Z39frXOPYVhTwsG6ZxP8mKsJ+GsQAchqJeVtlCT7MbpsEClTD8ZgvDNN9+85JJLzjjjDIIg1q5d+/HHH7e3tw/c4MiRI998882PfvQjgiDmz58/Z86crVu3Zr9frbMxFpvLwrcbYgxv6ATnmkybLLgLA6pWMM1pkG7CVCLFB2MuFK8RBCHL6hMNDQ3Tp09PP/Z4PC6Xq7GxsaysbOAGfr+fpr899aitrT19+vRIr5ZKpQ4ePBgMBtM/TpkypaioaKSNJTGZaGkSzVot+WFKk71fNdqSWm2/yPNxOqMzyp6vk84SU7z5eK6bZCiZv/+QIaYwFfxKLJvRl8nGmn7/2TaJLJaS7SeVbsi4mcwWm3+KvK8pQxCyLFtS8p+lvBiGCYVCgzagKGqUDQaKxWI//elPbbZvlyh66qmnLr744pE25jvaUn95YcItVx5/Zk9DieXwDqXbMUGpVErI7Cykr/uGUtdn3ceDuW6SoWT+/kOGUpIj0nFr15//t4kYu89C0+9/H3e2RSzo/vNOpRsybibKyfzwvwiCiEQyqmyiadpiGWMZcBmC0OPx9Pb29v/Y29s78HIwvUFfX9/ADWbNmjXSq1EU9cEHH/j9Iy5qOpDVWmV/6Ddj/iNVy9UkHN/S6nvoBqUbMkEZrocnxlOn1x+teuS/MHZCXliPMBfanz1e8IP/YfxjrxSm6fc//Mdm70yXZ/7NSjckK3K9/zKczsyZM2fv3r3px4cPHzabzVOnTh24waxZszo6OlpbW9M/7t27d86cOdnvVwcYPxntjItxnQ+rZxt5xk8iBUETDDLrE9souKrQQfgtGYLw+uuvb2pqeu655w4cOHD//fffcccdDMMQBPHII48888wzBEGUlZVdd91169atO3jw4BNPPCFJ0hVXXJH9fnXAZDHRXgfXqvOxlWwDZq8AzXDV0OHTOg/CJC8meZEqxVD6b8kQhAzDfPTRR7t27VqzZs38+fN/9atfpZ+vqKjw+Xzpxy+//HJVVdWPf/zj+vr6Dz74wGo1+mTn/ZxVNNuk86MufBpV2qAZ7hqa1XsQsk2Cs4LCUPp+8gTSWWed9e677w56ct26df2P3W73Cy9ouaolZ1xVlM5XBJUItkmYvgpBCNpAlTrEeCrWl3AU2pRuS65Emngn7osOoNWSJ91wVlJsk54H8HLBqI2xYsEz0AwDLBfKNgsYQTgQglBhtMeR5MQkJyrdkFxhT/Nu3BcFTdF9vUykWcAV4UAIQqWZCGcFyep3VqdwA+9CpQxoilvX9TKx3gQhETq+8TsBCELlOStpHU9vGMYVIWiNs5ISgrGUTsc1sbgcHAJBqDxnFaXXwtEkJyYjIl2GKm3QErPNTHsdEZ2Oa4qgg3AIBKHynBWkXocSss2CsxJV2qA9zipKr/dpIumjEgZAECqPLLKL8VSCTSrdEPmxqNIGbXJV0not5+ZaBcaPo/I7EIQqYCKcflKX92EiTZjGCTRJrx0Wsd4EYTbZ3RjO9B0IQlVg/BTXqsPTz0gLbsKAJul1XFOkRcD62EMhCFXB6ScjLXq7IkSVNmiYiWD8ZKRFb6enXGvUWTH2whpGgyBUBaaCiujuipBtElzVGDgBWuWq0mE3YaRVcKKDcAgEoSrQHkeCTSajuroPE2lGpQxomC67CSMtUQZXhEMgCNXBRDA+vQ2iYFEpA1rmqqIi+roiTESSqXiKLLIr3RDVQRCqBaOz0YQSwbVEcRMGtMtRaCNMRKwvoXRDZBNpiTorSIzrHQpBqBZOP6Wnnnm+I2ZzWayMRemGAEycs1JXF4Vcq8CgZHQ4CEK1YCooPRWORppRpQ2a56yk9DQhfqQ16ixHB+EwEIRqwXgd0Z54KqGTeX4jLTj3BM1zVuhqgC+uCEeCIFQLk8VEldj5YEzphsgDw5VAB5x+UjczjorxVDyUpEpRKTMMBKGKMOUk16aLu6MSwbWiUgY0z15gM5lN+qiX4duilNdhMqNUZhgIQhXRTRAKXTErbbHSqJQBzWMqKH2Uc3Nt6CAcEYJQRZhykgvo4pBrxaBd0AmnXiZa49qitA9H5fAQhCqimyvCCEYQgl4wehnXxLVFGVwRjgBBqCI2p9Vs0UOHRKRFQKUM6INO1s2WCC6IIBwRglBdGL8eLgqx8ifoBlmsh3Wzoz1xK2WxUui2Hx6CUF0Yn+aDECt/gq7oYt1sri3KoINwZAhCddFBNyFW/gSd0cG62eggHB2CUF1o7QchhtKDzuhg3WwE4egQhOpCexyxvoSmJ1qLoIMQ9EUH62ZzAQThaBCE6mKymKhSbU+0xrVFnX4ccqAfVKk9EU6KUa2enoqxVCKcJEswudqIEISqo+l6maQgJnmRLMYhB/phMpsor4MPavWo5ANRqgyTq40GQag6jI/U7iH3bXEajjjQF02fnuK+6JgQhKpD+0guoNVbo+iTB11iysmIZoOQD8Zor0PpVqgaglB1aK+D1+yMo1xblEEHIeiOpsc1cQEMIhwDglB1HIW2lCglIpqcyYJrxRUh6BBTTvKBKCEp3Y4J4YMxTLc9OgShGtFeUouFo1JKEjpitBeHHOiNlbJYaUu0O650Q8YtziYJSbK7MNPTaBCEasR4HVpcj0nojNsLrBYHPlSgQxrtJuQDWH1pbPjOUiPap8krQlTKgI4xforXYhAGY+ggHBOCUI0YnybrZbi2KFOOOWVAnzRaL8MFoigZHROCUI1oH8lpsGceV4SgY05tBiFujWYCQahGVspiIS2aW6GXaxUQhKBX5CR7gksmBVHphoyHRPDtMboMV4RjQBCqFOPTWL1MkhNTCclRaFO6IQC5YSJoH6mtPotoT9xKYz3esSEIVUpzhxyXvgODydVAvxgfyWmqig1D6TOEIFQpxqvFQw53YEDPGJ9DW4WjGEqfIQShStFaKxzlgxhKDzpHe0ltdVjwQZSMZgRBqFK0xyF0xqWUZipHuTZUyoDOMeUkH4hpqJybD8QYnJ5mAEGoUma72e62amZKp3RxGs49QdestMVsN2mlnFtKSUJ3nPJgcdCxIQjVi/Y6tDK/TLQnbqVQnAb6x5Rr5u5otCvuKLCabfiSHxveI/XS0NTbKE4Dg9BQOTcXQLd9phCE6sV4HVpZqp4PoDgNDIHRzrrZqJTJHIJQvbR2RYhDDvRPQ1eEKOTOHIJQvagyh9CtjcJRzGcIBkGnj0pRC0clrggzhiBUL7PV5CiwCZ1qLxxNJaVob4L24JAD/TNbTY5Cm9Ch9ls1kihFexJUKY7KjCAIVU0ThaNCe4wstpssmF0NDEET3YRCZ9xRZDNbcVRmBEGoarQW6mW4IDoIwUBonwaOSj4YZXBfNGMIQlXTRL0MSkbBUBgtTLTGYZbR8UAQqpomrgjRJw+GQvs00GGBo3JcEISqRnsc0Z6EykvUMJoeDIWcZE+wSTGWUroho8HYiXFBEKqayWIii1VdoibGUkleJIsxnyEYhclsojwOvl29R2UqKcV6E1QJjspMIQjVji5T9SHHB6JUmQPr8YKh0F5Vr5ImdMTISTYUcmcOQah2lJdUcxByQazzAoaj8io2vj1Gl+GoHAcEodrRZarumUefPBgQ43Pw7eq9IsSaaOOFIFQ7lReOok8eDIj2qnpMPU5PxwtBqHYqLxzlMd02GI+j0JaKp0RBpYWjfBC3RscHQah2JovJUWQTutQ442iST6WSkr3ApnRDAPLLRNBljmiHGo9KSZRivQmqFCWj44Ag1ADV3h2NtscxewUYE+0jhY6E0q0YBt8RI4tRMjo+CEINUG29TKwzgfkMwZhoryPWqcYgFFApM34IQg2gy1Q6gkLoSKBSBoyJ9pLRdjUGIerXJgBBqAGqvjWKc08wJMbnEFTZR8i3R+kyHJXjI1sQCoKwb9++pqamkTbo6urav39/S0uLXHs0DsrjiHarsXA02plAEIIx2ZxWk8kUZ5NKN2QwPhhDEI6XPEG4f//+qVOn3n///eeee+6jjz46dIObbrpp+vTpa9eunTt37jXXXBOPq/FMSrW+XRRbZYWj8VDCbDHZnFalGwKgDNJjU9tEa+mF6UksTD9O8gThww8/fN999/3jH//48ssvX3755W+++WbQBqtXrw4Gg59//vmpU6cOHz782muvybJf41DhUvVcIEaWoUQbjMtRalPbUSl0xMhiLEw/bjIEYXt7+6effnrnnXcSBOH3+5cvX/72228P2uayyy6z2+0EQbhcrtmzZ7e3t2e/X0NRYTch3x51lOByEIyLKrOpbYVevh33RSdChi+y5uZml8tVUlKS/nHy5Mmj9BSeOHHi008/ffLJJ0faQBTFnTt39r/a2Wef7fF4sm+k1tFlZPfhsNKt+A4+ECO9uCIE4yI99tAhlR2VKBmdkEyD8Omnnz506NCgJ+fPn//ggw/yPJ++2ksjSZLjuGFfpKen59prr3344Yfr6upG2lE8Hv/Nb37jcHx7UvPII49ccMEFI20sCILdbrdYLBn+KzTMLUbahEgkonQ7/oNt5YprKVU1yWg4jjOZcBNMMSlngmsTImxEPcuQhVsjhbMZgxyVGX7+aZo2m8e495lpEC5atGjq1KmDnvT7/QRBlJWV9fX1pVKp9M66u7u9Xu/QVwiFQsuXL7/00ksfe+yxUXZEUdTbb7+dfuUxWSwWgwQhXSMd7w0yFKOWCSMkItaVLKx2O51OpZtiXJIk4f1XkCRJVipsSzocRWqZZTDeFSiuLqCdhrgolPHzn2kQXnTRRSP9qqamprCwcM+ePelLt3/+858PPPDAoG04jvve975XV1f33HPPTbitRtZfOKqSDoBoT9xKWywOjEMFQ2N8Dj4YVUkQomR0wmT4InM4HGvXrl23bt3f/va3Rx99tLOz87rrriMIYteuXVVVVeltrrvuuqampnnz5m3cuPHVV1/97LPPst+v0aiqcBRdEQAEQdBlKlqPCSWjEyZP1d/69esnTZr00ksv+f3+Tz/9lCRJgiDKyspuvPHG9AaLFy+eM2dOQ0ND+sf+WhjIHF2WXgvUrXRDCIIg+CBWXwIgaK8jdHL4koj8Q8nohJkkSV3zlVRUVHz++ecZ9hEaqFiGIDq/DHUfDs+8rVLphhAEQRz7Py1FM5zUDIvL5VK6LcbFsizefwWxLEv0Wk++01b34OD6CUU0/a2DIIiq5UYps5fx848+Hs1Q163RAJbABiBor4PviEkpVVxO8O04KicIQagZlMcR7Y6rYcZRKSUJnXHKg0MOjM5iN9ud1liPKpahwCyjE4Yg1AyzVS1L1Ue74vYCq8WODw8AQfscaphfJl0yitPTicF3mZaoZIVeLhBjsDA9AEEQBEF7STUclQIWps8CglBL/l04qjA+GKURhAAEQRAE41XFFSFKRrOBINQSlZx78oEogz55AIIgCIL2kWqYEB9De7OBINQSlaxBweGQA/g3laybjZLRbCAItUQNh1wqKcV6E1Qp1p0AIAiCMFtNZLGN71D4Vg0XwK3RiUMQasm3haOdShaO8sEoVWpHnzxAP9pHKrtU/benpygZnSgEocakJ/lVsAHoigAYRPHJLoSOGDkJp6cThyDUGNpLcooecphlFGAQxkcqWziKmZ6yhCDUGNrrUPYmDBfAFSHAdyhezs0FYyjkzgaCUGNor8K12nwgSuOKEGAAqsSeiCTFWEqpBmBob5YQhBpDldpjoWQqocwhl4yKSUEki1AyCjCAiaBKlewmRM99lhCEGmMym6gSO9+uzCHHB2K010GgSx7gu2jlqtjEeCrOJslJNkX2rg8IQu1RsEOCD0YxyyjAUIxyfRZ8MEZ7HCYzzk8nDkGoPYxPsXoZVMoADIv2kVybYqenKBnNEoJQe2gvqdStUa4typQjCAEGY8pJrk1QZNfoIMweglB7aOVmu8dwJYBh2d1WwmSKs8n87xpDe7OHINQestie5MVkVMzzfmO9CbPNZHNa87xfAE1gfA6+TYEzVHRYZA9BqEEmgipzCHmvl+ECuC8KMCLaR3J5D8KkIIox0VGIktGsIAg1SZEpnbg2lIwCjEjJoxIVo9lBEGqSIoccH8DsFQAjwumpdiEINYkpV6BWG7dGAUZB+xxCZzzPy4XygSiNozJrCEJNYspJvi1K5PGISyWlaDcWPAMYkdlmdhTme7lQLoArQhkgCDXJSlvMdlOsL5G3PfLtMarEbraiLwJgREye62Wk9CBCnJ5mC0GoVUx5Xg85vg0dhABjYHx5HeMb7Y7bGKuVsuRtj3qFINQqpjyvPfNcAIN2AcZA+8h8Tn/ItWFNNHkgCLUqz4OWUCkDMCYFTk9xVMoBQahVTH7PPXFrFGBMZLE9KYhJIU+zPvGolJEJglCr6DJHtCeRnxV64+FkKiVh9gqAMZgIxkdyrXk6Q43g9FQmCEKtMllMVGmeVujlWgWnn8rDjgC0zllBRVrzsQyFGE8lwkmq1J6HfekeglDD8jaTRaQ16vTjxBNgbEx5nq4I+UCMwnq8MkEQahjtI/Mz2z3Xij55gIww/jwFIYbSywhBqGFOPxnJyyEXaRWYCtwaBRgb4yOF7ngqmfNpn7hWgcF9GpkgCDWMqaC41pxPtCZG0RUBkCmTxUSV2PNQ0R1pjSII5YIg1DAbY7E4zNHe3M5tyLUJtI9EVwRAhpx+Kue3aiSCDyAIZYMg1DbGT3ItuT3kcOIJMC6Mn+RyXDjKd8RsLquVxORq8kAQapvTn/Naba4NJaMA45CHehmMaJIXglDb8nHItQgMDjmAjDn9FBfIbed9pDXqrMDpqWwQhNrmrKAiLTm8IpREie+MM1jnBSBjFtJsc1qFzhxOdoHTU3khCLXNUWRLiVKcTebo9fn2GFlkM9vxOQEYh1yXjvsMAAASlklEQVQPbYq0oedeTviC0zxnOcXl7KIw0iw4K3HiCTA+zsoc3qqJ9sTNFpPdZc3R6xsQglDznBU5PPdkmwRXFYIQYHyclVSkKVdByLVGnZjgQlYIQs1j/FTuarUjzYKzks7RiwPolbOKirQKUionBTOR1ihmepIXglDznBVkJDdDCVNJSeiIMeWolAEYHytpsbusQmdOJrvgWgSMaJIXglDzqFJHgkvmYi1Qri1KldrNNnxIAMbNWUnn6O5opEXArVF54TtO+0yEsyInHRKRJt5ZhfuiABPhqqLYZl72l431JQiJcBRhlWw5IQj1wFVNs43yH3KolAGYsBzVy7CNgqsap6cyQxDqgauKYnNxRYixEwAT5fSTXDAm+3pMkSbeidNTuSEI9cBVTcsehGIsFetL0JhTBmBCzHYzVWLngzIXsrFNgqsaQSgzBKEe2N1Ws9UU7ZGzRC3SLDDlWH0JYOJkvzsqpSRUyuQCglAnnFUU2yjnIcfivihAdlxVFNss51HJB2OOQpuVwupLMkMQ6oSrio40yVkvE2niUSkDkA1nlcxVbGyTgA7CXEAQ6oRL7ivCcAPvnszI+IIARsP4HLG+RJKTbYwv28ijZDQXEIQ64ayiuEBUEuUpURM644TJhLFKANkwmU2uKjp8WraLwghGNOUGglAnLHYzWWznAvKUqIUbuIKpOPEEyJZ7Ch1ukCcIxVgq2h1nfJhcTX4IQv1wVVOsTOee7GneXYP7ogDZck9mwqc4WV6KbRIYP2myoJBbfghC/XBPZUIn5QnC0CnePQVXhADZcldTXCAqy7D68EmuYCpOT3MCQagfBVOZ0EmOyPqIS3BiIpzEUHqA7JntZrrMIUtFdwhBmDMIQv1wFNosdjPfEcvydcKnONdkGkPpAWThnsxk302YSkqRFsFVg/s0OYEg1JWCqUz4ZLYdEuEG3j0ZxxuAPNyT6dCpbIOQbeRpr8PiwDd2TuBt1RX3NDokQxByCEIAubin0GwDn+Vq9bgvmlMIQl35tpswC2IsxQdjGKsEIBeb02pzW/lAVn0W4ZM8gjB3EIS6QhbbTRaT0Dnx2bdDJzhXNY1V6QFkVDTD2XssMuE/l0SJbeJduE+TM/i+05ssLwp7j7FFM50ytgcAimY6+46xE/5ztkmgPA4ribm2cwVBqDfZBuHRSOEMBCGAnNxTGbZZEGOpif156CRXMAX3RXMIQag3hTOcvUfZifXMR3viqXiK8WIOJwA5WexmV+XEC9l6v2FxeppT8gRhKpXavXv3+++/39fXN8pmsVhs//793d3dsuwUhuUotDkKbBNbiaL3SKRopovAAEIAuRXOdPYenUg3YZIT+WCsYBquCHNIhiAURfGqq6666667Xn311ZkzZx46dGikLdevX79gwYLt27dnv1MYRfEsV883E+mQ6D2KDkKAnCiaMcFuwp4jbMF0xmzF+WkOyRCE27dvP3HixN69e//617/eddddTzzxxLCb7d27d+fOnXPnzs1+jzC64tnunn+Fx/tXkiiFT/E48QTIBcZHppePGO8f9vwrXDzblYsmQT8ZgnDr1q3XXHMNRVEEQaxatWrbtm3JZHLQNvF4/O67737ppZcsFhQ+5ZyrikpwYrRnfIdcuIGnSu02pzVHrQIwNBNRNMM13rujkij11XPFMxGEuSXDt15LS8v8+fPTj6uqqpLJZDAYrKioGLjNL37xixUrVtTV1Y35aslk8r333isuLk7/eP755w96qYFEURRF2VZ/1pOimc6uQyHfhcWZ/0nngb7CWc5xvZ94/5WF919Z433/C2cxgX/0eM4vyPxPQsc5stRupk34jx4qw/ffbDabTGPcWM4oCA8cOPDggw8OfX7z5s3V1dXxeNxm+3Yp8/SDWOw7cygcPHhw69atX3zxRSb7EkVx27Zt6etLgiDKyspKS0tH2jgWi0mShKvMoZy1js7Pw8XnZnqfU0pJXQfCM9aUD/q/G108Hh/X9iAvvP/KGu/7T0+x8VuibAdnL8j0CqTrUMg9g8L/8rAyfP9JkpQnCGtqatavXz/0+XREeb3erq6u9DOdnZ0EQfh8voGb/fKXv/T7/U8++SRBEM3NzVu2bHG5XN///veH3ZfD4di4caPf78+kYSaTyW63IwiHcpxFNr3bZTc5rFRGb07vEZYqcRT5x3GuShCEKIo0jdkuFIP3X1kTeP8nnVXAHY0XXuLOaGuJCB9rOePOKprGiKZhyPj5zygICwoKFi9ePNJvFy1a9N577z322GMEQXzyySdz584d1LjVq1c3NjamHzscjtLS0pKSkizaDGOzOMyFM52dX4Z8izK6O9r5ZchzzvhSEADGq/ScwtPbgv5LMvoCDJ3kLA4z40MK5pwMfYS33Xbbs88++9BDD82aNevxxx9//vnn088vWrTommuuefjhhy+//PL+jTdt2nTxxRdfdNFF2e8XRuc9v6jhvWAmQZiKp3q+YSd/z5uHVgEYWeE0Js4m+fYYXTb2wtfBPb1l5xfloVUgQ9VoUVHRnj177Hb7vn37Nm/efOONN6afv+eee4YG3r333jtv3rzsdwpjKqx1ivFUpHnskfXdh1lXNWVzoV4UIMdMRGmdu+ur0JgbJjmx9wjrmV+Yh0aBPN99VVVVTz/99KAnb7nllqFb3nbbbbLsEcZmIrznFQX39E6rHGNNpY59vaXn4HgDyIfSeYVHX2+uvKzUZB6tgqNjf1/xbFeGffyQJcw1qmeeBUVdX4dGn+o30iLw7bGSOnQQAuSDs4JyFNs6vxzjojC4p8eL+6L5giDUM7vLWjCN6dg32gSwTX/rqFhaigmcAPKm6nJP8987RpkZP3SCk1KEezKmecoTBKHOVS33NH3YkeCGH3YaaRG4QLRsAU48AfKnYCpjLxjxolBKSaf+Gqhe7sH093mDINQ5xkeWzClo/H/tw/4Wl4MAiqi+oqzpw+EvCgP/6LExVvRW5BOCUP+qr/D0HA5HmgaXj7bv7eU7YrgcBMg/92SaKrE3vj/4DDXOJps/6pz6fd+wfwU5giDUPytlqbnKW//nloEz3/d8wzb+3/bZP67G5SCAImbcUtnzDdv6aVf/M0lOrP9jc9l5RZRn7FGGICMMHTMEz/zCpCAeeOHUlJVeR5E90iw0f9Q560fVVCmONwBlWGnL7DU1B39zKimkCqbShEQcf6u1ZE5B9XKP0k0zHAShUZT/r0nuyfTxN1vNVhPjp85YXeWqGmN8IQDklKPQduaamsA/e5o/7IyFElOvLS+ehRWXFIAgNBBnBTX34WlKtwIA/oPyOKZcgx5BhaGPEAAADA1BCAAAhoYgBAAAQ9N2EHZ1dbEsq3QrjCsQCESjUaVbYVzNzc2iOPycQZBrqVSqqalJ6VYYVywWa2trk+vVtB2Ejz322Lvvvqt0K4zr9ttv379/v9KtMK7Fixf39PQo3QqDCoVCF154odKtMK6vv/562AWOJkbbQQgAAJAlBCEAABgaghAAAAzNJEkjromliKuvvvrgwYNmc0YJzbKszWYjSTLXrYJh9fX1MQxjs9mUbohBdXd3FxUVZXiwgLwkSeru7i4pKVG6IQaVTCZZli0qGnvNgO3bt59xxhmjb6O6IGRZtrOzU+lWAACAHlRUVNjt9tG3UV0QAgAA5BNuqgAAgKEhCAEAwNAQhAAAYGgIQgAAMDTNrEfY2dm5adOmzs7OK6+8csmSJUM3EEXxjTfeOHTo0MyZM++44w7U9Muos7Nz27ZtR44cKSwsvOGGG2prawdtkEqlfv/73/f/eOaZZ15wwQX5baOe7dixo76+Pv3YarWuXr166DZtbW2bN2/u6+tbuXIlpv6S18aNGwcWFQ79eNfX1+/YsaP/x5UrV3o8WGU+Kx0dHfv3729ubr744ounT5/e/3xLS8trr70WDoevvfba888/f+gfxuPxTZs2HT9+vK6u7pZbbslwcJE2rgh5nl+4cOGxY8eqq6tXrVr11ltvDd3m3nvvfemll2pra//4xz/+8Ic/zHsb9WzdunUffPCBz+fr6Oioq6vbtWvXoA2SyeSaNWuOHj166tSpU6dOdXd3K9JOvXr99dffeuut9Hvb0NAwdIO+vr4FCxa0trZWVFRcffXV77//fv4bqWMNDQ2n/u2BBx44dOjQoA127969YcOG/m1isZgi7dSTpUuXPvXUU48++uju3bv7n+zq6jr33HM7Ozt9Pt8VV1zx97//fegf3nzzzW+99VZtbe0LL7zw4IMPZro/SQs2b9587rnnplIpSZL+9Kc/nX322YM2CAQCDoejublZkqSenh6SJE+ePKlAQ3VKEIT+x2vWrFm9evWgDdJHPsuy+W2XUdx+++2//vWvR9ng+eefX7JkSfrxK6+8smjRory0y3D27dtHUVRvb++g519//fUVK1Yo0iS9EkVRkqTzzjvv9ddf73/y2WefveKKK9KPX3jhhf7PfL8jR47QNB0KhSRJamhooCiqq6srk91p44rws88+W7ZsmclkIghi2bJlBw8e7O3tHbjB7t27p02bVlFRQRBEUVHROeecs3PnTmXaqkcD5+6JRqNOp3PYzX73u9+9+OKLX375Zb7aZSB79uzZsGHDli1bEonE0N9+9tlnl112WfrxsmXLdu/ePexmkKVNmzZdf/31hYWFQ3/V0tKyYcOGzZs3Yz4QWQx7S3PQ53znzp2pVGrQBgsWLHC73QRB1NTUVFZWfv755xntLusG50MgECgtLU0/njRpktVqDQQCI21AEERZWZmMS1VBv927d2/duvUnP/nJ0F8tXbq0u7v7yJEjS5Ys2bBhQ/7bpmNVVVWlpaU9PT1PP/30ggULeJ4ftMHAz7/H40mlUsFgMO/N1DlBEN58881hO2gLCgrOPPPMUCi0devWmTNnDr13CrIY9DlPJBJdXV0DNwgGgwODwOPxZBgE2iiWsVqtyWQy/VgURVEUB02ZY7VaB65QmkgkxpxTB8br2LFj119//caNG6dNmzboV3a7/aOPPko/vvnmmy+99NK1a9cyDJP3NurTU0891f/gnHPO2bx587p16wZuMPAAST/A5192f/nLX4qKii666KKhv1q5cuXKlSvTj++5557//u//fuedd/LbOkMY83M+4SDQxhWh3+/vD/b0A5/PN2iD1tbW/h9bW1vLy8vz2ULdq6+vX7p06TPPPHPDDTeMvuXChQuTyWRLS0t+GmYoNpttwYIFQ+tlBh4gra2tNpsNk0HLbvPmzXfeeWe6g2YUF1xwwalTp/LTJKMZ9DmnaXrQbeoJB4E2gnDFihXbtm2LRqMEQbzzzjtLlixJX218/fXXjY2NBEEsXry4q6srvVr68ePHjx492n8rGbLX2Nh4+eWXr1+/ftCS0F988UX6YycIQv+T27dvp2m6pqYmz43Usf63l2XZHTt2zJ49myCIeDz+8ccfp8uUVqxYsXXr1nS/4JYtW6688kqLxaJgg/WnoaFh586dt956a/8zPM9//PHH6euS/v8gSZLef//9M888U5lW6t2KFSvefffd9DXfli1bVqxYkX7+iy++SAfk5ZdffuDAgfSZ4p49e3ieX7RoUUYvLUuFT64lk8nLLrts3rx5t95666RJk3bt2pV+fsmSJT//+c/Tj59//nmfz7d69erKysr+J0EWy5cvZxhm3r/dfffd6efnzp3729/+VpKkjRs3zp49+wc/+MHy5cvdbvcf/vAHRdurN8XFxStWrFi1apXP57vqqqsSiYQkSekj//Tp05IkRaPRCy+8cOHChatWrSopKfnqq6+UbrLePPbYY1deeeXAZ44ePUoQRHd3tyRJy5cvX7p06S233HLWWWfV1tY2NjYq1Ez9eOCBB+bNm8cwTE1Nzbx589Lf+TzPn3feeRdeeOFNN93k8XgOHz6c3vjss89+5ZVX0o8ff/zx6urq1atXe73el19+OcPdaWb1CVEUd+zY0dXVtXjxYq/Xm37y2LFjLper/+L38OHD6QH1c+fOVa6lOlRfX8+ybP+PLpcrPcT1X//6V2lpabrXev/+/Q0NDQUFBfPnz8doYnmdPn3666+/jsfj06dPr6urSz+ZSCS++uqrurq6dC9IIpH45JNP+vr6LrnkkoH1AiCL+vp6t9vd/81DEEQ0Gj148OC8efMsFktPT8/evXt7e3v9fv/ChQsxm0f2Tp482dfX1/9jbW1tuhY0fSMkHA4vXbp00qRJ6d8ePny4rKys/2O/f//++vr6OXPmzJo1K8PdaSYIAQAAckEbfYQAAAA5giAEAABDQxACAIChIQgBAMDQEIQAAGBoCEIAADA0BCEAABgaghAAAAwNQQgAAIaGIAQAAENDEAIAgKH9f7vKShSuo5uXAAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 10
}
],
"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",
"ψ = 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": 10
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.3"
},
"kernelspec": {
"name": "julia-1.7",
"display_name": "Julia 1.7.3",
"language": "julia"
}
},
"nbformat": 4
}