{ "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.143577066747 -0.42 8.0\n", " 2 -0.156035130922 -1.90 -1.10 1.0\n", " 3 -0.156768327093 -3.13 -1.56 2.0\n", " 4 -0.157045858392 -3.56 -2.32 2.0\n", " 5 -0.157055277264 -5.03 -2.93 2.0\n", " 6 -0.157056388742 -5.95 -3.72 1.0\n", " 7 -0.157056406007 -7.76 -4.34 2.0\n", " 8 -0.157056406917 -9.04 -5.95 2.0\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380295 \n AtomicLocal -0.3163466\n LocalNonlinearity 0.1212607 \n\n total -0.157056406917" }, "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{StaticArrays.SVector{3, Float64}}:\n [-0.05568221088373597, 0.0, 0.0]\n [0.05568170081166464, 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+gvaeTAAAgAElEQVR4nOzdd3iT1f4A8POO7KQZ3XvvltKWDQXKLpYh86KAInLhh6ggKigOFL2uy5UriiIguHAgiELZIEuGtEBLC917ps3eyTt+f4Tby2WUjjRJk/N5fHzScPLm5OTN+b7nvGcgNE0DCIIgCHJXqKMzAEEQBEGOBAMhBEEQ5NZgIIQgCILcGgyEEARBkFuDgRCCIAhyazAQQhAEQW4NBkIIgiDIrcFACEEQBLk1GAghCIIgtwYDIQRBEOTWnC4QFhQUWCyWTiamKAouEedAJEk6OgtuDZa/Y8Hydywblr/TBcLJkydLpdJOJjaZTBRF9Wp+oA7o9XpHZ8GtwfJ3LFj+jmXD8ne6QAhBEARB9gQDIQRBEOTWYCCEIAiC3BoMhBAEQZBbg4EQgiAIcmswEEIQBEFuDQZCCIIgyK3BQAhBEAS5NdzRGYB6l4UCJSq6VEVXaQBBATMFhEzgwwbRQqSfBGHACyEIchySBsVKukRFN+qBwgSYKEAQEMIDMUIkToRwYfVsL7CkXVOzAfxYQR1roP5spv25SJwIiRAAJgaYKKhQgwstoFhJlavp/p7I9DB0VhgSJkAcnWUIchetRvBrNbW3irokpf25SKIYCeACMQvoCEBSILcVlKqoCg2d6omMD0T/FoFEC+HPs3fBQOhqjjfQnxSRf7bQj4aiT8ei349Gxaz7p9QT4EILvbeaGvw72U+CrErCsoIR+IODoLt89NFHV65cufd5kiQxDOvSoVRmUKulZSbgxQY+bDCWheAoAAC0ANByR7IYACJpoDSBvSb6MwPg4iCIB3w57v7zHD169PLly3vjyIizLVodFBR0+fLlwMDAziQ2GAxMJrOr56KrOlxHr79Kai3gpX7o7AiU1+mLHAsFfq6k/lVIWSjw4SBsUlBnf24ajUYgEHQzu1CPwfK3j4yMjAkTJsTFxTk6I27t6tWrhYWFBw4caH/Ghuc/bBG6gmIlveoSWa0F7w5Ap4eiaBevGxkoeDwKfTwK/a2GWnWJ/JgHPhuORXm4+dUnBP1XZmbmiBEjHJ0Lt8bhcAoLC3vp4HCwRN9mocDb16iRB4mJQWjBDHxGWJej4J2mhaI3ZuBZwejQ34mNNyjSuToLIAiCegUMhH1YoYIe9BtxWUpdexRfmYTaZAgojoKVSejlaXhOLTUmh2jQwWAIQZCLs00g/O6776Kjo318fBYvXmwwGB6U7NChQwMGDDh48KBN3tTNbS+hxuQQzyWiORPxQJ6NuzEjBMiJyfiEIHTAfuJoPYyFEAS5MhsEwlu3bj3zzDO7du0qKSmpqanZsGHDfZOpVKo1a9bI5XKZTNbzN3VnegLMP01uLqLOZuOLYnqrTY8iYF1/9Kex+OJz5MYbcPdjCIJclg2q0a+++mr69OnDhw8Xi8Wvvfbajh077jsSdfXq1StXrvT29u75O7qzWi2dcZDAEHBpKh4n6vXxLCP9kItTse/LqafOkmYYDSEIckU2CITFxcUpKSnWxykpKVKpVKFQ3JXm1KlTFRUVTz31VM/fzp1daKGH/E7Mj0K/HoVx7DXgN5iHnJuCK0wg+yihsdjpTSEIguzGBrWpXC5vn8zh4eEBAGhra5NIJO0JdDrdc889t3//fqQT80F1Ol1QUFD7n99+++306dMflNit5hHur0NX5eJfDiHG+5u0Wnu/+64h4IVcfNQBau9Iizf7dotfp9N15juFegksf/ugKBfpDJFKpWlpafX19dY/t23b1tTU9MYbbzwo/ZEjR/bt2/fll1/aK4MPQZKk9o66r5PnP5fLRdGHNPlsEAglEolarbY+tj7w9PS8M8H69esHDhyoUqny8vJ0Ol11dXVNTU1oaOh9j8bj8QoLCzs5oR7DMDcJhB8XUh/foE48gqVImI7Kw45MsP4qOfk0eiILsw7PoWmaz+c7Kj8QLH/7eGg12ldQFNXW1mZ9bDAY3n777by8vA7ST5w48aWXXioqKkpMTLRLBh8Cw7A7T3gbnv82CITR0dFFRUXWx4WFhZ6ennc2BwEAFEXduHFj6dKlAIDq6updu3ZhGPbaa6/1/K3dAQ3Aq1fIA7X0halYkK1Hh3bV+jSMz6BG5ZAnJ2OhfNgWgSAHqKmpyc/Pj4yM3LNnz8iRI8eMGVNcXHzw4EGLxTJ16lRr0NLpdDk5OUVFRRwOZ8qUKfdGsr1796anp/v4+AAATp065ePjk5SU9PXXX8+YMUOn0509e3bOnDkIgjz++OOff/75p59+6oDPaUc2uNJZtGjRvn37rl+/rtfr33///UWLFlmbqxs2bDh8+DAAYOPGjbn/kZiYuH79ehgFO4mkwZJz5Jkm+mw27vAoaPViMroyCR2dQ1Zr4LQKCHKAgoKC5cuXr1ixQiwWoyiak5OTlZWFIAiXy83Kyjp79iwA4OrVq2fPng0JCaFpeuzYsbm5uXcdJCcnJzMz0/p4x44dp06dAgCsWrVKJpOVlZWtX7/e+k+ZmZk5OTn2+2wOYoMWYXJy8kcffZSdna3VarOyst58803r8yUlJcHBwXclDgwMhKsjdpKZAvP/IFVm+vhkvPMLh9rBigQUBWDMITJnNBIPv0zI/Vxto9dcIe3zXkN8kA3pd9/9UavV+/fvFwqFAICYmJgdO3aMGTMGAODr6/vee++NHDkyIyMjIyPDmpjFYn311VcDBgy48wg3btxYuHDhQ989JiamurparVZbx3+4KtvUr0uWLFmyZMldT3733Xf3pty3b59N3tHlGQgw6yTBRJHfJ+As57sHujwBJWnwyB/Ms1NoJ2mqQpDdhAmQNf3s9LP0Yt/nydjYWGsUVKvVZWVlGzZseO+99wAAGo1GKpUCAJqampYtW3br1i0Gg2EwGO7tGtVoNFwu96Hvbr0Jp9FoYCCE7E1HgKnHCH8usmskhjvrffpnE1GNgRx/GD3zCO7DcXRuIMiOJCwwLtCR138s1u3N1ayjBbdu3do+Rds6tOe1116Li4uzjtX/7LPPDh06dNcRvL292+e5sdlso9HY/k9Go5HNvh1+ZTIZiqJ3jX90Pc5ay7oxHQGmHCXCBcg3o5w3Clo9G0vMjUDGHybkJkdnBYLcEpvNzszM3L17t/g/zGYzAKClpSUiIgJBELPZ/MMPP9z7woEDBxYUFFgfp6SkHDt2jCAI658HDhzo37+/9XF+fn5KSkp7XHRVzl3Ruh+1BYw/RMSJkG0ZWE/2kbCb9WnY2ABk6jFCTzg6KxDklrZt23b48OG0tLRHH300KSnp/fffBwAsW7bs1VdfffTRRwcMGHDf2WizZ8+2DmYEACxdulQoFMbExGi12nHjxuXl5bWvlHnkyJFZs2bZ7bM4CuwadSIqM5h0hEiRIJ8Nx/pCELxt4xBsyTly2nEiZyLOhFdWENTLJk2aNHr06PY/w8LCLl26VF1d3dLSEh4e7uvrCwDIzs6+detWRUVFTEyMQCCwNhN9fHzq6uqsr8rMzDSZTAUFBf369WOxWHv37pVKpdHR0d99992QIUOsafR6/YEDBy5cuGDvT2h3sN5yFgoTGH+YGOKDfD6iL0VBAAACwBfDMQ6GLD4LdzCEoF7HYDDuGnuPIEh4ePiQIUOsUdDKx8dn6NChnp6eTCbTOuYFRdH2W4kIgmzevPny5ct3pscwzM/Pr/2ZK1euvPrqq9a5hq4NBkKnYI2CGX7Ix0P6WBS0wlHw0xisWku/dNlOY8ohCOqh4cOH3zXa//nnn7cORrUaNWqUmywQDQOh48lNYNxhItMf2TjY+eZJdBoHB7+Px4/U03DPJgjqo958802xWOzoXDgADIQOJjOBcYeIcQHIR305ClqJWeDIJOyTImp3BYyFkOtwtoFghw8fHj9+fC8d/ODBg8uXL+8gQXFx8YQJE+67117fBQOhIylMYNJhYlwg8sGgPh8FrYJ4yOFJ2AuXyGMNLvU7gdxWuZqu1TrXyZyUlLRy5creODJN02vWrFmxYkUHaeLi4ths9oEDB3ojA44CA6HDSA1gdA4xMQj50FWioFWCCPl5LL7gNJEvd67qA4K6qkkPJhwm/TjOdeOeIAi9Xg8AMJvNH3zwQUNDw/r169euXVtTU6PX6z/++OOXXnrp6tWr1sRtbW3btm1btWrVhg0bSktL2w9SW1v79ttvr1u3rry8fNOmTSqVCgBw5swZDoeTkJAAADh79qx1AdJdu3bV1tZqNJqNGzdaX7tgwYLPP//czp+6V8FA6BjNBpCZQzwahrwzwKWioNVIP+SzYVj2UdLZLqUhqPM0FvDIUWJxLCpiOTor/6uwsPCTTz4BAJhMprVr1y5evDgoKEilUo0fP37RokUAAIlEkpmZ2draCgC4ePFifX39oEGDmEzmiBEjysvLAQCtra2DBw82Go1xcXHLli175ZVX5HI5AODQoUPtEzOOHj1qXY/mk08+qaysVCqVr7/+uvWfRo8e/ccffxgMBgd8+N4B5xE6QIOOHnuIXBiNvtrfZS9EZoWjDTqQdYQ8PwUXO1k9AkEPZaHA7JPEIG9kXX/0yD3/aq4pVv7ymX1ywoxIEj26tIME77zzzoABA5566ikvL69x48ZZB4IeO3bs9OnTs2fPnjJlypQpU2iaViqV9fX1P/zww+uvv759+/ZRo0b94x//AACkpaUlJSVZD1VYWDh16tSHZsnb25vD4ZSXlycnJ9viIzoeDIT2Vq2hxx0ml8WjLya7bBS0ej4JrdPR048TR7Nwtgu2eyGXRQOw5BzJRJHPht//xMV9Q0RznrNPZlAOr+ME1jCGoqiXl1d7SPPx8ZHJZACAoqKip59+Wi6XCwSC5ubm7OxsAEBpaWlqaqo1ZUJCQvvKpXq9nsPp1MLBXC5Xo9F06wM5IxgI7apMRY8/TL7UD30mwcWjoNWHg7DHT5MLTpM/jekbK8ZBEADg9VyyWEWfmow/aFYvyuYyg6Ptm6kHwrD/RmvritsAAARBrAM7n3/++SeffNK6L/orr7xijY4ikah9xW2dTmdddwYA4O3tbe0jBQBwOBylUtl+ZIPB0B4jKYqSy+V3Tt7v69yiOnYSN+R05iHyjTR3iYIAABQBu0ZibUb6BTjRHuojvrhF/VxFH5iAc12imdAesZRK5c8//2x98pFHHtm9e3djYyMAoH0IDABg8ODB+fn51sdpaWknTpxob/b98ssv6enp1selpaUCgSAiIsJun6K3uUuN7HCXpfT4w8S/BqNPxbhXmbMw8Ot4/FQj/U840R5yer/XUBuuUUcmYd6ust3C6tWrn3766aysrIyMjLS0NOuT48aNW7p0aWpqalhYGEVRHA7HujfhnDlz2rehmDx58tSpU2NjY0tKSubPn79///7NmzdbX56TkzN79mwEcaFOHtrJBAYG1tfXdzKxXq8nCKJX82MTx+opn+/Mh2opR2fExtRqdSdT1mup0B8s35eTvZofd9P58oc6489myuc7c27r3b/TESNGnDt3ziFZui+CIAwGA03T1i7K9ueVSqXFYrE+1mq1RqPR+lgmk+Xn5xsMBoPBoNPp2tNTFEVRVFFRkVAopKjbn3r27Nl79+5tT6PX6+Pj43/99df2Z0iSTE5OLi4u7rXPd38HDhzIzs6+8xkbnv8u0fh3br9WU8v/JPeOxUf4udAFVBcF8pBDk7AxOYQnC5kY5L7lADmtcjU9+yT5zSg83cvZz08Mw6z3BREEuXNFtDuXCeXx/jvERiKRSCSSuw6yZMmSpKQki8XyxRdfrFu3rr159/777985WZ7D4bDZbJFI1P5MeXn5smXLYmNjbfqZHAwGwt7170Jq4w3q+GQ8Sezsv67eliBC9o7DZ5wgcibiA5y+roHcSr2OHneI/GAQ6j5XaY8//nheXh6CIN98882wYcPan4+IiHj++efvTPnkk0+GhIS0/xkTExMTE2O/jNoFDIS9hQbgravknkr63BQslO8uv66ODfdFdmTgU44Sp7PxWCEsE8gpqMwg+yj5bCI6P8qN7t+PHj36zk0NO/Dcc3aaKOJAbvTF25OBAHNPkn800uen4DAK3ik7BHl3IJZ1hGzUw0VnIMfTE+CRo8T4QGS1q8/rhToAv3vbazaA0TkECwPHsuCiKvfxVAy6NA6deJhUmBydFci9WZePifJAPuz7e79APQG7Rm3smox+9Dj5VCz6eiqcQf5Aa1LQNiP9yFHi+GScB89ByBEoGjx5hsRRZHtG39sNm6KoL7/88tq1a8HBwStWrLhzMAvUDbBFaEs/V1KTjhD/HIy+AaPgw3w4GIsXITOOEyY41R5yhGcvkg16+qcxGN4Ha8GlS5cePHgwKyurpqYmMzOTJOGvqEfg1bhtEBRYl0vuqaKPZ+H9JDAIPhwCwJcZ2N9OkY/9QfbRygjqu165Ql5ppU9M7uYquDfbSjdd+cLWmbq/fj4JK9KfvvOZ+vr6PXv21NXVCQSCadOmJSQkHD582LqIKNQ9MBDaQIsBzDtFsDBwZTruCW8KdhqGgO8zsenHiafOkrtGwcVIITv5x3XqYC19Ohv3YHTzCKHCoNWDO9rG3Yb4jLsX3S4pKYmMjBQIBAAABEHS0tJKSkpgIOwJGAh76mQj/cQZckks+noqCqvyrmKiYO9YfPJRYtmf5NYRfe9WDdTnfFxIfV1GnX6kR9esPAY3VhJlu0x1jV6vV6vV7X+q1WoPDw9HZcY1wA6p7rNQYF0u+cQZ8ptR2JtpMAp2EwcHBybgNxX08xfhfQ6od225SX1aRJ2cjPlzHZ2VnikvLz99+jQAoKqq6syZM5mZmY7OUd8GA2E3FSvpYb8T+TL62qP4mAAYA3uEzwA5E/HLUvqFSzAWQr3ly2LqwwLq5GQsiNfnf7D9+/ffsGHDsGHDBg4c+O6770ZFOax56hpg12iXkTT4pIh67zr5zgDs73HwSsI2hExwNAsff4hYfZncCCd1Qba2rZh69zp1ajIWJujzURAAwOPxTp48WV1d7e3tfeeyolD3wHq8a4oU9IgDxIEa6uJUHEZB2xIxwfHJ+NkmevVlEq46A9nQtmLqnevUqclYpIcrRMF2YWFhMAraBKzKO0tPgFevkJk5xKIY9OQjuIv9opyEiAmOZeHnm+nnLsBYCNnGpzepf+S7VBSMj49ftmyZo3PhUmAg7JRfqqjEvUStDhTMZPw9Dg6L6UViFjg+Gb8mo5eeJykYDKGe+ecNalMhdfoR14mCAICoqKj58+d3++U0TRuNRpvkxLrloU0O5VgwED7ElVZ61EHi3evUzpHYd6MxP46jM+QGPBjgyCS8XEUvOE1a4Lb2UHe9mUfuKKHOPOKC278UFhbm5+d3nObKlSulpaX3Pn/9+vW4uDibZCMhIeHatWvdfjlBEHv27KEox//IYSB8oJtKeuYJcsYJcmE0mjcdH+3var8lZ8ZngEOTcI0FzDxBGuFIUqiLaABWXiIP1NJns/HAvj9G9F4//vjjt99+23GaLVu23LnFrhMyGAxz5sxxhvXh4KjR+yiQ0/+4Tp1uol7qh303GuPAQnIENgb2jsOePENOOkL8Nh4XMh2dIaiPsFBg8VmySkufegQXueJpc+vWrVOnThEEsXbt2sjIyCVLlhAEsXXr1tzc3MDAwBUrVvj5+V2+fDkvL6+2tra1tTU1NXXu3Ln3PZTZbP7888+vXbsWGhq6YsUKb29v6/OnT5/ev3+/XC5PS0t7/vnnSZLctWvXX3/9RRBERkbGE088gaIPbERt3LgxMzNzz549Uql01qxZWVlZ1udLS0u3b98uk8lGjx49f/58BEG2bNkCAFi3bh2Kos8880xwcLCti6qzYIvwf5xooCcfJbKOkAO9kYq5jNXJKIyCDsRAwbejsRQJMuog0aR3dG6gvkBHgGnHCaUZHJvkmlEQACAUCv39/X19fdPT062bxc+fP//AgQOzZ88mSTI9PV2pVHp5eXl6egYFBaWnp4eHhz/oUDNnzjx16tScOXM0Gs3AgQN1Oh0A4Kuvvlq4cGFqaupjjz3W1tZGUZTRaCwtLZ06deqMGTN27NjxzjvvdJC9bdu2PfbYY4mJiRMnTnz66ad///13AEBFRcXQoUN9fHymTZu2adOmV155BQBg7aFNS0tLT0937PBXWM0DAIDSDL4to7YWUwgAq5LRX8ehLDiTzTmgCPj3UOy9fGr4ASJnIhYvcsFuLshWpAYw9TiRKEK2jujdZdzVVfrS3fW9+AZ3EEXzouYE3vlMQEBAfHy80WicPXs2AKCiouL3339vaGgQi8WTJ0/Ozc3dtWvXypUrw8LCkpKSrGnu68aNG2fPnm1oaODz+ZMnT7506dLu3buXLFny2muvff311+PHjwcATJo0CQDA5/M//PBDg8HQ0tKyZs2atWvXvvHGGx3kedmyZdbhPDqd7p///OfUqVM3b948c+bMF198EQAQFRWVmpr65ptvjhkzBgAwc+ZMBqO7q77aiFsHQgsFjjfQ35VTh+uoScHop8OwUf4IrGid0CspaDAPjMkhfhqLj/SDXxF0HyUq+pGj5IIo9I20Xh/XzQ9iJy0L6+U3uQ1lPiSkl5eXh4WFicVi659paWn3HSNzr7KyspiYGD6ff+cLlUplU1PTkCFD7kyp0+nmzZtXXFwcGhpqMBiampo6PnJKSor1Qf/+/a2ZKSsra18WPCEhAcOwmpqawMDABx7CvtwxEOoIcKKB2l9DH6ihYkXI/Cj002EMCdw1wrnNj0IDecick8RHg7EFUbBLH/ofp5voeaeI9wdhT0Tb49xAGSjb08Edr+3zFkQikVKpbH9eoVC03+rrmEgkUigUd74wICCAz+czmUyZTGbd3cJq586dGIZZQ9qVK1cmTJjQ8ZHb86NUKq0RWiwWtz+p1+uNRmN75HYG7lKhmEhwvpl+9zo1/jDh/73l05tUqidyfQb+5xT8/+JRGAX7hEx/5I9H8LeuUq9egVMMof/aUULNO0XsHoPbJwo6A29v79raWuvjfv36oSi6Z88eAEBdXd2+ffus41PuTHNfAwcOVKvVOTk5AIDKysqcnJxJkybhOD516tR33nnHOpiztbWVpmmDwYAgCACAJMlNmzY9NHtffvklQRAURX3xxRfWztWsrKxvvvnGGgs3b96clpbm6+vL5/M5HE7HmbQP12wRkjRo0NHlalCkoAsVdF4bfUtJJ4qRkX7I84lY5niE55qf2/XFi5DL0/AZx4mZJ8hvRmMCB99ZgByMoMBLf5GH6uiz2Xi00I36zOfOnbt3796wsLARI0Z89913P/zww5NPPvnWW29JpdK1a9eOGDECAPDUU08tXLgwLCxsxowZ//rXv9pfi6Iom80GAAgEgt27dy9evFggELS2tr7zzjtpaWkAgM8++2zRokXW7laj0Xjr1q0nnnji66+/jo+Ppyhq9uzZIpHIeig2m33f4aPR0dFJSUlmszksLOyLL74AAMybNy83Nzc2NlYsFrNYrN27dwMAEATZsGHDuHHjUBTdt29fe4eq/SHOti5AUFDQ5cuXO9l33Kgy/rsERxBEawFqM1BbgNRANxtAo572ZiNRHiBehCRLkFRPpL8n0r2tqKEOaDSaO/tP7MZMgecukOdb6P3jsSgXWjGkqxxV/k5CZgJzTxIMFPwwpncHiGZkZLz33nvW6OLMWltbJRIJhnW5ppNKpZ6enne90Lrroa+vL/KfgRONjY0SicQaRDsQFxe3c+fO1NRUrVbr5eV15z9ZLBa1Wu3p6dnVHAIADh48+Nz7Xwzd8JsfF1jX5bfh+d+3W0YYQosYAEWRIB7wYAAhE3izUX8uCOAicNinC2Oi4IsR2NZiasQBYnsGnh3ivrHQbV1to2efJGeFI/8YiMENna06eWvwXj4+Pvc+yeVyudz/2bYxICCg88dks9n3hkwGg9G9KGjlzwWTg5HeWIa4bwdCDwZ4uR+CYe5yYwC609I4NEWCzD1FXm5F1qfB2tCN7CihXs0ltwzDZobD377TefvttyMiInrjyBIW8njvDJSDpxHUhw3xQXKn45ek9PhDcMa9W9BawMLT5KZC6mw2DqOgc5ozZ46vr6+jc9E18EyC+jZvNjg6Cc8MQNP3W47WO9cNb8i28uX0wN8IJgYuT8Nj3WloDNTbYCCE+jwUAa+noj+Owf9+nlx1CS7S7YJoAD4upCYcJt5IRbdnYNy+fUsHcjowEEIuYqQfcv1RvFEPBv1G5Mth09B11OnoiYeJvVXU5an4vEhYZUG2B88qyHWIWeCnMdhL/dAJh4l3r1OE47c5g3rq6zJqwH5itD96JhsPE8DuUKhXwC4GyNUsiEIz/ZGnz5H7q6kdI7F+Elh79kn1OnrZebJBD45n4Y79EtPT0+fOnXvvZACaphG4OLG96PX64cOH99LB+/aEeoPBwGQyuzGBFLIJZ57QTQOwq5Rae4X8exy6rj/mksspOHP59wRFg63F1Jt55HOJ2JoUlOHofqsHLTOt0+kcu3mQu/Hz87tzaiOcUA9BD4EAsCgGnRSEPneRTN5LfD4cGxcIL977gHw5vew8iaPgdDae4By7bnE4nPtOjHPVCxE3BAMh5Mr8uWDPWCynjv77eXKgN7JxMBrEc4q6FbqXygzeyCN/rKTeHYAtju31rZQgqJ2jOx0gqPc9EowUzcTjhCD1V+If1ykD4egMQf+LosGOEir+F4uRBDdnMZ6GURCyLxgIIbfAwcFb6dhf0/DrMjr+F2J3BdzHyVmcbKTT9xPflFEHJuBbR2CecE80yO5g1yjkRsIFyM9jsfPN9OrL5MYb1AcD4Y1DR7ouo9deISvU4P2BKFwvDXIgGAghtzPCD7k0Dd9bRT1zgQzkgncGYMN8YTi0q2Ilvf4qdaaJei0V+3uc48eFQm4OnoCQO0IAmBWOFs3E50ejj58ms44QF1pgX6k9FCvphafJUTlEf0+kfC7jmQQYBSHHg+cg5L5wFDwVg5bMxmeEofNPk2MPEVVfx3UAACAASURBVMcbYDjsLVfb6L+dIkfnEHEipGwOY20KyoMdUpBzgGci5O6YKFgShy6KQXdXUC9cIpkYeDEZnRUOWyq2QQNwrJ7eeIMsVoJVyej2DAaf4eg8QdD/goEQggAAAEfBwmh0QTSaU0t/XEi+/Bf1TAK6OBb1vntdLaizdAT4rpz6pJBiYmBVEjovEl5bQE4KBkII+i8EgOwQJDsEz5fTm4uo2D2W7GB0aTw6HI6m6YoiBf1lMfV9OTXKH/10OJbpD0sPcmowEELQfaRIkO0Z2EeDsJ2l1JJzJALAU7Ho41GoH8fROXNiagv4uZLaWUrVaMFTMci1GXgwXMcH6gtgIISgBxKzwAvJ6AvJ6LlmemcplfCLZZgP8lgUOi0UDvT4LwsFjjXQu8upQ3XUmAB0bQo6ORjFYASE+g74a4agh8vwQzL8sM0E9ms19X059cyf5MQgdHY4khWMuu1u6QQF/mii91RR+6upGCHyWCT6yTAGXBcG6ovc9UcMQV3Hw8H8KHR+FCozgV+rqS+LqcXnyEx/dGooMjkY9XWPXlOVGRxroH6voQ/XUdFCZHY4mjsdD+HDBiDUh8FACEFd5skCT8eiT8eiChM4VEf9XkuvvmyJFCCTgpEJgehQHwR3reGRFA2uy+jjDfTReiq3jR7hh0wNQT8YhAdwYfyDXIEtA6HJZGKxOuoZsVgsDAacQwS5DjELPB6FPh4FCAr7s4U+1kC9cIksVdHDfJFR/miGHzLAC2H1zT2BSRoUyOlzzfSZJvpME+XDQSYEIS/2w0b7I27bGwy5Ktuc0fX19fPmzSsoKGAymR9//PH8+fPvSvDmm2/u3LmztbWVx+OtWLFi/fr1NnlfCHISOApG+SOj/LF3BwCFCZxtpk430SsvUsUqOkWCDPJGBnkj6V5IlBBx5jZUnY7Oa6NzW+lLUvpKKx3EQ0b4ITPDkU+HMfy5D385BPVRtgmEK1euTEpKOnPmTF5e3pgxYzIzMwMDA+9M4O3tffz48djY2OLi4pEjRyYmJs6ePdsmbw1BzkbMAtNC0WmhAACgtYDcNvqylN5TRb+SS8mNdD9PJFmMJEuQOBESL0IcOB9DbgIlKvqmgr6ppAvkdL6MRhGQ7oUM9EZWJ2ODfRAJHPkCuQeEpnu6uKJSqfT29i4rKwsLCwMAZGdnjxw58uWXX35Q+lmzZiUmJr711lv3/degoKDLly/fFUcfxGAwMJlMDOubfU99n0ajEQgEjs5FXyI3gXw5fUNO31TQt5T0LSVtJEGUBxIuQMIFIISPhPBBABcJ4AJv9sP7VDtT/gQFpEa6xQAadKBOR9dq6WotqNLQZSqapEGsEEkUI3EiJEWCJEsAvOfXJfD8dywblr8NWoQ1NTVMJtMaBQEACQkJFRUVD0osl8vPnz+/bNmyByWgaVqlUnG5tztieDwek8nseSYhyBlIWCDTH7lzpRWlGZSr6SoNXa0B5Wr6VCNo0lONeiA10CwMeLERCQuImEDIRDg44OGAgwP2fwKk2YwzmaT1sYUCWgswksBAAIWZVpmB0gzajLTGArzZwJeDBPFAEA8J4iFTQkC4AI3yQODqcRBkZYNAqFQqeTxe+58CgaCqquq+KS0Wy8KFCydOnDhu3LgHHU2v1w8dOhRFb4+627p16+TJkx+UGLYIHUun0yFOfc+rD8ABiGODODYA3nf/k9qCyExAaQZKM6K2AAMJDCRiIICJul3mOG1mgf88xugwNmCiNBcHQgYQMoGQQUtYQMJ8QJcPAbTa3vpQbgKe/47VyfLncrntAeVBbBAIvby8VCpV+58KhcLX1/feZCRJLliwAACwbdu2Do7G4/E63zWKYRgMhA5E0zSfz3d0LlwWH4CADhNoNBaBADbrHAae/45lw/K3wXSnsLAwHMeLioqsf167di0+Pv6uNBRFLVq0SC6X//LLL7CrE4IgCHIeNgiEPB7vscceW7duXVNT0w8//HDt2rV58+YBAAoKCrKzs61pli9ffvLkyWefffb8+fMnTpy4detWz98XgiAIgnrONtMnNm7cuGrVqmHDhvn7+//2228SiQQAQNM0QRDWBFqtNiEh4ZNPPrH+OXny5HtbjRAEQRBkfzaYPmFbcPpEHwKHjzsWLH/HguXvWDYsf9daEhGCIAiCugguGuiyFEZlYeutwtbiGnVdk7ZFZlBYKMJIGD2YAiHbw4frFSkOjxaHp/v19+SIHZ1ZCHIvWrPuaktBiaysQlndqG1RGdVqsxpHGSyMKWR5+PN9gz0CErxik70T/Hg+js6s64OB0NW06WXHq8+cq7tUq65P8o5L9IqbGj3Jn+fryZUwUQYbZ6tMapVJ06xrqVBU/1n/1+a87T5cr5HBwyZFjPHl3TOXDYIg21GbNMerz/xRc65SWdPPJyHeMzY7amKQIEDIEghZHhaKMJNmhVHVrG2pVtedr7u85epOIctjZPCQsaEjQ4XBjs6+y4L3CF1HXnP+r6U5+S1FI0OGjg4ZnurbD0cfXjIUTd2SlR6vOnOq5lysZ9S8+Blpfv06+Y7wHoljwfJ3rC6Vf7mi6oeb+y415g4NGDghYnSqTzIDe/hWPBRNF8tKz9ReOF59OlgQOC0ma3TIcBSBt7QAsOn5DwOhK7jQ8Ne3hT/rLYbZcVPHho3i4N2ZZG0hLSdqzv54cx8bZy9JWTDAv/9DXwIrYseC5e9YnSz/MkXltuvfViirZ8dNnRI1kcfozkYeBEWer7+0t+Sg3KCYnzhrQkQmhrh7vQcD4W0wEN5sK9ly9SsDYXwi+W8jgoagPV7wiQb02dqL2/O/9eF5P5O2OEIU2kFiWBE7Fix/x3po+bfq27Ze++ZqS/4TyX+bHDmegdrgVlS+tOjrGz+26mVLU58YETS45wfsu2AgvM2dA6HcoNhydWe+tHBxyvwJ4Zk9D4F3ImnyYPmxnQW7J4RnLur32IOamLAidixY/o7VQfmTNPlL8YHvi36ZHjN5XsKM7nXSdOCvxqufX9spYgtXDVwW4hFk24P3FTAQ3uaegZCi6f2lh76+8WN29IQFiXPYeG/tGqc0qb64uutqy401Q55N90u5NwGsiB0Llr9jPaj8K5TV713YJGaLVg5cGijw76V3p2jq19JD39z4KTtq/BPJf2Nibrd0JQyEt7lhIKxVN3x46RMUQV8cvCLEo1Ol1EN5zfkfXNo8JCB9edoi9v9e2MKK2LFg+TvWveVP0dT3Rb/sLTmwLHXRpIgxdsiD3KDYnLe9XFH58pDnkr3da7kuGAhvs1UgJE2UWWWx6EjSRJEmEgCAMlCMibLEDJaYgaBOsdMKDehfig98V7jnyX5/mxY92bZ9oR3TWfSbc7cVtRW/PvzFGElk+/OwInYsWP6OdVf5t+haN/y5kYUzXxnyvBfX0545OVd3aVPu1rGhGU+nzHeepqFJaTHKzJSFtlaqGBNFmSiDjzP5OM6zQesFBsLbuhEILVpC32IytJqNbSZDq9koM5sUFpqkmUKcwcdRJoqzMQAAaaYoM2WUmy0aghvAFsfwxfECj/DuDPeyCam+7b2Lmywk8eqwlQF8P4fk4WTNuc25X85PnDMzLhsBCIAVsaPB8nesO8v/XN3FjX9tmRv/6Nz4R+15kdpObdJs/GtLjbr+tWGrosQR9s+Alb7JKL+pUZRoNbUGnIWyvZgYC8VYtytV0kRZtIRFQ1AkzRYz2J5MtheT483ieDM53iyWiAG6UnIwEN7WQSCkCNqsspgUFqPCbJRZjG1mQ5vJKDUDDHC8WVwfFsebyfZisb2YbDED5z4wlFIEranRK0u1sgI1giEBGZ4+A0QIZtcT/XTtn5uufDE7btq8hBmOnULUpG158/wH/jzfl4c8y2NwYUXsWLD8Hcta/gRFbr3+9dnaC29lrInzjHZslo5Vnd5ydcdjibNmx01FuhRVeogGskJ14zmZsc3sleIhiuF7RPAw1gMrK9JEmeRmg8xsbDMbWs2GVpNBaiKMFMeLyfFmsj1v18wsMYMpYmDM+x8HBsLbNDJd80klgiA0RZNGijRThIEk9KRFQ1AWiunBsPZt3r7u8GJyvFjdb5LTQFmmbTjdZlJYouYE2qd1aCRMm/O2XW8pfH34aof/xqwspGVz3varLQUbRr7ihYphRexAMBA6lkajseDEG+c+4DO5rw5b5cF0iu+iSduy4c+NPAb31WErxWyRHd7RIDWV72kkzVTQGC/PZI9u30siTZSh1WRsMxvabvfVmZQWs9ICAGAIcAYPw7kYxsIYAjxyhj+AgbCdVqlTFxlQFAUIgnNQlIHiXIzBxRgCHOf01ggaWYG6cn+TZz+P8Kl+vXr7sEpZs/78h7Ge0asGLrP58Osesl54Lk16Iit2nKPz4r5gIHSsv2qufnj106nRExckzbFr8+thSJrcVfDjocoTrw5ded/x3jbUeFZWd6I1ZIK3/3DPXioD0kSZ1QShJwkDSZpImgLeqUIAA2E7R40aJYxkyTd1AEHiFgZ30PzviQPlR7df/2552qKJdhl71g3liqp1p98dE5axpP8CuOaTQ8BA6EC/lx3Zkf/9q8NWDg5Id3Re7u9ay413L3w8MSLzqX6P9cYyNDRFV/7apK7SJywOZYkfvlyczcFAeJsDp0/QFF2xr0lTo0/+v/AObjF2g86i33j5sxp1/ZsjXrbPBIlua5Q1/fP6FhzF3xj+Ip/Jc3R23A4MhA5BUOQnuV8WSIteHbAyxi/K0dnpiNKkevfCx0bC9Mbw1d5cLxsemSbpW7tqaZKOeyKklxoDDwX3I3Q8BEWiZgWIYvg3v6qlCJtdTJTKK5YcXsVn8j+f+JGTR0EAgIDJ/2jM+mCPgP87+mKtusHR2YGgXqc0qV44+ZrMIN8y8SN/nq+js/MQIpbww8w3hwSk//3I6osNuTY7Lg3Kfm4AACQ8HeqoKGhbrvAZHCg8248tYRR/XUtTNoiF+0pyXv5j/ZKUBS8M+j/nmQzUMQzBnk1fMi9h5nPHX/mr8aqjswNBvahCWb3syIspPkkbRr7KZXAcnZ1OQQDyeOKstzPWfnzl88+v7iQosufHrDncom8xxS4IdpI51j0HA2HPICBqTiBppOpOtPbkMBqz9rWz7x2pOrll4keZoSNslTu7mRw57p2Rr35w6ZOfb+13dF4gqFecq7u4+uTrf++/cHHK4w6ZKdgTyd7x27M21aobnj2+pknb0pNDteWrW6+rEp8OfdCshr7IdT6Jo6A4EvdEcPMFuabG0L0j3Gi9tfjQSn+ez2cTPnTUZPmeS/KO+3zSP49Xn3nv4iYLaXF0diDIZmhAf33jp8152z/MXD8mNMPR2ekmD5bgH6PXjQ0d9X9HX/yj5nz3DmJWExX7GmPnBzP4LrWpOwyENsDg45EzA0q+qyNNVJdeSNHU1zd+fOPc+6sGLnsmfbFNdmlxIB+u1+bx75tI8/Mn1skMCkdnB4JswEgY15/78K+mvC8m/vPO9QX7IgQgs+KmfJi5fkfB9x9c+sRIGLv2ehqU/dQQMFwiCOkb3cKdBwOhbXgme3iEc6tzutDn0KyTPnf81Xxp0fasTUMDB/Re3uyJjbPeHPHS0MABy46sLpaV9d4bUQStbzEpS7XW//TNpq5ehUB9EUXQBqlJVa6zfu+6BiNp7sXvvVknfebYGi6Ds2nsuxKOuPfeyJ5iJJHbsj6maXrJ4VUlsvLOv7DpopzQk0HjvHsvb47St5sgTiXiUf+898v8hop5/g+f/H608tSWqzvnJc6YEze9z91v6BgCkAVJcyJEYWtPv/1/qTaeB2mQmtry1bIClb7FxBIzby9OSAOz2mJUWHAOJgjleIRxRbH8znwLUJ9glJuVxVp1lV5TqzcpLCwxgym8vRS+RUMYWk0sMdMzWeDVT8i3aUvlWsuNt//85/zEWTNjp9jwsM6Ag7PXDn3+j5rza06/NSt26mOJMx86FZjQk7VHpMnLw11mgMyd4DxCW2o8J5Pf1CQtDesgjdKk2nh5S52m8bVhL0SJw+2VtV7R8TyeGlXdurP/GOSf9kz6Uz2fz2uUm2sOS1WlWq9UobX9fe8P0ig3a6oN6mqd4paWJmnPZA+vVKFHKNeZFv2wJdeeR6hvNrVeU8oK1BY9KY7jCyN5glAu14d197dJA22jUVagar2qYnsyw7J9+UE2CId7in/ffXPvG8NfTPVNflAaFyj/Vn3bexf/bSRMrw5bGSQI6CBl5f4miqCjZnWUxs7ghPrbnC0Q0iR99aPyiGl+4vj7fz1nai/8O3frxPAxT/V7jIE5YC0G23roiag16965sFFvMbyVsaYnyx42nGmrO9EakOEZOMqrk/OWDFJT23VV63UVZaZ9Bop8BojYnn1jRkrnuUBFfC+LlpDmKaVXlISB9Oov9E4R8oM5nbmUoUm65bKi9pjUK0UYPtWv2yvjGwnTR5c/rVXXvzPyFV+eTwcpXaP8aUD/WnJo140fnkie+2hM9n07qAyt5oJPKtLWRDvVGBkYCG9ztkAIAJAXaapzmlNfjLqrvaIwKjdd2Vqlql0z5NlErzhHZc+2OnMiUjT99Y0fcyqOrx/xcpJ3lz84aabKf2owyMzxT4awRN25dNA1GFuuKFqvqngBbL+hYs8kDztvHtJ7XKMivo0GyjJt80W5slQnSRL4DhILI3jdaMqTRqp0d71FT8Y9EcwUdLnWbtA0vX72vWhJxAuDlrMeNpfXlcq/XtP4waXNNE29POS5e5fyuLWzVhDKCRrjXHcHYSC8zQkDIQCgYHNlQIanV3+h9U8a0IcqTmy7/s0jkeOfSP5bX5kp3xmdPxEvNuR+cOmTxxNnzYqb0vnliUkTVfhFFdeXHTkrAMV7FL0ogpYXqpsvyvUtJt9BYr+hEoesjmhbrlEREzqy5Yqi6YIcZ6F+QyXeaSKM3bNBfDSoO9HaclmRvCK8SxdP1j0FF/V7bFp0VmfSu0b5t6No+reyQzsLfpgZO+WxxJntg9h1jcaibTUD1sX08DdoczAQ3uacgVB+U1N7RNr/hUgAQJWq9l9/fU5QltWDnunrdwTv1aUTsVknfePc+348nzVDnuMxHr6JFUXQN7fXsCWMqNmBNrzJZ5Cami7IpXlKjzCu/3BPcSy/795B7OsVsabW0PSnTF6okSQJ/Id5CkJtOdSl4UxbyyVF8ooIRid2XrNQxJfXvj5bd7FLewr29fK/rxZd679zt9Zrml4YuKy/bzIAoPT7eq4/O2iMLZcqtQkYCG9zzkAIaHD1o7KAbMlew/7j1acX9XtsStQkFxsaatXVE9FCWrZc23mpIfeNES/Ge8Z0lJQGxd/UAoDELgjqjVFqlJlqvaZq+lNOGEi/YRLfQeLOVJfOpo9WxKSZar2qar4gJ4yk/1CJ7yBx93cJ7VDNoRZFiTb5mfCO10Bp0ra8df4jT454zdDnurSnYB8t/844V3fp07ztiV6xSyIW1nwuH7gutqfN9F4AA+FtzhkIKZo6nvNn21VV08SaJf0XCFkejs5Rb+neiXiu7uK/rnwxJ27a3PhHH3R90HhW1nZdlfxMeG/fz9PUGpovyGU31OJ4vt9QSffuSzlKn6uIdY3G5ovy1msqYQTPb5jEDs3x0t31KAONmv3AsY4na85tzt02P3H2zLjsru4p2OfKv0tMpPn7ol80R82RXqHjFgxlO9mWqAAGwnZOGAjP11/efv1bL5bnzD/nJC+J5Ac63dljQ90+EaX6tnf//BeKoq8OXXnv7jC6RmPhF9UpKyPYEjvdTyUMpDRX2XxRTpPAd7DYZ6CoG+Ms7K+vVMSEkWy7pmq5rDCrCd/BYt/B4u6Ne+oG0kRd21gePsXPM/nu61GdRb/pytYSefnrw1dHiyO6cfC+Uv7dRujIv/5Rcn7MH3+p855MnpsVOa439jXsNhgIb3OqQJjXnL8j/3sTaVrSf8GQgAH1f7QZpKbouc6+lVJP9OREpGh69829vxT/tjxt8YTw0f99nqCv/6siaKyXT3r3p1t0m6ZG33xJIStQe4RzfQaKJIkezjZA4E5OXhHTFK0q07XkKhU3NaIYnu9gezQB76WpMdz6qqb/6iimx38vbq613Hj/4r8HB6QvT3uKjbO6eWTnLv+eazjdpm8yRc8LLJaVfXn9G6m+bVG/xzJDRjjJjR4YCG9zkkCY23T968KfVCb1E0lzM0MzrGeJRUvkvVc28I1Y19iv6756fiKWKSrfvfBxsCDghUH/Z51oWHNEapCa4hYG2yiP3UGaKVmBWpqr1NYbPJM9vFOFwiieEy6o4aQVMQ00tYbWa8q26yqWiOGdLvJJE/XSXcBOqj0q1TcZ454MAQAYCdO2/G/O1F58efCKQQFpPTmsk5a/7Vz9oCxqbqBH2O2hbXnN+TsLdmvM2oVJczNDRzx0MZreBgPhbY4NhCRNnq29uPvmXgtFzE+cNSZ05F0XSsW7akVxAr8hLrJE4b1sciJaSMuuGz/kVJxYnvbUKPHw6xvL+6+OslvXWcfMasJaoRtlFs9kD68UD2Ekz3mmITpXRUwDTa2+rUAty1ejDMSrv9A7Tcjx7mZjy7Yogr76fln0vMBKfuWHlzYnecc/O+DpLo2LuS/nKn9bU1fpy39uSFtz9xja3Kbr3xT+1KqXzY2fPilibLfb0z0HA+FtjgqEapMmp+L4r6U5/jzfuQmPDg0ccN/b7Ipbmtpj0pTn+/aK9R2w4YlYKq/44NIn44snJkfFJU7tzg2bXmWUm2X56rYCtaHNJI4TeCYKRLF8nOPgrghnqIgpM6Wq0MkKNfKbGpyDevUTevbz4AU43a3xuivNRUeqdibvWDVo2ZAA26xx7wzl33vKfmjgBrACR91/1kRRW/EPN3+9Ib05OXLc9Jisjpfg6SUwEN5m50BIA7pAevNA2dGLjVcygofOjM3u+B47TdG575QmLgnluugC0LatCJTVuvwdZZ8k//vRhMlz46c75xJ0ZjUhv6mRF6lVFTqeP1scyxfF8PkhHId0nDqsIqaBrtmoLNUqS7TqKj0/mCNJFHgmerC9nHGxCBrQRypPbbv27dKSpYljI4MG26zKduFASBqpKxtK0l95yJpqTdqWfaU5RypPJnrFTomaNCQw3Z6jaWAgvM1ugbBe03ii+szRyj9YGPORqAkTIzI72a9Sc0RKGsmI6f69nUOHsG1FcGNLlc9AEZ1g+TRve5Wy9tkBT9vqyr03UAStrtApSrTKUq1JYfEI5wojeR7hXH4wx259p/asiGmK1jeb1JV6VaVOVaHDWKgohi+O4Yti+E44w6xdiax8U+5WAMDKgUsDNIHFX9cOWBdjqy/IhQNh80W5slQb90RIZxKbSPOpmnM55cfqNU3jw0aNCx8VK4nq7RwCGAjb9XYgrFU3nK+/9EfNeZlBPjpkxMSIzK5+wUa5Of/flYPejHXCoRY9Z8MTUVNjKPm2Lv3VaGtB/dV0dXPuNj+e7/K0ReGiUJu8Re+xaAlVpV5doVNV6oxtZl4AWxDKFYRw+MEctoTZe+Mke7siNqss2nqjplavqTVoawwMD9wa74VRPCe5iduBVn3btvzvcpuuLem/cFLEGOvNi6Ivq736C30H2ea2vQsHwoLNlcHjvB+0ecCD1Gsaj1X9caL6LAKQ0aHDRwYNjfGM7OrszM6DgfC23giEZtJ8o/XW5ca8iw25BsI4ImjwqJBhKT6J3R4ilb+pIvQRX1E034aZdBI2PBFv7aoVRfH8R3i2P0NQ5G9lh74t/Hl40OBFyfO8uJ4dvNx5kCZKU2vQWoNHnYE0UbwANi+AzfVnc31ZXD+WDe8s2rYiJk2UQWrSt5j0TUZto1HXaAQ04Aex+SFcQQhHEMrtK4vv6Cz6H27u+63s8NToSY8nzOIy/rtym7JMV7mvMe3laJtUzq4aCM0qy7V/lg9aH9ftpnOJvPxM7YVzdZeMhHFwQPrggPR0v5Q7vwibgIHwNlsFQiNhKpaV5ktv5ksLb8lKI0Sh1i8vRmKDy5mGP9oMMrNT7eNlK7Y6EQ2tpoLNVQNfi0HvWQpLa9btvrn3QPnRRyLHz0uY0eeW6bHoSF2jQddo1Deb9M1GQ4sZwQDHm8X2ZrI9mWwJkyVmsMQMpgejGxMWu1f+NEVbNIRRbjEpLSa52SgzG9rMhlYzoSe5viyOL4vnx7JGbudv9t3FSJj2lx368ea+oYEDF/V7zOeetRoAANc/rgiZ4CNJtMF566qBsPGsTNdktMkc6Dp1w8XG3MuNeTfbSsKEIWm+yck+CUle8Xwmr+cHh4Hwtm4HQiNhqlLVlCuqSuUVxbKyWnV9pCg82Ts+xTepv0+Sba9cXLh31FYnYvmeRqYHHjLxgaMY2gzybwt/PlVzbkrUxNlx08RsYc/f1FHMGsIgNRllZqPMbJJbjHKzSWkxqwmcgzEFOMP6HxfDeRjOwXAOhrFQjIVibAxBgbU1iWCIdXKqVqvl8/kAAMpCUxYKAEAaKZqmCT1JWWjSSBJGijCQFi1B6EmLjrSoLWYNYdGRDB7GkjBZIgZbwmB7MtmeTI43iyVi9KHl5e5iJIy/lx354davKT6Ji5LnhQofOA+17bqq8Zys37M2GJnsqoGwYHNl8HgfcZwtO7HMpLmoreRay40brTeLZWXeXK94z+gYSWSkODxKFN69uGjD8u8D60j1kMqkbtXLpPrWJm1LvaapQdNUq66XGxShwuAocXi0ODIrYly0OLz3xiiyJUy2mKGq0IuibXAR5HosWqItX5X+SkdrcHtxJKsGLns8cdb3Rb8sPLB8fPjoufHTfXnOtTtaJzEFOFOACyPvPhksGsKsJSwawqwhCB1p0ZOGVjNpJEkTRZoowkgCChAGEgBAkzRpogAANE0jCAIAQBkIykABABgbRRAE52IoA8FYGM5GcS7GEjP4gRwGH2N4MJgCnMHHXOmaTG3W/Fpy6NfSg/19k/815u2H3lH27OdRBa4XwwAAIABJREFUndOiqTHYdrMLl2FWWQxSk80rKybGTPVNTvVNBgBQNFWlqi2WlZXIyk/VnKtU1nBxTqgwOFDgHyjwD+D7+XC9fLheIrbIbkvY9O0WoUwjP1xzEkVRvcVAUqSRNBoJk95i0Ji1KpNGZVIrjUoOg+PN8fTlefvyfAIF/kGCgFBhkB/Px57LItT/0WZ0xd5Rm1yRNZxp0zeZov/W2X4YuUHxc/FvORXHB/qnzo2fbp/xac7JVVskndSobd5T/NuJqrMjggfPS5h573ayD1L/R5ux1RQ1p6ddfy5Z/jbsF+08qb6tVlVfr2ls0DQ16Vqkujapvk1j1ghZHkKWhwfLw4PJ5zI4LIzFY3DZOOuJ5L8B2CJsR9GU1qxDEITL4GAMzAvzZOMsHoPLZ/CEbA8Ry0PEFrVvL+lAXike+ZsqI2f4u9KVuK20XFZEze7Cr07CES9LfXJB0pyD5cdeP/u+F0cyMzZ7ZMgwZ/iiITugaDq36dq+0oO32sqyo8bvyv7Uk9O1UaC+A0R575eFT/N34eUPu60tXxU83t6z461NwAH+/e98kqRJpVGtNKnUJo3KpDYQRhNh0hMGNmb7adl9u+7gM3hLUhY4fK3Rh7KOiVBX6e/tEHNzmho9TYH2xQw7j8fgzo2fPjtu6p/1f/1amrM5b/ukiDHZUROCBK7W7IbayQyKw5Unc8qP8Zm8R2MeeStjLQvrzhR+hgAXRvHarqt8B7vs8ofdY9YQ+hbb94t2D4ZgnhxxV69yuqdvB8I+RJIgUNzSwEB4l5bLCt/B4m6P0UARNCN4SEbwkHpN48HyY88efyWQ7zcpYuzokOE2GZYGOQMzab7QcOVI5cnC1uJRIcPWj3g51rOn/eG+g8X1J1thILyL4qZGFMN3ntV07QYGQjsRxwnKfmoIy3Z0PpwJaabaCtRpL9+9qm83BAkClqU+uaT/gsuNeUcqT225+lWaX8q4sJFDAgY4cFFgqCcIirzakn+q+tz5+ssxkshJEWPWj1hjq29THMcv/7lB32Li+sLT478UxRpJQh+boWQTMBDaiSCEY9ESJoWFJe5jc7N6T1u+ShjBu3OXuB7CEGxY4KBhgYN0Fv3Z2gs5Fcc/vLR5oH/qyOChQwIH8Bhd7oCF7M9MmnObr5+ru/Rn/V/BHgGjQ0Ys6b/Q5v1jCIr4DBRL/1KETfGz7ZH7LpqklaW6yBnueHMBBkJ7QYA4jq+4pfEbJnF0VpxF61VVL+1RxWNwsyLHZUWOU5s0f9ZfPlF9duNfW+I8o4cEDhgSMKDzwwshu2nVt1mXc7rWciNGEpkRPORBM+JtxSdNWLStJizbr+/OnrQtdZWe7c1kCNwxKLjjZ3YUcZyg9ZoSBkIrwkBqawziRZ1a1bfbPFgCa0Q0Esa85oKLDVd+Kf4dAchA/9R0v5Q0v359bqkaV2IgjNdbCq82519pvi43KAb6p2aGjlgz9Lme7xTYGVx/NsZCNbVwQuFtimKNpIuLi7oMGAjtRxzPr/ilkSLobiym5XpkBWpRLA+7Z021XsLG2cODBg0PGgQAqFHVXWm6fqzq9EeXP/Xleff3TernnZjsHd9XljPt09QmTWHbrQLpzXxpUZWyJt4rJt03Ze2Q52IkUXabPd3OM0XYVqCCgdBKflPb+em8LgYGQvvBORjXn6Wu0IliXXAB7q5qK1D7DhQ55K1DhcGhwuBZcVMomiqVV+RLi45Xn950ZSsLZyV5xcZ7xcZ7xkSJw+EoG5sgKLJSWV0sK7vZVnJTVtqml8V7xfTzTliW+mS8ZzSzW/MfbMWrn8etr2rDYe8oACalxaIlBMFuek0AA6FdieMF8lsaGAgJA6mp0sctfOCCkPaBImicZ3ScZ/Tc+OkAgDp1wy1Z6c220pPVZ6tUtf5832hxeLQ4IlIcHiEK69MLnNqTzqKvVFZXKKrLFJVl8soadX0g3y/WMyrRO25W3NQIUag9F3XqGC+AjWCItt7Ad9cA0E5xSyOO47vtBQEMhHYliuaX72lwdC4cT1aoFkbznG1dj2CPwGCPwAnhmQAAgiKrVDVl8soyReWfDVcqFdUogkaIQ0M8gkI9gkM8AoM8Any43vbvzXM2bQZ5vaaxTt1Qq26oVtZWq2q1Fl2YMCRSFBYjiZwcOS5S5NRta68Uj7Z8NQyEyjKd294gBDAQ2hk/mG1SWCxagsF365KXFai9+zt1AwtHsWhxRLT4v3sUyAyKalVttaquWlV7rv5ivbpRaVIH8H0D+P4BAl8/nq8fz8eX5+3D8xKxnPqjdY/OopfqWpt00hadtFkrbdQ2N2pbGjSNbJwVJAgMFQYFewSm+fYLF4X48rx7by9Wm/NMERbvqg3L9nV0RhyKBqpyXbgbzyRx6+rY/hAU8Qjnqip0XikuWFd2EmWhVBW6mMeDHJ2RrrGu9pTul9L+jJEwNWqbrf81aVuutRS06Nqk+lYjYfLhenlyxN5cLwlbJOGIJWyRiC0UsYQitlDIErBx2y+W2ENm0qw2aZQmtcKoVJpUCqNKblDIDAqZQd5mkEt1rSiC+nC9fHk+vjxvf75vvFdMAN8vUODf12dn8gPZgAZuPrNe32zE2Kg7T3GGgdDehNE8VblbB0JVuY4fxMHZzr5C7EOxcVaEKDTinn1/TKRZqmuVGRWt+jaFUSXTy6tUtQqjUmFQqkxqtVlD0bQHky9g8vlMPo/B5TE4XAaXz+RxcDYTY/KZPCbKYOEsJsZkYUwMQTn/2SBTwPyfu8s6vU6D6O58Rm/RkzQFADASJoIiCIowEEYLZTESJoPFYCLNeotBZ9H/5/96jVmrNms1Zg1JUUKWQMjyELGFErZYzBZKOOIIUaiELfbievpwvWy+vbjzEMXxFbc07hwIleU6UbRbD1yAgdDeRFH84ot1js6FI8lvacVxrnw3goUxrbcbH5TARJo1Jo3GrNWYtTqLQW/R6wmDxqw1EkalUdWgaTKRZjNpNpEmM2khacpgMQAAaEBrzf8T9iiKQtH/uc/KZXAxBAUAsHAWA8UxFOPiHAbKYOMsDs5m4kwBi+/H9+HgHB6Dy2Ny+QyegMn3YAk4ztdItRtJvKDxrCxwdC9O3ndyqnKdl3PfquhtMBDaGy+ATegIs8rCFLppR4SiWBPfy/PonRwLY7K4nj2ftuiS++HZnzCaV/J9HWminG30lp3QQF2pi3S53VK7xC2/eMdCgEckT1Whe3hKV2SQmmiC5vm5b/sDcjYYExWEcJVlWkdnxDG0DQYGH2f+f3t3HiVFefcLvHqvpbtnYaane3o2lgEElUEQRbyiIIpHiRh34hIxEUVyXHOORl/0NclRg2+Omlw1EtBjchMVIyeC90bjUZQEEEFlicCwDLN2zz7d1VXVW3XdP9pMxll7pqu7tu/nr56emq6Hpqu/VfX8nucx5Mxq/RCECiicxoROGDQIe49Gis5waaeoEAyh6Axn7xGDBmHoOFdg7A5CAkGoiIJaZ99xgwZhT3rcLoCaFJ/h6j3KKt0KZYROcgXTjL54J4JQAbTHIcZSsb6E0g3JNzGeYhv5wukIQlAXyuMwmU18MKZ0Q/JOIsKn+YIp2h4Dkz0EoRJMhKuaYht5pduRb+GTnLOCMmhJAqhb0Qxn7zHDXRTywaiNsRh8fg8CQagU92Qm3GC4IAyd4AoNfxMG1Ck9wFfpVuRbuIF3T8YhiSBUiHsybcAg7DuBbnlQqYKpTPgUL6UkpRuSV+EG3lVj9PuiBIJQKc5KSuiIibGU0g3Jn2RUFDpimN0Y1MnmtNoLbFxLVOmG5FX4NO+ejCBEECrEbDUxPjLSLCjdkPwJn+Rd1TQWJQbVKqxl+ox0dzTBJpO8SHuMO7dcP9mCkOf5pqamVGrESxxRFBsbG6NRY51wjcJlsLujfccjhbXojQD1KpjGhE4YaDRhqIF3T6YxqJeQKwhffPFFv9+/ZMmSmTNnHjlyZOgG+/btmzJlyrJly8rLy9944w1Zdqp17ho63GCg08/QCQxXAlUrmMqET/OSaJRuQraBw33RNBmCsLGx8Wc/+9muXbtOnDixatWq++67b+g2a9aseeihh+rr6z/88MO1a9d2d3dnv1+tc09h2NOCQTrnk7wY60k4K9BBCOplpS3UJLtxOixQKdNPhiB88803L7nkkjPOOIMgiLVr13788cft7e0DNzhy5Mg333zzox/9iCCI+fPnz5kzZ+vWrdnvV+tsjMXmsvDthhjDGzrBuSbTJgvuwoCqFUxzGqSbMJVI8cGYC8VrBEHIsvpEQ0PD9OnT0489Ho/L5WpsbCwrKxu4gd/vp+lvTz1qa2tPnz490qulUqmDBw8Gg8H0j1OmTCkqKhppY0lMJlqaRLNWS36Y0mTvV422pFbbL/J8nM7ojLLn66SzxBRvPp7rJhlK5u8/ZIgpTAW/Estm9GWysabff7ZNIoulZPtJpRsybiazxeafIu9ryhCELMuWlPxnKS+GYUKh0KANKIoaZYOBYrHYT3/6U5vt2yWKnnrqqYsvvnikjfmOttRfXphwy5XHn9nTUGI5vEPpdkxQKpUSMjsL6eu+odT1WffxYK6bZCiZv/+QoZTkiHTc2vXn/20ixu6z0PT738edbRELuv+8U+mGjJuJcjI//C+CICKRjCqbaJq2WMZYBlyGIPR4PL29vf0/9vb2DrwcTG/Q19c3cINZs2aN9GoURX3wwQd+/4iLmg5ktVbZH/rNmP9I1XI1Cce3tPoeukHphkxQhuvhifHU6fVHqx75L4ydkBfWI8yF9mePF/zgfxj/2CuFafr9D/+x2TvT5Zl/s9INyYpc778MpzNz5szZu3dv+vHhw4fNZvPUqVMHbjBr1qyOjo7W1tb0j3v37p0zZ072+9UBxk9GO+NiXOfD6tlGnvGTSEHQBIPM+sQ2Cq4qdBB+S4YgvP7665uamp577rkDBw7cf//9d9xxB8MwBEE88sgjzzzzDEEQZWVl11133bp16w4ePPjEE09IknTFFVdkv18dMFlMtNfBtep8bCXbgNkrQDNcNXT4tM6DMMmLSV6kSjGU/lsyBCHDMB999NGuXbvWrFkzf/78X/3qV+nnKyoqfD5f+vHLL79cVVX14x//uL6+/oMPPrBajT7ZeT9nFc026fyoC59GlTZohruGZvUehGyT4KygMJS+nzyBdNZZZ7377ruDnly3bl3/Y7fb/cILWq5qyRlXFaXzFUElgm0Spq9CEII2UKUOMZ6K9SUchTal25IrkSbeifuiA2i15Ek3nJUU26TnAbxcMGpjrFjwDDTDAMuFss0CRhAOhCBUGO1xJDkxyYlKNyRX2NO8G/dFQVN0Xy8TaRZwRTgQglBpJsJZQbL6ndUp3MC7UCkDmuLWdb1MrDdBSISOb/xOAIJQec5KWsfTG4ZxRQha46ykhGAspdNxTSwuB4dAECrPWUXptXA0yYnJiEiXoUobtMRsM9NeR0Sn45oi6CAcAkGoPGcFqdehhGyz4KxElTZoj7OK0ut9mkj6qIQBEITKI4vsYjyVYJNKN0R+LKq0QZtclbRey7m5VoHx46j8DgShCpgIp5/U5X2YSBOmcQJN0muHRaw3QZhNdjeGM30HglAVGD/Fterw9DPSgpswoEl6HdcUaRGwPvZQCEJVcPrJSIverghRpQ0aZiIYPxlp0dvpKdcadVaMvbCG0SAIVYGpoCK6uyJkmwRXNQZOgFa5qnTYTRhpFZzoIBwCQagKtMeRYJPJqK7uw0SaUSkDGqbLbsJIS5TBFeEQCEJ1MBGMT2+DKFhUyoCWuaqoiL6uCBORZCqeIovsSjdEdRCEasHobDShRHAtUdyEAe1yFNoIExHrSyjdENlEWqLOChLjeodCEKqF00/pqWee74jZXBYrY1G6IQAT56zU1UUh1yowKBkdDoJQLZgKSk+Fo5FmVGmD5jkrKT1NiB9pjTrL0UE4DAShWjBeR7QnnkroZJ7fSAvOPUHznBW6GuCLK8KRIAjVwmQxUSV2PhhTuiHywHAl0AGnn9TNjKNiPBUPJalSVMoMA0GoIkw5ybXp4u6oRHCtqJQBzbMX2Exmkz7qZfi2KOV1mMwolRkGglBFdBOEQlfMSlusNCplQPOYCkof5dxcGzoIR4QgVBGmnOQCujjkWjFoF3TCqZeJ1ri2KO3DUTk8BKGK6OaKMIIRhKAXjF7GNXFtUQZXhCNAEKqIzWk1W/TQIRFpEVApA/qgk3WzJYILIghHhCBUF8avh4tCrPwJukEW62Hd7GhP3EpZrBS67YeHIFQXxqf5IMTKn6Arulg3m2uLMuggHBmCUF100E2IlT9BZ3SwbjY6CEeHIFQXWvtBiKH0oDM6WDcbQTg6BKG60B5HrC+h6YnWIuggBH3RwbrZXABBOBoEobqYLCaqVNsTrXFtUacfhxzoB1VqT4STYlSrp6diLJUIJ8kSTK42IgSh6mi6XiYpiEleJItxyIF+mMwmyuvgg1o9KvlAlCrD5GqjQRCqDuMjtXvIfVuchiMO9EXTp6e4LzomBKHq0D6SC2j11ij65EGXmHIyotkg5IMx2utQuhWqhiBUHdrr4DU74yjXFmXQQQi6o+lxTVwAgwjHgCBUHUehLSVKiYgmZ7LgWnFFCDrElJN8IEpISrdjQvhgDNNtjw5BqEa0l9Ri4aiUkoSOGO3FIQd6Y6UsVtoS7Y4r3ZBxi7NJQpLsLsz0NBoEoRoxXocW12MSOuP2AqvFgQ8V6JBGuwn5AFZfGhu+s9SI9mnyihCVMqBjjJ/itRiEwRg6CMeEIFQjxqfJehmuLcqUY04Z0CeN1stwgShKRseEIFQj2kdyGuyZxxUh6JhTm0GIW6OZQBCqkZWyWEiL5lbo5VoFBCHoFTnJnuCSSUFUuiHjIRF8e4wuwxXhGBCEKsX4NFYvk+TEVEJyFNqUbghAbpgI2kdqq88i2hO30liPd2wIQpXS3CHHpe/AYHI10C/GR3KaqmLDUPoMIQhVivFq8ZDDHRjQM8bn0FbhKIbSZwhBqFK01gpH+SCG0oPO0V5SWx0WfBAloxlBEKoU7XEInXEppZnKUa4NlTKgc0w5yQdiGirn5gMxBqenGUAQqpTZbra7rZqZ0ildnIZzT9A1K20x201aKeeWUpLQHac8WBx0bAhC9aK9Dq3MLxPtiVspFKeB/jHlmrk7Gu2KOwqsZhu+5MeG90i9NDT1NorTwCA0VM7NBdBtnykEoXoxXodWlqrnAyhOA0NgtLNuNiplMocgVC+tXRHikAP909AVIQq5M4cgVC+qzCF0a6NwFPMZgkHQ6aNS1MJRiSvCjCEI1ctsNTkKbEKn2gtHU0kp2pugPTjkQP/MVpOj0CZ0qP1WjSRK0Z4EVYqjMiMIQlXTROGo0B4ji+0mC2ZXA0PQRDeh0Bl3FNnMVhyVGUEQqhqthXoZLogOQjAQ2qeBo5IPRhncF80YglDVNFEvg5JRMBRGCxOtcZhldDwQhKqmiStC9MmDodA+DXRY4KgcFwShqtEeR7QnofISNYymB0MhJ9kTbFKMpZRuyGgwdmJcEISqZrKYyGJVl6iJsVSSF8lizGcIRmEymyiPg29X71GZSkqx3gRVgqMyUwhCtaPLVH3I8YEoVebAerxgKLRX1aukCR0xcpINhdyZQxCqHeUl1RyEXBDrvIDhqLyKjW+P0WU4KscBQah2dJmqe+bRJw8GxPgcfLt6rwixJtp4IQjVTuWFo+iTBwOivaoeU4/T0/FCEKqdygtHeUy3DcbjKLSl4ilRUGnhKB/ErdHxQRCqnclichTZhC41zjia5FOppGQvsCndEID8MhF0mSPaocajUhKlWG+CKkXJ6DggCDVAtXdHo+1xzF4BxkT7SKEjoXQrhsF3xMhilIyOD4JQA1RbLxPrTGA+QzAm2uuIdaoxCAVUyowfglAD6DKVjqAQOhKolAFjor1ktF2NQYj6tQlAEGqAqm+N4twTDInxOQRV9hHy7VG6DEfl+MgWhIIg7Nu3r6mpaaQNurq69u/f39LSItcejYPyOKLdaiwcjXYmEIRgTDan1WQyxdmk0g0ZjA/GEITjJU8Q7t+/f+rUqffff/+555776KOPDt3gpptumj59+tq1a+fOnXvNNdfE42o8k1KtbxfFVlnhaDyUMFtMNqdV6YYAKIP02NQ20Vp6YXoSC9OPkzxB+PDDD993333/+Mc/vvzyy5dffvmbb74ZtMHq1auDweDnn39+6tSpw4cPv/baa7Ls1zhUuFQ9F4iRZSjRBuNylNrUdlQKHTGyGAvTj5sMQdje3v7pp5/eeeedBEH4/f7ly5e//fbbg7a57LLL7HY7QRAul2v27Nnt7e3Z79dQVNhNyLdHHSW4HATjospsaluhl2/HfdGJkOGLrLm52eVylZSUpH+cPHnyKD2FJ06c+PTTT5988smRNhBFcefOnf2vdvbZZ3s8nuwbqXV0Gdl9OKx0K76DD8RIL64IwbhIjz10SGVHJUpGJyTTIHz66acPHTo06Mn58+c/+OCDPM+nr/bSSJLkOG7YF+np6bn22msffvjhurq6kXYUj8d/85vfOBzfntQ88sgjF1xwwUgbC4Jgt9stFkuG/woNc4uRNiESiSjdjv9gW7niWkpVTTIajuNMJtwEU0zKmeDahAgbUc8yZOHWSOFsxiBHZYaff5qmzeYx7n1mGoSLFi2aOnXqoCf9fj9BEGVlZX19falUKr2z7u5ur9c79BVCodDy5csvvfTSxx57bJQdURT19ttvp195TBaLxSBBSNdIx3uDDMWoZcIIiYh1JQur3U6nU+mmGJckSXj/FSRJkpUK25IOR5FaZhmMdwWKqwtopyEuCmX8/GcahBdddNFIv6qpqSksLNyzZ0/60u2f//znAw88MGgbjuO+973v1dXVPffccxNuq5H1F46qpAMg2hO30haLA+NQwdAYn4MPRlUShCgZnTAZvsgcDsfatWvXrVv3t7/97dFHH+3s7LzuuusIgti1a1dVVVV6m+uuu66pqWnevHkbN2589dVXP/vss+z3azSqKhxFVwQAQRB0mYrWY0LJ6ITJU/W3fv36SZMmvfTSS36//9NPPyVJkiCIsrKyG2+8Mb3B4sWL58yZ09DQkP6xvxYGMkeXpdcCdSvdEIIgCD6I1ZcACNrrCJ0cviQi/1AyOmEmSVLXfCUVFRWff/55hn2EBiqWIYjOL0Pdh8Mzb6tUuiEEQRDH/k9L0QwnNcPicrmUbotxsSyL919BLMsSvdaT77TVPTi4fkIRTX/rIAiiarlRyuxl/Pyjj0cz1HVrNIAlsAEI2uvgO2JSShWXE3w7jsoJQhBqBuVxRLvjaphxVEpJQmec8uCQA6Oz2M12pzXWo4plKDDL6IQhCDXDbFXLUvXRrri9wGqx48MDQNA+hxrml0mXjOL0dGLwXaYlKlmhlwvEGCxMD0AQBEHQXlINR6WAhemzgCDUkn8XjiqMD0ZpBCEAQRAEwXhVcUWIktFsIAi1RCXnnnwgyqBPHoAgCIKgfaQaJsTH0N5sIAi1RCVrUHA45AD+TSXrZqNkNBsIQi1RwyGXSkqx3gRVinUnAAiCIMxWE1ls4zsUvlXDBXBrdOIQhFrybeFop5KFo3wwSpXa0ScP0I/2kcouVf/t6SlKRicKQagx6Ul+FWwAuiIABlF8sguhI0ZOwunpxCEINYb2kpyihxxmGQUYhPGRyhaOYqanLCEINYb2OpS9CcMFcEUI8B2Kl3NzwRgKubOBINQY2qtwrTYfiNK4IgQYgCqxJyJJMZZSqgEY2pslBKHGUKX2WCiZSihzyCWjYlIQySKUjAIMYCKoUiW7CdFznyUEocaYzCaqxM63K3PI8YEY7XUQ6JIH+C5auSo2MZ6Ks0lykk2RvesDglB7FOyQ4INRzDIKMBSjXJ8FH4zRHofJjPPTiUMQag/jU6xeBpUyAMOifSTXptjpKUpGs4Qg1B7aSyp1a5RrizLlCEKAwZhykmsTFNk1OgizhyDUHlq52e4xXAlgWHa3lTCZ4mwy/7vG0N7sIQi1hyy2J3kxGRXzvN9Yb8JsM9mc1jzvF0ATGJ+Db1PgDBUdFtlDEGqQiaDKHELe62W4AO6LAoyI9pFc3oMwKYhiTHQUomQ0KwhCTVJkSieuDSWjACNS8qhExWh2EISapMghxwcwewXAiHB6ql0IQk1iyhWo1catUYBR0D6H0BnP83KhfCBK46jMGoJQk5hykm+LEnk84lJJKdqNBc8ARmS2mR2F+V4ulAvgilAGCEJNstIWs90U60vkbY98e4wqsZut6IsAGBGT53oZKT2IEKen2UIQahVTntdDjm9DByHAGBhfXsf4RrvjNsZqpSx526NeIQi1iinPa888F8CgXYAx0D4yn9Mfcm1YE00eCEKtyvOgJVTKAIxJgdNTHJVyQBBqFZPfc0/cGgUYE1lsTwpiUsjTrE88KmVkgiDUKrrMEe1J5GeF3ng4mUpJmL0CYAwmgvGRXGuezlAjOD2VCYJQq0wWE1WapxV6uVbB6afysCMArXNWUJHWfCxDIcZTiXCSKrXnYV+6hyDUsLzNZBFpjTr9OPEEGBtTnqcrQj4Qo7Aer0wQhBpG+8j8zHbPtaJPHiAjjD9PQYih9DJCEGqY009G8nLIRVoFpgK3RgHGxvhIoTueSuZ82ieuVWBwn0YmCEINYyoorjXnE62JUXRFAGTKZDFRJfY8VHRHWqMIQrkgCDXMxlgsDnO0N7dzG3JtAu0j0RUBkCGnn8r5rRqJ4AMIQtkgCLWN8ZNcS24POZx4AowL4ye5HBeO8h0xm8tqJTG5mjwQhNrm9Oe8VptrQ8kowDjkoV4GI5rkhSDUtnwcci0Cg0MOIGNOP8UFctt5H2mNOitweiobBKG2OSuoSEsOrwglUeI74wzWeQHImIU025xWoTOHk13g9FReCEJtcxTZUqIUZ5M5en2+PUYW2cx2fE4AxiHXQ5somrtYAAASlElEQVQibei5lxO+4DTPWU5xObsojDQLzkqceAKMj7Myh7dqoj1xs8Vkd1lz9PoGhCDUPGdFDs892SbBVYUgBBgfZyUVacpVEHKtUScmuJAVglDzGD+Vu1rtSLPgrKRz9OIAeuWsoiKtgpTKScFMpDWKmZ7khSDUPGcFGcnNUMJUUhI6Ykw5KmUAxsdKWuwuq9CZk8kuuBYBI5rkhSDUPKrUkeCSuVgLlGuLUqV2sw0fEoBxc1bSObo7GmkRcGtUXviO0z4T4azISYdEpIl3VuG+KMBEuKootpmX/WVjfQlCIhxFWCVbTghCPXBV02yj/IccKmUAJixH9TJso+CqxumpzBCEeuCqothcXBFi7ATARDn9JBeMyb4eU6SJd+L0VG4IQj1wVdOyB6EYS8X6EjTmlAGYELPdTJXY+aDMhWxsk+CqRhDKDEGoB3a31Ww1RXvkLFGLNAtMOVZfApg42e+OSikJlTK5gCDUCWcVxTbKecixuC8KkB1XFcU2y3lU8sGYo9BmpbD6kswQhDrhqqIjTXLWy0SaeFTKAGTDWSVzFRvbJKCDMBcQhDrhkvuKMNzAuyczMr4ggNEwPkesL5HkZBvjyzbyKBnNBQShTjirKC4QlUR5StSEzjhhMmGsEkA2TGaTq4oOn5btojCCEU25gSDUCYvdTBbbuYA8JWrhBq5gKk48AbLlnkKHG+QJQjGWinbHGR8mV5MfglA/XNUUK9O5J3uad9fgvihAttyTmfApTpaXYpsExk+aLCjklh+CUD/cU5nQSXmCMHSKd0/BFSFAttzVFBeIyjKsPnySK5iK09OcQBDqR8FUJnSSI7I+4hKcmAgnMZQeIHtmu5kuc8hS0R1CEOYMglA/HIU2i93Md8SyfJ3wKc41mcZQegBZuCcz2XcTppJSpEVw1eA+TU4gCHWlYCoTPplth0S4gXdPxvEGIA/3ZDp0KtsgZBt52uuwOPCNnRN4W3XFPY0OyRCEHIIQQC7uKTTbwGe5Wj3ui+YUglBXvu0mzIIYS/HBGMYqAcjF5rTa3FY+kFWfRfgkjyDMHQShrpDFdpPFJHROfPbt0AnOVU1jVXoAGRXNcPYei0z4zyVRYpt4F+7T5Ay+7/Qmy4vC3mNs0UynjO0BgKKZzr5j7IT/nG0SKI/DSmKu7VxBEOpNtkF4NFI4A0EIICf3VIZtFsRYamJ/HjrJFUzBfdEcQhDqTeEMZ+9RdmI989GeeCqeYryYwwlATha72VU58UK23m9YnJ7mlDxBmEqldu/e/f777/f19Y2yWSwW279/f3d3tyw7hWE5Cm2OAtvEVqLoPRIpmukiMIAQQG6FM529RyfSTZjkRD4YK5iGK8IckiEIRVG86qqr7rrrrldffXXmzJmHDh0aacv169cvWLBg+/bt2e8URlE8y9XzzUQ6JHqPooMQICeKZkywm7DnCFswnTFbcX6aQzIE4fbt20+cOLF3796//vWvd9111xNPPDHsZnv37t25c+fcuXOz3yOMrni2u+df4fH+lSRK4VM8TjwBcoHxkenlI8b7hz3/ChfPduWiSdBPhiDcunXrNddcQ1EUQRCrVq3atm1bMpkctE08Hr/77rtfeukliwWFTznnqqISnBjtGd8hF27gqVK7zWnNUasADM1EFM1wjffuqCRKffVc8UwEYW7J8K3X0tIyf/789OOqqqpkMhkMBisqKgZu84tf/GLFihV1dXVjvloymXzvvfeKi4vTP55//vmDXmogURRFUbbVn/WkaKaz61DId2Fx5n/SeaCvcJZzXO8n3n9l4f1X1njf/8JZTOAfPZ7zCzL/k9Bxjiy1m2kT/qOHyvD9N5vNJtMYN5YzCsIDBw48+OCDQ5/fvHlzdXV1PB632b5dyjz9IBb7zhwKBw8e3Lp16xdffJHJvkRR3LZtW/r6kiCIsrKy0tLSkTaOxWKSJOEqcyhnraPz83DxuZne55RSUteB8Iw15YP+70YXj8fHtT3IC++/ssb7/tNTbPyWKNvB2QsyvQLpOhRyz6DwvzysDN9/kiTlCcKampr169cPfT4dUV6vt6urK/1MZ2cnQRA+n2/gZr/85S/9fv+TTz5JEERzc/OWLVtcLtf3v//9YfflcDg2btzo9/szaZjJZLLb7QjCoRxnkU3vdtlNDiuV0ZvTe4SlShxF/nGcqxIEIYoiTWO2C8Xg/VfWBN7/SWcVcEfjhZe4M9paIsLHWs64s4qmMaJpGDJ+/jMKwoKCgsWLF4/020WLFr333nuPPfYYQRCffPLJ3LlzBzVu9erVjY2N6ccOh6O0tLSkpCSLNsPYLA5z4Uxn55ch36KM7o52fhnynDO+FASA8So9p/D0tqD/koy+AEMnOYvDzPiQgjknQx/hbbfd9uyzzz700EOzZs16/PHHn3/++fTzixYtuuaaax5++OHLL7+8f+NNmzZdfPHFF110Ufb7hdF5zy9qeC+YSRCm4qmeb9jJ3/PmoVUARlY4jYmzSb49RpeNvfB1cE9v2flFeWgVyFA1WlRUtGfPHrvdvm/fvs2bN994443p5++5556hgXfvvffOmzcv+53CmAprnWI8FWkee2R992HWVU3ZXKgXBcgxE1Fa5+76KjTmhklO7D3CeuYX5qFRIM93X1VV1dNPPz3oyVtuuWXolrfddpsse4SxmQjveUXBPb3TKsdYU6ljX2/pOTjeAPKhdF7h0debKy8rNZlHq+Do2N9XPNuVYR8/ZAlzjeqZZ0FR19eh0af6jbQIfHuspA4dhAD54KygHMW2zi/HuCgM7unx4r5oviAI9czushZMYzr2jTYBbNPfOiqWlmICJ4C8qbrc0/z3jlFmxg+d4KQU4Z6MaZ7yBEGoc1XLPU0fdiS44YedRloELhAtW4ATT4D8KZjK2AtGvCiUUtKpvwaql3sw/X3eIAh1jvGRJXMKGv9f+7C/xeUggCKqryhr+nD4i8LAP3psjBW9FfmEINS/6is8PYfDkabB5aPte3v5jhguBwHyzz2Zpkrsje8PPkONs8nmjzqnft837F9BjiAI9c9KWWqu8tb/uWXgzPc937CN/7d99o+rcTkIoIgZt1T2fMO2ftrV/0ySE+v/2Fx2XhHlGXuUIcgIQ8cMwTO/MCmIB144NWWl11FkjzQLzR91zvpRNVWK4w1AGVbaMntNzcHfnEoKqYKpNCERx99qLZlTUL3co3TTDAdBaBTl/2uSezJ9/M1Ws9XE+KkzVle5qsYYXwgAOeUotJ25pibwz57mDztjocTUa8uLZ2HFJQUgCA3EWUHNfXia0q0AgP+gPI4p16BHUGHoIwQAAENDEAIAgKEhCAEAwNC0HYRdXV0syyrdCuMKBALRaFTpVhhXc3OzKA4/ZxDkWiqVampqUroVxhWLxdra2uR6NW0H4WOPPfbuu+8q3Qrjuv322/fv3690K4xr8eLFPT09SrfCoEKh0IUXXqh0K4zr66+/HnaBo4nRdhACAABkCUEIAACGhiAEAABDM0nSiGtiKeLqq68+ePCg2ZxRQrMsa7PZSJLMdatgWH19fQzD2Gw2pRtiUN3d3UVFRRkeLCAvSZK6u7tLSkqUbohBJZNJlmWLisZeM2D79u1nnHHG6NuoLghZlu3s7FS6FQAAoAcVFRV2u330bVQXhAAAAPmEmyoAAGBoCEIAADA0BCEAABgaghAAAAxNM+sRdnZ2btq0qbOz88orr1yyZMnQDURRfOONNw4dOjRz5sw77rgDNf0y6uzs3LZt25EjRwoLC2+44Yba2tpBG6RSqd///vf9P5555pkXXHBBftuoZzt27Kivr08/tlqtq1evHrpNW1vb5s2b+/r6Vq5ciam/5LVx48aBRYVDP9719fU7duzo/3HlypUeD1aZz0pHR8f+/fubm5svvvji6dOn9z/f0tLy2muvhcPha6+99vzzzx/6h/F4fNOmTcePH6+rq7vlllsyHFykjStCnucXLlx47Nix6urqVatWvfXWW0O3uffee1966aXa2to//vGPP/zhD/PeRj1bt27dBx984PP5Ojo66urqdu3aNWiDZDK5Zs2ao0ePnjp16tSpU93d3Yq0U69ef/31t956K/3eNjQ0DN2gr69vwYIFra2tFRUVV1999fvvv5//RupYQ0PDqX974IEHDh06NGiD3bt3b9iwoX+bWCymSDv1ZOnSpU899dSjjz66e/fu/ie7urrOPffczs5On893xRVX/P3vfx/6hzfffPNbb71VW1v7wgsvPPjgg5nuT9KCzZs3n3vuualUSpKkP/3pT2efffagDQKBgMPhaG5uliSpp6eHJMmTJ08q0FCdEgSh//GaNWtWr149aIP0kc+ybH7bZRS33377r3/961E2eP7555csWZJ+/MorryxatCgv7TKcffv2URTV29s76PnXX399xYoVijRJr0RRlCTpvPPOe/311/uffPbZZ6+44or04xdeeKH/M9/vyJEjNE2HQiFJkhoaGiiK6urqymR32rgi/Oyzz5YtW2YymQiCWLZs2cGDB3t7ewdusHv37mnTplVUVBAEUVRUdM455+zcuVOZturRwLl7otGo0+kcdrPf/e53L7744pdffpmvdhnInj17NmzYsGXLlkQiMfS3n3322WWXXZZ+vGzZst27dw+7GWRp06ZN119/fWFh4dBftbS0bNiwYfPmzZgPRBbD3tIc9DnfuXNnKpUatMGCBQvcbjdBEDU1NZWVlZ9//nlGu8u6wfkQCARKS0vTjydNmmS1WgOBwEgbEARRVlYm41JV0G/37t1bt279yU9+MvRXS5cu7e7uPnLkyJIlSzZs2JD/tulYVVVVaWlpT0/P008/vWDBAp7nB20w8PPv8XhSqVQwGMx7M3VOEIQ333xz2A7agoKCM888MxQKbd26debMmUPvnYIsBn3OE4lEV1fXwA2CweDAIPB4PBkGgTaKZaxWazKZTD8WRVEUxUFT5lit1oErlCYSiTHn1IHxOnbs2PXXX79x48Zp06YN+pXdbv/oo4/Sj2+++eZLL7107dq1DMPkvY369NRTT/U/OOecczZv3rxu3bqBGww8QNIP8PmX3V/+8peioqKLLrpo6K9Wrly5cuXK9ON77rnnv//7v9955538ts4QxvycTzgItHFF6Pf7+4M9/cDn8w3aoLW1tf/H1tbW8vLyfLZQ9+rr65cuXfrMM8/ccMMNo2+5cOHCZDLZ0tKSn4YZis1mW7BgwdB6mYEHSGtrq81mw2TQstu8efOdd96Z7qAZxQUXXHDq1Kn8NMloBn3OaZoedJt6wkGgjSBcsWLFtm3botEoQRDvvPPOkiVL0lcbX3/9dWNjI0EQixcv7urqSq+Wfvz48aNHj/bfSobsNTY2Xn755evXrx+0JPQXX3yR/tgJgtD/5Pbt22marqmpyXMjdaz/7WVZdseOHbNnzyYIIh6Pf/zxx+kypRUrVmzdujXdL7hly5Yrr7zSYrEo2GD9aWho2Llz56233tr/DM/zH3/8cfq6pP8/SJKk999//8wzz1SmlXq3YsWKd999N33Nt2XLlhUrVqSf/+KLL9IBefnllx84cCB9prhnzx6e5xctWpTRS8tS4ZNryWTysssumzdv3q233jpp0qRdu3aln1+yZMnPf/7z9OPnn3/e5/OtXr26srKy/0mQxfLlyxmGmfdvd999d/r5uXPn/va3v5UkaePGjbNnz/7BD36wfPlyt9v9hz/8QdH26k1xcfGKFStWrVrl8/muuuqqRCIhSVL6yD99+rQkSdFo9MILL1y4cOGqVatKSkq++uorpZusN4899tiVV1458JmjR48SBNHd3S1J0vLly5cuXXrLLbecddZZtbW1jY2NCjVTPx544IF58+YxDFNTUzNv3rz0dz7P8+edd96FF1540003eTyew4cPpzc+++yzX3nllfTjxx9/vLq6evXq1V6v9+WXX85wd5pZfUIUxR07dnR1dS1evNjr9aafPHbsmMvl6r/4PXz4cHpA/dy5c5VrqQ7V19ezLNv/o8vlSg9x/de//lVaWprutd6/f39DQ0NBQcH8+fMxmlhep0+f/vrrr+Px+PTp0+vq6tJPJhKJr776qq6uLt0LkkgkPvnkk76+vksuuWRgvQDIor6+3u1293/zEAQRjUYPHjw4b948i8XS09Ozd+/e3t5ev9+/cOFCzOaRvZMnT/b19fX/WFtbm64FTd8ICYfDS5cunTRpUvq3hw8fLisr6//Y79+/v76+fs6cObNmzcpwd5oJQgAAgFzQRh8hAABAjiAIAQDA0BCEAABgaAhCAAAwNAQhAAAYGoIQAAAMDUEIAACGhiAEAABDQxACAIChIQgBAMDQEIQAAGBo/x/jAkkkX8+z/wAAAABJRU5ErkJggg==", "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": 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 }