{ "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.143558520868 NaN 3.81e-01 0.80 8.0\n", " 2 -0.156032524694 -1.25e-02 7.92e-02 0.80 1.0\n", " 3 -0.156768415148 -7.36e-04 2.74e-02 0.80 2.0\n", " 4 -0.157045816374 -2.77e-04 4.82e-03 0.80 2.0\n", " 5 -0.157055207782 -9.39e-06 1.22e-03 0.80 2.0\n", " 6 -0.157056386565 -1.18e-06 1.97e-04 0.80 1.0\n", " 7 -0.157056406524 -2.00e-08 3.10e-05 0.80 2.0\n", " 8 -0.157056406917 -3.92e-10 1.57e-06 0.80 2.0\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown:\n Kinetic 0.0380295 \n AtomicLocal -0.3163466\n PowerNonlinearity 0.1212607 \n\n total -0.157056406917\n" }, "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.05568295449851397, 0.0, 0.0]\n [0.05568183401379799, 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+gvaeTAAAgAElEQVR4nOzdd3iT1f4A8POO7KQZ3XvvltKWDQXKLhZQ5kUBReTCD1FBVBAcKHpdlytXFAcguHAgiELZIBRkSAttaaF77zZ7J+/4/RFuRUbpSJM0OZ/HxycNJ29OTt6c73vOewZC0zSAIAiCIFeF2jsDEARBEGRPMBBCEARBLg0GQgiCIMilwUAIQRAEuTQYCCEIgiCXBgMhBEEQ5NJgIIQgCIJcGgyEEARBkEuDgRCCIAhyaTAQQhAEQS7N4QJhQUGB2WzuYmKKouAScXZEkqS9s+DSYPnbFyx/+7Ji+TtcIJw6dWpra2sXExuNRoqi+jQ/UCd0Op29s+DSYPnbFyx/+7Ji+TtcIIQgCIIgW4KBEIIgCHJpMBBCEARBLg0GQgiCIMilwUAIQRAEuTQYCCEIgiCXBgMhBEEQ5NJgIIQgCIJcGm7vDEB9y0yBEiVdqqSr1ICggIkCQibwYoNIITJAgjDghRAE2Q9Jg2IFXaKkG3VAbgRMFCAICOKBKCESI0K4sHq2FVjSzqlZD36ooI43UH80075cJEaEhAkAEwNMFFSowIUWUKygylX0QHfk4RB0dggSIkDsnWUIchVtBvBLNbWvirrUSvtykXgx4scFYhbQEoCkQE4bKFVSFWo62R2Z6I/+IwyJFMKfZ9+CgdDZnGigPyoi/2ihHwlGn4pGvxuLiln3TqkjwIUWel81NfQ3coAEWZ2AZQQi8AcHQXf44IMPrly5cvfzJEliGNatQylNoFZDS43Agw282GA8C8FRAABoAaDltmRRAITTQGEE+4z0J3rAxUEAD3hzXP3nOXbs2BUrVvTFkRFHW7Q6ICDg8uXL/v7+XUms1+uZTGZ3z0VndaSO3niV1JjBiwPQOWEor8sXOWYK/FRJ/aeQMlPg/SHYlICu/tzUarVAIOhhdqFeg+VvG2lpaZMmTYqJibF3Rlza1atXCwsLDx482PGMFc9/2CJ0BsUKevUlsloD3h6EPhyMot28bmSg4LEI9LEI9NcaavUl8kMe+GQkFuHm4lefEPSX9PT0UaNG2TsXLo3D4RQWFvbRweFgif7NTIE3r1GjDxGTA9CCmfjMkG5HwdvNCEavz8QzAtHhvxGbr1OkY3UWQBAE9QkYCPuxQjk95Fficit17RF8VQJqlSGgOApWJaCXZ+BZtdS4LKJBC4MhBEFOzjqB8Ntvv42MjPTy8lqyZIler79fssOHDw8aNOjQoUNWeVMXt6OEGpdFPBuPZk3G/XlW7sYMEyAnp+KTAtBBB4hj9TAWQhDkzKwQCG/evPn000/v3r27pKSkpqZm06ZN90ymVCrXrl0rk8mkUmnv39SV6Qiw4Ay5tYjKzsQXR/VVmx5FwIaB6I/j8SXnyM3X4e7HEAQ5LStUo19++eXDDz88cuRIsVj8yiuv7Ny5854jUdesWbNq1SpPT8/ev6Mrq9XQaYcIDAGXpuMxoj4fzzLaB7k4HfuunHoymzTBaAhBkDOyQiAsLi5OSkqyPE5KSmptbZXL5XekOX36dEVFxZNPPtn7t3NlF1roYb8RCyLQr8ZgHFsN+A3kIeem4XIjyDxGqM02elMIgiCbsUJtKpPJOiZzuLm5AQDa29slEklHAq1W++yzzx44cADpwnxQrVYbEBBwx5NVVVUeHh53J3apeYQH6tDVOfgXw4iJvkaNxtbvvnsYeD4HH3OQ2jfa7Mm+1eLXarVd+U6hPgLL3zYoykk6Q1pbW1NSUurr6y1/bt++vamp6bXXXrtf+qNHj+7fv/+LL76wVQYfgCRJzW11XxfPfy6Xi6IPaPJZIRBKJBKVSmV5bHng7u5+e4KNGzcOHjxYqVTm5uZqtdrq6uqamprg4OB7Ho3H4xUWFnZxQj2GYS4SCD8spD68Tp18CEuSMO2Vh53pYONVcuoZ9GQGZhmeQ9M0n8+3V34gWP628cBqtL+gKKq9vd3yWK/Xv/nmm7m5uZ2knzx58osvvlhUVBQfH2+TDD4AhmG3n/BWPP+tEAgjIyOLioosjwsLC93d3W9vDgIAKIq6fv36smXLAADV1dW7d+/GMOyVV17p/Vu7AhqA9VfIg7X0helYgLVHh3bXxhSMz6DGZJGnpmLBfNgWgSA7qKmpyc/PDw8P37t37+jRo8eNG1dcXHzo0CGz2Tx9+nRL0NJqtVlZWUVFRRwOZ9q0aXdHsn379qWmpnp5eQEATp8+7eXllZCQ8NVXX82cOVOr1WZnZ8+dOxdBkMcee+zTTz/9+OOP7fA5bcgKVzqLFy/ev39/Xl6eTqd79913Fy9ebGmubtq06ciRIwCAzZs35/xPfHz8xo0bYRTsIpIGS8+RZ5vo7Ezc7lHQ4oVEdFUCOjaLrFbDaRUQZAcFBQUrVqxYuXKlWCxGUTQrKysjIwNBEC6Xm5GRkZ2dDQC4evVqdnZ2UFAQTdPjx4/Pycm54yBZWVnp6emWxzt37jx9+jQAYPXq1VKptKysbOPGjZZ/Sk9Pz8rKst1nsxMrtAgTExM/+OCDzMxMjUaTkZHx+uuvW54vKSkJDAy8I7G/vz9cHbGLTBRY8DupNNEnpuJdXzjUBlbGoSgA4w6TWWORWPhlQq7naju99gppm/ca5oVsSr3z7o9KpTpw4IBQKAQAREVF7dy5c9y4cQAAb2/vd955Z/To0WlpaWlpaZbELBbryy+/HDRo0O1HuH79+qJFix747lFRUdXV1SqVyjL+w1lZp35dunTp0qVL73jy22+/vTvl/v37rfKOTk9PgNmnCCaK/DYJZznePdAVcShJg4d+Z2ZPox2kqQpBNhMiQNYOsNHP0oN9jyejo6MtUVClUpWVlW3atOmdd94BAKjV6tbWVgBAU1PT8uXLb968yWAw9Hr93V2jarWay+U+8N0tN+HUajUMhJCtaQkw/Tjhy0V2j8ZwR71P/0w8qtaTE4+gZx/CvTj2zg0E2ZCEBSb42/P6j8W6tbmaZbTg559/3jFF2zK055VXXomJibGM1f/kk08OHz58xxE8PT075rmx2WyDwdDxTwaDgc2+FX6lUimKoneMf3Q+jlrLujAtAaYdI0IFyNdjHDcKWjwTTcwLQyYeIWRGe2cFglwSm81OT0/fs2eP+H9MJhMAoKWlJSwsDEEQk8n0/fff3/3CwYMHFxQUWB4nJSUdP36cIAjLnwcPHhw4cKDlcX5+flJSUkdcdFaOXdG6HpUZTDxMxIiQ7WlYb/aRsJmNKdh4P2T6cUJH2DsrEOSStm/ffuTIkZSUlEceeSQhIeHdd98FACxfvnz9+vWPPPLIoEGD7jkbbc6cOZbBjACAZcuWCYXCqKgojUYzYcKE3NzcjpUyjx49Onv2bJt9FnuBXaMORGkCU44SSRLkk5FYfwiCt2wehi09R844QWRNxpnwygqC+tiUKVPGjh3b8WdISMilS5eqq6tbWlpCQ0O9vb0BAJmZmTdv3qyoqIiKihIIBJZmopeXV11dneVV6enpRqOxoKBgwIABLBZr3759ra2tkZGR33777bBhwyxpdDrdwYMHL1y4YOtPaHOw3nIUciOYeIQY5oV8Oqo/RUEAAALAZyMxDoYsyYY7GEJQn2MwGHeMvUcQJDQ0dNiwYZYoaOHl5TV8+HB3d3cmk2kZ84KiaMetRARBtm7devny5dvTYxjm4+PT8cyVK1fWr19vmWvo3GAgdAiWKJjmg3w4rJ9FQQscBT+Ow6o19IuXbTSmHIKgXho5cuQdo/2fe+45y2BUizFjxrjIAtEwENqfzAgmHCHSfZHNQx1vnkSXcXDw20T8aD0N92yCoH7q9ddfF4vF9s6FHcBAaGdSI5hwmJjgh3zQn6OghZgFjk7BPiqi9lTAWAg5D0cbCHbkyJGJEyf20cEPHTq0YsWKThIUFxdPmjTpnnvt9V8wENqT3AimHCEm+CPvDen3UdAigIccmYI9f4k83uBUvxPIZZWr6FqNY53MCQkJq1at6osj0zS9du3alStXdpImJiaGzWYfPHiwLzJgLzAQ2k2rHozNIiYHIO87SxS0iBMhP43HF54h8mWOVX1AUHc16cCkI6QPx7Fu3BMEodPpAAAmk+m9995raGjYuHHjunXrampqdDrdhx9++OKLL169etWSuL29ffv27atXr960aVNpaWnHQWpra998880NGzaUl5dv2bJFqVQCAM6ePcvhcOLi4gAA2dnZlgVId+/eXVtbq1arN2/ebHntwoULP/30Uxt/6j4FA6F9NOtBehbxSAjy1iCnioIWo32QT0ZgmcdIR7uUhqCuU5vBQ8eIJdGoiGXvrPxdYWHhRx99BAAwGo3r1q1bsmRJQECAUqmcOHHi4sWLAQASiSQ9Pb2trQ0AcPHixfr6+iFDhjCZzFGjRpWXlwMA2trahg4dajAYYmJili9f/vLLL8tkMgDA4cOHOyZmHDt2zLIezUcffVRZWalQKF599VXLP40dO/b333/X6/V2+PB9A84jtIMGLT3+MLkoEl0/0GkvRGaHog1akHGUPD8NFztYPQJBD2SmwJxTxBBPZMNA9Ohd/2qqKVb8/IltcsIMSxA9sqyTBG+99dagQYOefPJJDw+PCRMmWAaCHj9+/MyZM3PmzJk2bdq0adNomlYoFPX19d9///2rr766Y8eOMWPG/Otf/wIApKSkJCQkWA5VWFg4ffr0B2bJ09OTw+GUl5cnJiZa4yPaHwyEtlatpiccIZfHoi8kOm0UtHguAa3T0g+fII5l4GwnbPdCTosGYOk5kokin4y894mLeweJ5j5rm8ygHF7nCSxhDEVRDw+PjpDm5eUllUoBAEVFRU899ZRMJhMIBM3NzZmZmQCA0tLS5ORkS8q4uLiOlUt1Oh2H06WFg7lcrlqt7tEHckQwENpUmZKeeIR8cQD6dJyTR0GL94dgj50hF54hfxzXP1aMgyAAwKs5ZLGSPj0Vv9+sXpTNZQZG2jZT94Vhf0Vry4rbAAAEQSwDO5977rknnnjCsi/6yy+/bImOIpGoY8VtrVZrWXcGAODp6WnpIwUAcDgchULRcWS9Xt8RIymKkslkt0/e7+9cojp2ENdldPph8rUUV4mCAAAUAbtHY+0G+nk40R7qJz67Sf1URR+chHOdopnQEbEUCsVPP/1kefKhhx7as2dPY2MjAKBjCAwAYOjQofn5+ZbHKSkpJ0+e7Gj2/fzzz6mpqZbHpaWlAoEgLCzMZp+ir7lKjWx3l1vpiUeI/wxFn4xyrTJnYeCXifjpRvrfcKI95PB+q6E2XaOOTsE8nWW7hTVr1jz11FMZGRlpaWkpKSmWJydMmLBs2bLk5OSQkBCKojgcjmVvwrlz53ZsQzF16tTp06dHR0eXlJQsWLDgwIEDW7dutbw8Kytrzpw5COJEnTy0g/H396+vr+9iYp1ORxBEn+bHKo7XU17fmg7XUvbOiJWpVKoupqzXUMHfm78rJ/s0P66m6+UPdcUfzZTXt6actjt/p6NGjTp37pxdsnRPBEHo9Xqapi1dlB3PKxQKs9lseazRaAwGg+WxVCrNz8/X6/V6vV6r1XakpyiKoqiioiKhUEhRtz71nDlz9u3b15FGp9PFxsb+8ssvHc+QJJmYmFhcXNxnn+/eDh48mJmZefszVjz/naLx79h+qaZW/EHuG4+P8nGiC6hu8uchh6dg47IIdxYyOcB1ywFyWOUqes4p8usxeKqHo5+fGIZZ7gsiCHL7imi3LxPK4/01xEYikUgkkjsOsnTp0oSEBLPZ/Nlnn23YsKGjeffuu+/ePlmew+Gw2WyRSNTxTHl5+fLly6Ojo636mewMBsK+9d9CavN16sRUPEHs6L+uvhYnQvZNwGeeJLIm44Mcvq6BXEq9lp5wmHxvCOo6V2mPPfZYbm4ugiBff/31iBEjOp4PCwt77rnnbk/5xBNPBAUFdfwZFRUVFRVlu4zaBAyEfYUG4I2r5N5K+tw0LJjvKr+uzo30Rnam4dOOEWcy8WghLBPIIShNIPMY+Uw8uiDChe7fjx079vZNDTvx7LM2mihiRy70xduSngDzTpG/N9Lnp+EwCt4uMwh5ezCWcZRs1MFFZyD70xHgoWPERH9kjbPP64U6Ab9762vWg7FZBAsDxzPgoir38GQUuiwGnXyElBvtnRXItVmWj4lwQ97v/3u/QL0Bu0at7JqUfuQE+WQ0+moynEF+X2uT0HYD/dAx4sRUnAfPQcgeKBo8cZbEUWRHWv/bDZuiqC+++OLatWuBgYErV668fTAL1AOwRWhNP1VSU44S/x6Kvgaj4IO8PxSLFSEzTxBGONUesodnLpINOvrHcRjeD2vBZcuWHTp0KCMjo6amJj09nSThr6hX4NW4dRAU2JBD7q2iT2TgAyQwCD4YAsAXadg/TpOP/k7208oI6r9evkJeaaNPTu3hKrg32ku3XPnM2pm6twFecStTn7r9mfr6+r1799bV1QkEghkzZsTFxR05csSyiCjUMzAQWkGLHsw/TbAwcOVh3B3eFOwyDAHfpWMPnyCezCZ3j4GLkUI28q886lAtfSYTd2P08AjBwoA1Qzvbxt2K+Iw7F90uKSkJDw8XCAQAAARBUlJSSkpKYCDsDRgIe+tUI/34WXJpNPpqMgqr8u5iomDfeHzqMWL5H+Tno/rfrRqo3/mwkPqqjDrzUK+uWXkMbrQkwnqZ6h6dTqdSqTr+VKlUbm5u9sqMc4AdUj1npsCGHPLxs+TXY7DXU2AU7CEODg5Owm/I6ecuwvscUN/adoP6uIg6NRXz5do7K71TXl5+5swZAEBVVdXZs2fT09PtnaP+DQbCHipW0CN+I/Kl9LVH8HF+MAb2Cp8Bsibjl1vp5y/BWAj1lS+KqfcLqFNTsQBev//BDhw4cNOmTSNGjBg8ePDbb78dEWG35qlzgF2j3UbS4KMi6p088q1B2D9j4JWEdQiZ4FgGPvEwseYyuRlO6oKsbXsx9XYedXoqFiLo91EQAMDj8U6dOlVdXe3p6Xn7sqJQz8B6vHuK5PSog8TBGuridBxGQesSMcGJqXh2E73mMglXnYGsaHsx9VYedXoqFu7mDFGwQ0hICIyCVgGr8q7SEWD9FTI9i1gchZ56CHeyX5SDEDHB8Qz8fDP97AUYCyHr+PgG9a98p4qCsbGxy5cvt3cunAoMhF3ycxUVv4+o1YKCWYx/xsBhMX1IzAInpuLXpPSy8yQFgyHUO/++Tm0ppM485DxREAAQERGxYMGCHr+cpmmDwWCVnFi2PLTKoewLBsIHuNJGjzlEvJ1H7RqNfTsW8+HYO0MuwI0Bjk7By5X0wjOkGW5rD/XU67nkzhLq7ENOuP1LYWFhfn5+52muXLlSWlp69/N5eXkxMTFWyUZcXNy1a9d6/HKCIPbu3UtR9v+Rw0B4XzcU9KyT5MyT5KJINPdhfKyvs/2WHBmfAQ5PwdVmMOskaYAjSaFuogFYdYk8WEtnZ+L+/X+M6N1++OGHb775pvM027Ztu32LXQek1+vnzp3rCOvDwVGj91Ago/+VR51pol4cgH07FuPAQrIHNgb2TcCeOEtOOUr8OhEXMu2dIaifMFNgSTZZpaFPP4SLnPG0uXnz5unTpwmCWLduXXh4+NKlSwmC+Pzzz3Nycvz9/VeuXOnj43P58uXc3Nza2tq2trbk5OR58+bd81Amk+nTTz+9du1acHDwypUrPT09Lc+fOXPmwIEDMpksJSXlueeeI0ly9+7df/75J0EQaWlpjz/+OIretxG1efPm9PT0vXv3tra2zp49OyMjw/J8aWnpjh07pFLp2LFjFyxYgCDItm3bAAAbNmxAUfTpp58ODAy0dlF1FWwR/s3JBnrqMSLjKDnYE6mYx1iTiMIoaEcMFHwzFkuSIGMOEU06e+cG6g+0BJhxglCYwPEpzhkFAQBCodDX19fb2zs1NdWyWfyCBQsOHjw4Z84ckiRTU1MVCoWHh4e7u3tAQEBqampoaOj9DjVr1qzTp0/PnTtXrVYPHjxYq9UCAL788stFixYlJyc/+uij7e3tFEUZDIbS0tLp06fPnDlz586db731VifZ2759+6OPPhofHz958uSnnnrqt99+AwBUVFQMHz7cy8trxowZW7ZsefnllwEAlh7alJSU1NRU+w5/hdU8AAAoTOCbMurzYgoBYHUi+ssElAVnsjkGFAH/HY69k0+NPEhkTcZiRU7YzQVZS6seTD9BxIuQz0f17TLuqipd6Z76PnyD24gieRFz/W9/xs/PLzY21mAwzJkzBwBQUVHx22+/NTQ0iMXiqVOn5uTk7N69e9WqVSEhIQkJCZY093T9+vXs7OyGhgY+nz916tRLly7t2bNn6dKlr7zyyldffTVx4kQAwJQpUwAAfD7//fff1+v1LS0ta9euXbdu3WuvvdZJnpcvX24ZzqPVav/9739Pnz5969ats2bNeuGFFwAAERERycnJr7/++rhx4wAAs2bNYjB6uuqrlbh0IDRT4EQD/W05daSOmhKIfjwCG+OLwIrWAb2chAbywLgs4sfx+Ggf+BVB91CipB86Ri6MQF9L6fNx3fwAdsLykD5+k1tQ5gNCenl5eUhIiFgstvyZkpJyzzEydysrK4uKiuLz+be/UKFQNDU1DRs27PaUWq12/vz5xcXFwcHBer2+qamp8yMnJSVZHgwcONCSmbKyso5lwePi4jAMq6mp8ff3v+8hbMsVA6GWACcbqAM19MEaKlqELIhAPx7BkMBdIxzbggjUn4fMPUV8MBRbGAG79KG/OdNEzz9NvDsEezzSFucGykDZ7nbueO2YtyASiRQKRcfzcrm841Zf50QikVwuv/2Ffn5+fD6fyWRKpVLL7hYWu3btwjDMEtKuXLkyadKkzo/ckR+FQmGJ0GKxuONJnU5nMBg6IrcjcJUKxUiC883023nUxCOE73fmj29Qye5I3kz8j2n4/8WiMAr2C+m+yO8P4W9cpdZfgVMMob/sLKHmnyb2jMNtEwUdgaenZ21treXxgAEDUBTdu3cvAKCurm7//v2W8Sm3p7mnwYMHq1SqrKwsAEBlZWVWVtaUKVNwHJ8+ffpbb71lGczZ1tZG07Rer0cQBABAkuSWLVsemL0vvviCIAiKoj777DNL52pGRsbXX39tiYVbt25NSUnx9vbm8/kcDqfzTNqGc7YISRo0aOlyFSiS04VyOredvqmg48XIaB/kuXgsfSLCc87P7fxiRcjlGfjME8Ssk+TXYzGBne8sQHZGUODFP8nDdXR2Jh4pdKE+83nz5u3bty8kJGTUqFHffvvt999//8QTT7zxxhutra3r1q0bNWoUAODJJ59ctGhRSEjIzJkz//Of/3S8FkVRNpsNABAIBHv27FmyZIlAIGhra3vrrbdSUlIAAJ988snixYst3a0Gg+HmzZuPP/74V199FRsbS1HUnDlzRCKR5VBsNvuew0cjIyMTEhJMJlNISMhnn30GAJg/f35OTk50dLRYLGaxWHv27AEAIAiyadOmCRMmoCi6f//+jg5V20McbV2AgICAy5cvd7HvuFFp+G8JjiCIxgxUJqAyg1Y93awHjTrak41EuIFYEZIoQZLdkYHuSM+2ooY6oVarb+8/sRkTBZ69QJ5voQ9MxCKcaMWQ7rJX+TsIqRHMO0UwUPD9uL4dIJqWlvbOO+9Yoosja2trk0gkGNbtmq61tdXd3f2OF1p2PfT29kb+N3CisbFRIpFYgmgnYmJidu3alZycrNFoPDw8bv8ns9msUqnc3d27m0MAwKFDh55997Phm3714QLLuvxWPP/7d8sIQ2gRA6AoEsADbgwgZAJPNurLBX5cBA77dGJMFHw2Cvu8mBp1kNiRhmcGuW4sdFlX2+k5p8jZoci/BmNwQ2eLLt4avJuXl9fdT3K5XC73b9s2+vn5df2YbDb77pDJYDB6FgUtfLlgaiDSF8sQ9+9A6MYALw1AMMxVbgxAt1sWgyZJkHmnycttyMYUWBu6kJ0l1PocctsIbFYo/O07nDfffDMsLKwvjixhIY/1zUA5eBpB/dgwLyTnYfxSKz3xMJxx7xI0ZrDoDLmlkMrOxGEUdExz58719va2dy66B55JUP/myQbHpuDpfmjqAfOxese64Q1ZV76MHvwrwcTA5Rl4tCsNjYH6GgyEUL+HIuDVZPSHcfg/z5OrL8FFup0QDcCHhdSkI8RryeiONIzbv2/pQA4HBkLISYz2QfIewRt1YMivRL4MNg2dR52WnnyE2FdFXZ6Ozw+HVRZkffCsgpyHmAV+HIe9OACddIR4O48i7L/NGdRbX5VRgw4QY33Rs5l4iAB2h0J9AnYxQM5mYQSa7os8dY48UE3tHI0NkMDas1+q19LLz5MNOnAiA7fvl5iamjpv3ry7JwPQNI3AxYltRafTjRw5so8O3r8n1Ov1eiaT2YMJpJBVOPKEbhqA3aXUuivkP2PQDQMxp1xOwZHLvzcoGnxeTL2eSz4bj61NQhn27re63zLTWq3WvpsHuRofH5/bpzbCCfUQ9AAIAIuj0CkB6LMXycR9xKcjsQn+8OK9H8iX0cvPkzgKzmTicY6x6xaHw7nnxDhnvRBxQTAQQs7Mlwv2jsey6uh/nicHeyKbh6IBPIeoW6G7KU3gtVzyh0rq7UHYkug+30oJgjrYu9MBgvreQ4FI0Sw8RgiSfyH+lUfpCXtnCPo7igY7S6jYn80GEtyYzXgKRkHItmAghFwCBwdvpGJ/zsDzpHTsz8SeCriPk6M41UinHiC+LqMOTsI/H4W5wz3RIJuDXaOQCwkVID+Nx84302suk5uvU+8NhjcO7SlPSq+7QlaowLuDUbheGmRHMBBCLmeUD3JpBr6vinr6AunPBW8NwkZ4w3BoU8UKeuNV6mwT9Uoy9s8Y+48LhVwcPAEhV4QAMDsULZqFL4hEHztDZhwlLrTAvlJbKFbQi86QY7KIge5I+TzG03EwCkL2B89ByHXhKHgyCi2Zg88MQRecIccfJk40wMdCTxEAACAASURBVHDYV6620/84TY7NImJESNlcxroklAc7pCDHAM9EyNUxUbA0Bl0che6poJ6/RDIx8EIiOjsUtlSsgwbgeD29+TpZrACrE9EdaQw+w955gqC/g4EQggAAAEfBokh0YSSaVUt/WEi+9Cf1dBy6JBr1vHNdLairtAT4tpz6qJBiYmB1Ajo/HF5bQA4KBkII+gsCQGYQkhmE58vorUVU9F5zZiC6LBYdCUfTdEeRnP6imPqunBrji348Ekv3haUHOTQYCCHoHpIkyI407IMh2K5Sauk5EgHgyWj0sQjUh2PvnDkwlRn8VEntKqVqNODJKOTaTDwQruMD9QcwEELQfYlZ4PlE9PlE9FwzvauUivvZPMILeTQCnREMB3r8xUyB4w30nnLqcB01zg9dl4RODUQxGAGh/gP+miHowdJ8kDQfbCuB/VJNfVdOPf0HOTkAnROKZASiLrtbOkGB35vovVXUgWoqSog8Go5+NIIB14WB+iNX/RFDUPfxcLAgAl0QgUqN4Jdq6otiask5Mt0XnR6MTA1EvV2j11RpAscbqN9q6CN1VKQQmROK5jyMB/FhAxDqx2AghKBuc2eBp6LRp6JRuREcrqN+q6XXXDaHC5Apgcgkf3S4F4I71/BIigZ5UvpEA32snsppp0f5INOD0PeG4H5cGP8gZ2DNQGg0GlmsznpGzGYzgwHnEEHOQ8wCj0Wgj0UAgsL+aKGPN1DPXyJLlfQIb2SML5rmgwzyQFj9c09gkgYFMvpcM322iT7bRHlxkEkByAsDsLG+iMv2BkPOyjpndH19/fz58wsKCphM5ocffrhgwYI7Erz++uu7du1qa2vj8XgrV67cuHGjVd4XghwEjoIxvsgYX+ztQUBuBNnN1JkmetVFqlhJJ0mQIZ7IEE8k1QOJECKO3Iaq09K57XROG32plb7SRgfwkFE+yKxQ5OMRDF/ug18OQf2UdQLhqlWrEhISzp49m5ubO27cuPT0dH9//9sTeHp6njhxIjo6uri4ePTo0fHx8XPmzLHKW0OQoxGzwIxgdEYwAABozCCnnb7cSu+tol/OoWQGeoA7kihGEiVIjAiJFSF2nI8hM4ISJX1DTt9Q0AUyOl9KowhI9UAGeyJrErGhXogEjnyBXANC071dXFGhUHh6epaVlYWEhAAAMjMzR48e/dJLL90v/ezZs+Pj49944417/mtAQMDly5fviKP3o9frmUwmhvXPvqf+T61WCwQCe+eiP5EZQb6Mvi6jb8jpmwr6poI2kCDCDQkVIKECEMRHgvjAj4v4cYEn+8F9ql0pf4ICrQa6RQ8atKBOS9dq6GoNqFLTZUqapEG0EIkXIzEiJEmCJEoAvOfXLfD8ty8rlr8VWoQ1NTVMJtMSBQEAcXFxFRUV90ssk8nOnz+/fPny+yWgaVqpVHK5f+uIEYlEiEN3KUFQl0hYIN0XuX2lFYUJlKvoKjVdrQblKvp0I2jSUY060KqnWRjwYCMSFhAxgZCJcHDAwwEHB+z/BUiTCWcySctjMwU0ZmAggZ4AchOtNAGFCbQbaLUZeLKBNwcJ4IEAHhLAQ6YFgVABGuGGwNXjIMjCCoFQoVDweLyOPwUCQVVV1T1Tms3mRYsWTZ48ecKECfc7mk6nGz58OIr+bdRdXl6eu7v73Ylhi9C+tFotvEDpJRyAGDaIYQPgeec/qcyI1AgUJqAwISoz0JNATyJ6AhipW2WO0yYW+N9jjA5hAyZKc3EgZAAhEwgZtIQFJMz7dPkQQKPpqw/lIuD5b19dLH8ul3tHQLmbFQKhh4eHUqns+FMul3t7e9+djCTJhQsXAgC2b9/eydF4PF7Xu0YxDIOB0I5omubz+fbOhdPiA+DXaQK12iwQwGad3cDz376sWP5WmO4UEhKC43hRUZHlz2vXrsXGxt6RhqKoxYsXy2Syn3/+mclk9v5NIQiCIMgqrBAIeTzeo48+umHDhqampu+///7atWvz588HABQUFGRmZlrSrFix4tSpU88888z58+dPnjx58+bN3r8vBEEQBPWedaZPbN68efXq1SNGjPD19f31118lEgkAgKZpgiAsCTQaTVxc3EcffWT5c+rUqXe3GiEIgiDI9qwwfcK64PSJfgQOH7cvWP72BcvfvqxY/s61JCIEQRAEdRNcNNBpyQ2KwrabhW3FNaq6Jk2LVC83U4SBMLgxBUK2mxfXI1wcGikOTfUZ6M4R2zuzEORaNCbt1ZaCEmlZhaK6UdOiNKhUJhWOMlgYU8hy8+V7B7r5xXlEJ3rG+fC87J1Z5wcDobNp10lPVJ89V3epVlWf4BkT7xEzPXKKL8/bnSthogw2zlYaVUqjulnbUiGv/qP+z625O7y4HqMDR0wJG+fNu2suGwRB1qMyqk9Un/295lylomaAV1yse3RmxOQAgZ+QJRCy3MwUYSJNcoOyWdNSrao7X3d529VdQpbb6MBh44NHBwsD7Z19pwXvETqP3Ob8X0qz8luKRgcNHxs0Mtl7AI4+uGQomropLT1RdfZ0zblo94j5sTNTfAZ08R3hPRL7guVvX90q/3J51fc39l9qzBnuN3hS2Nhkr0QG9uCteCiaLpaWnq29cKL6TKDAf0ZUxtigkSgCb2kBYNXzHwZCZ3Ch4c9vCn/SmfVzYqaPDxnDwXsyydpMmk/WZP9wYz8bZy9NWjjId+ADXwIrYvuC5W9fXSz/Mnnl9rxvKhTVc2KmT4uYzGP0ZCMPgiLP11/aV3JIppcviJ89KSwdQ1y93oOB8BYYCG+0l2y7+qWeMDye+I9RAcPQXi/4RAM6u/bijvxvvHieT6csCRMFd5IYVsT2Bcvfvh5Y/m269s+vfX21Jf/xxH9MDZ/IQK1wKyq/teir6z+06aTLkh8fFTC09wfsv2AgvMWVA6FML992dVd+a+GSpAWTQtN7HwJvR9LkofLjuwr2TApNXzzg0fs1MWFFbF+w/O2rk/InafLn4oPfFf38cNTU+XEze9ZJ04k/G69+em2XiC1cPXh5kFuAdQ/eX8BAeItrBkKKpg+UHv7q+g+ZkZMWxs9l4321a5zCqPzs6u6rLdfXDnsm1Sfp7gSwIrYvWP72db/yr1BUv3Nhi5gtWjV4mb/At4/enaKpX0oPf339x8yIiY8n/oOJudzSlTAQ3uKCgbBW1fD+pY9QBH1h6Mogty6VUi/lNue/d2nrML/UFSmL2X+/sIUVsX3B8revu8ufoqnvin7eV3JwefLiKWHjbJAHmV6+NXdHubzypWHPJnq61nJdMBDeYq1ASBopk9Js1pKkkSKNJAAAZaAYE2WJGSwxA0EdYqcVGtA/Fx/8tnDvEwP+MSNyqnX7QjunNeu25mwvai9+bdSLkeKwjudhRWxfsPzt647yb9G2bfpjMwtnvjzsOQ/uPbaN6zvn6i5tyfl8fHDaU0kLHKdpaFSYDVITZaYtlSrGRFEmyuDjTD6O86zQeoGB8JYeBEKzhtC1GPVtJkO7Ud9mMkhNRrmZJmmmEGfwcZSJ4mwMAECaKMpEGWQms5rg+rHFUXxxrMAttCfDvayiVdf+zsUtZpJYP2KVH9/HLnk4VXNua84XCxPmzoqeZnkGVsT2Bcvfvm4v/3N1Fzf/uW1e7CPzYh+x5UVqB5VRvfnPbTWq+ldGrI647WrVxnRNBtkNtbxEo67V4yyU7cHEWCjGulWpkkbKrCHMaoIiabaYwXZnsj2YHE8Wx5PJ8WSxRAzQnZKDgfCWTgIhRdAmpdkoNxvkJoPUbGg36duNhlYTwADHk8X1YnE8mWwPFtuDyRYzcO59QylF0OoanaJUIy1QIRjil+buNUiEYDY90c/U/rHlymdzYmbMj5tp3ylETZqW18+958v3fmnYMzwGF1bE9gXL374s5U9Q5Od5X2XXXngjbW2Me6R9s3S86sy2qzsfjZ89J2Y60q2o0ks0kBaqGs9JDe0mjyQ3URTfLYyHse5bWZFGyigz6aUmQ7tJ32bStxn1rUbCQHE8mBxPJtv9Vs3MEjOYIgbGvPdxYCC8RS3VNp9SIAhCUzRpoEgTRehJQkea1QRlpphuDEvf5q3rDg8mx4PV8yY5DRRlmoYz7Ua5OWKuv21ahwbCuDV3e15L4asj19j9N2ZhJs1bc3dcbSnYNPplD1QMK2I7goHQvtRqtRknXjv3Hp/JXT9itRvTIb6LJk3Lpj828xjc9SNWidkiG7yjvtVYvreRNFEB4zzcE916fC+JNFL6NqOh3aRvv9VXZ1SYTQozAIAhwBk8DOdiGAtjCPDwmb4ABsIOGoVWVaRHURQgCM5BUQaKczEGF2MIcJzTVyNopAWqygNN7gPcQqf79OntwypFzcbz70e7R64evNzqw697yXLhuSzh8YzoCfbOi+uCgdC+/qy5+v7Vj6dHTl6YMNemza8HIWlyd8EPhytPrh++6p7jva2oMVtad7ItaJKn70j3PioD0kiZVAShIwk9SRpJmgKeyUIAA2EHe40aJQxkydd1AEFiFgV20vzvjYPlx3bkfbsiZfFkm4w964FyedWGM2+PC0lbOnAhXPPJLmAgtKPfyo7uzP9u/YhVQ/1S7Z2Xe7vWcv3tCx9ODkt/csCjfbEMDU3Rlb80qap0cUuCWeIHLxdndTAQ3mLH6RM0RVfsb1LX6BL/L7STW4w9oDXrNl/+pEZV//qol2wzQaLHGqVN/87bhqP4ayNf4DN59s6Oy4GB0C4Iivwo54uC1qL1g1ZF+UTYOzudURiVb1/40EAYXxu5xpPrYcUj0yR9c3ctTdIxjwf1UWPggeB+hPaHoEjEbD9RFP/Gl7UUYbWLiVJZxdIjq/lM/qeTP3DwKAgAEDD5H4zbGOjm93/HXqhVNdg7OxDU5xRG5fOnXpHqZdsmf+DL87Z3dh5AxBK+n/76ML/Ufx5dc7Ehx2rHpUHZTw0AgLingu0VBa3LGT6DHYVm+rAljOKvamnKCrFwf0nWS79vXJq08Pkh/+c4k4E6hyHYM6lL58fNevbEy382XrV3diCoD1UoqpcffSHJK2HT6PVcBsfe2ekSBCCPxc9+M23dh1c+/fTqLoIie3/MmiMtuhZj9MJAB5lj3XswEPYOAiLm+pMGqu5kW28OozZpXsl+52jVqW2TP0gPHmWt3NnM1PAJb41e/96lj366ecDeeYGgPnGu7uKaU6/+c+CiJUmP2WWmYG8kesbuyNhSq2p45sTaJk1Lbw7Vnq9qy1PGPxV8v1kN/ZHzfBJ7QXEk5vHA5gsydY2+Z0e43nZzyeFVvjyvTya9b6/J8r2X4Bnz6ZR/n6g++87FLWbSbO/sQJDV0ID+6vqPW3N3vJ++cVxwmr2z00NuLMG/xm4YHzzm/4698HvN+Z4dxKQiKvY3Ri8IZPCdalN3GAitgMHHw2f5lXxbRxqpbr2Qoqmvrv/w2rl3Vw9e/nTqEqvs0mJHXlyPrRPfNZKm505ukOrl9s4OBFmBgTBsPPf+n025n03+d5Qk3N7Z6RUEILNjpr2fvnFnwXfvXfrIQBi693oalP3Y4DdSIgjqH93CXQcDoXW4J7q5hXKrs7rR59CsbX32xPr81qIdGVuG+w/qu7zZEhtnvT7qxeH+g5YfXVMsLeu7N6IIWtdiVJRqLP/pmo3dvQqB+iOKoPWtRmW51vK9axsMpKkPv/dmbevTx9dyGZwt49+WcMR990a2FCUJ357xIU3TS4+sLpGWd/2FTRdlhI4MmODZd3mzl/7dBHEoYY/45r5b5jNczPN98OT3Y5Wnt13dNT9+5tyYh/vd/YbOIQBZmDA3TBSy7syb/5ds5XmQ+lZje75KWqDUtRhZYuatxQlpYFKZDXIzzsEEwRy3EK4omt+VbwHqFwwyk6JYo6rSqWt1RrmZJWYwhbeWwjerCX2bkSVmuicKPAYI+VZtqVxruf7mH/9eED+7Y3Fdp8HB2euGP/d7zfm1Z96YHT390fhZD5wKTOjI2qOtiStCnWaAzO3gPEJrajwnld1QJywL6SSNwqjcfHlbnbrxlRHPR4hDbZW1PtH5PJ4aZd2G7H8N8U15OvXJ3s/nNchMNUdalaUaj2Shpf199w/SIDOpq/Wqaq38poYmafdEN49koVsw15EW/bAm555HqGs2tl1TSAtUZh0pjuELw3mCYC7Xi3Xnt0kDTaNBWqBsu6pkuzNDMr35AVYIh3uLf9tzY99rI19I9k68XxonKP82Xfs7F/9rIIzrR6wKEPh1krLyQBNF0BGzO0tjY3BC/S2OFghpkr76QXnYDB9x7L2/nrO1F/6b8/nk0HFPDniUgdlhLQbreuCJqDFp37qwWWfWv5G2tjfLHjacba872eaX5u4/xqOL85b0rcb2PGVbnpIy0V6DRV6DRGz3/jEjpeucoCK+m1lDtOYqWq8oCD3pMVDomSTkB3K6cilDk3TLZXnt8VaPJGHodJ8er4xvIIwfXP64VlX/1uiXvXlenaR0jvKnAf1LyeHd179/PHHeI1GZ9+yg0reZCj6qSFkb6VBjZGAgvMXRAiEAQFakrs5qTn4h4o72ityg2HLl8ypl7dphz8R7xNgre9bVlRORoumvrv+QVXFi46iXEjy7/cFJE1X+Y4Neaop9Iogl6smlg7bB0HJF3nZVyfNj+wwXuye42XjzkL7jHBXxLTRQlGmaL8oUpVpJgsB7iFgYxutBU540UKV76s06MubxQKag27V2g7rp1ex3IiVhzw9ZwXrQXF5nKv96deN7l7bSNPXSsGfvXsrj5q5aQTAnYJxj3R2EgfAWBwyEAICCrZV+ae4eA4WWP2lAH644uT3v64fCJz6e+I/+MlO+K7p+Il5syHnv0kePxc+eHTOt68sTk0aq8LMqrjc7fLYfivcqelEELStUNV+U6VqM3kPEPsMldlkd0bqcoyImtGTLFXnTBRnOQn2GSzxTRBi7d4P4aFB3sq3lsjxxZWi3Lp4sewouHvDojMiMrqR3jvLvQNH0r2WHdxV8Pyt62qPxszoGsWsbDUXbawZtiOrlb9DqYCC8xTEDoeyGuvZo68DnwwEAVcra//z5KUGZ1wx5ur/fEbxbt07EZm3ra+fe9eF5rR32LI/x4E2sKIK+saOGLWFEzPG34k0+faux6YKsNVfhFsL1Hekujub33zuI/b0iVtfqm/6QygrVkgSB7wh3QbA1h7o0nG1vuSRPXBnG6MLOa2aK+OLaV9l1F7u1p2B/L/97atG2/Tfn83p10/ODlw/0TgQAlH5Xz/VlB4yz5lKlVgED4S2OGQgBDa5+UOaXKdmnP3Ci+sziAY9Oi5jiZENDLbp7IppJ87Zruy415Lw26oVY96jOktKg+OtaAJDohQF9MUqNMlFt15RNf8gIPekzQuI9RNyV6tLR9NOKmDRRbVeVzRdkhIH0HS7xHiLu+S6hnao53CIv0SQ+Hdr5GihNmpY3zn/gzhGvHf5st/YU7Kfl3xXn6i59nLsj3iN6adiimk9lgzdE97aZ3gdgILzFMQMhRVMnsv5ov6psmlyzdOBCIcvN3jnqKz07Ec/VXfzPlc/mxsyYF/vI/a4PGrOl7XnKxKdD+/p+nrpW33xBJr2uEsfyfYZLenZfyl76XUWsbTQ0X5S1XVMKw3g+IyQ2aI6X7qlHGWjEnPuOdTxVc25rzvYF8XNmxWR2d0/Bflf+3WIkTd8V/aw+Zgr3CJ6wcDjbwbZEBTAQdnDAQHi+/vKOvG88WO6z/pibuDSc7+9wZ48V9fhEbNW1v/3Hf1AUXT981d27w2gbDYWfVSetCmNLbHQ/ldCTrTmK5osymgTeQ8Veg0U9GGdhe/2lIiYMZPs1ZctluUlFeA8Vew8V92zcUw+QRura5vLQaT7uiXdej2rNui1XPi+Rlb86ck2kOKwHB+8v5d9jhJb8818l58f9/qcq94nEeRnhE/piX8Meg4HwFocKhLnN+TvzvzOSxqUDFw7zG1T/e7u+1Rg5z9G3UuqN3pyIFE3vubHv5+JfV6QsmRQ69q/nCTrvPxUB4z28Uns+3aLH1DW65ktyaYHKLZTrNVgkiXdztAECt3PwipimaGWZtiVHIb+hFkXxvIfaogl4N3WN/uaXNQPXRDDd/rq4udZy/d2L/x3ql7oi5Uk2zurhkR27/Huv4Uy7rskYOd+/WFr2Rd7Xrbr2xQMeTQ8a5SA3emAgvMVBAmFOU95XhT8qjarHE+alB6dZzhKzhsh9p2zwa9HOsV/XPfX+RCyTV7594cNAgd/zQ/7PMtGw5mirvtUYsyjQSnnsCdJESQtUrTkKTb3ePdHNM1kojOA54IIaDloR00Bdq2+7pmjPU7JEDM9UkVeKqI/uAnZR7bFWXZMh5okgAICBMG7P//ps7cWXhq4c4pfSm8M6aPlbz9X3yiLm+buF3Braltucv6tgj9qkWZQwLz141AMXo+lrMBDeYt9ASNJkdu3FPTf2mSliQfzsccGj77hQKt5dK4oR+AxzkiUK72aVE9FMmndf/z6r4uSKlCfHiEfmbS4fuCbCZl1nnTOpCEuFbpCa3RPdPJLchOE8x5mG6FgVMQ3Utbr2ApU0X4UyEI+BQs8UIcezh40t66II+uq7ZZHz/Sv5le9f2prgGfvMoKe6NS7mnhyr/K1NVaUr/6khZe2dY2hzmvK+LvyxTSedF/vwlLDxPW5P9x4MhLfYKxCqjOqsihO/lGb58rznxT0y3H/QPW+zy2+qa4+3Jj3Xv1es74QVT8RSWcV7lz6aWDw5MSImfnpPbtj0KYPMJM1XtReo9O1GcYzAPV4giubjHDt3RThCRUyZKGWFVlqolt1Q4xzUY4DQfYAbz8/hbo3XXWkuOlq1K3Hn6iHLh/lZZ417Ryj/vlP2fQPXj+U/5t6zJorai7+/8cv11htTwyc8HJXR+RI8fQQGwltsHAhpQBe03jhYduxi45W0wOGzojM7v8dOU3TOW6XxS4O5TroAtHUrAkW1Nn9n2UeJ/30kbuq82Icdcwk6k4qQ3VDLilTKCi3Ply2O5oui+Pwgjl06Tu1WEdNA22xQlGoUJRpVlY4fyJHEC9zj3dgejrhYBA3oo5Wnt1/7ZlnJsvjx4QFDrVZlO3EgJA3UlU0lqS8/YE21Jk3L/tKso5Wn4j2ip0VMGeafasvRNDAQ3mKzQFivbjxZffZY5e8sjPlQxKTJYeld7FepOdpKGsiwh337Ood2Yd2K4Pq2Kq/BIjrO/HHujipF7TODnrLWlXtfoAhaVaGVl2gUpRqj3OwWyhWG89xCufxAjs36Tm1ZEdMUrWs2qip1ykqtskKLsVBRFF8cxRdF8R1whlmHEmn5lpzPAQCrBi/zU/sXf1U7aEOUtb4gJw6EzRdlilJNzONBXUlsJE2na85llR+vVzdNDBkzIXRMtCSir3MIYCDs0NeBsFbVcL7+0u8156V62digUZPD0rv7BRtkpvz/Vg55PdoBh1r0nhVPRHWNvuSbutT1kZaC+rPp6tac7T487xUpi0NFwVZ5i75j1hDKSp2qQqus1BraTTw/tiCYKwji8AM5bAmz78ZJ9nVFbFKaNfUGda1OXavX1OgZbrgl3gsjeA5yE7cTbbr27fnf5jRdWzpw0ZSwcZabF0VfVHsMFHoPsc5teycOhAVbKwMneN5v84D7qVc3Hq/6/WR1NgKQscEjRwcMj3IP7+7szK6DgfCWvgiEJtJ0ve3m5cbciw05esIwKmDomKARSV7xPR4ilb+lIvghb1Ek34qZdBBWPBFv7q4VRfB8R7l3PENQ5K9lh78p/GlkwNDFifM9uO6dvNxxkEZKXavXWIJHnZ40Ujw/Ns+PzfVlc71ZXB+WFe8sWrciJo2UvtWoazHqmgyaRoO20QBowA9g84O4giCOIJjbXxbf0Zp139/Y/2vZkemRUx6Lm81l/LVym6JMW7m/MeWlSKtUzs4aCE1K87V/lw/ZGNPjpnOJrPxs7YVzdZcMhGGoX+pQv9RUn6TbvwirgIHwFmsFQgNhLJaW5rfeyG8tvCktDRMFW768KIkVLmcafm/XS00OtY+XtVjrRNS3GQu2Vg1+JQq9ayksjUm758a+g+XHHgqfOD9uZr9bpsesJbWNem2jQdds1DUb9C0mBAMcTxbbk8l2Z7IlTJaYwRIzmG6MHkxY7Fn50xRtVhMGmdmoMBtlJoPUpG836dtMhI7kerM43iyeD8sSuR2/2XcHA2E8UHb4hxv7h/sPXjzgUa+71moAAOR9WBE0yUsSb4Xz1lkDYWO2VNtksMoc6DpVw8XGnMuNuTfaS0KEQSneiYlecQkesXwmr/cHh4Hwlh4HQgNhrFLWlMurSmUVxdKyWlV9uCg00TM2yTthoFeCda9cnLh31FonYvneRqYbHjT5vqMY2vWybwp/Ol1zblrE5DkxM8RsYe/f1F5MakLfajRITQapySgzG2Qmo8JsUhE4B2MKcIblPy6G8zCcg+EcDGOhGAvF2RhAgaU1iWCIZXKqRqPh8/kAAMpMU2YKAEAaKZqiCT1JmWjSQBIGitCTZg1B6EizljSrzCY1YdaSDB7GkjBZIgZbwmC7M9nuTI4niyVi9KPl5e5gIAy/lR39/uYvSV7xixPnBwvvOw+1PU/ZeE464BkrjEx21kBYsLUycKKXOMaanVgm0lTUXnKt5fr1thvF0jJPrkese2SUJDxcHBohCu1ZXLRi+feDdaR6SWlUtemkrbq2Jk1LvbqpQd1Uq6qX6eXBwsAIcWikODwjbEKkOLTvxiiyJUy2mKGs0IkirXAR5HzMGqI9X5n6cmdrcHtwJKsHL38sfvZ3RT8vOrhiYujYebEPe/Mca3e0LmIKcKYAF4bfeTKY1YRJQ5jVhElNEFrSrCP1bSbSQJIGijRRhJ4ENCD0JACAJmnSSAEAaJpGEAQAgDIQlIECADAWiqAIzkFRJoqxMJyN4lyMJWbw/TkMPsZwYzAFOIOPOdM1mcqk/qXk8C+lhwZ6J/5n3JsPvKPsPsCtOqtFEPnfPgAAIABJREFUXaO37mYXTsOkNOtbjVavrJgYM9k7Mdk7EQBA0VSVsrZYWlYiLT9dc65SUcPFOcHCQH+Br7/A14/v48X18OJ6iNgimy1h079bhFK17EjNKRRFdWY9SZEG0mAgjDqzXm3SKI1qpVGlMCg4DI4nx92b5+nN8/IX+AYI/IKFAT48L1sui1D/e7vBGXtHrXJF1nC2XddkjPxHV/thZHr5T8W/ZlWcGOybPC/2YduMT3NMztoi6aJGTfPe4l9PVmWPChw6P27W3dvJ3k/97+2GNmPE3N52/Tll+VuxX7TrWnXttcr6enVjg7qpSdvSqm1v1bWrTWohy03IcnNjubkx+VwGh4WxeAwuG2c9nvgPAFuEHSia0pi0CIJwGRyMgXlg7mycxWNw+QyekO0mYrmJ2KKO7SXtyCPJLX9LZfhMX2e6EreWlsvyiDnd+NVJOOLlyU8sTJh7qPz4q9nvenAks6IzRweNcIQvGrIBiqZzmq7tLz10s70sM2Li7syP3TndGwXqPUiU+25Z6AxfJ17+sMfa85WBE209O97SBBzkO/D2J0maVBhUCqNSZVQrjSo9YTASRh2hZ2PWn5bdv+sOPoO3NGmh3dcafSDLmAhVle7uDjEXp67R0RToWMyw63gM7rzYh+fETP+j/s9fSrO25u6YEjYuM2JSgMDZmt1QB6lefqTyVFb5cT6T90jUQ2+krWNhPZnCzxDgwghee57Se6jTLn/YMyY1oWuxfr9oz2AI5s4Rd/cqp2f6dyDsRyRxAvlNNQyEd2i5LPceKu7xGA0UQdMCh6UFDqtXNx4qP/7MiZf9+T5TwsaPDRpplWFpkCMwkaYLDVeOVp4qbCseEzRi46iXot172x/uPVRcf6oNBsI7yG+oRVF8x1lN12ZgILQRcYyg7MeGkEx758ORkCaqvUCV8tKdq/r2QIDAb3nyE0sHLrzcmHu08vS2q1+m+CRNCBk9zG+QHRcFhnqDoMirLfmnq8+dr78cJQmfEjZu46i11vo2xTH88p8adC1Grjc8Pf4iL1ZL4vrZDCWrgIHQRgRBHLOGMMrNLHE/m5vVd9rzlcIw3u27xPUShmAj/IeM8B+iNeuyay9kVZx4/9LWwb7JowOHD/MfxGN0uwMWsj0TacppzjtXd+mP+j8D3fzGBo1aOnCR1fvHEBTxGixu/VMeMs3Hukfuv2iSVpRqw2e64s0FGAhtBQHiGL78ptpnhMTeWXEUbVeVfbRHFY/BzQifkBE+QWVU/1F/+WR19uY/t8W4Rw7zHzTMb1DXhxdCNtOma7cs53St5XqUJDwtcNj9ZsRbi1eKsGh7TUimT/+dPWldqiod25PJELhiUHDFz2wv4hhB2zUFDIQWhJ7U1OjFi7u0qm+PubEElohoIAy5zQUXG678XPwbApDBvsmpPkkpPgP63VI1zkRPGPJaCq82519pzpPp5YN9k9ODR60d/mzvdwrsCq4vG2Oh6lo4ofAWebFa0s3FRZ0GDIS2I47lV/zcSBF0DxbTcj7SApUomofdtaZaH2Hj7JEBQ0YGDAEA1CjrrjTlHa8688Hlj715ngO9EwZ4xid6xvaX5Uz7NZVRXdh+s6D1Rn5rUZWiJtYjKtU7ad2wZ6MkETabPd3BPUnYXqCEgdBCdkPT9em8TgYGQtvBORjXl6Wq0IqinXAB7u5qL1B5DxbZ5a2DhYHBwsDZMdMomiqVVeS3Fp2oPrPlyucsnJXgER3rER3rHhUhDoWjbKyCoMhKRXWxtOxGe8kNaWm7ThrrETXAM2558hOx7pHMHs1/sBaPAW43v6wNhb2jABgVZrOGEAS66DUBDIQ2JY4VyG6qYSAk9KS6Shez6L4LQtoGiqAx7pEx7pHzYh8GANSpGm5KS2+0l56qzq5S1vryvSPFoZHisHBxaJgopF8vcGpLWrOuUlFdIa8uk1eWySprVPX+fJ9o94h4z5jZMdPDRMG2XNSpczw/NoIhmno931UDQAf5TbU4hu+yFwQwENqUKJJfvrfB3rmwP2mhShjJc7R1PQLd/APd/CeFpgMACIqsUtaUySrL5JV/NFyplFejCBomDg5yCwh2Cwxy8w9w8/Pietq+N8/RtOtl9erGOlVDraqhWlFbrazVmLUhwqBwUUiUJHxq+IRwkUO3rT2S3NrzVTAQKsq0LnuDEMBAaGP8QLZRbjZrCAbfpUteWqDyHOjQDSwcxSLFYZHiv/YokOrl1craamVdtbL2XP3FelWjwqjy43v78X39BN4+PG8fnpc3z9OL5yFiOfRH6xmtWdeqbWvStrZoW5s1rY2a5kZNS4O6kY2zAgT+wcKAQDf/FO8BoaIgb55n3+3FanXuScLi3bUhmd72zohd0UBZrg114ZkkLl0d2x6CIm6hXGWF1iPJCevKLqLMlLJCG/VYgL0z0j2W1Z5SfZI6njEQxkZNs+W/Jk3LtZaCFm17q67NQBi9uB7uHLEn10PCFkk4YglbJGILRSyhiC0UsgRs3PqLJfaSiTSpjGqFUSU3KBRGpdyglOnlUr1cqpe162Wt2jYUQb24Ht48L2+epy/fO9Yjyo/v4y/w7e+zM/n+bEADF59Zr2s2YGzUlac4w0Boa8JInrLcpQOhslzLD+DgbEdfIfaB2DgrTBQcdte+P0bS1Kptkxrkbbp2uUEp1cmqlLVyg0KuVyiNKpVJTdG0G5MvYPL5TD6PweUxOFwGl8/kcXA2E2PymTwmymDhLCbGZGFMDEE5/9sgU8D8291lrU6rRrS3P6Mz60iaAgAYCCNBEQRF6AmDmTIbCKPerDeSJp1ZrzXr/vd/ndqkUZk0apOapCghSyBkuYnYQglbLGYLJRxxmChYwhZ7cN29uB5W317ccYhi+PKbalcOhIpyrSjSpQcuwEBoa6IIfvHFOnvnwp5kNzXiGGe+G8HCmJbbjfdLYCRNaqNabdKoTRqtWa8z63SEXm3SGAiDwqBsUDcZSZOJNBlJo4k0kzSlN+sBADSgNaa/hT2KolD0b/dZuQwuhqAAABbOYqA4hmJcnMNAGWycxcHZTJwpYPF9+F4cnMNjcHlMLp/BEzD5biwBx/EaqTYjiRU0Zkv9x/bh5H0HpyzXejj2rYq+BgOhrfH82ISWMCnNTKGLdkTIi9WxfTyP3sGxMCaL6977aYtOuR+e7QkjeSXf1ZFGytFGb9kIDVSV2nCn2y21W1zyi7cvBLiF85QV2gendEb6ViNN0Dwf121/QI4GY6KCIK6iTGPvjNiHpkHP4ONMl1xZrQMMhHbw/+3deZwU5f0n8Oqzru6eu4/pOYABZgQVUBSNiscaxVeCLsT8dInHSlQM8goe+BKjUXJ6YDDGrL7UiOZYk6grCRBXRc2iQhAHCQMiMAPD3D33THdXVV/VtX+0knHOnunqruvz/qunp6broenqT1U93+d58meygw0GDcL+I+GC05zaKSoEQyg4zdH/hUGDcLCeyzN2ByGBIFRE3izHQL1Bg7AvNW4XQE0KT3P2Hwkp3QplDB7n8mYaffFOBKECGDcpRpPRgbjSDck1MZYMNfH5sxGEoC60mzSZTXwgqnRDck4igif5vBnaHgOTOQShEkyEs5IONfFKtyPXgsc5Rxlt0JIEULeCakf/UcNdFPKBiI21GHx+DwJBqBTXdDbYaLggHGzg8g1/EwbUKTXAV+lW5FqwkXdNxyGJIFSIazpjwCAcaEC3PKhUXhUbPMFLSUnphuRUsJF3TjP6fVECQagURzktdEXFaFLphuROIiIKXVHMbgzqZHNY7Xk2rjWidENyKniSd01HECIIFWK2mlgfFW4RlG5I7gSP885KBosSg2rlz2IHjHR3NB5KJHiRcRt3brlTZAtCnuebm5uTyTEvcURRbGpqikSMdcI1DqfB7o4O1IfzZ6E3AtQrbyY72GCg0YSDjbxrOoNBvYRcQfib3/zG7/dfdtllNTU1X3zxxcgNamtrZ8yY8c1vfrO0tPQPf/iDLDvVOtc0JthooNPPwQYMVwJVy6tigyd5STRKN2GokcN90RQZgrCpqelHP/rR7t27GxoaVqxYsXbt2pHbrFq16t577z127Ni77767evXq3t7ezPerda4ZbOikYJDO+QQvRvvijjJ0EIJ6WRkLXWQ3TocFKmVOkSEI//KXv1x66aWnnXYaQRCrV6/+4IMPOjs7h27wxRdfHD58+NZbbyUIYuHChfPmzduyZUvm+9U6G2uxOS18pyHG8A42cM7pjMmCuzCgankzHQbpJkzGk3wg6kTxGkEQsqw+0djYOHv27NRjt9vtdDqbmpo8Hs/QDfx+P8N8eeoxa9askydPjvVqyWSyrq4uEAgMfXLevHlW6yhNlcREvLVZNGu15IctSfTvb7IltNp+kedjTFpnlH3/TjiKTbGW+mw3yVDSf/8hTWx+MrBf9FQPpLOxpt//ULtEFUqJzuNKN2TSTGaLzT9D3teUIQhDoVBx8X+W8mJZdnBwcNgGNE2Ps8FQ0Wj0vvvus9m+tkTR1q1b8/PzR27Md7Un/8/TU2+64vjT+xqLLYf+n9LtmKJkMimkdxYy0PtfJc4Pe+sDE28KaUv//Yc0JSUy3HVjz5//l4mYuM9C0+//AHemRczr/fNHSjdk0ky0g/2fPyYIIhxOq7KJYRiLZYJlwGUIQrfb3d/ff+rH/v7+oZeDqQ0GBgaGbjBnzpyxXo2m6XfeecfvH3NR06Gs1gr7vc9M+I9ULWezUP96m+/e/1K6IVOU5np4Yix58uEjFet/jLET8sJ6hNnQ+Xh93vd+xfonXilM0+9/8E8t3hqne+H/ULohGZHr/ZfhdGbevHl79+5NPT506JDZbK6qqhq6wZw5c7q6utra2lI/7t27d968eZnvVwdYPxXpjokxnQ+rDzXxrJ9CCoImGGTWp1CT4KxAB+GXZAjC7373u83NzU8++eSBAwfuuuuuW265hWVZgiDWr1//2GOPEQTh8XiuvfbaNWvW1NXVPfLII5IkXXXVVZnvVwdMFhPjJbk2nY+tDDVi9grQDOc0JnhS50GY4MUEL9IlGEr/JRmCkGXZ9957b/fu3atWrVq4cOETTzyRer6srMzn86UeP/fccxUVFbfddtuxY8feeeedUStfjMlRwYSadX7UBU+iShs0wzWNCek9CEPNgqOMxlD6U+QJpDPOOOPNN98c9uSaNWtOPXa5XE8/reWqlqxxVtA6XxFUIkLNwuwVCELQBrqEFGPJ6ECczLdNvLU2hZt5B+6LDqHVkifdcJTToWY9D+DlAhEba8WCZ6AZBlguNNQiYAThUAhChTFuMsGJCU5UuiHZEjrJu3BfFDRF9/Uy4RYBV4RDIQiVZiIcZVRIv7M6BRt5JyplQFNcuq6XifbHCYnQ8Y3fKUAQKs9Rzuh4esMgrghBaxzltBCIJnU6rimEy8EREITKc1TQei0cTXBiIiwyHlRpg5aYbWbGS4Z1Oq4pjA7CERCEynOUUXodShhqERzlqNIG7XFU0Hq9TxNOHZUwBIJQeVSBXYwl46GE0g2RXwhV2qBNznJGr+XcXJvA+nFUfg2CUAVMhMNP6fI+TLgZ0ziBJum1wyLaHyfMJrsLw5m+BkGoCqyf5tp0ePoZbsVNGNAkvY5rCrcKWB97JAShKjj8VLhVb1eEqNIGDTMRrJ8Kt+rt9JRrizjKJl5Yw2gQhKrAltFh3V0RhpoFZyUGToBWOSt02E0YbhMc6CAcAUGoCoybjIcSiYiu7sOEW1ApAxqmy27CcGuExRXhCAhCdTARrE9vgyhCqJQBLXNW0GF9XRHGw4lkLEkV2JVuiOogCNWC1dloQongWiO4CQPaRebbCBMRHYgr3RDZhFsjjjIK43pHQhCqhcNP66lnnu+K2pwWK2tRuiEAU+co19VFIdcmsCgZHQ2CUC3YMlpPhaPhFlRpg+Y5ymk9TYgfbos4StFBOAoEoVqwXjLSF0vGdTLPb7gV556geY4yXQ3wxRXhWBCEamGymOhiOx+IKt0QeWC4EuiAw0/pZsZRMZaMDSboElTKjAJBqCJsKcW16+LuqERwbaiUAc2z59lMZpM+6mX49gjtJU1mlMqMAkGoIroJQqEnamUsVgaVMqB5bBmtj3Jurh0dhGNCEKoIW0pxHbo45NowaBd0wqGXida49gjjw1E5OgShiujmijCMEYSgF6xexjVx7REWV4RjQBCqiM1hNVv00CERbhVQKQP6oJN1syWCCyAIx4QgVBfWr4eLQqz8CbpBFeph3exIX8xKW6w0uu1HhyBUF9an+SDEyp+gK7pYN5trj7DoIBwbglBddNBNiJU/QWd0sG42OgjHhyBUF0b7QYih9KAzOlg3G0E4PgShujBuMjoQ1/REa2F0EIK+6GDdbK4DQTgeBKG6mCwmukTbE61x7RGHH4cc6AddYo8HE2JEq6enYjQZDyaoYkyuNiYEoepoul4mIYgJXqQKcciBfpjMJtpL8gGtHpV8R4T2YHK18SAIVYf1Udo95L4sTsMRB/qi6dNT3BedEIJQdRgfxXVo9dYo+uRBl9hSKqzZIOQDUcZLKt0KVUMQqg7jJXnNzjjKtUdYdBCC7mh6XBPXgUGEE0AQqg6Zb0uKUjysyZksuDZcEYIOsaUU3xEhJKXbMSV8IIrptseHIFQjxktpsXBUSkpCV5Tx4pADvbHSFitjifTGlG7IpMVCCUKS7E7M9DQeBKEasV5Si+sxCd0xe57VQuJDBTqk0W5CvgOrL00M31lqxPg0eUWIShnQMdZP81oMwkAUHYQTQhCqEevTZL0M1x5hSzGnDOiTRutluI4ISkYnhCBUI8ZHcRrsmccVIeiYQ5tBiFuj6UAQqpGVtlgoi+ZW6OXaBAQh6BVVZI9ziYQgKt2QyZAIvjPKeHBFOAEEoUqxPo3VyyQ4MRmXyHyb0g0ByA4TwfgobfVZRPpiVgbr8U4MQahSmjvkuNQdGEyuBvrF+ihOU1VsGEqfJgShSrFeLR5yuAMDesb6SG0VjmIofZoQhCrFaK1wlA9gKD3oHOOltNVhwQdQMpoWBKFKMW5S6I5JSc1UjnLtqJQBnWNLKb4jqqFybr4jyuL0NA0IQpUy2812l1UzUzqlitNw7gm6ZmUsZrtJK+XcUlISemO0G4uDTgxBqF6Ml9TK/DKRvpiVRnEa6B9bqpm7o5GeGJlnNdvwJT8xvEfqpaGpt1GcBgahoXJurgPd9ulCEKoX6yW1slQ934HiNDAEVjvrZqNSJn0IQvXS2hUhDjnQPw1dEaKQO30IQvWiPaTQq43CUcxnCAbBpI5KUQtHJa4I04YgVC+z1UTm2YRutReOJhNSpD/OuHHIgf6ZrSYy3yZ0qf1WjSRKkb44XYKjMi0IQlXTROGo0BmlCu0mC2ZXA0PQRDeh0B0jC2xmK47KtCAIVY3RQr0MF0AHIRgI49PAUckHIizui6YNQahqmqiXQckoGAqrhYnWOMwyOhkIQlXTxBUh+uTBUBifBjoscFROCoJQ1Rg3GemLq7xEDaPpwVCoIns8lBCjSaUbMh6MnZgUBKGqmSwmqlDVJWpiNJngRaoQ8xmCUZjMJtpN8p3qPSqTCSnaH6eLcVSmC0GodoxH1Ycc3xGhPSTW4wVDYbyqXiVN6IpSRTYUcqcPQah2tJdScxByAazzAoaj8io2vjPKeHBUTgKCUO0Yj6p75tEnDwbE+ki+U71XhFgTbbIQhGqn8sJR9MmDATFeVY+px+npZCEI1U7lhaM8ptsG4yHzbclYUhRUWjjKB3BrdHIQhGpnspjIApvQo8YZRxN8MpmQ7Hk2pRsCkFsmgvGQkS41HpWSKEX743QJSkYnAUGoAaq9OxrpjGH2CjAmxkcJXXGlWzEKvitKFaJkdHIQhBqg2nqZaHcc8xmCMTFeMtqtxiAUUCkzeQhCDWA8Kh1BIXTFUSkDxsR4qUinGoMQ9WtTgCDUAFXfGsW5JxgS6yMFVfYR8p0RxoOjcnJkC0JBEGpra5ubm8faoKenZ9++fa2trXLt0ThoNxnpVWPhaKQ7jiAEY7I5rCaTKRZKKN2Q4fhAFEE4WfIE4b59+6qqqu66665zzjnngQceGLnB9ddfP3v27NWrVy9YsGDZsmWxmBrPpFTry0WxVVY4GhuMmy0mm8OqdEMAlEG5bWqbaC21MD2FheknSZ4gXLdu3dq1az/++OPPPvvsueeeO3z48LANVq5cGQgEPvnkkxMnThw6dOjll1+WZb/GocKl6rmOKOVBiTYYF1liU9tRKXRFqUIsTD9pMgRhZ2fnzp07v//97xME4ff7lyxZ8tprrw3b5oorrrDb7QRBOJ3OuXPndnZ2Zr5fQ1FhNyHfGSGLcTkIxkV7bGpboZfvxH3RqZDhi6ylpcXpdBYXF6d+nD59+jg9hQ0NDTt37tywYcNYG4ii+NFHH516tZTFixenctSwGA/VeyiodCu+hu+IUl5D/6eAwVFu++BBlR2VKBmdknSD8NFHHz148OCwJxcuXHjPPffwPD80pSiK4jhu1Bfp6+v7zne+s27duvnz54+1o1gs9swzz5Dk105q5syZ43K5Rm4sCILdbrdYLGn+KzTMJYbbhXA4rHQ7/iPUxhXOolXVJKPhOM5kwk0wxSQdca5dCIfC6lmGLNgWzp/LGuSoTPPzzzCM2TzBvc90g/CCCy6oqqoa9qTf7ycIwuPxDAwMJJPJ1M56e3u9Xu/IVxgcHFyyZMnll1/+4IMPjrMjmqZfe+211CtPyGKxGCQImWlSfX+ApVm1TBghEdGeRH6ly+FwKN0U45IkCe+/giRJstJBW4IkC9Qyy2Csp6OwMo9xGOKiUMbPf7pBuHjx4rF+NW3atPz8/D179nzjG98gCGLXrl133333sG04jrv66qvnz5//5JNPTrmtRnaqcFQlHQCRvpiVsVhIjEMFQ2N9JB+IqCQIUTI6ZTJ8kZEkuXr16jVr1rz99tsPPPBAd3f3tddeSxDE7t27KyoqUttce+21zc3NZ5999osvvvjCCy98+OGHme/XaFRVOIquCACCIBiPitZjQsnolMlT9ffwww8XFRU9++yzfr9/586dFEURBOHxeK677rrUBhdffPG8efMaGxtTPw6rhYF0MJ7UWqCj9JXmHh/A6ksABOMlB4+PXhKReygZnTKTJKlrvpKysrJPPvkkzT5CAxXLEET3Z4O9h4I1N5Ur3RCCIIij/7u1oNpBV1ucTqfSbTGuUCiE919BoVCI6Lcef6N9/j3D6ycU0fx2F0EQFUvcSjckR2T8/KOPRzPUdWu0A0tgAxCMl+S7olJSFZcTfCeOyilCEGoG7SYjvTE1zDgqJSWhO0a7cciB0VnsZrvDGu1TxTIUmGV0yhCEmmG2qmWp+khPzJ5ntdjx4QEgGB+phvllUiWjOD2dGnyXaYlKVujlOqIsFqYHIAiCIBgvpYajUsDC9BlAEGrJV4WjCuMDEQZBCEAQBEGwXlVcEaJkNBMIQi1Rybkn3xFh0ScPQBAEQTA+Sg0T4mNobyYQhFqikjUoOBxyAF9RybrZKBnNBIJQS9RwyCUTUrQ/Tpdg3QkAgiAIs9VEFdr4LoVv1XAduDU6dQhCLfmycLRbycJRPhChS+zokwc4hfFRyi5V/+XpKUpGpwpBqDGpSX4VbAC6IgCGUXyyC6ErShXh9HTqEIQaw3gpTtFDDrOMAgzD+ihlC0cx01OGEIQaw3hJZW/CcB24IgT4GsXLublAFIXcmUAQagzjVbhWm++IMLgiBBiCLrbHwwkxmlSqARjamyEEocbQJfboYCIZV+aQS0TEhCBSBSgZBRjCRNAlSnYTouc+QwhCjTGZTXSxne9U5pDjO6KMlyTQJQ/wdYxyVWxiLBkLJagimyJ71wcEofYo2CHBByKYZRRgJFa5Pgs+EGXcpMmM89OpQxBqD+tTrF4GlTIAo2J8FNeu2OkpSkYzhCDUHsZLKXVrlGuPsKUIQoDh2FKKaxcU2TU6CDOHINQeRrnZ7jFcCWBUdpeVMJlioUTud42hvZlDEGoPVWhP8GIiIuZ4v9H+uNlmsjmsOd4vgCawPpJvV+AMFR0WmUMQapCJoD2kkPN6Ga4D90UBxsT4KC7nQZgQRDEqkvkoGc0IglCTFJnSiWtHySjAmJQ8KlExmhkEoSYpcsjxHZi9AmBMOD3VLgShJrGlCtRq49YowDgYHyl0x3K8XCjfEWFwVGYMQahJbCnFt0eIHB5xyYQU6cWCZwBjMtvMZH6ulwvlOnBFKAMEoSZZGYvZbooOxHO2R74zShfbzVb0RQCMic1xvYyUGkSI09NMIQi1ii3N6SHHt6ODEGACrC+nY3wjvTEba7XSlpztUa8QhFrFlua0Z57rwKBdgAkwPiqX0x9y7VgTTR4IQq3K8aAlVMoATEiB01MclXJAEGoVm9tzT9waBZgQVWhPCGJCyNGsTzwqZWSCINQqxkNG+uK5WaE3FkwkkxJmrwCYgIlgfRTXlqMz1DBOT2WCINQqk8VEl+RohV6uTXD46RzsCEDrHGV0uC0Xy1CIsWQ8mKBL7DnYl+4hCDUsZzNZhNsiDj9OPAEmxpbm6IqQ74jSWI9XJghCDWN8VG5mu+fa0CcPkBbWn6MgxFB6GSEINczhp8I5OeTCbQJbhlujABNjfZTQG0smsj7tE9cmsLhPIxMEoYaxZTTXlvWJ1sQIuiIA0mWymOhiew4qusNtEQShXBCEGmZjLRbSHOnP7tyGXLvA+Ch0RQCkyeGns36rRiL4DgShbBCE2sb6Ka41u4ccTjwBJoX1U1yWC0f5rqjNabVSmFxNHghCbXP4s16rzbWjZBRgEnJQL4MRTfJCEGpbLg65VoHFIQeQNoef5jqy23kfbos4ynB6KhsEobY5yuhwaxavCCVR4rtjLNZ5AUibhTLbHFahO4uTXeD0VF4IQm0jC2xJUYqFEll6fb4zShXYzHZ8TgAmIdtDm8Lt6LmXE77gNM9RSnOv0c9iAAARqklEQVRZuygMtwiOcpx4AkyOozyLt2oifTGzxWR3WrP0+gaEINQ8R1kWzz1DzYKzAkEIMDmOcjrcnK0g5NoiDkxwISsEoeaxfjp7tdrhFsFRzmTpxQH0ylFBh9sEKZmVgplwWwQzPckLQah5jjIqnJ2hhMmEJHRF2VJUygBMjpWy2J1WoTsrk11wrQJGNMkLQah5dAkZ5xLZWAuUa4/QJXazDR8SgElzlDNZujsabhVwa1Re+I7TPhPhKMtKh0S4mXdU4L4owFQ4K+hQCy/7y0YH4oREkAVYJVtOCEI9cFYyoSb5DzlUygBMWZbqZUJNgrMSp6cyQxDqgbOCDmXjihBjJwCmyuGnuEBU9vWYws28A6enckMQ6oGzkpE9CMVoMjoQZzCnDMCUmO1mutjOB2QuZAs1C85KBKHMEIR6YHdZzVZTpE/OErVwi8CWYvUlgKmT/e6olJRQKZMNCEKdcFTQoSY5D7kQ7osCZMZZQYda5Dwq+UCUzLdZaay+JDMEoU44K5hws5z1MuFmHpUyAJlwVMhcxRZqFtBBmA0IQp1wyn1FGGzkXdNZGV8QwGhYHxkdiCc42cb4hpp4lIxmA4JQJxwVNNcRkUR5StSE7hhhMmGsEkAmTGaTs4IJnpTtojCMEU3ZgSDUCYvdTBXauQ55StSCjVxeFU48ATLlmsEEG+UJQjGajPTGWB8mV5MfglA/nJV0SKZzz9BJ3jUN90UBMuWazgZPcLK8VKhZYP2UyYJCbvkhCPXDVcUOHpcnCAdP8K4ZuCIEyJSrkuY6IrIMqw8e5/KqcHqaFQhC/cirYgePc0TGR1ycE+PBBIbSA2TObDczHlKWiu5BBGHWIAj1g8y3Wexmviua4esET3DO6QyG0gPIwjWdzbybMJmQwq2Ccxru02QFglBX8qrY4PFMOySCjbxrOo43AHm4pjODJzINwlATz3hJC4lv7KzA26orrpnMoAxByCEIAeTimsGEGvkMV6vHfdGsQhDqypfdhBkQo0k+EMVYJQC52BxWm8vKd2TUZxE8ziMIswdBqCtUod1kMQndU599e7CBc1YyWJUeQEYF1Y7+o+Ep/7kkSqFm3on7NFmD7zu9yfCisP9oqKDGIWN7AKCgxjFwNDTlPw81C7SbtFKYaztbEIR6k2kQHgnnVyMIAeTkqmJDLYIYTU7tzwePc3kzcF80ixCEepNf7eg/Eppaz3ykL5aMJVkv5nACkJPFbnaWT72Qrf9wCKenWSVPECaTyX/961//+Mc/BgYGxtksGo3u27evt7dXlp3CqMh8G5lnm9pKFP1fhAtqnAQGEALILb/G0X9kKt2ECU7kA9G8mbgizCIZglAUxW9/+9u33377Cy+8UFNTc/DgwbG2fPjhh88999zt27dnvlMYR+EcZ9/hqXRI9B9BByFAVhRUT7GbsO+LUN5s1mzF+WkWyRCE27dvb2ho2Lt379///vfbb7/9kUceGXWzvXv3fvTRRwsWLMh8jzC+wrmuvs+Dk/0rSZSCJ3iceAJkA+ujUstHTPYP+z4PFs51ZqNJcIoMQbhly5Zly5bRNE0QxIoVK7Zt25ZIJIZtE4vF7rjjjmeffdZiQeFT1jkr6DgnRvomd8gFG3m6xG5zWLPUKgBDMxEF1c7J3h2VRGngGFdYgyDMLhm+9VpbWxcuXJh6XFFRkUgkAoFAWVnZ0G1+/vOfL126dP78+RO+WiKR2Lp1a2Fh4dAnly5dSpKjzAEtiqIoyrb6s54U1Dh6Dg76LiyceNOvdB8YyJ/jmNT7ifdfWXj/lTXZ9z9/DtvxcZ/7vLz0/2SwnqNK7GbGhP/okdJ8/81ms8k0wY3ltILwwIED99xzz8jnN2/eXFlZGYvFbLYvlzJPPYhGvzaHQl1d3ZYtWz799NN09iWK4rZt21LXl6dccsklo24cjUYlScJV5kiOWWT3J8HCc9K9zyklpZ4DwepVpcP+78YXi8UmtT3IC++/sib7/jMzbPzrkVAXZ89L9wqk5+Cgq5rG//Ko0nz/KYqSJwinTZv28MMPj3y+pKSEIAiv19vT05N6pru7myAIn883dLNf/OIXfr9/w4YNBEG0tLS8/vrrTqdz+fLlo+6LJMkXX3zR7/en0zCTyWS32xGEI5FnUM1v9thNpJVO683p/yJEF5MF/kmcqxIEIYoiw2C2C8Xg/VfWFN7/ojPyuCOx/EtdaW0tEcGjrad9v4JhMKJpFDJ+/tMKwry8vIsvvnis315wwQVbt2598MEHCYL45z//uWDBgmGNW7lyZVNTU+oxSZIlJSXFxcUZtBkmZiHN+TWO7s8GfRekdXe0+7NB91mTS0EAmKySs/JPbgv4L03rC3DwOGchzawPKZh1MvQR3nTTTY8//vi99947Z86chx566Ne//nXq+QsuuGDZsmXr1q278sorT2380ksvXXLJJYsXL858vzA+73kFjVsD6QRhMpbsOxyafrU3B60CMLL8mWwslOA7o4xn4oWvA3v6PecV5KBVIEPVaEFBwZ49e+x2e21t7ebNm6+77rrU8z/4wQ9GBt6dd9559tlnZ75TmFD+LIcYS4ZbJh5Z33so5KykbU7UiwJkmYkome/q2T844YYJTuz/IuRemJ+DRoE8330VFRWPPvrosCdvuOGGkVvedNNNsuwRJmYivIsKAnv6Z5ZPsKZSV21/yVk43gByoeTs/COvtJRfUWIyj1fB0bVvoHCuM80+fsgQ5hrVM/e5BT3/Hhx/qt9wq8B3Rovno4MQIBccZTRZaOv+bIKLwsCePi/ui+YKglDP7E5r3ky2q3a8CWCb3+4q+28lmMAJIGcqrnS37OgaZ2b8wQZOShKu6ZjmKUcQhDpXscTd/G5XnBt92Gm4VeA6Ip5zceIJkDt5Vaw9b8yLQikpnfh7R+USN6a/zxkEoc6xPqp4Xl7T/+0c9be4HARQROVVnuZ3R78o7Pi4z8Za0VuRSwhC/au8yt13KBhuHl4+2rm3n++K4nIQIPdc0xm62N70j+FnqLFQouW97qrlvlH/CrIEQah/Vtoy7dveY39uHTrzfd/hUNNbnXNvq8TlIIAiqm8o7zscatvZc+qZBCce+1OLZ1EB7Z54lCHICEPHDMG9MD8hiAeePjHjv3vJAnu4RWh5r3vOrZV0CY43AGVYGcvcVdPqnjmREJJ5VQwhEfV/bSuel1e5xK100wwHQWgUpRcVuaYz9X9pM1tNrJ8+bWWFs2KC8YUAkFVkvu30VdM6dvW1vNsdHYxXfae0cA5WXFIAgtBAHGX0gnUzlW4FAPwH7SZnLEOPoMLQRwgAAIaGIAQAAENDEAIAgKEhCAEAwNC0HYRPPfXUjh07lG6Fcd1///2ff/650q0wrhtvvHFwcOIFfSAbwuHwihUrlG6FcR05cmTdunVyvZq2q0aPHDni9WI5WcXs379/YGC8Gb0hq3bt2hWLxSbeDrIgHo9//PHHSrfCuAYHBz/77DO5Xk3bV4QAAAAZQhACAIChmSRpzDWxFHHNNdfU1dWZzWkldCgUstlsFEVlu1UwqoGBAZZlbTab0g0xqN7e3oKCgjQPFpCXJEm9vb3FxcVKN8SgEolEKBQqKJh4zYDt27efdtpp42+juiAMhULd3d1KtwIAAPSgrKzMbrePv43qghAAACCXcFMFAAAMDUEIAACGhiAEAABDQxACAIChaSYIN2/e7PP5nE7nsmXLRp3NpK6u7uyzz6Zp+swzz9y7d2/uW6hjL7/88vz580mS9Hg899xzTyKRGLZBPB6vGuJnP/uZIu3Uq/vuu+/Ue3v66aePus0zzzzjdrtdLtf111/PcVyOW6hvs2fPHvrxfuyxx4Zt8Oabbw7d4MCBA4q0U082btx41VVXzZw588033xz6/KZNm0pKSlwu1w033CAIwsg/3LVr19y5c2maPueccw4fPpzu/iQtOHLkiNPp3LdvnyAIy5cvv/POO4dtkEwm58yZs3HjRlEUX3jhhcrKykQioUhTden555//+OOPI5FIQ0PD7NmzN27cOGyDaDRKEMTBgwePHz9+/Pjxnp4eRdqpVzfffPODDz6Yem9PnDgxcoPa2tqCgoLPP/+c47grrrhi/fr1uW+kjh3/ysGDB2ma/vDDD4dt8Morr1x22WWnNotEIoq0U08ef/zxV199de7cua+88sqpJ3ft2lVcXHz06NFwOHzppZc+8sgjw/4qFouVlpZu3rxZFMVf/vKXCxYsSHN32gjCBx54YMWKFanHn376qcvlisfjQzfYs2dPfn5+LBaTJCmZTPr9/nfffVeBhhrA/fff/73vfW/Yk6kgDIVCijRJ926++eZNmzaNs8Gdd955xx13pB6///77brc7J+0ynN///vezZ89OJpPDnn/llVeWLl2qSJP0bdGiRUOD8NZbb127dm3q8VtvvVVeXj5s+23btk2bNi31OBKJOJ3O/fv3p7Mjbdwara+vP3VH6PTTTw8Gg4FAYOgGx44dq6mpSU1xYjKZ5s6de+zYMQUaqneJRGLHjh3nn3/+qL8988wzq6qqbrnllmH/O5C5jRs3+ny+iy666O233x7522EHSFdXF1alyIbNmzevXLnSZDKN/NVHH33k8/kWLFiwadOmZDKZ+7YZwbDPeUtLy7C7o/X19WeccUbqMUmSM2fOTDMItBGE/f39Docj9ZiiKJvN1tfXN2wDlmVP/ehyuYZtALK4//777Xb7qlWrhj1vsVi2bt26e/fut956q7+/f/ny5Yo0T69uvfXWnTt37t+//6abblq2bNnISfeHHiBOp5MgCHz+ZXfixIndu3ffeOONI3+1aNGiHTt21NXVbdy48amnnvrtb3+b++YZwcjPeX9//7ANphYE2gjCoqKiU2e4HMfF4/GSkpKhGxQXFweDwVM/DgwMDNsAMveTn/xkx44d27dvt1qHr95lsViWLl3q9Xqrq6tfeumlPXv2tLa2KtJIXbrwwgtnzZrl9Xpvu+22b33rW3/729+GbTD0AEmVkuHzL7vf/e53V111VWlp6chf1dTULFy4sKSk5PLLL1+/fv3rr7+e++YZwbDPuclkKioqGrbB1IJAG0FYU1NTV1eXelxXV1dQUOB2u4duUF1dfeTIkVRPVTKZPHToUHV1tQIN1a9f/epXr7766jvvvDPskzdS6saRhKn7smPU+3LDDhC/33/qxBlkIYrin/70p5UrV0645aj/QSCLYZ/z6dOnkyQ5coPUl48gCPX19ekGQeb9mTlw4sQJh8Px3nvv9fb2Llmy5N577009/+Mf//i1115LPT7rrLMeeuihUCj0xBNPzJo1SxRF5dqrN08//XReXt62bdtqa2tra2uPHj2aev7uu+9+6623JEn69NNP33///c7Ozvr6+uXLly9atEjR9urN888/39jY2NnZuXnzZoqiamtrJUnq6+u7+uqru7q6JEmqq6vLy8vbtWtXd3f3RRddtGHDBqWbrDfbtm3zeDypcryUtra2q6++OlUg9te//vXgwYO9vb0ffPBBeXn5U089pVxLdaKhoaG2tvb000/fsGFDbW3t4OCg9FV19CeffNLZ2Xneeec9+uijqY3Xrl379ttvS5KUSCSmTZu2adOmUCi0fv36888/P83daSMIJUl644035s6d6/F4Vq5cyXFc6sk1a9a8/PLLqcf19fWXX355UVHRxRdffOjQIcUaqkd33HHH2UOcKlC8+eab33jjDUmSdu7cee6557rd7lSxTFtbm6Lt1ZtrrrmmoqLC4/EsXrw4deYhSVJPT8+iRYs6OjpSP/7xj3+srq72+XyrV69G+b7sHnrooccff3zoM01NTYsWLUp9Qf/0pz+tqakpKiqaN2/epk2bcBaeubvvvnvod87u3btTz2/evHn27NmlpaU//OEPT52X3HDDDVu2bEk9/ve//33hhRcWFRVdeeWVo441GhVWnwAAAEPTRh8hAABAliAIAQDA0BCEAABgaAhCAAAwNAQhAAAYGoIQAAAMDUEIAACGhiAEAABDQxACAIChIQgBAMDQEIQAAGBoCEIAADC0/w/KNlSDHuR97gAAAABJRU5ErkJggg==", "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 kpoint, all G components, first eigenvector\n", "ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n", "ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\n", "\n", "using Plots\n", "x = a * vec(first.(DFTK.r_vectors(basis)))\n", "p = plot(x, real.(ψ), label=\"real(ψ)\")\n", "plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n", "plot!(p, x, ρ, label=\"ρ\")\n", "plot!(p, x, tot_local_pot, label=\"tot local pot\")" ], "metadata": {}, "execution_count": 10 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" }, "kernelspec": { "name": "julia-1.6", "display_name": "Julia 1.6.2", "language": "julia" } }, "nbformat": 4 }