{
"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": [],
"cell_type": "code",
"source": [
"x1 = 0.2\n",
"x2 = 0.8\n",
"nucleus = ElementGaussian(1.0, 0.5)\n",
"atoms = [nucleus => [[x1, 0, 0], [x2, 0, 0]]];"
],
"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",
" PowerNonlinearity(C, α),\n",
"]\n",
"model = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,\n",
" spin_polarization=:spinless); # use \"spinless electrons\""
],
"metadata": {},
"execution_count": 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 `nucleus` we have to specify\n",
"a starting density and we choose to start from a zero density."
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy Eₙ-Eₙ₋₁ ρout-ρin α Diag\n",
"--- --------------- --------- -------- ---- ----\n",
" 1 -0.143587409681 NaN 3.81e-01 0.80 7.0\n",
" 2 -0.156036336865 -1.24e-02 7.92e-02 0.80 1.0\n",
" 3 -0.156769577585 -7.33e-04 2.73e-02 0.80 2.0\n",
" 4 -0.157044839994 -2.75e-04 4.93e-03 0.80 2.0\n",
" 5 -0.157048932300 -4.09e-06 3.02e-03 0.80 1.0\n",
" 6 -0.157056386160 -7.45e-06 1.82e-04 0.80 1.0\n",
" 7 -0.157056406546 -2.04e-08 3.35e-05 0.80 2.0\n",
" 8 -0.157056406852 -3.06e-10 1.27e-05 0.80 1.0\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380286 \n AtomicLocal -0.3163446\n PowerNonlinearity 0.1212596 \n\n total -0.157056406852"
},
"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×1 Matrix{StaticArrays.SVector{3, Float64}}:\n [-0.055663898206394974, 0.0, 0.0]\n [0.05566243989578465, 0.0, 0.0]"
},
"metadata": {},
"execution_count": 8
}
],
"cell_type": "code",
"source": [
"hcat(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+gvaeTAAAgAElEQVR4nOzdd3xTVd8A8HPvzR5N0nQPundLaUvZq+xC2UsQUBAeeRQHLhBUUPBxojyiOAABB6iIoFD2LMjsoKWlm+7d7J3c8f4R3srDKB1pkibn++GPNJzcnJzcnN89556BUBQFIAiCIMhZobbOAARBEATZEgyEEARBkFODgRCCIAhyajAQQhAEQU4NBkIIgiDIqcFACEEQBDk1GAghCIIgpwYDIQRBEOTUYCCEIAiCnBoMhBAEQZBTs7tAmJeXZzKZOpiYJEm4RJwNEQRh6yw4NVj+tgXL37YsWP52FwgnTZrU3NzcwcQGg4EkyR7ND9QOrVZr6yw4NVj+tgXL37YsWP52FwghCIIgyJpgIIQgCIKcGgyEEARBkFODgRCCIAhyajAQQhAEQU4NBkIIgiDIqcFACEEQBDk1GAghCIIgp0azdQagnmUiQbGCKlFQFSqAk8BIAgEDeLBAmADp64rQ4YUQBNkOQYEiOVWsoOq1QGYADBQgCOjDBeECJFKIcGD1bC2wpB1Tow78Uk6erCP/bqS8OUikEAnmAwYGGCgoV4LLTaBITpYpqX5iZHogOjsQCeQjts4yBDmLFj04WEkeqCCvNlPeHCRGhPhwgIgJNDggSJDZAkoUZLmKShAj43zRJ4KRMAH8efYsGAgdzak66osC4u8makYAuiwC/XkUKmI+PKUWB5ebqAOV5MC/iL6uyKpYLNUfgT84CLrPJ598cuPGjQefJwgCw7BOHUphBNVqSmIAbizgwQJjmAgNBQCAJgCa7kkWDkAIBeQGcMBAfaUDHBrw4wJPtrP/PEeNGvXcc8/1xJERe1u02s/P79q1a76+vh1JrNPpGAxGZ89FR3WshtqQTahN4PW+6JxglNvhixwTCX67Q36WT5pI8PEAbKJfR39uKpWKz+d3MbtQt8Hyt47hw4ePHz8+MjLS1hlxatnZ2fn5+YcPH257xoLnP2wROoIiObXqKlGpBu/3R6cHoGgnrxvpKHgyFH0yFP2zilx1lficC74aioW6OPnVJwT9IyUlZdiwYbbOhVNjs9n5+fk9dHA4WKJ3M5HgvRxyxBF8gh+aN5M2M7DTUfBe0wLQWzNpqf7o4L/wzbdIwr46CyAIgnoEDIS9WL6MGvAnfq2ZzJlBezkWtcgQUBoKXo5Fr02jpVeTo9PxOg0MhhAEOTjLBMKffvopLCzMw8PjmWee0el0j0p29OjR/v37HzlyxCJv6uR2FJOj0/EXY9D0CTRfroW7MYP5yOlJtPF+aP9D+IlaGAshCHJkFgiEhYWFzz///O7du4uLi6uqqjZu3PjQZAqFYvXq1TKZTCKRdP9NnZkWBwvPE1sLyIw02pLwnmrTowhY1w/9dQztmYvE5ltw92MIghyWBarR77//fvr06UOHDhWJRG+99dbOnTsfOhL11VdfXbVqlZubW/ff0ZlVq6nhR3AMAVen0iKFPT6eZYQXcmUq9nMZuTSDMMJoCEGQI7JAICwqKoqPjzc/jo+Pb25ulslk96U5e/bsnTt3lixZ0v23c2ZXmqlBf+FPhqB7RmJsaw349eciF6fQZAaQdgJXmaz0phAEQVZjgdpUKpW2TeZwcXEBALS2trq6urYl0Gg0L7744qFDh5AOzAfVaDR+fn5tf/7444/Tp09/VGKnmkd4qAZ9JYv2zQB8vI9Brbb2u+8eBF7JpI08TB4YYXJn3W3xazSajnynUA+B5W8dJOkgnSHNzc2JiYm1tbXmP7dv397Q0PDOO+88Kv3x48f/+OOP7777zloZfAyCINT31H0dPP85HA6KPqbJZ4FA6OrqqlQqzY/ND8Ri8b0JNmzYkJycrFAosrKyNBpNZWVlVVVVQEDAQ4/G5XLz8/M7OKEewzAnCYSf55Of3yJPpmL9xAxb5WFnCtiQTaSeQ89MwszDcyiK4vF4tsoPBMvfOh5bjfYWJEm2traaH+t0uvfeey8rK6ud9BMmTHj99dcLCgpiYmKsksHHwDDs3hPegue/BQJhWFhYQUGB+XF+fr5YLL63OQgAIEny1q1bzz77LACgsrJy9+7dGIa99dZb3X9rZ0ABsPYGcbiaujwV87P06NDO2pCI8ejkyHTidCoGlyeFIJuoqqrKzc0NCQnZv3//iBEjRo8eXVRUdOTIEZPJNHXqVHPQ0mg06enpBQUFbDZ7ypQpD0ayAwcOJCUleXh4AADOnj3r4eERGxu7Z8+emTNnajSajIyMuXPnIgjy5JNPfv31119++aUNPqcVWeBKZ8mSJX/88UdOTo5Wq/3ggw+WLFlibq5u3Ljx2LFjAIDNmzdn/r+YmJgNGzbAKNhBBAWWXyQuNFAZaTSbR0Gz1+LQVbHoqHSiQgWnVUCQDeTl5T333HMrV64UiUQoiqanp6empiIIwuFwUlNTMzIyAADZ2dkZGRl9+vShKGrMmDGZmZn3HSQ9PT0lJcX8eOfOnWfPngUArFq1SiKRlJaWbtiwwfxfKSkp6enp1vtsNmKBFmFcXNwnn3wydepUlUo1adKk9evXm58vKSnx9/e/L3FkZOR9HafQoxhJsPAcoTBSpybROr5wqBU8H42iCBhzlDgyComGS11Czie7lVp9g7DOew3yQDYm3X/3R6lUHjp0SCAQAADCw8N37tw5evRoAICnp+cHH3wwYsSI4cOHDx8+3JyYyWR+//33/fv3v/cIt27dWrx48WPfPTw8vLKyUqlUmsd/OCrL1K/Lly9fvnz5fU/++OOPD6bcs2ePRd7R4elwMPsMzkCRv8bTmPZ3D/TfUShOgrRzjIwplJ00VSHIagL5yOq+VvpZurEe8mRERIQ5CiqVytLS0o0bN37wwQcAAJVK1dzcDABoaGhYsWJFYWEhnU7X6XQPdo2qVCoOh/PYdzffhFOpVDAQQtamwcG0k7gXB9k9AqPZ6336F2JQlY4Ydwy9MJnmwbZ1biDIilyZYKyvLa//mMy7m6uZRwt+++237u7u5mfMQ3veeuutyMhI81j9r7766ujRo/cdwd3dvW2eG4vF0uv1bf+l1+tZrLvhVyKRoCjq8N149lrLOjENDqaexL04yJ6R9hsFzV6IwOcFI+OO4VKDrbMCQU6JxWKlpKTs3btX9P+MRiMAoKmpKTg4GEEQo9G4b9++B1+YnJycl5dnfhwfH3/y5Ekcx81/Hj58uF+/fubHubm58fHxbXHRUdl3Ret8lCYw/hgeLkB+GIlhvaHHcUMiNsYHmXoS1+K2zgoEOaXt27cfO3YsMTFxxowZsbGxH374IQBgxYoVa9eunTFjRv/+/R86G23OnDnmwYwAgGeffVYgEISHh6vV6rFjx2ZlZbWtlHn8+PHZs2db7bPYCuwatSMKI5h4HI93RbYN7RVB8K7Ng7DlF4lpp/D0CTQGvLKCoB42ceLEUaNGtf0ZGBh49erVysrKpqamoKAgT09PAEBaWlphYWF5eXl4eDifzzc3Ez08PGpqasyvSklJMRgMeXl5ffv2ZTKZBw4caG5uDgsL++mnnwYNGmROo9VqDx8+fPnyZWt/QquD9Za9kBvBuGP4QA/k62G9KQoCABAAvhmKsTHkmQy4gyEE9Tg6nX7fzuwIggQFBQ0aNMgcBc08PDwGDx4sFosZDIZ5zAuKom23EhEE2bp167Vr1+5Nj2GYl5dX2zM3btxYu3atea6hY4OB0C7IDGDcUXy4F7JlUC+LgmY0FPw6GqtUU69fs9KYcgiCumno0KH3jfZ/6aWXzINRzUaOHLl06VKr58sGYCC0PakBjD2Gj/JGNg+0v3kSHcamgb/G0Y7XUnDPJgjqpdavXy8SiWydCxuAgdDGJAYw9ig+1gf5pDdHQTMRExyfiH1RQO4th7EQchz2NhDs2LFj48aN66GDHzly5LnnnmsnQVFR0fjx4x+6117vBQOhLckMYOIxfKwv8tGAXh8Fzfy4yLGJ2CtXiZN1DvU7gZxWmZKqVtvXyRwbG/vyyy/3xJEpilq9evXKlSvbSRMZGclisQ4fPtwTGbAVGAhtplkHRqXjE/yQjx0lCppFC5HfxtAWncdzpfZVfUBQZzVowfhjhBfbvm7c4ziu1WoBAEaj8aOPPqqrq9uwYcOaNWuqqqq0Wu3nn3/++uuvZ2dnmxO3trZu37591apVGzduLCkpaTtIdXX1e++9t27durKysi1btigUCgDAhQsX2Gx2dHQ0ACAjI8O8AOnu3burq6tVKtXmzZvNr120aNHXX39t5U/do2AgtI1GHUhJx2cEIpv6O1QUNBvhhXw1BEs7QdjbpTQEdZzKBCafwJ+JQIVMW2flf+Xn53/xxRcAAIPBsGbNmmeeecbPz0+hUIwbN868+bmrq2tKSkpLSwsA4MqVK7W1tQMGDGAwGMOGDSsrKwMAtLS0DBw4UK/XR0ZGrlix4s0335RKpQCAo0ePtk3MOHHihHk9mi+++OLOnTtyufztt982/9eoUaPOnTun0+ls8OF7BpxHaAN1GmrMUWJxGLq2n8NeiMwOQus0IPU4cWkKTWRn9QgEPZaJBHPO4APckXX90OMP/K+xqkj++1fWyQkjOFY449l2EmzatKl///5Lly51c3MbO3aseSDoyZMnz58/P2fOnClTpkyZMoWiKLlcXltbu2/fvrfffnvHjh0jR478z3/+AwBITEyMjY01Hyo/P3/q1KmPzZK7uzubzS4rK4uLi7PER7Q9GAitrVJFjT1GrIhCX4tz2Cho9lIsWqOhpp/CT6TSWA7Y7oUcFgXA8osEA0W+GvrwE5fm2Uc490XrZAZlc9tPYA5jKIq6ubm1hTQPDw+JRAIAKCgoWLZsmVQq5fP5jY2NaWlpAICSkpKEhARzyujo6LaVS7VaLZvdoYWDORyOSqXq0geyRzAQWlWpghp3jHi9L/p8tINHQbOPB2BPnicWnSd+HY2h9nWfBYIe6e1MokhBnZ1Ee9SsXpTFYfiHWTdTj4Rh/0Rr84rbAAAEQcwDO1966aWnn37avC/6m2++aY6OQqGwbcVtjUZjXncGAODu7m7uIwUAsNlsuVzedmSdTtcWI0mSlEql907e7+2cojq2E7ekVMpR4p1EZ4mCAAAUAbtHYBI9teoqnGgP9Q7fFJL7K6gj42kch2gmtEUsuVz+22+/mZ+cPHny3r176+vrAQBtQ2AAAAMHDszNzTU/TkxMPH36dFuz7/fff09KSjI/Likp4fP5wcHBVvsUPc1ZamSbu9ZMjTuGfzYQXRruXGXOxMDBcbRzDdSncKI9ZPf+rCI35pDHJmIP3QWwN3r11VeXLVuWmpo6fPjwxMRE85Njx4599tlnExISAgMDSZJks9nmvQnnzp3btg3FpEmTpk6dGhERUVxcvHDhwkOHDm3dutX88vT09Dlz5iCIA3XyUHbG19e3tra2g4m1Wi2O4z2aH4s4WUt6/GQ8Wk3aOiMWplQqO5iyVk0G7DP9VEr0aH6cTcfLH+qIvxtJj5+MWS33/06HDRt28eJFm2TpoXAc1+l0FEWZuyjbnpfL5SaTyfxYrVbr9XrzY4lEkpubq9PpdDqdRqNpS0+SJEmSBQUFAoGAJO9+6jlz5hw4cKAtjVarjYqKOnjwYNszBEHExcUVFRX12Od7uMOHD6elpd37jAXPf4do/Nu3Q1Xkvy8RB8bQhnk50AVUJ/lykVOp2Kh0wo2FTPBz3nKA7FaZkppzhvhhJC3Rzd7PTwzDzPcFEQS5d0W0e5cJ5XL/GWLj6urq6up630GWL18eGxtrMpm++eabdevWtTXvPvzww3sny7PZbBaLJRQK254pKytbsWJFRESERT+TjcFA2LP+m09uvkWemkSLFdn7r6unhQmQ38ZgM0/j6RNo/e2+roGcSq2GGnuU+GgA6jxXaU8++WRWVhaCID/88MOQIUPang8ODn7ppZfuTfn000/36dOn7c/w8PDw8HDrZdQqYCDsKRQA72YT++9QF6dgATxn+XW1b6gnsnM4bcoJ/HwaLUIAywSyCwojSDtBvBCDLgx1ovv3o0aNundTw3a8+KKVJorYkBN98dakw8G8M8S5eurSFBqMgvdK64O8n4ylHifqtXDRGcj2tDiYfAIf54u86ujzeqF2wO/e8hp1YFQ6zsTAyVS4qMpDLA1Hn41EJxwjpAZbZwVybublY0JdkI97/94vUHfArlELy5FQM04RSyPQtxPgDPJHWh2PSgxU2gn81CQaF56DkC2QFHj6AkFHkR3De99u2CRJfvfddzk5Of7+/itXrrx3MAvUBbBFaEm/3SEnHsc/HYi+A6Pg43w0AIsSIjNP4QY41R6yhZWXiXot9ctojNYLa8Fnn332yJEjqampVVVVKSkpBAF/Rd0Cr8YtAyfBukxifwV1KpXW1xUGwcdDAPhuODb/LDH/HPFb76yMoN5rzQ0is5U6M6mLq+Debi3ZcuMbS2fq4fp6RK9MWnbvM7W1tfv376+pqeHz+dOmTYuOjj527Jh5EVGoa2AgtIAmHZh/Fmdi4MZ0mhjeFOwwDAE/pWDTT+FLM4jdI+FipJCV/OcmmV5NnU+j8eldPEKAwO/Vge1t425BPPr9i24XFxeHhITw+XwAAIIgiYmJxcXFMBB2BwyE3XWmnnrqArE8An07AYVVeWcxUHBgDG3SCXzF38S3w3rfrRqo1/k8n9xTSp6f3K1rVi6dE+EaarlMdY5Wq1UqlW1/KpVKFxcXW2XGMcAOqa4zkWBdJvHUBeKHkdj6RBgFu4hNA4fH027LqJeuwPscUM/adpv8soA8Mwnz5tg6K91TVlZ2/vx5AEBFRcWFCxdSUlJsnaPeDQbCLiqSU0P+wnMlVM4M2mgfGAO7hUcH6RNo15qpV+AmFVCP+a6I/DiPPDMJ8+P2+h9sv379Nm7cOGTIkOTk5Pfffz801GbNU8cAu0Y7jaDAFwXkBzeJTf2xf0XCKwnLEDDAiVTa+GP4q9eIzXBSF2Rp24vI/9wkz07GAvm9PgoCALhc7pkzZyorK93d3e9dVhTqGliPd85tOTXsMP5XFXllKg1GQcsSMsDJVFpGA/XKVQKuOgNZ0HdF5Kab5JlJWLBDRME2gYGBMApaBKzKO0qLg7U3iFFH8KfD0bOTaSEuDvWLshNCBjg1iXa5mXrhMoyFkGVsLSA/yCXPTcYc5jcbFRW1YsUKW+fCocBA2CG/V5AxB/AqNcibRX82Eg6L6UHmduFNCfXsJYKEwRDqnk9vkf8tIM9Pdqi2YGho6MKFC7v8coqi9Hq9RXJi3vLQIoeyLRgIH+NGCzXyCP7+TXLXCOznFMyLbesMOQEXOjg+kVamoBadJ0xwW3uoq9ZnETuLyQuTHXD7l/z8/Nzc3PbT3Lhxo6Sk5MHnb968GRkZaZFsREdH5+TkdPnlOI7v37+fJG3/I4eB8JEK5dTsM8TM08TiMDRrOm2Ut6P9luwZjw6OTqSpTGDWaUIPR5JCnUQB8NIV4nA1lZFG8+39Y0Qf9Msvv/z444/tp9m2bdu9W+zaIZ1ON3fuXHtYHw6OGn2IW1Lq/ZvkuQby9b7YjyMxNiwkW2Bh4MBYbEkGMeEY/ud4mpBh6wxBvYSJBEsziEo1dW4yTeCIp01hYeHZs2dxHF+zZk1ISMjy5ctxHP/2228zMzN9fX1Xrlzp5eV17dq1rKys6urqlpaWhISEefPmPfRQRqPx66+/zsnJCQgIWLlypbu7u/n58+fPHzp0SCqVJiYmvvTSSwRB7N69+/r16ziODx8+/KmnnkLRRzaiNm/enJKSsn///ubm5tmzZ6emppqfLykp2bFjh0QiGTVq1MKFCxEE2bZtGwBg3bp1KIo+//zz/v7+li6qjoItwv9xuo6adAKfeJxIdkfK59Jfi0NhFLQhOgp+GIkluCEjj+Bw/0KoIzQ4mHoSV5rAyYmOGQUBAAKBwNvb29PTMykpybxZ/MKFCw8fPjxnzhyCIJKSkuRyuZubm1gs9vPzS0pKCgoKetShZs2adfbs2blz56pUquTkZI1GAwD4/vvvFy9enJCQsGDBgtbWVpIk9Xp9SUnJ1KlTZ86cuXPnzk2bNrWTve3bty9YsCAmJmbChAnLli3766+/AADl5eWDBw/28PCYNm3ali1b3nzzTQCAuYc2MTExKSnJtsNfYTUPAAByI/ixlPy2iEQAWBWHHhyLMuFMNvuAImDLIOzDXHLYYSJ9AhYldMBuLshSmnVgykk8zhX5ZmjPLuOurNCW7K3twTe4hzCMGzrX995nfHx8oqKi9Hr9nDlzAADl5eV//fVXXV2dSCSaNGlSZmbm7t27X3755cDAwNjYWHOah7p161ZGRkZdXR2Px5s0adLVq1f37t27fPnyt956a8+ePePGjQMATJw4EQDA4/E+/vhjnU7X1NS0evXqNWvWvPPOO+3kecWKFebhPBqN5tNPP506derWrVtnzZr12muvAQBCQ0MTEhLWr18/evRoAMCsWbPo9K6u+mohTh0IcRKcrKN+KiOP1ZAT/dEvh2AjvRFY0dqhNfGoHxekpOO/jqaNhDdroYcpVlCTTxCLw6yxFSjPjxW7IrCH3+QulPGYkF5WVhYYGCgSicx/JiYmPnSMzINKS0vDw8N5PN69L5TL5Q0NDYMGDbo3pUajmT9/flFRUUBAgE6na2hoaP/I8fHx5gf9+vUzZ6a0tLRtWfDo6GgMw6qqqnx9fR95COtyxkCowcGZOvJgFXW4iowQIgtD0S+H0F3hrhH2bWEo6stF5p3FPxmILQqFXfrQ/zjfQM0/i384AHsqzBrnBkpHWWIbd7y2zVsQCoVyubzteZlM1narr31CoVAmk937Qh8fHx6Px2AwJBKJeXcLs127dmEYZg5pN27cGD9+fPtHbsuPXC43R2iRSNT2pFar1ev1bZHbHjhLhWIgwKVG6v2b5LhjuPfPpq23yQQxcnMm7e8ptH9HoTAK9gop3si5ybR3s8m1N+AUQ+gfO4vJ+WfxvaNp1omC9sDd3b26utr8uG/fviiK7t+/HwBQU1Pzxx9/mMen3JvmoZKTk5VKZXp6OgDgzp076enpEydOpNFoU6dO3bRpk3kwZ0tLC0VROp0OQRAAAEEQW7ZseWz2vvvuOxzHSZL85ptvzJ2rqampP/zwgzkWbt26NTEx0dPTk8fjsdns9jNpHY7ZIiQoUKehypTgtpzKl1JZrdRtORUjQkZ4IS/FYCnjEK5jfm7HFyVErk2jzTyFzzpN/DAK6/J+cpBjwEnw+nXiaA2VkUYLEzhRn/m8efMOHDgQGBg4bNiwn376ad++fU8//fS7777b3Ny8Zs2aYcOGAQCWLl26ePHiwMDAmTNnfvbZZ22vRVGUxWIBAPh8/t69e5955hk+n9/S0rJp06bExEQAwFdffbVkyRJzd6tery8sLHzqqaf27NkTFRVFkuScOXOEQqH5UCwW66HDR8PCwmJjY41GY2Bg4DfffAMAmD9/fmZmZkREhEgkYjKZe/fuBQAgCLJx48axY8eiKPrHH3+0dahaH2Jv6wL4+fldu3atg33H9Qr9f4tpCIJoTEBpAkojaNZTDVpQr6XcWUioC4gSIrEiJMENSRAjXduKGmqHSqW6t//EaowkeOkKkdFAHRqHOVX1dx9blb+dkBjAvDM4HQX7Rvfs7Jrhw4d/8MEH5uhiz1paWlxdXTGs0zVdc3OzWCy+74XmXQ89PT2R/x84UV9f7+rqag6i7YiMjNy1a1dCQoJarXZzc7v3v0wmk1KpFIvFnc0hAODIkSMvffjNkE1/erLBpwMxYNHzv3e3jDCEEtIBiiK+HODCAAIG8GCh3hzgw0HgsE8HxkDB10Ox74rI4Ufw7cOxKX2cpUMMapPdSs05Q8wOQv6TjMENnc06eGvwQR4eHg8+yeFwOJz/2bbRx8en48dksVgPhkw6nd61KGjmxQHj/ZCeuC3SuwOhCx280RfBMFgPOqN/RaLxrsjcs8TVZuq9JFgbOpGdxeTaTGLbEGxWEPzt25333nsvODi4J47sykR6aKAcPI2gXmygB5I5nXajhRp7FM64dwpqE1h8ntiST2ak0WAUtE9z58719PS0dS46B55JUO/mzgLHJ9LG+KD9D+HHa2EsdGQ3JVT/QzgDA9en0SKc+N4wZHEwEEK9HoqAtxLQX0fTVlwiVl2Fi3Q7IAqAz26RE47jGxLRHcPh8r+QhcFACDmI4V5IzgxavRYM+BPPlcKmoeOo0VATjuEHq8hrU2lPhMAqC7I8eFZBjkPEBL+Oxl7vi44/hr9/k8Rtv80Z1F17Ssn+h/BR3uj5ybRAB9pcF7IrsIsBcjSLQtEUb2TZReJQJblzBNbXFdaevVKthlpxiajTglOpNNt+iUlJSfPmzXtwMgBFUQhcnNhatFrt0KFDe+jgvXtCvU6nYzAYXZhAClmEPU/opgDYXUKuuUH8KxJd1w9zyOUU7Ln8u4OkwLdF5Pos4sUYbHU8Srd1v9WjlpnWaDS23TzI2Xh5ed07tRFOqIegx0AAWBKOTvRDX7xCxB3Avx6KjfWFF++9QK6UWnGJoKHgfBot2j523WKz2Q+dGOeoFyJOCAZCyJF5c8D+MVh6DfWvS0SyO7J5IOrHtYu6FXqQ3Ag2ZBO/lJOb+mPPRPT4VkoQ1MbWnQ4Q1PMm+yMFs2jRQiThIP7+TVKH2zpD0P8iKbCjmIz+3aTHQcFs+jIYBSHrgoEQcgpsGlifiN6YRsuTUlG/43vL4T5O9uJMPZV0CP+pjDwynvbNMEwM90SDrA52jUJOJJCP/Doau9RIvXqN2HyL/CgZ3ji0pZsSas0NolwJPkxG4XppkA3BQAg5nWFeyNVptAMV5POXCV8O2NQfG+IJw6FVFcmpDdnkhQbyrQTsX5G2HxcKOTl4AkLOCAFgdhBaMIu2MAx98jyRehz/uwn2lVpDkZxafJ4YmY73EyNl8+jPR8MoCNkePAch50VDwdJwtHgObVYQuvg8MTodP1kHw2FPyW6l5oxy3X8AACAASURBVJ0lRqXjUSKkdC59TTzKhR1SkH2AZyLk7BgoWBaBPh2G7rtDvnaNoCHgtb7onCDYUrEMCoATtdRnt4giOXglDv1+BB3GP8jewFMSggAAgIaCRaHowlD0aA31+S3ijevk89HosgjU/f51taCOUpvAT2XkFwUkEwOrYtH5IfDaArJTMBBC0D8QACb7I5P9aXlSamsBGbHfNNkffTYSHeYFR9N0Qr6M+q6I3FtGjvJBtw3FRnnD0oPsGgyEEPQQfV2R7cOxjwdgu0vJf10iAABLI9CFoagX29Y5s2NKE/jtDvl9MVmjAUvDkZszaXAdH6hXgIEQgh5JxASrYtFVseilRmpXCRn9u2mwB7IgBJ0eCAd6/MNEgpN11N4y8mgNOdoHXdsPS/VHMBgBod4D/poh6PGGeSHDvLCtOHawktxbTq68TEzwQ+cEIan+KMdZf0MmEpxvoPZXkIcqyXABsiAE/WIIHa4LA/VGzvojhqDO49DAk6Hok6GoxAAOVpLbi8mlGUSKDzo1AJnk7yy9pgojOFlH/llFHashIwTInGD0rem0PjzYAIR6MRgIIajTxEywLAJdFoHKjeBoDflnFfXaNVMIH5ngh4zzRYd4Ig42PJKkQI6EOllHnawls1qp4V7I1AD0kwF0b87jXwtB9s+SgdBgMDCZ7fWMGI1GBoNhwXeEINsSMsCCEHRBCMBJ7O8m6kQt+do1okRBDfZERnmjwzyRZHeE2Tv3BCYokCuhLjVR5xuoCw2kJxsZ54u81hdL8UactjcYclSWOaNra2vnz5+fl5fHYDA+++yzRYsW3Zdg/fr1u3btamlp4XA4L7zwwoYNGyzyvhBkJ2goGOmNjPTG/pMMZAaQ0Uieb6BWXSUL5VS8GBnojiS7I0luSKgLYs87DFWrqaxWKrOVutpMZbZQflxkmBcyOwj5aghs/EGOzDKB8OWXX46Njb1w4UJWVtbo0aNHjx7t6+t7bwJ3d/dTp05FREQUFRWNGDEiJiZmzpw5FnlrCLI3IiaYFoBOCwAAAA0OMluoay3UgUpqXSYp0VNxrkicKxInQiKFSJQQsWGAkRhAsZy6LacKZNQtKXVTQtFRkOSGJLujr8WhAz0QVzjyBXIOCEV1d3FFuVzu7u5eWloaGBgIAEhLSxsxYsQbb7zxqPSzZ8+OiYl59913H/q/fn5+165duy+OPopOp2MwGBjWO/ueej+VSsXn822di95EZgB5UipPShXIqCIFdVtG6QgQ5oIE8ZEgPujDQ/y5wJuD+HKBB+vxfaodKX8TCVr0VKMW1GmpGjWo0VCValChokoVFEmBCCESI0QihUi8GIkT2TIq90bw/LctC5a/BVqEVVVVDAbDHAUBANHR0eXl5Y9KLJVKL126tGLFikcloChKoVBwOHd/kVwuF95WhByGiGnuQf2ne1RhBKVKqkJFVapAmZI6Vw8adGS9BjTrKToKxEzEjQUEDCBgIGwM8OiAiYG2W3RGI43BIMyPDQTQ4kCHAx0B5EZKbgByI2jRUxoceLAQLw7w4QA/LuLHRab0AUF8NNQFgavHQZCZBQKhXC7ncrltf/L5/IqKioemNJlMixcvnjBhwtixYx91NK1WO3jwYBS9O+ru22+/nTRp0qMSwxahbWk0GgSx41tevQEGQCQLRLIAcL//v1Q4IjEAmQEoTIjSBHQ40BKIgQA64m6Z0ygjE9x9LKJTHDZgYRQLA0I6cKFTAgYQM4GI8YguHxyo1T31oZwEPP9tq4Plz+Fw2gLKo1ggELq5uSkUirY/ZTKZp6fng8kIgjAPotm+fXs7R+NyuR3vGsUwDAZCG6Ioisfj2ToXDosHgHe7CVQqE58Pm3U2A89/27Jg+VtgulNgYCCNRisoKDD/mZOTExUVdV8akiSXLFkilUp///132NUJQRAE2Q8LBEIul7tgwYJ169Y1NDTs27cvJydn/vz5AIC8vLy0tDRzmueee+7MmTMvvPDCpUuXTp8+XVhY2P33hSAIgqDus8z0ic2bN69atWrIkCHe3t5//vmnq6srAICiKBzHzQnUanV0dPQXX3xh/nPSpEkPthohCIIgyPosMH3CsuD0iV4EDh+3LVj+tgXL37YsWP6OtSQiBEEQBHUSXDTQYcn08vyWwvyWoiplTYO6SaKTmUhcj+tdGHwBy8WD4xYiCgoTBSV59ROzRbbOLAQ5F7VRk92UVywpLZdX1qubFHql0qikoXQmxhAwXbx5nv4uPtFuEXHu0V5cD1tn1vHBQOhoWrWSU5UXLtZcrVbWxrpHxrhFTg2b6M31FHNcGSidRWMpDEqFQdWoaSqXVf5de31r1g4PjtsI/yETg0d7ch+YywZBkOUoDapTlRfOVV28I6/q6xEdJY5IC53gx/cRMPkCpouJxI2EUaZXNKqbKpU1l2qubcveJWC6jPAfNCZgRIDA39bZd1jwHqHjyGrMPViSnttUMKLP4FF9hiZ49qWhjy8ZkiILJSWnKi6crboYIQ6dHzUz0atvB98R3iOxLVj+ttWp8i+TVey7/cfV+szBPsnjg0cleMTRMfpjX0VSVJGk5EL15VOV5/35vtPCU0f1GYoi8JYWABY9/2EgdASX667/mP+b1qSbEzl1TOBINq0rk6xNhOl0VcYvt/9g0VjL4hcmeyc89iWwIrYtWP621cHyL5GW78j9qVxeOSdy6pTQCVx6V1Z0xUniUu3VA8VHpDrZwpjZ44NTMMTZ6z0YCO+CgfB2a/G27O91uP6puCeG+Q1Cu73gEwWojOorO3J/dOe4PZ/0TIgwsJ3EsCK2LVj+tvXY8m/Wtn6bsyenKe+puCcmhYyjoxa4FZXbXLDn1i8tWsmzCU8N8xvY/QP2XjAQ3uXMgVCqk23L3pXbnP9M/MLxQSndD4H3IijiSNnJXXn7xgWNWtp3waOamLAiti1Y/rbVTvkTFLG/6K+9BQdmhE9+InpG1zpp2nG9PvvrnF1ClmBV8oo+Ln6WPXhvAQPhXc4ZCEmKOlRydM+tX9LCxi+Kmcui9dSucQqD8uvsXdlNt1YPeiHJK/7BBLAiti1Y/rb1qPIvl1d+cHmLiC1clbzCh+fVQ+9OUuTBkqM/3Po1LXTcU3FPMDCnW7oSBsK7nDAQVivrPr76BYqgrw1c2celQ6XUTVmNuR9d3TrIJ+m5xCWs/72whRWxbcHyt60Hy5+kyJ8Lfj9QfHhFwpKJwaOtkAepTrY1a0eZ7M4bg16Mc3eu5bpgILzLUoGQMJBGhcmkIQgDSRgIAABKRzEGyhTRmSI6gtrFTisUoH4vOvxT/v6n+z4xLWySZftC26cxabdmbi9oLXp76GvhriFtz8OK2LZg+dvWfeXfpGnZ+PdmJo3x5qCX3Dhia+bkYs3VLZnfjgkYvix+of00DQ1yk15iJE2UuVLFGCjKQOk8GoNHo3Et0HqBgfCuLgRCkxrXNhl0LUZ9q0HXYtRLjAaZiSIohoBG59FQBkpjYQAAwkiSRlIvNZpUOMeHJQrniaL4LkE228C7Wdv6wZUtJgJfO+Tlnutsad+ZqotbM79bGDN3VmQaAhAAK2Jbg+VvW/eW/8WaK5uvb5sXNWNe1AxrXqS2URpUm69vq1LWvjVkVago2PoZMNM26KW3VbJitapaR2OiLDcGxkQx5t1KlTCQJjVuUuEkQbFEdJaYwXJjsN2ZbHcG253JFNJBZ0oOBsK72gmEJE4ZFSaDzKSXGfUSk77VqGs16JuNAANsdybHg8l2Z7DcmCw3BktEp3EeGUpJnFJVaeUlakmeEsEQn+Fij/5CBLPqiX6++u8tN76ZEzltfvRM204halA3rb/0kRfXY/WgF7l0DqyIbQuWv22Zyx8niW9zdmfUXHl3+OpIcZhts3Sy4vy27J0LYmbPiZyKdCqqdBMFJPnK+osSfavRLd5FGM5zCeZizEdWVoSBNEiNOolR32rUtRh1LQZdswHXk2w3BtudwRLfrZmZIjpDSMcYDz8ODIR3qSWahjNyBEEAReE6kjSRuJYwaQmTGieNJMOFbu7bvHvd4cZguzG73iSngLxUXXe+1SA3hc7xtU7rUI8btmZtv9mU//bQV23+GzMzEaYvs3dmNtzcNHKtGyqCFbENwUBoWyqVykjD11/8iM/gvjnkZReGXXwXDeqmjX9v5tI5a4e8LGIJrfCOumZD2e/1hJ70G+0m7uvS5XtJhIHUtRj0rUZd692+OoPcZJCbEAQwXGg0Dkbj0DAGynChBc/wBjAQtlHLNcoCHYqiAAAaB0PpKI2N0TgYg09rp5HXTZI85Z1DDeK+LkFTvXr09mGFvGrDpY8jxGGrkldYfPh1N52qPP9V1s5nY59KjRhr67w4LxgIbet6VfbH2V9OC5u4MHaOVZtfj0NQxO68X47eOb128MsPHe9tQfUZkprTLX3Gu3sPFfdQGRAG0qg04VoC1xKEkaQI4J4oADAQtrHVqFFcTxT/UAMQJHKxfzvN/+44XHZix82fnktcMsEqY8+6oExWse78+6MDhy/vtwiu+WQTMBDa0F+lx3fm/rxuyKoBPom2zsvD5TTdev/y5xOCU5b2XdATy9BQJHXnYIPyjjZ6WQBT9Pjl4iwOBsK7bDh9giKp8gMN6mpt7HNBNLYlM6AxaTdf+6pKWbt+2BvWmSDRZfWShk9vbqOhtHeGvsZjcG2dHacDA6FN4CTxReZ3ec0Fa/u/HO4VauvstEduULx/+XM9bnhn6KvuHDcLHpkiqMLd1RRBRS7ug7Fscx0M9yO0PQRFQuf4CMJ4t3dWk7jFLiZKpOXLj63iMXhfT/jEzqMgAIDP4H0yekMfF98VJ16rVtbZOjsQ1OPkBsUrZ96S6KTbJnzizfW0dXYeQ8gUfJyyfpBP0r+Ov3qlLtNix6VA2W/1AIDoZQG2ioKW5QifwYaCpnixXOlFe6op0gKx8I/i9DfObVgev+iVAf+2n8lA7cMQbGXSsiejZ7146s3r9dm2zg4E9aAyWcWzx16N94jdOGIth862dXY6BAHIkzGz3xu+5vMbX3+dvQsnie4fs+pYk6ZJH7HI307mWHcfDITdg4DQub6Enqw53dKdw6iM6rcyPjhecWbbhE9SAoZZKndWkxoydtOItR9d/eK3wkO2zgsE9YiMmiuvnX1nRcLTz8Q/aZOZgt0R5x61I3VLtbLuhVOrG9RN3TlUa66y5aYiZlnAo2Y19EaO80lsBaUhkU/5N16Wqqp0XTvCrZbCZ46+7M31+Gr8x7aaLN99se6RX0/89FTlhQ+ubDERJltnB4IshgLUnlu/fJm14+OUDb3xOtXMhcn/z6h1YwJG/vvEa+eqLnXtIEYlXv5HfcRCfzrPoTZ1h4HQAug8Wsgsn+KfaggD2akXkhS559Yv71z8cFXyiueTnrHILi025MFx2zruQwNhfPH0WolOZuvsQJAF6HH9hosfX2/I/mbCp/euL9gbIQCZHTnl45QNO/N+/ujqF3pc37nXU6D01zqfoa78Pr2jW7jjYCC0DHGci0sQp+poJ/ocGjXNL55am9tcsCN1y2Df/j2XN2ti0Zjrh70+xDd5xfFXiySlPfdGJE5pmwzyErX5n7bR0NmrEKg3InFK12xQlGnM37umTt+j33ujpvn5k6s5dPaWMe+7skU990bWFO4asj31c4qilh9bVSwp6/gLG65IcS3uN9a95/JmK727CWJXgmd4Z31Y6jVIxPF+/OT3E3fObsveNT9m5tzI6b3ufkP7EIAsip0bIgpcc/49i6/Br2s2tN5UtN5S6poMTFcGU0AHCAAUMCpNepmJxsb4AWyXQI4wgsftwLcA9Qp6qVFepFZWaFXVWoPMxBTRGYK7S+GbVLiu1cgU0sSxLuJ4gWVbKtmNeRsvb14YM2dWRJoFD2sP2DTWmsEvnau6tPr8u7Mjpi6ImfXYqcC4lqg+3hz3XJDDDJC5F5xHaEn1FyWyQlXMvwLbSSM3KDZf21ajqn9ryCuhoiBrZa1HtD+Pp0pRsy7jPwO8E59PWtr9+bx6ibHqWLOiVO2WIHDr68IP5Dz4g9RLjapKnbJSIytUUwQljnNxSxC4BHDsadEPS3LseYTaBn3LTYUkT2nSEqJIniCEyw/gcDyY93+bFNDU61tvKVuy5CwxIzDNk+dngXC4v+ivvbcPrB/6Wj/PuEelcYDyb9G2fnDlv3rcsHbIy358n3ZS3jnUQOJU6Oz20lgZnFB/l70FQoqgsj8pC57mJYp6+NdzofryfzO/nRA0emnfBXTMBmsxWNZjT0S1UbPp8matSffu8NXdWfaw7kJrzekWn+Fi35FuHVzKR9dsaM1VtuQoCAPhmSzySBayxL1jRkrHOUBF/CCTGm/Okjdfl+N6wr2fwC1ewPNnd+RShiKppquy6pPNbvGCoKleXV4ZX48bPrn2ZbWydtOINz25Hu2kdIzypwB1sPjo7lv7noqbNyM87aEdVLoWY94X5Ymrw+xqjAwMhHfZWyAEAEgLVJXpjQmvhd7XXpHp5VtufFuhqF496IUYt0hbZc+yOnIikhS159Yv6eWnNgx7I9a90x+cMJJlv9bpJMaop/swhV25dNDU65uuy1qyFVwfltcgkTjOxcqbh/Qcx6iI76KAvETdeFUqL9G4xvI9B4gEwdwuNOUJPVmyt9akJSKf8mfwO11r16ka3s74IMw1+JUBzzEfN5fXkcq/VlX/0dWtFEW+MejFB5fyKNxVze/D9htjX3cHYSC8yw4DIQAgb+sdn+Fit34C858UoI6Wn95+84fJIeOeinuit8yU74iOn4hX6jI/uvrFkzGzZ0dO6fjyxISBzP+mguPJCpntg9K6Fb0ogpLcUjZekWkb9R4DRN5DXG2yOqJlOUZFjGuIpuuyhstSGhvzGixyTxB2d7ESCtScbmm6JotbGdSpiyfznoJL+i6YFpbakfSOUf5tSIr6s/Torrx9syKmLIiZ1TaIXVOvL9he1X9deDd/gxYHA+Fd9hkIpbdV1ceb+70SAgCoUFR/dv1rnDS9OuD53n5H8EGdOhEbNc3vXPywbS/Dx6Yncer2jiqWKz10jq8Fb/LpWgyNl6XNmXJ+IMd7qFgUweu9dxB7e0WsqtY1/C2R5qvEsS5elh6UX3ehtemqLG5lML0DO6+ZSPy7nD2d3VOwt5f/QzVpWv6b+W2tquGV5BXm+6MlP9dyvFl+oy25VKlFwEB4l30GQkCB7E9KfdJcD+gOnao8v6TvgimhEx1saKhZZ09EE2HalrPral3mO8NeixKHt5eUAkU/VAOARCzy64lRaqSJbMlRNPwtxbWE1xBXzwGijlSX9qaXVsSEkWzJVjReluJ6wnuwq+dAUQ9tmlZ1tElWrI57Pqj9NVAa1E3vXvpEzBatHvxip/YU7KXl3xEXa65+mbUjxi1iefDiqq+lyesi7HBNURgI77LPQEhS5Kn0v1uzFQ0Tqpb3WyRgutg6Rz2layfixZorn934Zm7ktHlRMx51fVCfIWm9qYh7Pqin7+epq3UNl6WSW0pRFM9rsGvX7kvZSq+riDX1+sYr0pYchSCY6zXE1QrN8ZK9tSgdDZ3zyLGOZ6oubs3cvjBmzqzItM7uKdjryr9TDITx54LfVSeMIW4BYxcNZtnZlqgABsI2dhgI/669vv3mD25M8ay/58YtD+H52t3ZY0FdPhGbta3v//0ZiqJrB7/84O4wmnp9/jeV8S8Hs1ytdD8V1xHNmfLGKzKKoDwHijz6CxkudjQ67lF6S0WM64nWHEXTNZlRhXsOEHkOFHVt3FMXEAYyZ3NZ0BQvcdz916Mak3bLjW+LpWVvD301TBTchYP3lvLvMlxDXP9P8aXR564rs56Om5caMrYn9jXsMhgI77KrQJjVmPt93s86k355v8WDffvXnmvVNRvC5tn7Vkrd0Z0TkaSovbcP/F7053OJz4wPGvXP8zh187NyvzFuHkldn27RZaoqXdM1aWuu0iWI45EsdI1xsbcBAvey84qYIilFqabphlxWqBKG8zwHimxyR1ZVpSv8vqrfq6H3XtzkNN368Mp/B/okPZe4lEVjdvHI9l3+3Vd3vlXbYAib71skKf3u5g/N2tYlfRek9BlmJzd6YCC8y04CYWbDzd23flEYFE/FPTE6YIT5LDGp8awPSpPfieihLeztQfdPxFLZnfcvf+7P93llwL/NEw2rjjfrmg2Ri/0tlMeuII1ka56yOUuurtGJ41zcEwSCUK4dLqhhpxUxBVTVutab8pYcBVNE90gSuicKe+guYAdVn2jWNugjn+4DANDjhu25P1yovvLGwJXd3FzeTsvfcrI/Kg2d5+sSeHdoW1Zj7q68vSqjenHsvJSAYY9djKanwUB4l20DIUERGdVX9t4+YCLxJ2NmjQkYcd+ZUbS7WhjJ9xrkIEsUPsgiJ6KJMO2+tS+9/PRziUtHiobe3FzW79VQq3Wdtc+oxFty5K03FXqJSRzn4hbvIgjh2s80RPuqiCmgqta25ikluUqUjrj1E7gnCtnudjFZiMSp7A9Lw+b7lnPLP7n2ZZxH9MqkZzo1Luah7Kv8LU1ZoS37rS5x9f1jaDMbbv6Q/2uLVjIvavrE4DFdbk93HwyEd9kqECoNqvTyUwdL0r25nvOiZwz27f/Q2+yyQlX1yeb4l3r3ivXtsOCJWCIt/+jq1rGF4/uGR8ZM6coNmx6llxolucrWPKWu1SCK4Itj+cIIHo1t464Ie6iISSMpL9NIC1TSAiWNg7n1FYj7unB97O7WePWNxtvHK3bF7XxlwL8H+iRZ5Jj2UP49p3RfHceH6Tvy4bMmClqL9t0+eKv59qSQsdPDU9tfgqeHwEB4l5UDIQWovObbh0tPXKm/Mdx/8KyItPbvsVMklbmpJGZ5QEeW4e6NLFsRKKo0N3eUfhH33+nRqU9EzbDPJeiMSlx6WyUtUCrKNVxvljCCJwzn8fuwbdJxarOKmAKaBr28VC0vVisrtDx/tmsMXxzjwnKzi/bffShAHS8/813OjytKVsSMDvEbZLEq24EDIaEnb2wsTnrzMWuqNaib/ihJP37nTIxbxJTQiYN8k6w5mgYGwrusFghrVfWnKy+cuHOOiTEmh46fEJzSwX6VquPNhJ4Inu7d0zm0CctWBLe2VXgkC6lo05dZOyrk1S/0XzbIx353pyJxSnlHIytWy0vUBonJJZjjEsJ1CeLw/dlW6zu1ZkVMkZS20aC8o1GUaxXlGhobFYbxzNcB9nwXvFhStiXzWwQgLyc/663yKdpT3X9duKW+IAcOhI1XpPISdeRTfTqS2EAYz1ZdTC87WatqGBc4cmzQyAjX0J7OIYCBsE1PB8JqZd2l2qvnqi5JdNJRfYZNCE7p7Beslxpz/3tnwPoIOxxq0X0WPBFVVbriH2uS1oaZC+p6Q/bWzO1eXM/nEpcECQMs8hY9x6QhlOUaRblGWaHVtRi4Pix+Hw4/gM3zY7PEjJ4bJ9nTFbFBbtLU6VXVWlWVTlWtZQro/ECOIIQrCOXayU3cdrRoW7fn/pTZkPOvfk9NCE4x37wo+K7SrZ/Ac4Blbts7cCDM23rHf6z7ozYPeJRaVf3JinOnKzMQgIwKGDrCb3C4OKSzszM7DgbCu3oiEBoJ462Wwmv1WVfqMnW4fpjfwJF9hsR7xHR5iFTulvKAyZ7CMJ4FM2knLHgiFu6uFoZyvYeJ257BSeLP0qM/5v821G/gkrj5bhxxOy+3H4SBVNfoVFVaVbVOXasjdCTXl8XxZnF9WBxPJseTacHxk5atiAkDqWs2aBsNmka9pk6vqdcDAHj+bL4/mx/A5vfh0HrJ4jsak3bf7T/+LD02NWzik9GzOfR/Vm6Tl2ru/FGf+EaYRSpnRw2ERoUp59OyARsiu9x0LpaWXai+fLHmqh7XD/RJGuiTlOQVf+8XYREwEN5lqUBoIIyFkpK85oKbTfmFkpJgYYD5ywt3tcDlTN25Vp3EaFf7eFmKpU5EXYshb2tF8lvh6ANLYamNmr23DxwuOzE5ZNwT0TOETEH3386acA2hrtdp6vXaRoO20aBrMgAMsN2ZbDcGS8xgiRlMEZ0pojNc6F2YsNi18qdIyqTC9VKTQWbSS416iVEvMeqaDbie5Hgw2Z5MrheT68Pi+rAYAntv9t1Hj+sPlhz9tfDgYN/kJX0XeDywVgMA4Obn5X3Ge7jGWOC8ddRAWJ8h0TToLTIHukZZd6U+81p91u3W4iBBQIJXXJx7VJx7dEdWG34sGAjv6nIg1OOGCkVVmayiRFpeJCmtVtYGCwPj3KP6ecbFe8RY5Ev6570ct3fUUidi2W91DAG9z4RHjmJo1Ul/yt9/pipjSuiEOZHTRKxeFg7vZVLhuhaDTmLUS0wGqVEvNRlkRqMSp7ExOo/GcKHReTQ6F6NxMBoHo7EwjIliLBRjYQgKzONUERQxL/yoVqt5PB4AgDSSJE4BAHA9CUgK1xGEkSQMJKEjcR1h0uC4ljBpCKPSZFLhJg1B52JMVwZTSGe5MlhiOsuNwXZjMoX0XrS83H30uP7P0uO/FB6M94hZEjc/QPDIeaituYr6DEnfFywwMtlRA2He1jv+4zxEkZbsxDISxoLW4pymW3nNBcXSMg+OW6Q4LNw1JEQUFCoM4jG4XTimBcu/F6wj1U0Kg7JFK2nWtjSom2pVDXWqhmplrVQnCxD4h4qCwkQhqcFjQ0WBPbc7EsuVwRLRFeVaYVhXvmyHZ1LjrXnKpDfbW4Pbje36cvKzC2Jm/Vzw++LDz40LGjUvaron1752R+sgOp9G59Ncgu8/GUxq3KTCjSrcpMJNWgLXEPpWI64jCANJGEhCT1AkwHUEAIAiKMJAAgAoikIQBACAMlBzg5LGQgGC0DgYSkcwFkZjoTQ2xhTReX5sOhdjuNDpfBqDT+u9Ae9BSqPqYPHRgyVH+nnGfTb6vcfeURbHuVQeaVJVafkBlrzedRhGhUnXbLB4ZcXAoA0KTAAAIABJREFUGAmecQmecQAAgiIqFTVFktJiSdnZqot35FUcGjtA4O/L9/ble/vwvDw4bh4cNyFLaLUlbHp3i1Cikh6rOoOiqNakI0hCT+j1uEFr0qmMaoVBpTAo5Xo5m852Z4s9ue6eXA9fvrcf3ydA4OfF9bDmsgi151r1jtg7apErsrrzrdpGQ9gTHe2HkepkvxX9mV5+Ktk7YV7k9AixNcan2SdHbZF0UL26cX/Rn6crMob5D5wfPevB7WQfxVLLHzpk+VuwX7TjmrWt1YraWlV9naqhQdPUrGlt1raqjCoB00XAdHFhurgweBw6m4kxuXQOi8Z8Ku4JAFuEbUiKVBs1CIJw6GyMjrlhYhaNyaVzeHSugOUiZLoIWcK27SVtyC3eJXfLnZCZ3o7XO9p9TddkoXM78atzZYtWJDy9KHZuetnJdy5+KGa7zopIG9FniD180ZAVkBSV2ZjzR/GRwtbStNBxu9O+FLM7NwrUs78w68PS4One9jzxw1ZacxX+46w9O97cBOzv3e/eJwmKkOuVcoNCaVApDEodrjfgBi2uY2GWn5bdu+sOHp27PH6RzdcafSyWK4MpoisrtIIQ2Dv6P1RVWooCbYsZdhyXzpkbNX125NS/a68fLEnfmrVjYvDoySHj/DvcLIB6HYlOduzOmfSyk3wGb3r4pHeHr2F26Y4GnU8ThHJbbyo8Bzrs8oddY1Th2ibL94t2DYZgYraos1c5XdO7A2Ev4hrNlxWqYCC8T9M1medAUZdvWaEIOtx/0HD/QXWqhsNlJ148vdaX5zUxeMyoPkO7dvsdskNGwni57sbxO2fyW4pG9Rm6Ydgb3e8P9xwoqj3TAgPhfWS3VcJwnv2spms1MBBaiSiSX/prXWCarfNhTwgj2ZqnTHzj/lV9u8CX770i4enl/RZdq886fufstuzvE73ixwaOGOTT34aLAkPdgZNEdlPu2cqLl2qvhbuGTAwevWHYakt9m6JIXtn+em2TgeMJT49/yIpUrtEOu5F4O2AgtBJ+H7ZJjRtkJqaol83N6jmtuQpBMNeCW+BiCDbEd8AQ3wEakzaj+nJ6+amPr25N9k4Y7j9okE9/2EbsFQyEMavx5sWaq3/XXvd38RnVZ9jyfost3j+GoIhHf2HzdVngFC/LHrn3oghKXqIJmeloY/o6AgZCa0GAKJInK1R5DXG1dVbsRUu2oof2qOLSOakhY1NDxioNqr9rr52tuvjZ9a8jxKGDfZMH+fTv+PBCyGpatK3m5Zxymm5FiEOH+Q181Ix4S/FIFBRsrwpM83KkySTdoazQstwZdL4zBgVn/My2Iorkt+TIYSA0w3WEukonWtKhVX27zIXJN0dEPW7Iasy9Unfj96K/EIAkeyckecUneMX1uqVqHInWpMttLshuzL3ReFOqkyV7J4wOGL568Ivd3ymwIzjeLIyJqqp1/AALL/3VS8mKVK6dXFzUYcBAaD2iKF757/UkTnVhMS3HI8lTCiO42ANrqvUQFo051G/AUL8BAIAqRc2NhpsnK85/ev0rd45bgmdsX/eYOPeo3rKcaa+mNKjyWwvzmm/nNudXyKuj3SISPfuuGfRiuGuo1WZPtxHHC1rzFDAQmklvq608fdB+wEBoPTQ2xvFmKss1wggHXIC7s1rzlJ7JQpu8dYDAP0DgPztyCkmRpdI7N5vzT1We33LjWyaNGesWEeUWESkOCxMFw1E2FoGTRLm8okhSWthacltS0qqVRLmF93WPXpGwJFocbttdJ936uhR+Xx0Ee0cBMMhNJjXO7+Ok1wQwEFqVKIovLVTBQIjrCFWFNnLxIxeEtA4UQSPEoRHi0HlR0wEANcq6QknJ7daSs1UZd+TV3jzPMFFQqCg4VBQULAzs1QucWpPGpL0jryyXVZbK7pRK71Qpa315XhHi0Bj3yNmRU4OFAdZc1Kl9XB8WgiHqWh3P30kDQBtZoUoUwXPaCwIYCK1KGMYr219n61zYniRfKQjj2tu6Hv4uvv4uvuODUgAAOElUKKrKpHdKZRVX6m6UyypRBA0S9gkQ+Ae4+Pdx8fVz8fHguFu/N8/etOqktar6GmVdtaK2QlFdpahRmzSBgj4hwsBw15BJIWNDhEH23LZ2i3dpzVXCQCgv1TjtDUIAA6GV8fxZBpnJpMbpPKcueUme0r2fXTewaCgWJgoOEwWn/v8zEp2sUlFdpaitUtZcqr1ao6yTG5TePE8fnqcPz9uL5+HN9fDkenhw3RxyAI7aqGnWtjZqmhvVzY2apnp1U4O6sU7VwKIx/fi+AQI/fxffJO/4QEEfT657z+3FanHieEHR7urANE9bZ8SmKKAo0wQ58UwSp66OrQ9BEZcgjqJc4xbvgHVlB5EmUlGuCX/Sz9YZ6Rzzak9JXvFtzxgIY72qoV7dVK9ubNQ05TblN2lamrQtetzgzhG7sV3dOW6uLKErW+TKEgpZAiFTIGQJBEw+i2b5xRK7yUgYlQaV3KCU6eVyg0KmV0h1MolOJtFJW3XSZk0LiqAeHDcvnocn18OL6xHtFmHeKMCye5ZZH8+XBSjg5DPrtY16jIk68xRnGAitTRDGVZQ5dSBUlGl4fmway95XiH0sJsYIEgY8uO+PgTC2aFslOlmLtlWql0u00gpFtUwvl+nkCoNSaVSRFMVncPkMHp/B49K5XDqbQ+fwGFwWjcnEmDwGl4bS2DQWA6MzMSaKoG1be3NobAz9p9w0Wo0K0dz71iqj2vxAh+sJkjCRJj1uMBImA2HQmXR6wqAz6dUmjc6k15g0apNWbVSrjBqVUUWQpIDJFzBdhCyBK0skYglc2aJgYYCY7Spmu3pw3Cy+vbj9EEbyZIUqZw6E8jKNnawvaiswEFqbIJTXeKXG1rmwJWmhWhTpyHcjmBjDj+/jx3/kCh0GwqgyqlVGtdqoVhu1WpNWi+vURo0O1ysMyjpVw/8HMKOBMJIUqTXpzC/U4jqCJNqOQ5Ikiv7PfVY+4+44LBaNSUNpdJTOojHNAZVNZzMxBp/J8+J5sGlsHoPDpXN5/x+P2fbXSLUa1yh+fYbEd1QPTt63c05+aQ5gILQ+ng8L1+BGhYkhcNKOCFmRKqqH59HbOSbGYLJd3djdXVrBIffDsz5BGLf45xrCQNrb6C0roYDyjibE4XZL7RSn/OJtCwEuIVxFuebxKR2RrtlA4RTXy3nbH5C9wRgovw9HXqq2dUZsQ12no/NoDKdcWa0NDIQ2IAzlKsqcNBDKitSiKH7vGVQIOQVRFE9W6KSBUFGqEYQ5+8xmGAhtQBDGk5c6aSCUFqr+r707D3OqvOMFfrKfJckszGSSySwMu6AyrIp4RUEEH6Vi3alLxV5RpI9K7fNoaWlrF7V4+6DtVSsFfGxvq2LlqeC91VpFaQFRBIEKDMswazL7TJKzZDk594/Y6ThrZnKSs30/f2XCmeQ15sz3nPf9ve9bMM3oZx2oTeEFrq6TYaVboYyes2zeJENXyhAIQkXQHocYTUa740o3JNfEWDJcx+VPQRCCulAeh8ls4oJRpRuScxIROs/lTdD2HJjMIQiVYCJclVS4jlO6HbkWOss6yyiDliSAuhVMdXadMtxNIRcUbIzF4Ot7EAhCpbirmFCt4YKw5wybb/hOGFCn1ARfpVuRa6Fazl2FUxJBqBB3FW3AIOw+g2F5UKm8iUzoHCclJaUbklOhWs413uj9ogSCUCnOcopvjYrRpNINyZ2EIPKtUaxuDOpkc1rteTa2UVC6ITkVOs+5qxCECEKFmK0mxkdGGnilG5I7obOcq5LGpsSgWvmTmW4j9Y7Gw4kEJ9Ie464t10u2IOQ4rr6+Ppkc8hZHFMW6ujpBMNYF1zBcBusd7T4dMfh6hqByeZOYnjMGmk3YU8u5q2hM6iXkCsLnn3/e7/cvXrx42rRpJ06cGHjAZ599NmHChKVLl5aWlr766quyvKnWucfToVoDXX72nMF0JVC1vIlM6DwniUYZJgzXsugXTZEhCOvq6n7wgx/s27fvzJkzq1atevjhhwces2bNmu9973s1NTXvvffe2rVrOzo6Mn9frXNPYMLneYMMzic4MdoZd5ZhgBDUy0pbqHF24wxYoFKmlwxB+Nprr1111VUXXHABQRBr1679xz/+0dLS0veAEydOfPnll9/5zncIgpg7d+7MmTN37tyZ+ftqnY2x2FwWrsUQc3h7zrCuKtpkQS8MqFreJKdBhgmT8SQXjLpQvEYQhCy7T9TW1k6ZMiX12OPxuN3uurq6kpKSvgf4/X6a/urSY/LkyefPnx/q1ZLJ5NGjR4PBYOrHCRMmFBQUDHWwJCbijfWiWaslP0xxoutwnS2h1faLHBej07qi7DyScBaZYg2ns90kQ0n/84c0MfnJ4GGxZGp3Ogdr+vMPN0tkoZRoOat0Q0bNZLbY/BPkfU0ZgjAcDhcV/XcrL4Zhenp6+h1AUdQwB/QVjUa///3v22xfbVH05JNPXnnllUMdzLU2J//y3Jhbrjzuws7aIsvxPUq3Y4ySySSf3lVId8etxa6PO04Hs90kQ0n/84c0JSVHpPWu9j//bxMx8piFpj//bvZii5jX8ee9Sjdk1EyUk/n2jwiCiETSqmyiadpiGWEbcBmC0OPxdHV19f7Y1dXV93YwdUB3d3ffA6ZPnz7Uq1EU9e677/r9/nTe2mqtsH/vNyP+R6qWq54/vaPJ971blW7IGKW5H54YS57feLLi8R9h7oS8sB9hNrQ8czrvW/+L8Y+8U5imP//QHxu801yeuXco3ZCMyPX5y3A5M3PmzIMHD6YeHz9+3Gw2T5w4se8B06dPb21tbWpqSv148ODBmTNnZv6+OsD4SaEtJsZ0Pq0+XMcxfhIpCJpgkFWfwnW8qwIDhF+RIQhvueWW+vr6Z5999siRI4888si9997LMAxBEI8//vjTTz9NEERJScnNN9+8bt26L774YuPGjZIkXXvttZm/rw6YLCba62CbdD63MlyL1StAM1zj6dB5nQdhghMTnEgVYyr9V2QIQoZh3n///QMHDjzwwAPz5s371a9+lXq+rKzM5/OlHr/44ouVlZX333//mTNn3n33XavV6Iud93JW0OF6nZ91ofOo0gbNcI+nw3oPwnA97yyjMJW+lzyBdNFFF7355pv9nly3bl3vY7fbvXnzZlneS2dcFZTOdwSViHA9P2UVghC0gSp2iLFktDvuyLcp3ZZsidRzTvSL9qHVkifdcJZT4Xo9T+Blg4KNsWLDM9AMA2wXGm7gMYOwLwShwmiPI8GKCVZUuiHZEj7PudEvCpqi+3qZSAOPO8K+EIRKMxHOMjKs31WdQrWcC5UyoCluXdfLRLvihETouON3DBCEynOW0zpe3jCEO0LQGmc5xQejSZ3OawrjdnAABKHynBWUXgtHE6yYiIh0Caq0QUvMNjPtdUR0Oq8pggHCARCEynOWkXqdShhu4J3lqNIG7XFWUHrtp4mkzkroA0GoPLLALsaS8XBC6YbIL4wqbdAmVzmt13Jutoln/DgrvwZBqAImwuknddkPE6nHMk6gSXodsIh2xQmzye7GdKavQRCqAuOn2CYdXn5GGtEJA5qk13lNkUYe+2MPhCBUBaefjDTq7Y4QVdqgYSaC8ZORRr1dnrJNgrNs5I01jAZBqApMGRXR3R1huJ53VWLiBGiVq0KHw4SRJt6JAcIBEISqQHsc8XAiIeiqHybSgEoZ0DBdDhNGGgUGd4QDIAjVwUQwPr1NogijUga0zFVBRfR1RxiPJJKxJFlgV7ohqoMgVAtGZ7MJJYJtFNAJA9rlyLcRJiLaHVe6IbKJNArOMhLzegdCEKqF00/paWSea43aXBYrY1G6IQBj5yzX1U0h28QzKBkdDIJQLZgySk+Fo5EGVGmD5jnLKT0tiB9pEpylGCAcBIJQLRivQ+iMJeM6Wec30ohrT9A8Z5muJvjijnAoCEK1MFlMVJGdC0aVbog8MF0JdMDpJ3Wz4qgYS8Z6ElQxKmUGgSBUEaaUZJt10TsqEWwTKmVA8+x5NpPZpI96Ga5ZoLwOkxmlMoNAEKqIboKQb49aaYuVRqUMaB5TRumjnJttxgDhkBCEKsKUkmxAF6dcEybtgk449bLQGtss0D6clYNDEKqIbu4II5hBCHrB6GVeE9ssMLgjHAKCUEVsTqvZoocBiUgjj0oZ0Aed7JstEWwQQTgkBKG6MH493BRi50/QDbJQD/tmC50xK2WxUhi2HxyCUF0Yn+aDEDt/gq7oYt9stllgMEA4NAShuuhgmBA7f4LO6GDfbAwQDg9BqC609oMQU+lBZ3SwbzaCcHgIQnWhPY5od1zTC61FMEAI+qKDfbPZAIJwOAhCdTFZTFSxthdaY5sFpx+nHOgHVWyPhxKioNXLUzGajIcSZBEWVxsSglB1NF0vk+DFBCeShTjlQD9MZhPldXBBrZ6VXECgSrC42nAQhKrD+EjtnnJfFafhjAN90fTlKfpFR4QgVB3aR7IBrXaNYkwedIkpJSOaDUIuGKW9DqVboWoIQtWhvQ5OsyuOss0CgwFC0B1Nz2tiA5hEOAIEoeo48m1JUYpHNLmSBduEO0LQIaaU5AICISndjjHhglEstz08BKEa0V5Si4WjUlLiW6O0F6cc6I2Vslhpi9ARU7ohoxYLJwhJsruw0tNwEIRqxHgdWtyPiW+L2fOsFge+VKBDGh0m5ALYfWlk+JulRrRPk3eEqJQBHWP8FKfFIAxGMUA4IgShGjE+TdbLsM0CU4o1ZUCfNFovwwYElIyOCEGoRrSPZDU4Mo87QtAxpzaDEF2j6UAQqpGVslhIi+Z26GWbeAQh6BU5zh5nEwleVLohoyERXEuULsEd4QgQhCrF+DRWL5NgxWRccuTblG4IQHaYCNpHamvMQuiMWWnsxzsyBKFKae6UY1M9MFhcDfSL8ZGspqrYMJU+TQhClWK0NpWQDQiMDz0woGeMz6GtwlEuiHm9aUEQqhStta5RnHKge7SX1NpZKdC4PE0DglClaI+Db4tJSc1UjrLNqJQBnWNKSS4Q1VA5NxeIMrg8TQOCUKXMdrPdbdXMkk6p4jRMVwJds9IWs92klXJuKSnxHTHKg81BR4YgVC/a69DKMKHQGbNSKE4D/WNKNdM7KrTHHHlWsw1/5EeGz0i9NLT0NorTwCA0VM7NYtg+bQhC9WK8Dq1sVc8FsM8LGAKjnX2zOSyuljYEoXpp7Y4Qpxzon4buCFHInT4EoXpRJQ6+QxuFo1jPEAyCTp2VohbOyiDuCNOFIFQvs9XkyLPybWovHE0mJKErTntwyoH+ma0mR76Nb1V7V40kSkJnnCrGWZkWBKGqaaJ3lG+JkoV2kwWrq4EhaGKYkG+LOQpsZivOyrQgCFWN1kK9DBvEACEYCO3TwFnJBQUG/aJpQxCqmibuCFEyCobCaGGhNTaIs3IUEISqpok7QozJg6HQPg2sdIGzclQQhKpGexxCZ1zlJWqYTQ+GQo6zx8MJMZpUuiHDwdyJUUEQqprJYiILVV2iJkaTCU4kC7GeIRiFyWyiPA6uRb1nZTIhRbviVBHOynQhCNWOLlH1KccFBKrEgf14wVBor0PN0+r51ig5zoZC7vQhCNWO8pJqDkI2iH1ewHBUXsXGtUTpEpyVo4AgVDu6RNUj8xiTBwNifA6uRb13hNgTbbQQhGqn8sJRjMmDAdFeVc+px+XpaCEI1U7lhaMcltsG43Hk25KxpMirtHCUC6JrdHQQhGpnspgcBTa+XY0rjia4ZDIh2fNsSjcEILdMBF3iEFrVeFZKohTtilPFKBkdBQShBqi2d1RoiWH1CjAm2kfyrXGlWzEIrjVKFqJkdHQQhBqg2nqZaFsc6xmCMdFeR7RNjUHIo1Jm9BCEGkCXqHQGBd8aR6UMGBPtJYUWNQYh6tfGAEGoAaruGsW1JxgS43Pwqhwj5FoEugRn5ejIFoQ8z3/22Wf19fVDHdDe3n7o0KGGhga53tE4KI9D6FBj4ajQFkcQgjHZnFaTyRQLJZRuSH9cMIogHC15gvDQoUMTJ0585JFH5s2b98QTTww84Pbbb58yZcratWtnz5594403xmJqvJJSra82xVZZ4WisJ262mGxOq9INAVAG6bGprasmtTE9iY3pR0meIHzssccefvjhf/7zn4cPH37ppZe+/PLLfgesXr06GAx+8skn586dO378+Pbt22V5X+Ogvaqrl2EDUbIEJdpgXI5im9rOSr41ShZiY/pRkyEIW1paPvroo/vuu48giNLS0mXLlr3xxhv9jrnmmmvsdjtBEC6Xa8aMGS0tLZm/r6GocJiQaxEcRbgdBOOiSmxq26GXa0G/6FjI8IesoaHB5XIVFRWlfqyqqhpmpPDMmTMfffTRT37yk6EOEEVx7969va928cUXezyezBupdXQJ2XE8pHQrvoYLREkv7gjBuEiPveeYys5KlIyOSbpB+NRTTx07dqzfk3Pnzl2/fj3Hcam7vRSSJFmWHfRFOjs7b7rppscee6y6unqoN4rFYr/5zW8cjq8uah5//PHLLrtsqIN5nrfb7RaLJc3/Cg1zi5FmPhKJKN2O/wo3c4WTSVU1yWhYljWZ0AmmmKQzzjbzkXBEPduQhZoi+TMYg5yVaX7/aZo2m0fo+0w3CBcuXDhx4sR+T/r9foIgSkpKuru7k8lk6s06Ojq8Xu/AV+jp6Vm+fPnVV1+9YcOGYd6Ioqg33ngj9cojslgsBglCerx0uivIUIxaFoyQiGhbPL/S43Q6lW6KcUmShM9fQZIkWamQLeFwFKhllcFYe6CwMo92GuKmUMbvf7pBeMUVVwz1T1VVVfn5+QcOHEjduv3rX/969NFH+x3Dsuw3vvGN6urqZ599dsxtNbLewlGVDAAInTErbbE4MA8VDI3xObigoJIgRMnomMnwh8xut69du3bdunV/+9vfnnjiiba2tptvvpkgiH379lVUVKSOufnmm+vr6+fMmbNly5aXX375448/zvx9jUZVhaMYigAgCIIuUdF+TCgZHTN5qv42btw4bty4F154we/3f/TRRyRJEgRRUlJy2223pQ5YtGjRzJkza2trUz/21sJA+uiS1F6gbqUbQhAEwQWx+xIAQXsdPWcHL4nIPZSMjplJktS1XklZWdknn3yS5hihgYplCKLt856O46Fpd5cr3RCCIIhT/6exYKqTmmpxuVxKt8W4wuEwPn8FhcNhost69s3m6vX96ycUUf+3VoIgKpYbpcxexu8/xng0Q11dowFsgQ1A0F4H1xqVkqq4neBacFaOEYJQMyiPQ+iIqWHFUSkp8W0xyoNTDozOYjfbndZopyq2ocAqo2OGINQMs1UtW9UL7TF7ntVix5cHgKB9DjWsL5MqGcXl6djgb5mWqGSHXjYQZbAxPQBBEARBe0k1nJU8NqbPAIJQS/5TOKowLijQCEIAgiAIgvGq4o4QJaOZQBBqiUquPbmAwGBMHoAgCIKgfaQaFsTH1N5MIAi1RCV7ULA45QD+I7VvdjKhcBUbSkYzgSDUEjVsVZ9MSNGuOFWMfScACIIgzFYTWWjj2xTuqmED6BodOwShlnxVONqmZOEoFxSoYjvG5AF60T6SU3SY8KvLU5SMjhWCUGNSi/wq2AAMRQD0wyi92AXfGiXH4fJ07BCEGkN7SVbRUw6rjAL0Q/tIZQtHsdJThhCEGkN7Hcp2wrAB3BECfA3tJTlF96Bgg1EUcmcCQagxtFfhWm0uINC4IwTogyqyx9mEGE0q1QBM7c0QglBjqGJ7tCeRjCtzyiUEMcGLZAFKRgH6MBFUsZLDhBi5zxCCUGNMZhNVZOdalDnluECU9joIDMkDfB2tXBWbGEvGwglynE2Rd9cHBKH2KLi+DBcUsMoowECMV7F6GS4YpT0OkxnXp2OHINQexqdYvQwqZQAGRfsUq5fhgigZzRSCUHtoL6lU1yjbLDClCEKA/phSkm3mFXlrDBBmDkGoPbRSq91LmK4EMDi720qYTLFQIvdvjam9mUMQag9ZaE9wYkIQc/y+0e642WayOa05fl8ATVBqzAIDFplDEGqQiaBKHHzO62XYAPpFAYZE+0i2OddBmOBFMSo68lEymhEEoSYxSizphAFCgGEodlb6SMxoyhCCUJMUOeW4AFavABgSU6pcEEJmEISaxJSSbHPOu0ZxygEMjfY6+LZYjrcL5QICjX6ajCEINYkpJblmgcjhGZdMSEJXHDt/AgzFbDM78nO9XSgbwOWpDBCEmmSlLWa7Kdodz9k7ci1RChueAQyLyXG9jJSaRIjL00whCLWKKc3pKcc1owcGYASML6dzfIWOmI2xWilLzt5RrxCEWpXjkXk2gEm7ACOgS8lcTiVkm7EnmjwQhFqV40lLqJQBGFGOu0YxtVcuCEKtYnw5vfbE3AmAEZGF9oQgJrgcrfrEoVJGJghCraJLHEJnPDc79MZCiWRSwuoVACMw5fSmMNKMy1N5IAi1ymQxUcU52qGXbeKdfioHbwSgdc4yKtKUi20oxFgyHkpQxfYcvJfuIQg1LGfry0SaBKcfF54AI2NKSbYpF2clF4hS2I9XJghCDaN9JJeTThi2SWAQhABpYPxkpDEnZyUGCOWDINQwp5+M5OTaM9LEM+gaBUgD4yOFzlgykfVln9gmHpenckEQahhTRrFNWV9oTRQwFAGQLpPFRBXZc1DRHUE/jXwQhBpmYywWh1noyu7ahmwzT/tIDEUApMnpp7LeVSMRXABBKBsEobYxfpLN8oAELjwBRoXxk2yWC0e51qjNZbWSWFxNHghCbXP6s16rzTajZBRgFBh/1gtHMaNJXghCbcvFKdeIShmAUXD6KTaQ3cH7SJPgLMPlqWwQhNrmLKMijVm8I5REiWuLMdjnBSBtFtJsc1r5tiwudoHLU3khCLXNUWBLilIsnMjS63MtUbLAZrbjewIwCtme2hRpxsi9nPAHTvOcpRSbtZv1HK2+AAASgUlEQVTCSAPvLMeFJ8DoOMuz2FUjdMbMFpPdZc3S6xsQglDznGVZvPYM1/OuCgQhwOg4y6lIfbaCkG0SnGU4K+WEINQ8xk9lr1Y70sA7y+ksvTiAXjkrqEgTLyWzUjATaRIYBKGsEISa5yzL1tqGyYTEt0aZUlTKAIyOlbTYXVa+LSuLXbCNPGY0yQtBqHlUsSPOJhK8/HuBss0CVWw32/AlARg1Zzmdpd7RSCOPrlF54W+c9pkIZ1lWBiQi9ZyzAv2iAGPhqqDCDZzsLxvtjhMS4SjALtlyQhDqgauSDtfJf8qhUgZgzLJULxOu412VuDyVGYJQD1wVVDgbd4SYOwEwVk4/yQajsu/HFKnnnLg8lRuCUA9clbTsQShGk9HuOI01ZQDGxGw3U0V2LihzIVu4nndVIghlhiDUA7vbaraahE45S9QiDTxTit2XAMZO9t5RKSmhUiYbEIQ64aygwnVynnJh9IsCZMZVQYUb5DwruWDUkW+zUth9SWYIQp1wVdCRejnrZSL1HCplADLhrKDD5+U8K8P1PAYIswFBqBMuue8IQ7Wcu4qR8QUBjIbxOWI9iQQr2xzfcB2HktFsQBDqhLOCYgOCJMpTosa3xQiTCXOVADJhMpuclVRIvpvCCGY0ZQeCUCcsdjNZaGcD8pSohWrZvIm48ATIlLuKDtXKE4RiNCl0xBgfFleTH4JQP1yVlFwDEuHznHs8+kUBMuWuYkLnWFleKlzPM37SZEEht/wQhPrhnsj0nJUnCHvOce4JuCMEyJS7kmIDgizT6kNn2byJuDzNCgShfuRNZHrOskTGZ1ycFeOhBKbSA2TObDfTJQ5ZKrp7EIRZgyDUD0e+zWI3c63RDF8ndI51VdGYSg8gC3cV03Mu0yBMJqRII+8aj36arEAQ6kreRCZ0NtMBiVAt567C+QYgD1nqZcJ1HO11WBz4i50V+Fh1xT2J7pEhCFkEIYBc3BPocC2X4W716BfNKgShrnw1TJgBMZrkglHMVQKQi81ptbmtXCCjMYvQWQ5BmD0IQl0hC+0mi4lvG/vq2z1nWFcljV3pAWRUMNXZdSoy5l+XRClcz7nQT5M1+HunNxneFHadChdMc8rYHgAomObsPhUe86+H63nK47CSWGs7WxCEepNpEJ6M5E9FEALIyT2RCTfwYjQ5tl/vOcvmTUC/aBYhCPUmf6qz62R4bCPzQmcsGUsyXqzhBCAni93sKh97IVvXl2FcnmaVPEGYTCb379//zjvvdHd3D3NYNBo9dOhQR0eHLG8Kg3Lk2xx5trHtRNF1IlIwzUVgAiGA3PKnObtOjmWYMMGKXDCaNwl3hFkkQxCKonj99dfff//9L7/88rRp044dOzbUkRs3bpw/f/7u3bszf1MYRuF0V+eXYxmQ6DqJAUKArCiYOsZhws4T4bwpjNmK69MskiEId+/efebMmYMHD/71r3+9//77f/zjHw962MGDB/fu3Ttr1qzM3xGGVzjD3fnv0Gh/SxKl0DkOF54A2cD4yNT2EaP9xc5/hwpnuLLRJOglQxDu3LnzxhtvpCiKIIhVq1bt2rUrkUj0OyYWiz3wwAMvvPCCxYLCp6xzVVBxVhQ6R3fKhWo5qthuc1qz1CoAQzMRBVNdo+0dlUSpu4YtnIYgzC4Z/uo1NjbOnTs39biioiKRSASDwbKysr7H/PznP1+xYkV1dfWIr5ZIJN5+++3CwsLUj5deemm/l+pLFEVRlG33Zz0pmOZsP9bju7ww/V9p+6I7f7pzVJ8nPn9l4fNX1mg///zpTOCfnZ5L89L/lZ7TLFlsN9Mm/I8eKM3P32w2m0wjdCynFYRHjhxZv379wOe3b99eWVkZi8Vstq+2Mk89iEa/tobC0aNHd+7c+emnn6bzXqIo7tq1K3V/SRBESUlJcXHxUAdHo1FJknCXOZBzsqPtk1DhvHT7OaWk1P5FaOqa0n7/74YXi8VGdTzIC5+/skb7+dMTbNwOIdzK2vPSvQNpP9bjnkrh//Kg0vz8SZKUJwirqqp+8pOfDHw+FVFer7e9vT31TFtbG0EQPp+v72G/+MUv/H5/6hUaGhp27Njhcrm++c1vDvpeDodjy5Ytfr8/nYaZTCa73Y4gHMhxEVn/Vrvd5LBSaX04XSfCVJGjwD+Ka1WCIERRpGmsdqEYfP7KGsPnP+6iPPZkLP8qd1pHS0ToVOMF91XQNGY0DULG739aQZiXl3fFFVcM9a8LFy58++23N2zYQBDEhx9+OGvWrH6NW716dV1dXeqxw+EoLi4uKirKoM0wMovDnD/N2fZ5j29hWr2jbZ/3eGaPLgUBYLSKZ+ef3xX0X5XWH8Ces6zFYWZ8SMGsk2GM8O67737mmWfWr18/ffr0H/3oR5s3b049v3DhwhtvvPGxxx5btmxZ78Fbt2698sorh4lVkIv30oLat4PpBGEyluz8Mlz1DW8OWgVgZPmTmFg4wbVE6ZKRN74OHugqubQgB60CGapGCwoKDhw4QJLk4cOHt2/fftttt6Wef/DBBwcG3kMPPTRnzpzM3xRGlD/ZKcaSkYaRZ9Z3HA+7KimbC/WiAFlmIoqr3e2He0Y8MMGKXSfCnrn5OWgUyPO3r6Ki4pe//GW/J++8886BR959992yvCOMzER4LykIHuiaVD7Cnkqtn3UVz8b5BpALxXPyT77SUH5Nsck8XAVH66HuwhmuNMf4IUNYa1TPPPML2o/0DL/Ub6SR51qiRdUYIATIBWcZ5Si0tX0+wk1h8ECnF/2iuYIg1DO7y5o3iWn9bLgFYOv/1lq2pBgLOAHkTMUyT8PfW4dZGb/nDCslCXcVlnnKEQShzlUs99S/1xpnB592Gmnk2YBQMh8XngC5kzeRsecNeVMoJaVzfw1ULvdg+fucQRDqHOMji2bm1f2/lkH/FbeDAIqovLak/r3BbwoD/+y0MVaMVuQSglD/Kq/1dB4PRer7l4+2HOziWqO4HQTIPXcVTRXZ697pf4UaCyca3m+b+E3foL8FWYIg1D8rZRl/vbfmz419V77v/DJc939bZvzPStwOAihi6p3lnV+Gmz5q730mwYo1f2wouaSA8ow8yxBkhKljhuCZm5/gxS+eOzdhpddRYI808A3vt03/TiVVjPMNQBlW2jJjzfijvzmX4JN5E2lCIk6/3lQ0M69yuUfpphkOgtAoSv/HOHcVffq1JrPVxPipC1ZXuCpGmF8IAFnlyLdduGZ84F+dDe+1RXviE28qLZyOHZcUgCA0EGcZNeuxSUq3AgD+i/I4JtyIEUGFYYwQAAAMDUEIAACGhiAEAABD03YQtre3h8NhpVthXIFAQBAEpVthXA0NDaI4+JpBkG3JZLK+vl7pVhhXNBptbm6W69W0HYQbNmx46623lG6Fcd1zzz2HDh1SuhXGtWjRos7OTqVbYVA9PT2XX3650q0wriNHjgy6wdHYaDsIAQAAMoQgBAAAQ0MQAgCAoZkkacg9sRRxww03HD161GxOK6HD4bDNZiNJMtutgkF1d3czDGOz2ZRuiEF1dHQUFBSkebKAvCRJ6ujoKCoqUrohBpVIJMLhcEHByHsG7N69+4ILLhj+GNUFYTgcbmtrU7oVAACgB2VlZXa7ffhjVBeEAAAAuYROFQAAMDQEIQAAGBqCEAAADA1BCAAAhqaZ/Qjb2tq2bt3a1tZ23XXXLV68eOABoii++uqrx44dmzZt2r333ouafhm1tbXt2rXrxIkT+fn5t9566+TJk/sdkEwmf//73/f+eOGFF1522WW5baOe7dmzp6amJvXYarWuXr164DHNzc3btm3r7u5euXIllv6S15YtW/oWFQ78etfU1OzZs6f3x5UrV3o82GU+I62trYcOHWpoaLjyyiunTJnS+3xjY+P27dtDodBNN9106aWXDvzFWCy2devW06dPV1dX33nnnWlOLtLGHSHHcQsWLDh16lRlZeWqVatef/31gcc89NBDL7zwwuTJk//4xz9++9vfznkb9WzdunXvvvuuz+drbW2trq7et29fvwMSicSaNWtOnjx57ty5c+fOdXR0KNJOvXrllVdef/311GdbW1s78IDu7u758+c3NTWVlZXdcMMN77zzTu4bqWO1tbXn/uPRRx89duxYvwP279+/adOm3mOi0agi7dSTJUuWPPnkk0888cT+/ft7n2xvb583b15bW5vP57v22mv//ve/D/zFO+644/XXX588efJzzz23fv36dN9P0oJt27bNmzcvmUxKkvSnP/3p4osv7ndAIBBwOBwNDQ2SJHV2dpIkefbsWQUaqlM8z/c+XrNmzerVq/sdkDrzw+FwbttlFPfcc8+vf/3rYQ7YvHnz4sWLU49feumlhQsX5qRdhvPZZ59RFNXV1dXv+VdeeWXFihWKNEmvRFGUJOmSSy555ZVXep985plnrr322tTj5557rvc73+vEiRM0Tff09EiSVFtbS1FUe3t7Om+njTvCjz/+eOnSpSaTiSCIpUuXHj16tKurq+8B+/fvnzRpUllZGUEQBQUFs2fP3rt3rzJt1aO+a/cIguB0Ogc97He/+93zzz//+eef56pdBnLgwIFNmzbt2LEjHo8P/NePP/74mmuuST1eunTp/v37Bz0MMrR169ZbbrklPz9/4D81NjZu2rRp27ZtWA9EFoN2afb7nu/duzeZTPY7YP78+W63myCI8ePHl5eXf/LJJ2m9XcYNzoVAIFBcXJx6PG7cOKvVGggEhjqAIIiSkhIZt6qCXvv379+5c+d3v/vdgf+0ZMmSjo6OEydOLF68eNOmTblvm45VVFQUFxd3dnY+9dRT8+fP5ziu3wF9v/8ejyeZTAaDwZw3U+d4nn/ttdcGHaDNy8u78MILe3p6du7cOW3atIF9pyCLft/zeDze3t7e94BgMNg3CDweT5pBoI1iGavVmkgkUo9FURRFsd+SOVarte8OpfF4fMQ1dWC0Tp06dcstt2zZsmXSpEn9/slut7///vupx3fcccfVV1+9du1ahmFy3kZ9evLJJ3sfzJ49e9u2bevWret7QN8TJPUA33/Z/eUvfykoKLjiiisG/tPKlStXrlyZevzggw/+9Kc/ffPNN3PbOkMY8Xs+5iDQxh2h3+/vDfbUA5/P1++Apqam3h+bmppKS0tz2ULdq6mpWbJkydNPP33rrbcOf+SCBQsSiURjY2NuGmYoNptt/vz5A+tl+p4gTU1NNpsNi0HLbtu2bffdd19qgGYYl1122blz53LTJKPp9z2nabpfN/WYg0AbQbhixYpdu3YJgkAQxJtvvrl48eLU3caRI0fq6uoIgli0aFF7e3tqt/TTp0+fPHmytysZMldXV7ds2bKNGzf22xL6008/TX3teJ7vfXL37t00TY8fPz7HjdSx3o83HA7v2bNnxowZBEHEYrEPPvggVaa0YsWKnTt3psYFd+zYcd1111ksFgUbrD+1tbV79+696667ep/hOO6DDz5I3Zf0/g+SJOmdd9658MILlWml3q1YseKtt95K3fPt2LFjxYoVqec//fTTVEAuW7bsiy++SF0pHjhwgOO4hQsXpvXSslT4ZFsikbjmmmvmzJlz1113jRs3bt++fannFy9e/LOf/Sz1ePPmzT6fb/Xq1eXl5b1PgiyWL1/OMMyc/3jggQdSz8+aNeu3v/2tJElbtmyZMWPGt771reXLl7vd7j/84Q+KtldvCgsLV6xYsWrVKp/Pd/3118fjcUmSUmf++fPnJUkSBOHyyy9fsGDBqlWrioqKDh8+rHST9WbDhg3XXXdd32dOnjxJEERHR4ckScuXL1+yZMmdd9550UUXTZ48ua6uTqFm6sejjz46Z84chmHGjx8/Z86c1N98juMuueSSyy+//Pbbb/d4PMePH08dfPHFF7/00kupxz/84Q8rKytXr17t9XpffPHFNN9OM7tPiKK4Z8+e9vb2RYsWeb3e1JOnTp1yuVy9N7/Hjx9PTaifNWuWci3VoZqamnA43Pujy+VKTXH997//XVxcnBq1PnToUG1tbV5e3ty5czGbWF7nz58/cuRILBabMmVKdXV16sl4PH748OHq6urUKEg8Hv/www+7u7uvuuqqvvUCIIuamhq32937l4cgCEEQjh49OmfOHIvF0tnZefDgwa6uLr/fv2DBAqzmkbmzZ892d3f3/jh58uRULWiqIyQUCi1ZsmTcuHGpfz1+/HhJSUnv1/7QoUM1NTUzZ86cPn16mm+nmSAEAADIBm2MEQIAAGQJghAAAAwNQQgAAIaGIAQAAENDEAIAgKEhCAEAwNAQhAAAYGgIQgAAMDQEIQAAGBqCEAAADA1BCAAAhvb/AZaMSmWLViOEAAAAAElFTkSuQmCC",
"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.0"
},
"kernelspec": {
"name": "julia-1.7",
"display_name": "Julia 1.7.0",
"language": "julia"
}
},
"nbformat": 4
}