{
"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 a custom\n",
"potential"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"struct ElementCustomPotential <: DFTK.Element\n",
" pot_real::Function # Real potential\n",
" pot_fourier::Function # Fourier potential\n",
"end"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"We need to extend two methods to access the real and Fourier forms of\n",
"the potential during the computations performed by DFTK"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function DFTK.local_potential_fourier(el::ElementCustomPotential, q::Real)\n",
" return el.pot_fourier(q)\n",
"end\n",
"function DFTK.local_potential_real(el::ElementCustomPotential, r::Real)\n",
" return el.pot_real(r)\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;"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"We define the width of the Gaussian potential generated by one nucleus"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"L = 0.5;"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"We set the potential in its real and Fourier forms"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"pot_real(x) = exp(-(x/L)^2)\n",
"pot_fourier(q::T) where {T <: Real} = exp(- (q*L)^2 / 4);"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"And finally we build the elements and set their positions in the `atoms`\n",
"array. Note that in this example `pot_real` is not required as all applications\n",
"of local potentials are done in the Fourier space."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"nucleus = ElementCustomPotential(pot_real, pot_fourier)\n",
"atoms = [nucleus => [x1*[1,0,0], x2*[1,0,0]]];"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"Setup the 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": 9
},
{
"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.293945996322 NaN 4.29e-01 7.0 \n",
" 2 +0.441608793303 1.48e-01 5.35e-01 2.0 \n",
" 3 +0.260909838321 -1.81e-01 1.72e-01 2.0 \n",
" 4 +0.243050951794 -1.79e-02 2.67e-02 2.0 \n",
" 5 +0.242664372678 -3.87e-04 6.06e-03 1.0 \n",
" 6 +0.242635928800 -2.84e-05 2.15e-03 1.0 \n",
" 7 +0.242633958912 -1.97e-06 9.62e-04 2.0 \n",
" 8 +0.242633645891 -3.13e-07 7.31e-04 2.0 \n",
" 9 +0.242633377825 -2.68e-07 4.57e-05 2.0 \n",
" 10 +0.242633376334 -1.49e-09 2.99e-06 1.0 \n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown:\n Kinetic 0.0304518 \n AtomicLocal 0.0972469 \n PowerNonlinearity 0.1149346 \n\n total 0.242633376334 \n"
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"Ecut = 500\n",
"basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\n",
"ρ = zeros(complex(eltype(basis)), basis.fft_size)\n",
"scfres = self_consistent_field(basis, tol=1e-8, ρ=from_fourier(basis, ρ))\n",
"scfres.energies"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"Computing the forces can then be done as usual:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "2×1 Array{StaticArrays.SArray{Tuple{3},Float64,1,3},2}:\n [-0.03868686200624871, 0.0, 0.0]\n [0.038689393066859326, 0.0, 0.0]"
},
"metadata": {},
"execution_count": 11
}
],
"cell_type": "code",
"source": [
"hcat(compute_forces(scfres)...)"
],
"metadata": {},
"execution_count": 11
},
{
"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": 12
},
{
"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+gvaeTAAAgAElEQVR4nOzdd5xU5dk38OuUOWf6zO5sB7aw9N5ERAULiAYFRbGXV9AHS2yJMZYnMWhiRI2iMZbEgEbRR8UeS6woiIACoiCdhYVdYOv0Of1+/zhkXWCBLXPKzF7fD3+ws1Pu3dl7fufuFCEEEEIIoe6KtroACCGEkJUwCBFCCHVrGIQIIYS6NQxChBBC3RoGIUIIoW4NgxAhhFC3hkGIEEKoW8MgRAgh1K1hECKEEOrWMAgRQgh1a10Nwm3btiUSifbfX9O0Lr4iSi9CCG6zZzdYTexGVVWri4AOkt53pKtBOHv27O+++679908mk/ixayuyLMuybHUp0M8IIclk0upSoIPgO2I36X1HsGsUIYRQt4ZBiBBCqFvDIEQIIdStYRAihBDq1jAIEUIIdWsYhAghhLo1DEKEEELdGgYhQgihbg2DMMPEq1NEwx0JEMps8RpBU7Ai2wVrdQFQB6TqpR+fqnLmcX1mlvjK3FYXByHUYaqg7fpw/97lTT0m5pWfXWh1cRAABmEmIbDttZqyqYUOD7txYXXJhLyep+VZXSaEupF169Z98803XXkGOaHULml0F/HBft53n2ko+THE5zjSVbzslpeXd8EFFxj05BiEGWPvN01EJSUnhYACX5lr3eM7ekwMUQxldbkQ6i4++eSTl156ady4cZ1+BqISKIJmhqrZDqQEmrfvpVmswsfW2Ni4bds2DMLuTgzL1R/VDb2xAigAAGeIcxXwzZviuYN9VhcNoW5k0qRJjzzyiNWl6HbWrl07a9Ys454fJ8tkhvrV4bwRAXch33JLwZhg3XdhC4uEEELZAYMwMzRvjucO9La+JW+EP7wlrqTwmDSEEOoSDMIMoEpafHfKX+lpfSPrZIL9vA3rolaVCiGEsgMGYQaIbE34ytwMd+ibhb2jCCHUdRiEGaB5czzY33v47TkDvKk6UYrg+fIIIdR5GIQZILw5ltNWEFIM5S93R3emzC8SQghlDQxCuxOaJFXQPMXONr/rLXXFq5MmFwkhhLIJBqHdhTfFgwO8cIRFt75SV6waW4QIIdR5GIR217w53ma/qM5b6orvwW24EUIdE4/HW3YGWLNmzVdffdXm3RobG1966SUTy2UNDEK7i+1M+nt7jvRd1snwQUdyn2hmkRBCmS4ajf7ud78DAFVVZ8+eXVBQ0ObdQqHQk08+uXbtWnNLZzYMQluT44qmEj54tG15vdg7ihACIIQ0Njaq6oFNNuLxeMv/W25JpQ79rPjggw/y8/MHDBhwpKe97rrrHn744bSX1lYwCG0tUSN4StqeJtPCV+rG+TIIdWe33Xbb9ddfP2zYsFGjRv3www/ffvvtiBEjjjvuuF69ev3pT38CgC1btgwYMGDs2LH9+vWbPHlyY2Njy2MXLVo0ffp0ALjwwgtfe+01AMjNzY1Gox9++OGZZ54JANOmTXvrrbeSyWz+kMFNt20tUdueIHTt+6bJnPIghA7xfSOpNTcjTi+heOagW5LJ5IcffrhixYqSkpJYLDZ48OCFCxeefvrpkUhk/PjxJ5544ogRI77++utQKKRp2o033vjQQw/dcsst+mO//vrrW2+9FQDi8bgoigDQ3NxMCJFlORaLAUBubm5JScl33303YcIEU39OE2EQ2lpirxDoc8QBQp2nxCk0SqqoMTy27xEy23/2kC/3aWa+4vEF7CFBCAAXXnhhSUkJAHz99dcOh4MQ8umnnwLA8OHDP/vss1NOOWXlypVPP/10PB6vr6/ftm2b/ihCyN69e480QNiisLCwpqYm/T+JbWAQ2lqiViiZEDr6fSiGchc543tSgcpjRCZCKO1+O5z+7XDrr0Hz8g4c011XV6coip6CAFBaWjpq1KjXXnvtf//3f++4446ePXt6PJ4PPvhA/y5FUTzPS5Kkf0nIz/PPCSEUdWDZliAILpfLpJ/EChiE9kVUkqqXWh+9dCS+Mle8GoMQIQQDBw4UBOHee+9tHV2XX375Lbfccs011wDAd9991/r+ffv23blz54ABA3r06LFr166W26uqqnr06AEAhJDq6uq+ffua9RNYAIPQvpJ1ojPXQTuOfbHp7elq3hQ3oUgIIZs77rjjTjrppBkzZtxyyy0Oh2PVqlXjxo0bOHDgokWLhg4dunnz5ueffz4/P7/l/pMnT162bNmZZ5559dVXT58+PRQKAcDChQvnzZv3/PPPA8DGjRudTuegQYOs+olMgEFoX+2ZKaNzF/C1XzUe+34IoWw0ZcqU3Nzcli9fe+21559//rXXXpNlefDgwYMGDRo/fjzDMH/729/69u37f//3f6tWrfL5fLfffjsAXHPNNeecc8599903fvz4zz777MUXXwSA7du3f/DBByNHjgSAV1555dprr23pJs1KGIT2lWx3ELoK+VS9CASOtBMbQiiLzZgxo/WXDMPMnj179uzZrW+88847W/4/ZswYALj//vsBoH///ieeeOIbb7wxc+bMYcOGPfzww4888sgf//jHQCAAAPF4/M0331y+fLkZP4Z1MAjtK14rlJx8jJkyOoanGZ4WI/LRl94jhNDhnn32WUVR2vyW0+lcuXKl13vEXR6zAwahfbW/RQgA7kI+VSdiECKEOorjOI7jWr5sPXeUZdmsT0HAnWVsqz2bq7XmKuBTdbjjKEIIdRgGoU0l9gpHOoOwTa4CPolBiBBCHYdBaFOpOslVcOwVhC1c+XyqTjKuPAgh21q1atXGjRvT+IRLly5tvR9pa+vXr9++fXsaX8sOMAhtSmiQnHncse/3X/oYoXHlQQjZ1j//+c/33nsvXc9WVVV1/fXX+/3+I93h8ssvbz2OmAVwsoxNpRpEf293++/PBx1KUsUdRxHqhp599tk0Ptsjjzwye/Zsh6PtCQpDhgzhef6TTz4544wz0vii1sIPTZvqaIsQKHDmc6l6bBQi1O3MnTt3wYIF//rXv+66666LLrooNzd3zJgxO3fu/O1vf1tQUNC/f/9vvvkGAO68886Kiori4uJx48atWLECAPbs2XPGGWcUFRWddNJJDz744L333qsoyssvv3z++eenUqkxY8aIovjZZ5/deOONAHDvvfcuWrQIAC644AJ905msgS1CWyIgNMvOUEeCEMBdwKfqRG/PbN4bFyG7UcMNmpAw8xUdBT2BPuj4iZqaGqfTSVHUk08++Z///OeVV1657rrrTjzxxPvvv3///v1PPfXUrbfeunLlyrPPPnvu3Lk8z7/33nsXXXTRjh07Zs2aNWLEiI8++qiqqmrChAmTJk1at26d0+ksLS1NJBKrV6/WNK25uXnr1q0AsHPnTn2V/dixYx988EEzf2SjYRDakRiWWRfDcB1rr7sKcL4MQmaLf/W28NMqM18x/5cP0d5gm9+aPHny+PHjAWD69OmffPLJrFmzAOC88877zW9+AwBjx479+OOPq6qqkslkU1PT5s2bP/3008WLF9M0XVlZeckll9TX1+/evbuoqOjoBSgpKamtrZVl+UjdpxkHg9COUg2Sq0P9ogAA4MrnG9dHjSgPQuhIAtOuCUy7xupSHKBvmQ0APM+37D7qdDpTqZQgCCeccMLgwYPHjh2bk5PjcDiqq6s5jmuZFFNQUFBfX88wjKa1fbxiy8FMiqLQNE3T2TOylj0/STYRGjs4QAgA+o6jOHEUIdSWNWvWSJL00ksv3XzzzRdffHEsFsvPz6dpWu/2BIAffvgBACoqKvbs2UMIcbvdXq933759Lc/QcoTvnj17ysrKGOaw04EzFgahHQkNUkcHCAHAnc8JDRJk1axmhFB6FBUV7dmzZ8mSJdu3b58zZw7LsgzD3HLLLVdeeeW777774IMPfvnllwAwePBgjuO2bNlCUdTMmTNvvvnmTZs2xWKxJ598cu3ataeffjoArFix4tRTT7X6B0on7Bq1I6FBzBsR6OijaI5m3YwYlvmcLOm4Rwi1x9ixYwsLCwkhLf2chYWFkyZN0v/PcdzMmTN79+79t7/97d5772VZ9pe//GVubm4wGHzggQcWLFjw7rvvDh069Oqrr04mkxRFXX311a+++urvf//7p5566qmnnlq0aFFVVdVPP/30zTff6MOHr7766hNPPGHZT2sE0jUTJkxYsmRJ++8fi8U0Tevii2a9NQ9vje1OduKB6/66I7wt3qGHiKIoimInXgsZRNO0WCxmdSnQQaLRKCHk4Ycf/vWvf211WdKpsbFRURRCSHNz8+DBg998801CSF1dXf/+/ROJhH6f119/ffLkyS0PWbp06ZlnnmlyOdesWTNixIjWt+jvSLpg16gddW6MEACcuQ6hCSeOIoTaZenSpeXl5WPGjOnfv//kyZPPPfdcAMjPz3///fdVVW3zIeXl5frhvdkEu0ZtR4oqtINmnZ0ZiHbmcmKTnPYiIYSy0vTp088888xIJBIKhVpPfqmsrGz5/1lnnXXiiSe2fNmzZ09Ti2gKDELbERrETqyd0PG5jsg2U9f2IoQyGs/z+lzQI/F4PB6Px7TyWAK7Rm0n1dHN1Vpx5nJiM7YIEUKoAzAIbafDu4y24gxxQiOOESKEUAdgENpO57aV0XEBVo4rRMW1hAgh1F4YhLYjNEnOUAeO5G2Noiku4MDeUYQQaj8MQtsRm2RnbudXxOMKCoS6oZUrV5522mkDBgy49tpro1Hcc7hjcNaovWiypgqqw9v594XP5QRcQYGQWWRVVkjbS+4M4mR5CqjWt+zevXv69OkLFy4cOXLkHXfcMWvWrMWLF5tZpEyHQWgvYrPM5zgO/iPvGGeOQ8QWIUJmeWTV35ZULzfzFV+e9mzIldP6lldeeeWss84666yzAOCxxx4rKCiIRqMt262hY8IgtBcxLPM5nZwpo+NDXPPGeLrKgxA6urtOuPWuE261tgy7du0qLy/X/x8Khfx+f01NDQZh++EYob0ITV3dMtuZy2GLEKFuRVGUPXv26P+PRqPRaLS4uNjaImUWDEJ7EcMyH+xiEOJkGYS6nTfeeGPVqlWiKN5zzz1nnnlmMNj2EfaoTRiE9iI2S3wXpowCAOd3KClNU3ApIULdyFVXXTV37tzKysqampqFCxdaXZwMg2OE9iI2y3ywS2OEQAEfZMUmyVXQycWICKGMU1hY+Nhjj1ldikyFLUJ7EZu7tIhQhysoEEKo/bBFaCcEpIjMBboahDhfBqFu5dZbb3W5XFaXIoNhENqIFJVZN0OzXVhFCAA4XwahbmbgwIFWFyGz0QBw33339ezZs2fPnn/4wx8IIQAQjUYvu+yy/Pz8AQMGvPnmm1YXsrsQm7u6iFDH5+BhTAgh1F7s66+/vnDhwqVLl9I0ffrppw8YMODiiy++5557otHo9u3bv/322/POO++4447r1auX1UXNfkJzVxcR6vggK4YxCBFCx9bQ0OB2u91ud+ceXlNTU1RU1Pp0+0xEP/fcczfddFNFRUVZWdlNN9303HPPSZL0r3/969577/X7/aeffvqkSZNeeOEFq8vZLYjNsjMdQcgFHRiECHUrixcvrq2tPcodXn311bq6usNvv+qqq7rS7de/f/+jv+6RrFy5ctWqVZ1+3fSiN23aNHz4cP2L4cOHb968uba2NhqNHnKjdSXsRsSwlJ4WYcAhxxSi4VJChLqL+++/f8uWLUe5w+9+97sdO3aYVp5jeuONN95++22rS3EA29TU5PP59C/8fn9DQ0NjY6PT6XQ4DnwiBwKB+vr6Iz0+HA6fdtppFHXQ/I6bbrrp/vvvb/P+iUSCEHLI/ZEuUZdylTni8TTsFMq46PD+qMN37P4KSZIAgOPSMDaJ0oIQkkwmrS4FOkgikaAoShRFqwvStpdeemnPnj1/+ctfXnnllZkzZ06aNOmjjz568cUXZVmeNm3a5ZdfvnDhwrq6unnz5hUUFFx66aUTJ048/EkWL178+uuv0zR98cUXT58+HQAkSfrrX/+6fPlyAJg5c+aMGTPmz5+/cuVKVVUnTJhw0003tcREC1EUb7755quvvnr+/PkURV1//fUTJkwAgJ07dz788MPV1dWjR4/+zW9+s3nz5o8//piiqMbGxsGDB998883H/Bk1TWv92ai/I+355bjdbpo+xkJBNjc3t+Xwqmg0GgqFQqGQIAiSJOkfjpFIJC8v70iPDwaDn3zyif6jtqBp+igv7PF4MAjbpMb3BYp9Hq+z60/lzOFYifN6jz2jGoPQbvQrRa/Xa3VB0M8IIV6vl+dtuknF6NGjg8HgKaecMnz48H79+r333nv/8z//88wzz3i93htvvLGxsXHSpEler/e0004bOHBgZWXl4c+wcOHC+++//+mnn1ZVdc6cOclk8pJLLpkxYwbDMHfccYeqqlVVVfpnxe233w4Af/jDHxoaGh544IFDnkdRlL///e/bt2+/7777du3add555y1ZsqRnz57HH3/8bbfdduWVVz7++OPTpk175ZVX+vXrxzDMzJkzQ6FQe35GmqZbVwr9Hen8r+xgbN++fTds2DBp0iQAWL9+fd++fYuLiz0ez08//TRixAgA2LBhw+DBg4/yFAzDsCwuw0gDscs7brfggg4xIvsAlxYhZKxtr9U0/GDqQbijftuX8x30kTtw4ECv1zt69OhTTjkFAGbNmnXPPfforbpHH330xhtvvOWWW9xu93HHHTdu3Lg2n3P+/Pl//vOfp0yZAgD333///Pnz+/fv/8033+zevVufR6O3du64447a2tp9+/ZdeeWVDzzwwOFBqPvjH/84bty48ePHr169+plnnhk4cOCwYcPuvPNOAHjuuedKSkr27t1bXl7OsqwePZZjr7766rlz51500UUURT355JN33303z/OXXnrpn/70pxdffHHdunUfffTRvHnzrC5n9lMFjRDCutIz+YoPOiScL4OQ8SqmF5dPLTLzFVn3MT4lqqqqhgwZov9/6NCh1dXViqJ06CE7duzYunVr//79W88mbW5unjp1qqZpFRUVALB///4jPVvLusZBgwa9/vrrPM+3tKbcbnfv3r1tNVoJAOyll176448/6r+C2bNnX3755QAwb9682bNnFxQU5ObmPvPMM/qPjQwlNEtpWUSo4wMOXEqIkAkYngYb9Ji2Hm/Ky8trmdhRV1cXDAZZlj36gNQhD8nLyysoKDgk6hYuXFhRUbFo0SIAWLZs2ccff3ykZ2toaAgEAgBQX1+fl5eXl5f3008/tXy3vr4+Pz8fAPRl63ZAUxT14IMPNjQ0NDQ0zJs3Tx/by8nJefPNN6PR6M6dOy+99FKrC9ktiGlaRKjjcxxiBIMQoe6isLBww4YNerSce+65TzzxRDKZVBRl3rx55513nn6H9evXH+nh55577l/+8hdJkgRBePTRR88777xx48YRQp566in9Dtu2bdM0LRwOa5omCMJDDz10lML85S9/IYQ0NjYuWLBg+vTp06ZNe++99/QsfOGFFyiKGj16dFFR0ZYtW2TZFh9TuOm2XUiRrp5E2BouJUSoW/nd7363aNGiPn36LFiw4LbbbhsyZEifPn1KS0sJIXpozZ079x//+EdlZaXepGuRl5fn8Xjmzp0bCATKy8t79+7dq1eve+65x+Vyvfvuu4sWLSopKSkpKXniiSdmz54djUZ79OgxePDgCRMmlJeXA0CPHj0OnyPSs2fPioqKgQMHnnvuueeff/6QIUOefPLJX/ziFz169HjsscfeeOMNl8t1+eWXK4oyaNCgq6++2qxf0hFRXWycTpw48b777mtzMm6b4vE4zhpt066P6igaSs8oSMuziWH5hyd2HPf7/se8J84atRtCSCKRwFmjthKLxXw+3yOPPLJv375HHnnE6uK0CyFEVdUOzWRUVZWiqEPm/EuS5HA4Wj60RVE8yuxZ/U9XEASHw6Fp2iGvLsvy4Ssu2mPt2rWzZs1au3Ztyy36O9KJp2oTtgjtQgrLfJfPnWjB+Vk5jmvqEeq+KIrq6Hx+hmEOX/nGcVzrpks715DQNH34q3cuBU2AQWgXYjoOYGpB0RTrYaXoMaaKIYRQGvE8/9prr9k28I4Eg9AupHA6xwgBV1AghEzHsuzMmTOPuZOL3WRYcbOYFFG4QDr3JeBxvgxCCLUDBqEtqKJGtLStptfhYUwIIdQeGIS2IEVkLq39ooArKBBCqH0wCG1BDMt8WvtFAYAPcjhGiBBCx4RBaAtSREnjlFEdH2RxcxmEEDomDEJbENO6rYyOD+J2owghdGx4fJItSGHZXZyGYwhbc/hYJakSlVAM7uODUBr06dPn0Ucf/eyzz9r8riZrQFN0uqubKmiMs7u3WFKplMtl4KFyGIS2IEWVYP80vxcUTTm8rBRV0riXN0Ld2bnnntu7d29VVdv87o6394aG+gOVnvS+6Ia/7+x7Sc9DDiDshlofCJV23f2XaxNiOJ1HT7TgAg4pYsgzI9Q9DRs27Ejfor7w9xlX4u2Z5oYLPcDfu6zYV2ZgDKDu3uK2ifTur9aCD7Ai7rKGkCkMq8UOMYy12FgYhNYjKlFTKudNf+tcbxGm/WkRQocgKlGSqsOTzj0xdFiLTYBBaD0xIjt8LBgwo4ULsFiFEDKBFFU4H0vR6a/GGIQmwCC0nhRR0ngAU2tcwCFGsFMFIcNJxvSLAi4INgUGofWM2F9Nx+O1JEKmECNK2jeH0nEBh4SXswbDILSeGEn//mo6LsBiFULIBFLUuBYhbhpsOAxC60nh9O+vpsPRBYTMIUXk9B6j1uJALSZGPDc6AIPQekbsr6ZjOBpoUIS21/8ihNJFNGC7YB3NUjRPy0msxQbCILSeFJE5v1E7G/A4wICQ8aSIzPmN2rmCDzjwJBlDYRBaT4ooBk2WAewdRcgUUkTmg0ZdznIBhxTFWmwgDEKrEZCisnEbCXIBnHuNkOGkiGJki5DFdVCGwiC0mJJUaZ6mHUa9Edg1ipDRFEGlaIrhjarFnB9bhMbCILSYGJV5wy4kQa9C2CJEyEhSRDFoyqiO8+M6KGNhEFrM8CqESwkRMphx28rocIzQaBiEFjN0shnoW9djixAhI4mGLSLUYYvQaBiEFpOiinFrJwCAC2LXKELGMm67YB22CI2GQWgx43Zm0jm8jJxQiYb7UiBkFKO7Rh1eRkliLTYQBqHFjG4RUjTl8LJyDPtVEDKKFDG8FrMerMUGwiC0mEGnWrfG4SIkhIwkRmRDu0YBZ70ZDIPQYkZfSwIexoSQwYw7Sa0F78dZbwbCILQUATlueBDiOfUIGYdoRE6oDi9j6KtwflaKYovQKBiEVpLiCutmKJoy9FV4PKceIcPIMcXhMbwW48RRQ2EQWkkyfmgBcH8mhIwkGrlpfgtsERoKg9BKJgwQAg6zI2Qkky5ncaTfSBiEVjJ6EaEOW4QIGcfoFVA6bBEaCoPQStgiRCjTSVFjd0nUcbhXopEwCK1kTouQdTGaSlRJM/qFEOqGzLmcdbgZTdQ0BTeXMQQGoZVEU6oQAHA+VsZ+FYQMYM7lLFA4xmEgDEIrmVSFsF8FIcOY0yIEPIPCSBiEVpIiCm9aFcIWIUIGEM27nGWxRWgQDELLEJWogurwmhGEPFYhhAygKUQTNYfb2G1ldJzfgS1Cg2AQWkaKKg4fC8buR3GAw+/AFiFCaXdgyqgptZgLYL+OUTAILWPO8iMdjysoEDKAFFUcZtVizo9r6o2CQWgZow/zbA3nmyFkBCki8wGzghBbhIbBILSMOetwdbimHiEjSBHFvFqMl7OGwSC0jJldo9ipgpARpKjCmdkixMtZY2AQWkaKmFeFGJ4GGlQBN5dBKJ3M7NdhnQzRiCpiLU4/DELLiFGZ85lUhQD7VRAygJmXswDgwAXBxsAgtIzJVYgPsHg8L0LpZWaLEAB4PJ7XGBiElpGiMm9iFcIWIUJpJ0ZkMy9ncZc1g2AQWkNTiCZprCkbUui4AIvzZRBKI1XSiAas08RajJezxsAgtIYUMW9DCh2Hm8sglFZmLiLU4abBBsEgtIaZs6512KmCUHqZuYhQh/tuGwSD0Bomj7EDAIfD7AillRQ1dYAQsF/HMBiE1jB5yihgixChdLOmRYi12AAYhNYwc1sZ3YEWITHzNRHKZhbUYpwsYwwMQmuY3zVKsxTN0UpSNfNFEcpi5tdihqcBcIuo9MMgtIYUMftaEgD4gEPEy0mE0kQ0fYADcLDfGPTbb79NHWzNmjUPPPBA61t2795tdTmzjRQ17wymFjjAgFAaSRGZN78W+3GLqPSjzz33XPJfL730UkVFxciRIwHguuuua7m9V69eVpcz21jSIsQzKBBKIylm3qm8LbBFaISDukYXLFhwzTXXUNSBZd6KgtcdhtAkTVMJ6zJvQwodHuyJULooKZVmKIYze3SJx+nfBvj5Xayqqvrqq6+uuOIK/cvnn3/e7/cXFRXdf//9hBxxrqGmadXV1ZsOVl9fb3jBM5n5q+l1XABbhAilhxQ2exGhzuHHNfXp9/MbuXDhwilTpui9oBdffPENN9wQDAbXrFlz9tlnl5SUzJ49u83Hx+Pxu+++2+Vytb7xiiuuuPXWW9u8fyKR0DStpdHZPcX3CqyXjsViJr+u6pATjalDXleSJADgOM7kwqAjIYQkk8mjXH0i88Xj8UNv2S/QHsr8WqxxSrJJMP917ebwd+RI3G43wxyj++1AEGqa9sILLzz22GP6l71799b/M2rUqDlz5nzwwQdHCkK/3z9//vyJEye2s0wURXk8nm4ehIKiuXJ4n89n8utShWxDMn7I62IQ2g0hhKZpr9drdUHQQQ6pOElJ8eS6zK/FpJBuXp0y/3VtKI2/hANdox9//HEikZg6derh94hEIh6PJ12vh+DA8iMrukaD2DWKUHrgAEc2OfBGLliw4KqrruJ5Xv/yoYceGjVqVF5e3tdff/3ss8++88471pUwC0lRxWHuOlydw8vICZVohKK7dYscoa6TIrK70Gn+6+IBFEZgAZyilQIAACAASURBVEBRFJ7nr7322pZbVVV94IEH4vF4r1693nnnnUmTJllXwiwkRRRviQVViKIph4eRY4r5SxgRyjJSVAn2s6BFSDtomqWUlGr+tPMsxgIAy7Ivvvhi61vvuuuuu+66y6IiZT9LVtPruIBDimAQItRV5u+v1kLfQB+DMI1wizULSBFrJl4D7rKGUJqIYWvGCAHX1BsAg9ACUtTs01ta4GFMCKUBATluweZQOi7AijhfJq0wCM2mCCpQB3aRNx9OOUOo66SY4nAzVk064/wOvJxNLwxCs0kRxfyNelvglDOEuk6KWDbMDwB8ADeXSTMMQrNJEcvG2AFHFxBKB6sWEer0KW9WvXpWwiA0m2TFGWYtuAArhrEKIdQl1rYI8RiZtMMgNJuFaycAgPdjixChrhKjls2UgQOTZfByNp0wCM0mWnESYQvWw2gy0WTNqgIglAUsOZK3Bedn5bhCNNyWPW0wCM1mbRUCnC+DUJdJEWu2C9ZRNOVwM3JctaoA2QeD0GxS1LLV9DrOz+IAA0JdIUUt3p4J10GlFwah2Szf4YwLOLBFiFBXWDtZBgC4APbrpBMGobn0DSl8VrYIedyWAqEu0BSiiprDbeVWn9giTC8MQlPJcYVxMhRj5SlIuC0FQl0hRWWHjwVLjzLj/Xg5m04YhKYSIzIftPjkB7yWRKgrpLA9ajF2jaYPBqGpJEvXTuj4ICuGMQgR6iTR6mF+wMvZdMMgNJW1q+l1uD8TQl1hixYhHiOTVhiEprJ2fzXdge1GcTEuQp0iWneeaAs8WDS9MAhNJVq647aOZinGycgJvJxEqDOsPUBGx7oZTdRwi6h0wSA0lRRVeKuvJUFfQYHDhAh1ihiROau7RoECzo/zZdIGg9BUlq/D1eEwIUKdJoUt3iVRh2vq0wiD0FR2CcKgAxchIdQZ+p4YVs/9Bpw4mlYYhOaxw4YUOj7gkLBrFKGOk2IK67J4TwwdbhGVRhiE5pGiMud3WLshhY4LYhVCqDOksA0GCAEAcIwwnTAIzWOH1fQ6HscIEeoUMWrlAUytcQE8RiZtMAjNI9rnWjLgwFmjCHWCFJb5HPvUYrycTQ8MQvPYYUMKHR/EYXaEOkO0wSJCHR/Ekf60wSA0jx02pNAxPA0UKAKecI1Qx0hhW0z8BtwiKq0wCM0j2eDoiRY4TIhQJ4j2WAEFLVtExbEWpwEGoXnEsF06VUC/nMR+FYQ6SIrIfNAW/ToAwOOC4DTBIDSPZIedmf6LD+J8GYQ6TIoolm8X3ILDvRLTBIPQJEQjUswuyycA514j1HFKSqUYiuHt8rHJBx0SThxNB7u8o1lPiioOL0vRNlhODwD6MS44RohQR4i2mSmj4wLYNZoeGIQmsdXQAgBwuIICoQ6SIoqtajGug0oXDEKT2GqmDOC1JEIdZ8cWIY4RpgMGoUlsNVMGAPgAi7NGEeoQW62AAgA+iCP96YFBaBLJBmfTt+bwsCqecI1QR0gRu2w0quNwpD9NMAhNItpmf7UDKOCCuFchQh0gNMt8Dmd1KX7GcDTNUkoCt4jqKgxCk4hhe02WAQBnjkMMS1aXAqGMITbbZcftFrimPi0wCE0iRRRbjRGCfk59M1YhhNrLbmOEgPNl0gSD0BTkv6fy2gmfw2EQItROSkKlaButptfhCoq0sNebmq2kuMK4GJq1y2p6He6yhlD7CbY5ibA1LsDifJmuwyA0gxSWbbWIUIdBiFD72edI3tZ43D0/HTAIzWC7KaMAAMDn4BghQu0lNkt2DEK8nE0HDEIzSBGZs9mUUcAqhFBH2PNyFvdKTAsMQjOIEXvtr6ZjeFyEhFB7ic12DEK8nE0LDEIzSDbborAFH8SdKRBqF9Fmq+l1DE8DBYqAl7NdgkFoBtF+y490fA6OtCPULvatxQGHhIP9XYNBaAaxyY7zzQCAz3HI2CJE6FiIRuSYwgVsN9IP+qw3vJztGgxC4xGQona9lsQBBoTaQYlptjpYuzXcGaPrMAgNJ0Zk1s1QjC2rUJCTcN9thI5Fiij27NQBAD7HIWAQdg0GoeHEsMzn2m6MXcfnODAIETom2d5BiC3CLsIgNJzYZNN+UcCuUYTaRwortq3FzhyH2IzHyHQJBqHhxLDkzLVpFeICrJJQiUqsLghCtiZHVdsGIY4Rdh0GoeHsuQ5XR9EU62HkGC5CQuhopLBiw0WEOi7AynEFL2e7AoPQcPZch9sC514jdEx2bhFSNOXwsVIUB/s7D4PQcDY81bo1LsjimnqEjk4K23eyDADwQYfQhMOEnYdBaDi7B2GAxYmjCB2FKmpEJaybsbogR+TMxWHCLsEgNJaSUgkhrMu+VYgLOkQMQoSOTGySbHh6TGscTv/uGgxCY9l8gBAAnLkOETtVEDoyoUnmcmwdhM5cXErYJRiExhLDsm3XTuj4XIfYhC1ChI5IaJIc9m4R8kFcStglGITGsvNqeh2nH0CBU68ROgKhUeJz7R2EOEbYNawgCLW1tS1fFxQUeL1eAIhEIt9++20oFBo5cqR1xct4Yliy80wZAKBZinHRUtSmJyYiZDmxSfaV8FaX4mj4IHaNdgm7evXqU089tVevXvrXjz766PTp07/77rupU6eOGTNmy5YtY8eOXbRokbWlzFxis+wpcVldimNw5nJCEwYhQm0TmqRQ0GN1KY6G4WmKpeSE6vDYd16enbEA0KdPn59++qn1rXffffett9561113RSKRQYMGffXVVxMmTLCohJlNsPfaCR2X6xAaJX+F2+qCIGRHQpNk88ky8N8dRx0eu1922xMNAKqq/vjjjzt37tQ0DQCi0einn356+eWXA0AgEJg2bdqbb75pcTEzls0XEer4HBZX4yLUJiWhUhTFOO0+nQJ3HO0KFgDq6upmzZq1Z8+e4uLixYsXi6JI03SPHj30e5SWlq5evfpIj5dl+Ysvvqirq2t9Y79+/YYOHdrm/TVN0zSNoux4OF/aEZUoCYX10voVhj1pmsYF2WS1ZOdCdiuEEL2aWF0QBACQbBCcIc7+7wgXZIWmblSL2/+O0PSxL2LYUaNG1dXVORwOWZZnz5594403/vnPf2ZZtuXBPM+nUqkjPV6SpM8++2zdunWtb5w6dWrfvn3bvL8gCAzDdJMgFJtkh48VJdHqghyNJEmUl6QaRUEQrC4LAgAghAiCwLJ274vrJmL7E44ALYqiw2Hrrh3GRyUbUt2nFrf/HXE6ncfMQtblOtCn7HA45syZ84tf/KKoqEgUxVgs5vP5AKC+vr64uPhIj/d4PPfdd9/EiRPbWXpN09xudzcJQqkm7gxxbretx95YlmWKHbXhsM3L2X0QQggh+HbYRFM86S5wuVwum78j3kKlfk83qsWqqqbxhz0oJ7du3VpQUFBYWNi7d+8lS5boNy5ZsuSEE05I1+t1K0Kj7AzZelsZHRdg5Rge44JQG4QmyZmbAbXYGeKERhzp7yT24YcfTqVSlZWV27Ztmz9//iOPPEJR1O23337zzTdHo9Fvv/22pqbm4osvtrqcGUlolDIiCCma4vwOMZwZsY2QmcQmKTTYZ3Upjs2FQdgF7MSJExcvXvzxxx8XFha+//7748ePB4Drr78+FAp9+OGHhYWFy5cv93hsvYbGtoRGKTTEb3Up2sUZcmRKbCNkJr1fRwG7ZwzjpCmWkuOKw4ujyx3Gjh07duzYsYd/48ILL7zwwgvNL1A2yaBo4XM5oQnnXiN0MAJiWOaDDkW0exDCf3tHMQg7we6LYzKa0Cg58zIjCHGAAaHDSVGZddE0lxmfk85crMWdlBlvcCZSBJWoJFN2PMLDmBA6nNAkZ8RMGZ0zxAmN2K/TGRiERhEaMqY5CAe2G8UgROggQqPEZ8joBuhBiLW4UzAIjZJBA4QAwOO1JEKHEZokm58n2hoOcHQaBqFRhMZM6lThvKwmaarQXfZnQqg9hIZMupx15XFCAwZhZ2AQGiWDZsoAAFDgzONSDbbeDQ4hk6XqJVe+rU8ibI0LsHJc0RTcGaPDMAiNIjRmxoYULVz5fKoOLycR+lmqXnTlZ0wtpmiKC+Kst87AIDRKZo0RAoArH1uECP1MSaigQWYty8Nhws7BIDQE0YgUkTNomB0AXPm8UI9VCKEDUvWisyCTrmUBV1B0FgahIcRm2eFjKSaTDtlw5nOpemwRInRAZg0Q6nAFRedgEBpCaJRcGTRTBgBwjBChg6XqxYyrxdg12jkYhIbIuAFCAHB4GKBBjitWFwQhW0g1ZGaLEIOw4zAIDSE0ynxGTRnVufK5FA4TIgQAmTZlVIeHMXUOBqEhMrFFCHrvKA4TIgQABIQGyZlpQcg4aZql5Bj263QMBqEhUnWiqyDDOlUAd6ZA6L+kqExzNOvMjE3zW3Pi5WzHYRAagEAqAyfLAFYhhP4rE6eM6twFfBJnvXUQBmH6Cc2Sw80wfOb9bnGMECFdJg4Q6lwFuA6qwzLvw9r+UnVSJvaLgj5G2CABblWIur1Ug5SpQZjPp/ZjEHYMBmH6ZegAIQAwPM06aTGCO1Og7i5VJ7ryMrIWuwpwgKPDMAjTL1UnujNtZ6YWrnwee0cRStVnbIswjxObZTyDokMwCNMvVS9m6DA7ALjyOQEvJ1H3RjQiNmXUMWqtUAzF5+Bqwo7BIEy/5P5M7RoFAFchn8QBBtS9CQ0SF3DQjkz9eHQVcKk6rMUdkKnvtG2poqYIGh/MpHMnWvMUO5P7BKtLgZCVEvtET7HT6lJ0nquAxyDskEw6aisjpOpEVz4HmXTsxEHcRc5ELQZhZ8gaNAjQIJBGEZpF0iRCWIKwRKISxGWIyRCVSVIB/Z+ogqBCSiUAkJBB0g55Mh7goClLQQ4oCmiAAEcBQA4PDAV+Djws5WHB64AcHgIcFeQgyEEOT+XykOeEfCfly9RLMisl9wru4kzt1AEAVz4fq0paXYpMgkGYZqk60Z2x/aIAwPlZAJBjisOHfxuHUgnsTZJdcdiTILVJ2B0n+1OwJ0HqBNifIjEZ8njIc1IhJ+RwVC4PQR6CHNXbB14H+DnwO2gXCx4WXCw4GeBpcLMUALhZ4FttYEIISSQSXq+39UuHJSAEVAJRmRACYQkUDWIyxBWSkCGuQFiEsESq4xCWoEnQmiWoS0GDQFQC+U6qxAOFLqrYBcVuqpcXStxUqRfKvJQH3+S2JPcJoWEBq0vRee4Cbv+qZqtLkUmwHqRZsj5TFxG2cBc5E/uEoM977Ltmr5gMWyNka5Rsj8L2KKmKkZ1xqE2SPCdV5oWeHqrEDaVeakw+9HDT+S4ocFJ5RvalBf87byPPeUhvwzE6H1IK1AlkbxLqUqQ2CXuTZNk+qE1q1XHYFScuBsp9VIWP6u2DSj/Vx0/1DUAvT8Z2aKRJYq9YekYG12LsGu0oDMI0S+0XQ0P9VpeiS9zFzuReMdi3GwVhWIL1TWRDmGxoJhvDZFMYGkXSz0/1DVCVfhhXQF3ahy73QqmXyrj5Ey4WyrxUmRfajMx6AapiZGeMbI/Bqnry8nZtc5jEZBgQpAYGqUE51OAgDMmlKnxU98lGTSFiWM7oy1mHlwUAOa7o/0HHhL+mNEvVi66MXUSo8xTz8d0pq0thrJoE+a6BrG0kaxvhhybSKJDBOdSQHGpQDnV2Kd0/AKXebvHRn++EfCc1Nv+gnzUqw+Yw+SlMfmomz24iG5qhSSRDc6kRIWpkiBqdRw3JybwLgvZL7RddIY5iMvv9dxXwqXoJg7Cd8NeUVgSEDDzM8xDuIuf+VWGrS5FmCQW+rSff1JGVdeTbeqIQMjqPGhWiruxDjQjRvf3dIvbaye+A4/Kp41qlY1iCdY3k+0aydB95fL22M06G5lBjC6jj86nxhVSZN6t+eYkMnymjc+fzyf2iv8JtdUEyAwZhOgnNEpuZ22235i7mk/sEIMccfrK7JhGW7tOW7CVf7yc/NZPhIWpcAXVZH+rxE+gs++w2WpCDicXUxOIDv7S4DGsbycp6sriK/HqlSlPUSYXUhCLqlBJqYDDjLymS+0R3UQavndC5CnEpYQdgEKZTcq/ozuTlRzrWybAuRmiWnLmZ18ebUuDLfeSzGu2zWrI9SsYXUhOL6UePp8fkUxl4tJxNeR1wchF1chEFQwGA2R4lX+8nX+4lj67XYjI5rYQ+rZg6o2emthQTe4Xi8blWl6Kr3EXOvcsarS5FxsAgTKdEreApyfgghP/Ol8mgINwcIe9Xk4/2aCvqyKg86vQS+snx9Nh8is3sxnlmqPRTlX7qyr4AANVx8lkt+byW/G61GuSoKT2pqaX0xCKKz5yrkOQ+wV2U8V2jnhJnogYXBLcXBmE6JWqFTJ8yqvMU84m9Qu5gn9UFORqVwPL95N1d2tu7iKDC1F7ULwfRb06ivbiE3DqlXurqftTV/YAAs66RfLSH3LdG/bGJTOpBTy+jppbSIXtHjCpoSkLNoEvAI+GDDk0jOHG0nfB3lE6JWqF0SoHVpUgDd5GzeVPc6lK0TSXweS15o0p7a5fWw01NL6NfP50eEcrIXrgsRgGMCFEjQtSdw+lGEd6v1t6tJjctl8cWUOeX0zPK6QKX1UVsS0JvDmbFX5On2JmoFYL9utE6qE7DIEwbTdbEiJyhR7ccwl3srPmiwepSHIQALN9PXt6uLa7Syr3UBRX0N9PY3r6s+MTKdiEeruxLX9kXUgrz4R7tjSpy17fycfnUJZX0+RW0304t+OReIQuG+XUYhO2HQZg2yb2iOz/jlx/p3IW80ChpCqFZ63+cHTHywhbtpW3ExcJlfTD/MpiLhRnl9IxySCnM+7u1l7eTX62Qp/Skr+xLT+lpi6qTqBUyervt1jwlzsj2hNWlyAwYhGkTrxU8Jbbs7uk4mqVcBXyiRvCVWfYTCSq8UaU9t1nb0Ewu7UMvnsSMxP7PbOFi4YIK+oIKaBaZV3do961Vr10KV/WlrhlAW3uVE6tOFYwOWliANPKUOGuX4sTRdsEgTJtEreAusfdMgI7wlrpi1UlLgnBblDy9UXtxqzYqj/rlIPqcMprDyZ9ZKoeH6wbS1w2kfwqTf27Wxr2jjAxR1w+izymlzW8gagpJ7Rc9PbKkRegu4lP1ElGJLdra9oYfMGmTyKIWIQD4erlM3miNAHy0h5z1kXLiewpHw8rp7EdnsudXYAp2C4OC1F+OZ3Zf4riyL/3ID1rvV5UH12mN5q4IT9QIrgI+c8/jPQTtoPkcRxKX1bcDtgjThEByr+DJ/J2ZWvhKXXvMmi8jqPDCVm3+es3Nws2D6bcm07j4vXviGbisD31ZH/r7RvLEBq3va/KFFfSvhtL9Ama0aWLVSW9p9lzLgr6aMItGPY1jdhDWC9Si3RpFUfoZbG6W8rDgc0AuDyEnZav5Yx0ihmXaQWXTkh13kVOKykpKZV0GhlKzCE9t1J7coI4toJ89iZlQhH04CABgRIhaMIH5c4p5eqN68r+Vkwrp3w6nD9kcPO3iu1OBSo+hL2EyT4kzWSvAaKvL0VmiCo0iaRIhJkNChmaJAMAZPehAuufmm/3BnVBhTSMAEP2Q7qSixRWIydAkQqNABBUKXVRPDxS7qVIvlHup3j7oG6Aq/ZTN+8cStYKnR1ZdSwIF3h6u+J6UQecx1aVg/nr175u0aWX051PZgUGMQHSoQhf8YRRzxzBmwRbtos/VPn743xFMy5anaRerTvU4Nc+gJ7eEpyQzNlrbm4QtEbI9RnbGyK44VMfJ/hTUJEhKhTwn5PKUzwFeB+RwFACcVAgZH4TlHvLMiTR1hI15RRXqBbI7ceAc8KoY+aSGbIvC7gTp5aGG5FBDc2FEiBqdR9nt7NCs7H/wlrri1ekPwnoBHv5BXbBZu7iSXn0em6E7UiLTuFn45SB6zgB60Tbtf5apxW6YOyr9cagIqhSR3YXZM7oBdt1oLanA2kbyfSP5vpGsbyYbw4RnoK+f6hugyr3UaSXQy0MXu6HYTQXNWpVtr648noGeHqqnBw45+EDWYGuUbGgm3zeSf2zSrm8gFAVj8+nxhdRJhdSYfOvbi4maVGhYwOJCpJuvl6t+bSSNTxiR4OEf1Gc2apdU0utmsD1sdjWD7MxBw//rR1/Rl355u3btMrXUA38awxxfkLY/ofhuwdPDRdFZ9Tepb7QmRRXOb/FH/c4YWbqffL2PrKgj26JkcA6lN2mu6ksPyqFyrb78sFcQHomDhkFBalCQmllx4JbqOFlZT5btIzd9o22NkBMKqFNL6Ck9qREha06Bie5KlU8rsuKVDeQtdVW9uy8tTyWo8NcN2iM/queU0mtnsHZr0KNMwVBwRR/6kt70C1u1mZ+pY/KpB8bQA9LRrx63aLGQ0Xyl7tiupCV7INcL8EmN9mkN+WIvEVVyUiF9UhE1uz89PGR90+UQmRGEhyv1UqXeA7kYkeDLvdpnteSSz7WIRH7Riz6njDqjB+0264cTm2WikizYqPcQzlxOU7p6OakReHm79r/faaPzqC+nsmn5zELdHEvD7P70ZX3ov/2kTXxfmVFO3zuKKepaisWqU/kjs61TBwD8ZS6Tg3BdE3l3F/l3tbYlQk4toSeVUHcMS8/FinEyNQhbC3AwrYyeVgYAUBUj7+8mT/2kXbVEPa2EvqCCOqfM8M0MY7uS/vLsPAna28sVq06GhnSyFi3fT25doTIULDqVObHQ1jXBUKqoyTFFTqpKUlUFVUmpiqBpkqbJRBFUIKCKKlF/vr+iKCzb3PIl46QpmmI4mmIp1kUzHMM4adbFsC6adbOsh3F4mCzr02sPJwO/HkrP6kf/6Xt16Bvyr4cytw2hO33eU3x3qmJ6tnXqAICv3F39cb3Rr0IAVtWRxVXa4p2EpWB6GTVvLHNiIZUpazKzIQhbq/BRvxxE/XIQHZbg3V3a/+3QblyuTulJX1ZJndXLqJWy0Z1JX3kWdqoAgL/cHa3qTBDWJslvVmpL95EHx9KXVHaDD2kCUkwRGiWxWRLDshRRxGZZiitSWJbjCsVQDg/LehjWzbBOhnUzjJNmHDTrpvlcB0UBw9Mt238QQkRRdDp/nnulChrRiKoHZ1ITw7Ka0hRBVZKqklTlpCrHFYeHdXhZLsByfpbP4fgAy+dwzlwHn8vZYcNY4+Tw8MjxzPUD6d+s0gYtVh4dR08v63A9F5qkrOzUAQBfqTtRkzJuf5nNEbJom7ZoG+EYmFlBvTOZGZabeX9v2RaELYKcvuE93STC61Xawz9q/7NMvawPPbs/PSjdjfTYrlT52YXpfU6bCPT17HijtkMPkTWYv1576Ad1zgD6HyezpnVQm0kVtdR+MbFfSNVJqToxVS8JDSLrYvgQ58xx8EGHM48L9PFwPpYLOBxepkOXYISQRCLh9XZksi4BOaHIMUWMKFJMFpvl+O5Uww9RsUkWmyWHl3Xmca583lXAuQud7iKeD2bsit0jqPRTb05iPqslNy1X/75Je/wEpo+/A9U8sjURMGaZkOUYJ+3M5RK1grdXOi/W4zL83w5t4RZtR5Rc2od+/XRmVF7m5V+LbPyUOlguD3MG0HMG0Fsj5Pmt2uQP1D5+mDOQviBNe3dpCknuFXxp/SOzD18vl9AoKwmV9bSry2nZPnL912qpF76Zxnbok8jmxLAc35NK7BESe4VEjSDFFXcB7y7kXYV8/qigK59z5nGMhRMAKHB4WYeXdRcf+h2iETEsCw1Sql5K7hebN8aT+wRNIu5i3lPi9PRweXs6PcXO7NiO8vQSat0M9vEN2gnvKjcNZn47rL09peGtiWDfrFpK35qv3B3blUxXEP7YRJ7eqL26Q5tQRN85nD6rJ81mSP/nUWR/ELboG6D+NIb5wyh4r1p7eqN2+0p1zgDmuoF0Ydf+PBJ7Uq4CPls3xKQYyt/bHd6WyBt+jN7RsAS/Wal+tIc8No6+oCLjfxuapMV2p6JVyXh1KladJBp4e7m8PZz5o4LlZztdeVymnN1K0ZQzl3PmcsF+P9+oJNVErZCoFaI7ErVfNYpNkrvY6St1+crc/go3n5PB7UUHDbcPpS/qTd3yjTb8TeXZk9qx4pBAZGu87BfZcKR2m3xlrvCWRPFJXXoSjcC71doT67XNEZgzkP4hu5Y/daMg1Dlo/UQ0ekMz+esGbdBi+YIK+tdd2MwwuitrBwh1gb6e8Nb40YNwcZV2yzfaeeXUhgvYzN0nTxW16I5EZHsisj2R3Cu6S5z+Mnf+6EDv84ozOhsOx7qZQB9PoM+BNpAqaYk9qdiuVMO6yI539tIMFaj0+Cs9gUpPhh403ctDvTmJeXeXdsUS9axe1LyxzFGWZif3izRHZ+UAoc5f7t79Sefny6QUeGGr9uh6LcTDLYPp8yuyZlvyn3W7IGwxOId65iTmj2OYpzZqE/6tTCii7x5Bj+j4iXexnanQEJ8RJbSJYF/vvuXVR/ru/hTc8LW6KUxeP50Zn4HzQolG4rtTzZvi4S3xRK3gK3X5Kz3lZxf5ytzZPcekNYaj/b09/t6eHgAAIDRI+tXA7o/rgIJgP29Of2+wn5d1Z9hW6NPK6FNK6DtXqUPfUJ4+kTm7tO03NLw1nsX9ogDgyueVlCrHFIevYx/4MRme3qjNX68en08vnJDNE7+7bxDq8pzw+5H07UPpv2/SzvlYHRmi7h1Fj+7IqG9sVzJbZ8roPMVOJaWJzfLhraJXtmu3rVBn9aNfPpXt9LR1SygptXljvOmnWPPmGB9w5Azwlk4p8Fe4s/Bat+OceZwzjys8PgcAUnVieEuibnV422u17iI+d7AvZ5Avg3YT9DvgqROZC3uTa5aqr+2gnhjfRtMwvCWeNYfxto0CX5k7urMDqwmjMjy+XvvrBnVKT/rTX7Bpn2BoN909CHVuFm4dQl8/kH5us3buJ+roPOr+0fTQdkwCFptlTSXOUNZ2qgAAUBDo4wlvixcel9NyW4MA13+tbgyT96ewHbpusJYcUxp+jDb+EI1VJwOVDTS2GAAAIABJREFUntzB/oppRZbvPmVnrgLeVcAXn5RLVBLZnmjaENv4z2qgIDTUnzfM7ytzZ8RY6SnF1LoZ7F3fqsPeUP5xMjOl58+FJhqJ7kj2vaiHhcUzgb/MFWtfECYVeHyDNn+9emZPenl2TXk7CvwI+BnPwI2D6Nn96Wc3amd8qEzuQc8dTVf4jvZ3EN4SN+hwBlsJ9vNGtiRagvDf1WTOMvWyPtSLp7AZcXCgklAb1kXqv48kaoTcQb7i8bmDZpVm6/wmg1AMFeznDfbz9j6vOFErNP4Y3fparSqoecP9+SOD9j/Gz8PCEycw55aR2UvVM3tSjxzPeFgAgPhugc9xZNMBam0K9PMecx2UrMFzm7U/rtVOLqKWns2acwakTWT5298JTgZuGULP7k8/ul477m3l//Wj7xnB5BxhT9jmTfHcwdk8QKgL9vVUf7QfCMQV+PVK9dMa8uppzEm2PztQU0jT+mjdd+HojmTOIG+PiXnB/t7uM/JnHE+J01PiLJ1SkNwvNnwf2bxoNwDkjw4WjAnafMrJaSXUuhnszcvVUW8p/5rIHF9Ahbdk+QChzlfqEprlowwTvrVTu/NbrcIH/57CjOz4VIlMh0HYNq8Dfj+SnjOA/sMadcBi+X9HMNcNPHT8iGgkvC3ee8ZhS7eyjjPEObzsyu8TV2zjTy6kvp/B+uw9iTJeI+xf0VT/fcTbw1VwXLD/Fb0YHtt/6ecu5EunFJROKYhVp+pXh9fN3+Eu4guPz8kbHrDtBYffAc9PZN7cqZ37iXLDIGbqukifblCFKZoK9vU0b44XjDl0NHRNA7lthRqR4K/jmTN62PRdMxoG4dEUuuDpE5lfDqJ/tUJ9eqM2/4SD/lDiu1N8wMF1cCJWJlIJbCj0b/wo/OeLi8+38RpBVdLq10T2LW9Skmrh2ODIX/fJvi1U7MlX6vKVuiqmFTVtiO1b2bzj7b0Fo4LFJ+a6Cqw+X+cIZpTT4wqoX3+QGtWo5ua5LTiawXTB/r5DgrAuBXd/p36wW7tvNDOrH90NNkI8ouz/EO+6wTnUf85i/11NbvxaHZpLPXo8Xe6jAKB5UzxnQPYPEO5JkCuWqHku36+EXePKSqwuTttS9dLeZY11a8KB3p7yqYXBft6MmMSRZSiGCg3zh4b5xWZ534qmH5+qchXwJSeFcof4bLgneImbetAf/7bSd/x7yuPjmIsr7XuFlxY5A7y7PtgPBIACRYO//aT96Xv1//WjN810ZO7a33TBIGyvs0upyT3YR9drx72j3DKY+c0wunlzvOysbF44AQBv7dSu/1q9ZQjz22GuH6od0R3JllXYNhHemqj9siFWnSo8PgebgDbB5zjKziosPaOg8cdozZKGqnf3FZ+cW3R8LuO0V9g0fh+ZcmGPD73spV+o/6khfz2B8Wbvnw8fdDg8TLwm9T3rvGG5WuyCr87GY9EOYAGgqqpqyZIlkUhk5MiREydOBICNGzeuX7++5U5Tp051u7PzmKEO4Rm4azh9WSV1yzfauNfFv9aKgd5Z+2tJKfCrleonNeSdyax+DnjeiED99xGbBCHRSMP30Zov6jWFlEzMG/D/Sm07KNVtUQyVNyKQNyIQ352q+bLh2083F47N6TExzybrVZL7RFXUfGXuURSsPo+9ebk6+m3llVMze/Poo3P28T73UfRRj+PRrNgHMY3Y9evXT5w48eyzz87Ly3vssccmTZr0z3/+86233lqwYMGoUaP0O5122mkYhC1KvdRbk5n/fBb73u16Zqn2l+OZLu5WakPrm8kln6vDcqk15/28ZVreCP+6+TsqZxRb282lKaTu2+Y9nzdwAUfZLwpzBviwF9TmvL1c/S/vJTbLNV82rHloa94wf8/T8y1ffdvwfSRveED/4/Gw8M8JzKs7tLP+o9w5nLl1iP16cruGADy/RXt1v+u6SNNPVxVmccO3c9iKiorq6mqPxwMAN9xwQ58+ff74xz8CwOTJk59++mmri2dfvesSA07x73bCsDflB8Yws/pnT815eqN272r14eOZq/oedM3ozOWcuY7ItkSwnzUjo5pC9n3TVPNFg6fE2e/Snv4KvDjLJHyOo/e5xaWTC2qXNa57fEfOAG+vSfkWzqap/z7S/9KerW+5qDd9fD512RL10xpt4QS2IFsucDdHyJxlakqBZ873iY/vdSoqODJh/a+JaI/Ho6cgAHg8HoqiKIoCgN27d7/44ouff/65oiiWltCONFlr3BAtHumfN5b59Cz2uc3aKf9WNoWJ1eXqqiYRZnyqPrdZW3YOe0gK6grH5uxd1mR+wYhK9i5r+u5PWyJbEwNnlQ66pgxTMEOxHqZ0SsGYu/u5Cvgfnqza8vIeoUEyvxjhLXGKpg4/majcR305lR0Zoka9rXxck/E1WlRh7hrt5PeUCyrob6axIwuZ4ABvw7qo1eWynYM66++8884LLrigqKjI6XSqqvrFF1+sWLGCZdklS5bk5ua2+XhBEJ577rlPPvmk9Y0nnHDCpEmT2ry/KIosy+pZm7mafoh5evKEU0VR7eeBz8+Av2+hTn5Pu64//GYIyaxdNyVJAgBCyFf7qVnLqAvKyQsnEo5WRbGNOweGuXZ+tD+yJ+Y061ACopHGtbGaTxvdRXyfK4o8PZwAILZZuGyhn1DvcGR17xUFBSf7Q2O9+5c1f//49pxBnpLTQ1zAvLHD3Z/XFZ0cFKW2/5B+NxQmFFCzvyQXlJP7RhKOBlEUOc7WewUcbnkddcMKqp8fVk7VStyaLAEA5Az37PuqKWdkxl9Htv8d4TjumInz81X/Aw888O233+rdob/61a8+/PDDBQsW/Pjjj6FQaN68eUd7CpqmDtOe8mWuxu9joZE/Lz2iKbiuP1l1NvmhmTr+fXp5XYb9+LIGv19LXbWUevoE8uBocpStx2gHXTAusG9psxnFItC8Pr7h8eqGNdE+lxT3vapET0GUNRieLjk9NPTXZayHXf/Ert0fNihJ1YTXTe4VU/ul3OFH6+GfWEhWnUOq4tTEj+jNERMKlU5hCX65grpiKfWHEdprp2glrVIv0NedqpPEZtm60tnRgUuwxx577Pnnn1+yZEkoFGr9bYZhzjzzzGXLlh3p8U6nc9asWfpc0/aQZZnn+YxOSjmhxncKA68qO+RE8goe3pkCb+3Urlym/aIX9eBxR9yYzVa2RKmrl1E9vPS685n8dgRNr4mFq/+8pffZTEePdOmQaFWy6r19mqz1PrekOyzWbI0QoigKz2fCX0868DxUTnP3OrWg+j916x+r7nFKqOTkkKHHgOz6ur5kQp7TfYw/92Ie3j4D/rFJm/Qx3DXYfetIPiM+tl6v0m5boU0vozZcwATaajLljwxE1qd6TcrsaiVJUhrrCA0A//jHP5544omPP/64pOTAcmlCfu4cX7ZsWWVlZbpeLws0fB/JHeRjjtBuOq+c3nAB66Bh8Bvyy9s1k8vWIQTgiQ3a6f+hru5L3pncrhQEAIeHyR8VqF3aaFCpUvXSxoXVm1/aXXxi7shf9eluKdhtcT62zwUlw26qiFenVj+4tW51GIwZoRPDctOmWNEJOce+KwAAXDuAXnYOu6iKOesjpSZh61HDqhiZ+h/lvjXa66czfxvfdgoCQMHoYP3qsLlFszt27dq1c+bMGTdu3B133KHfNHfu3BtuuKGysjInJ2flypW7du169tlnrS2lrdStDpeeUXCUO/gd8OR45oo+9Jxl6vNbtL+NZ/rabx/3qhi5ZqkqqvDlWaSyg9uG///27js8iutcGPiZ2d67tqsLSahLgBBqgDAGC0JJbMd5cuNrJ44T248T34s/f7lOnGLf2DeJE/vGcflcE7fYuAWMiUHCRoCEAIEKEgh1rbRd2+vMzsz3xzqyQgxIaKVdac/v4Y/d2V3Ni1Yz78w57zlHWy/venpYXSOL7YCwiJ8YP2S1nXNr18tz/00PxwUmIY6ClffvqZ6RwMg+s/HoVMYOlSgrxuNWx/9uVVVK57RsygoRcrgBe2aYV/5R5MlK2rezE24EHkaC3/eQT/YQ/1lE+8+ia9xOC1K5FEn5JoJ83XKpi503uk6ne+edd2ZuksvlTz/9dHt7u9vtvv/++xsbGzkc+Pv6gt8UCjsw8YprH5yVKciZnfRn+sjq/ZG789GflNC4CTGMGJAU+FMf+eg54v+U0B4oRAl8zjV7bBlTvU469L4x/47UmIREEZTx+NREs11eKix/KIfBW1LlRlCsCTO4Jfdn2jrdA29P8rTs9G0qToyKs1yXfO5Bf9mD2XP9IB0FPytDt6Uidxwl3hkmn6+maXmJcqHWbKTuayVyhMipHfSrLxv3BQSkVIjNbc7sm+GJ/QvIzFbQ61BfX/+rX/1q9n2EPp8vOkhjPjuNo/7XDXwdR7tBPvuPmAJgTztx3EL9Zg16a2acryV7HNQPThA0BLxUS4uuNxatGp1rRRxFUOeeHErdrJCXiuYVEAXs3e7Rjy1cFSt9m4qrTJaOsaugKMrv9/P5sE0YkBHK2DI1+bldUSZK3ZxCn98VEoGR534zmPWN6+l19nq9AoEAAICT4PEu8ple4pFy2j35cZ6oesRLPXiKPGen/rAW/VraHM4tuJ/oeOJS+Z5spmipFidPfyMxsdg3KSELNvg3C0AAnUNDaAiNidK5NAaPxuDTGQI6S8yI/kvMuUKC1rBrwJ99y9wWs1ZzwZsbaMfM1I/aiD/1kb+rpK1RxOG/58XBo+eIPw+Qj1XQvpc33+MXoSE5t2ovvDomyuFf9w2ce8g/+rGZIkDOrdoEmbkNSigoHdFtlCsrJeOfWjueuKRdL9fUyq57ReWxAxZRFm+evc4MFDxSht6cgfzwBPGXAfLZatqqeEzJ5sbA/3QRL/aTPy6kvbGeNtf1sRk8mnK1ZOIze+bOBF2CKhIkwg487MYxN455IxE/gQcIEiOJMJnzTW3MpxRe7ETIlDKyblYjFBIJERRBERgZ/R+GHJhnNIC58ZADj/gjLCmTq2JxU1g8DZunZXPkrERIjYZmm6ZWdn0r29WqkDM76a8NkF9vIqqVyKMV6KJ1HJIUeG2A/NkZ8kYd0rObEav5MgRpHOUaSd9LowV3p891oXrfRHD879aAOZy2NUVRLk6ELxdKWAweLWu3WlMrG/vE0vHEgH6TQlkpQWhz+6MxHZ9y9HlL/yM2dX/5YuSzRvpfBsgdh4jNOuTXq2jqxRqYFybACxfJX3cSjXq0cxf9ultotRvkZ/9nQN+gWNDy71mKBAjfZNBvDAXM4YAlHLSGKYJiy5gsMYMppDOEDLaMydfTaCyUxkLpnNh3nSz2rwBlIDwd5+pNo2SECtmxgCUctIRsne7RAxbcHxGkcgVpHGEGT5jBjcsiqyEH5uzzZj58/RdQKALuXIF+MxN9upes3h9pTEV/WopmCRcwCVAA7BsjH+kgBQzwt82xv3RN26okwmTvC3PIhb7xoKHZ5h0P6hvk+XekzvV0BiUtjoKZd7veNx4c+9Q6ccSua5CnrJbMsqLKeGzK2DJVdE9GDM+hCAC356C709FfdxJF7+N35aEPFtOkC9m0j5Hg5X7y8U6yVIY03UQvlMzr2GEK6CkV4smj9vRtqlhFOAcUCJhD7qGAdyzgGQvgvghPw+Fr2YI0rnKNhKNgMviLmpuWRh9hJEB4xwLesaB72O8zBLlqtmQFX5zLF6RxFm0C6EtvTrBkzLQtV6sXnT03Bp46T/6pj9ioQR8sRitinaJICnw0Rj7RRWIEeHQVuj31ipcO19dH+CUKDP/N5B7yZ39DI0i74lUxRVCOPq+xZSrkwLT1clWVZEEHii1psI/wmjyjgYkmm9cQVK+TKislV2koi4QIw6e2qfOeonsyWJLrb0+7eo+UwU89do78YJT8QT5638rYz8LvxsDzF8g/9pGlUvCLiphd0YZd+LknB8v+I3s+v5m57tF50ee65HMP+OhcuiiLK0jnCtK43JQ5t/nFto9waSTCmcgI5R0NOPt9rn5fyIlJcgWyQoEkX7Cgt4lTPZ7Rjy1l/5l13f0TX8mHg5f7yT+cJ9VccHceeksmOv/KUnsIvDFI/qmPTOGAB4vRHWnXuFKYbyIEAFDAesY1dtAizOCqqqR8PWf6u8C8keiXNdXl5qrZqiqpvESYgGu0JhSYCGcpYAmbjk3Zu9w8DVtWJBJmcrkq1hd/XRQIWMPOi96JI3ZZgSBtq3KeDYCzOe0Oe6nf95BvD5G70tG7ctHo4mXzdNpGvXCR/GCU3J6K7ilCi6QxPnYmmm2uS/7CH6QvYPcEBXwTwakej6PPi7lxcZ5AkssXr+DPc/BVsifCmTA37ujzTp33ekcCwiyuvFQkKxTGPCPivsi53w3m35F6lTue+SAo8ImB/H8XyeNm6iY9enMmskmDznWdFDcGDhrI90apI0Zyeyr6g3y0anbHYQwSIQAAABIjJ1umnH1evzHE4NMBAiIBAiBAmM4VZvEUZSK4ZO4swUQ4J2SEcl7wOvq83tFAyIkzuDSEhkQCBJ1LE2ZwtevlPE0MZuab/WnXHgIv95OvXCKZKLgtC92RhhTMsQ2TAqDbQX0wQu4doTASfC8XvWMFukBrvVEk1fPMiLxMpKmVXfvdc/vRwDsetHe67N0elIHKioTSAoEgNWZteDARfoVIiHD2em2dbs9QQLyCp6gQS/IFsRmRTYELr45zVay0mxZ8MXpbCHwwSr43Qp60UmUypF6NlEiRMhmSykf+tR0RJ8GIl7rgok5aqeMWqnuKqlejO9ORr6ejV5pR4ivFKhFOowgq7MIBADQ2jcGlwUKYuYKJ8LpFi+8okqKx0Nh2Ms31tEsBcNxMvTdC/m2MoiGgRoXUKJESGZIjRL5y5kV7CPQ6qc4pqt1GHTGSAgayIw25JRNdrVjwoWYhO9b1v8PF92XEakmsgDls7XDZz7lRBiIvFclLRFxV7PtOYSK8mkiQmOr2WDtcflNIXixKWSUWpnOv/1xMgcH3jQFjqOjejMUs6whEwHELdcJMdjlAt4Oa9FNSFkjhfJEOAxHgwihnGGh5SK4IrFGg1UqkSonwruvAj3kihOYJJsIENJ/Tbq+TOmGhWi3UeSd1yU2hCJCxvkyHzjAwBiguHeSLkRIZskqObNQgafxFvX60nHYaPrUV/iCdLb/+8wDmjdjOumxnXLifUJSJFBXimNyLX8nSHke40OgcmrJSoqyUhF247ax7cK+RxEhFhSilQjzX6x2KpAbfMYYcWMHd6Ytc3Milg81aZLP2iwo3CgBzANhCFE4CAACHDsRMIGcjMe2vhCBoQRRIkAIJ8v28L546wsAZplz/mNBJwgJqDsKJ65lYuVpCEaDn2ZGCu9PnOqkFESanejy2sy7vWFBWJMzYoRZl8ZZcO9ByS4TTWGKGbqNct1HuN4asHa6e50YZfJqiVCQrFnIU1/6mA+bw8IcmgICCu9JiWyBzHRAA1Fyg5i61Py4Igv6FlAWkibeOhWqtBGUgPc+OpG1Vqiol18xkRJh0XvTZO12ufr8wi6tcI8m/I3XploIv20Q4jadhZ2hUGdtU7mG/vct9/rlRlIFI8gWiTJ4gnXtZ5RJFUN6xgPWse6rbo79BoV4nhQPdIAhKBikVYq6SNfyR2XRiSrdeLs7lX9bPSkYo30TQOxJw9vu84wFhOldeIsq+RbsQI9wX2fJPhF9AgCiLJ8riZe0CPmPIddFrOeUc3DsJEIQlZjAEdIqgIiEiZMXYcqZkpaDi/+bQuUv+24UgCJo9vo5TfF+Gvdtj63ANfWBiCukMPh1lokSYxFw45otwlSxhOlddLc2/IzUuE5sskKRJhNMQwNey+f9Y6xz3RsJuHPdGEBpCY9PYMiZc+gCCoGQmLxbKi4UUSQXM4UiAIHESZaIsMYMpYizXxdGSLxH+M4aAnghT7UEQBCUUBEUWtOwzoSyfe1sIgiAIug4wEUIQBEFJDSZCCIIgKKnBRAhBEAQlNZgIIQiCoKQGCyYhaLGFIiGcjAAA/HiApMh/bAzjJA4AoCgqGAxysS+WOhEwv5x0lEPn0FEaAIDP5CFLbhorCEpUMBFC0HXCCdyNeT1hryfs8WI+D+bzYj4f5vfjgQAeCEZCPswfjAQxAg/gwWAkGCGJaOZj01kMlAEA4DK4NOSLVpnpjQAAkiRR9IvtXsw3vccAHiQoYnojl8GhITQeg8ukMTh0DpfB4dDZHDqby+DymTw+g8dn8vhMnpApELD4IpZQxBLyGAuylBgELWkwEULQFXnCXnvQYfHbnCGXLWB3htz2oMMVcjlDbkfIiRO4kCUQsoRCJl/IEgiYfAGTz2fyJGw1l8Hl0NkCJp9DZzNpzGiKoqP0aOq6+k5nv/pENC/6MD9O4sFIyI8FQkQoiIcCkaAX8/kxvzVg92E+D+bzYF53yOPBvDiBi1hCMVss50jFbKGMI5VxJDKOVM6RKbgyGUcaveOEoKQCEyEEAWfIbfKZTT6LyWcx+S0Wv80asFv8VibKlHGlSq5CyhEruLI0ka5CVSxmi0UsoZQt5jN58Q2by+CAf247vSacwF1hdzSRu0Iee9Ax6TV3WXungg6L3+4KuUQsoZKXouIplLwUNV8Z/afiKWGChJYxmAih5BKKhMY9kwbP5KjbMOE1TniNk14TDaVp+eroST9Xml2vX5fCU6h4CjZ9uc2swaAxFFy5giv/yldJinQEnSa/1eK3mv22fsfg5+MnTD6LLTCl4Mp0Ao1OqEkValOFujShTs6N9ZrmEBQnMBFCy1mEJEbd48OusWHX6Kh7fNRtcIZcOoE2VajVC7XVujU6gUYrUM/ppmoZQxFUzpXJubIiRf7M7RGSMPstE17juGdyxDV+dLx1zDMRjoTTRPpMcVqGKC1TnJYlSRexhPGKHILmAyZCaFkJRcIDzuFLjqEB5/CAY9jgmVDzlZni9CxJxvbsLRniVBVPiSKw3nJu6ChNJ9DoBJq1mlXTG72Yb9Q9PuIaH3aNtRhah1yjHDonR5KRI83MkWTlybKvdN8JQYkGJkJoaSMpatQ93mfv77P3X5wamPSZ00X6FdKsAnnuzpytGeI0Fo0Z7xiXJwGTX6RYWaRYOb3F5LMMOocvOYYPDB36/ennAAB50px8eU6BPC9PlgMLVqGEBRMhtPSEIqE++6VuW+9528U+e7+MI8mX566UrdixYmuWOAOWdcRLtJO1Vl8VfWoN2C9ODfTZ+//c89cB57CKl1KoyC9S5BenFKh4KfENFYJmgokQWhpCkVCP7UKnpafTen7QOZojyShSrNydu+2R6j1CliDe0UFfIYUrT+HK6/RVAACCIgYdIz22CycmTj137jUGyihVFpamFJarimFShOIOJkIocZEUeWHq0hlT1xlz54BjaIU0q0xZ/L2Sf1spz4UNnksLDaHlyrJzZdnfyNsOABj3THZbz3eYu17sep1FY1aoSlapSitUJfCaBooLmAihhGMPOk5OnjltOnfW3J3Ck69Slf5b4c3FigI2nRXv0KDYSBVqU4Xabdk3AgBG3YYOc+ehkc9+2/5MqlC3RlNWqanIl61AETgTMrRIYCKEEgJJUf2OgdaJ022Tpy0B22p1WbVuzY9X3y1hi+MdGrSw0kX6dJH+67nbIyRx3tZ3ynTu96eeswWmKjXl67Rr1mjKYZUNtNAWOxEGIsEJhwn5R/06A2Ww6SwWjcWkMeBYriSEE/g5S8+xiZMnJtoFTP467Zofrf5+gTwP3g0kITpKK1UWlSqLvl/6HVvA3jZ55u/DR37b/ky+fEWtbm21bg0cj5GEQpEwRmI+zE+QRCASBABkiFKZse4ZWexEaA86nup8YfopTuKhSDhEhDEC82F+Dp3NY3AFTL6AJYjOYiVhi6UcsZwjk3OlCo5MypEscsDQQsAIrN3YcXS87aTxTIY4rUZX+czm/9HwVfGOC0oUCq78azlbvpazJRQJnTadOz7R/nL3mzqBuk6/bn1qtZqvjHeAUAyEImFrwGYLTNkDU46QayrocIU87rDHHfZ4MZ8X8/lxPwNlMGlMAZOPImi0beCXtQ/F/A8AoShqPp+vr6//1a9+VV9fP8v3+3w+Ho+HXGFEczAS8mE+H+b3YD5XyO0Ku50hlyPosgWm7MEpW8Dux4MqXoqKn6Llq7UCdXR+EBUvBd5AXDcMwwAATOZi1J7gBH7S2PHZ+PF2Y0euNLs+dV2tbi28uLnM7CfdTioERXRazh8db20xtKXw5BtSazam1SgXq+LU6/UKBLCQ5/p5MK/BMznunpj0mSe9JpPPYvZbAnhQyVPIuTIFVy5li2UcqYglELNEIpZQwOJHV1C50rk9tt9IYvURRheRuUoDSJjAor/BCY9p0mdsN3YYPJOusFsv0KaLUrMk6Vni9CxJhgyeWxMJSZEd5q6m0aOtE6ezJOkb02p/tOr7cDouaE5oCK1CVVKhKvnx6ru7rL1Hxo7ddfA/9EJNQ1r9hrQaCVsU7wChL2EENuIaH3SNDDlHR93jI+7xcCSsF2pThVqdQFOtW6Phq1S8lMS5CE6sRHhNLBoz2rUONF9uDEXCY27DiHt8yDV6+sKHA85hBkrPlWbnyrLzZSvyZDnwnBsv/Y7BwyOfN48dU3IVm9Lrv196O7xGgeYJRdAyZVGZsuhHq+7uMHc2jR59ufuNAnne5oz1Nbq1sLQ4LiIkMeQauTg1cHFqoH9qcMJr1At12ZL0THH6Ou3qdJE+wadoX2KJ8Cux6azoEKXpLRa/9ZJj6KJj8J0LH12cGpByJAXy3CJFfpFiZapIB5f2Xmj2wNSh0c8/HT6CEfjmjA1/vOFxnUBz7Y9B0FzQUVqlpqJSUxGKhE9MtB8a+fyp0y/U6Cq3ZG4sURbCw3yheTFfj+3CeduF87YLlxxDar4yX7YiX7ZiZ85NmeI0Bo0R7wDnYDkkwn+l5KUoeSnRqZ5Iihpzj5+3X+yy9r3R+17Bi8T8AAAUCUlEQVQoEipJKSxJKSxXFaeL9PGOdFnBCOz4RPvB4eaL9oH61HV7Ku8rVOTB8xG00Nh0VkN6XUN6nTPkOjx69I8dL/kw/5bMhi2ZG2FZTWz5MP85S0+ntafT0mvymfNkOcUpK28v+uZSn0s2sYplFoEtYD9nOd9p6Tlr6cYIrEJVulpdtlpdmrTj1WJSLDPoHD4w1NQ82pIjzdya2VCrr4Izv1w3WCwzfwPO4YNDzc2jLVmS9JuyNtXpq+ZZcJ/MxTIRkjhvv3DadO6MqdPgmSyQ55WriktSCnNlWTQkbvP6xvYbSbpEOJPJZzlj7jxj6jxr7lbyFJWa8kpNRaEiP6lqUOeTCP14oGn06IHBw+6w56asTVsyG5Q8RawDTDowEcYKTkZOTLQfGDp8cWpgU3pdY9bmbEnG9f2oJEyEFr+t3dhxynT2rLk7VahbrS5bpS5dKc9loAnRjggTYeyRFNlr7283drQbO8x+62p1WbWuslJdzmfy4h3agru+RNhj6/t48NDxifbV6rJtWZvLVSVwkb9YgYkw5qwB+ydDhz8ZapawRduyNzek1XEZnDn9hCRJhCRFXZwaODHZ3jZx2h50VGoq1moqVqlLE7DeECbChRWd6LJ18lSn5XyeLKdGt7Zat2YZ3+jMKRG6w55Ph498PHQYALAt64YbMzcm4BGy1MFEuEBIijpjOvfx0KGz5u46fdW27M0r5bmz/OzyToQYgXWYu49PnDwxcUrCFq3TranWrsmTrUjkq9vlPI4wEcg50m3Zm7dlbw4TWHRKiz/3/FXJU9Tq19bpq9KSsr6GAtRZc/fHg4dOmc7W6CofrLx35nKsELQkoAiyRlO+RlPuDLn+Pnzk121P0VH69uzNN2SsFzKXbZK7Cj8eOGnsOGZoO2U8myPNqtFVfrvg5uQsL4J3hNdGUmS3tbfFcPKYoY3D4NTpq+r0VSukWfGOKzaufkdoD0wdHG7+ZKiJy+BGTxlLujZsSYB3hIuDAlSXtffA4KHWydNrNasas24oUxVdqch5Od0ResLeExPtLYa2LmtvccrKWn1Vja5yyTXtwKbRuKEAdXFqoGW87fPxEwCA+tR1dfqqfPmKJT1C4CsTIU7gJyZPfTLUdMF+aUNaTWP2DbnS7Cv8ACjGYCJcZF7Md3jk6IGhQ348sDWz4StrvpZBInSGXMcMJ48aWi/YL61Wl9Xpq9ZqVy3d61qYCBPCoHO4xdB2dLzVjwdr9WvrU9cVK1YuxXLTyxJh/9Tg30eONI+2ZIrTbsraVJ9aDQdCLDKYCOOl3zF4cKi5eawlR5K5JbOhTl81PU/N0k2E1oC9ZbytxdA67Bqr1FTUp65boy5fBvPvwESYWMY9Ey3jbUcNrVa/vUZfWaNbW6EqjvkqIQsnmginMGfzaMuh0c8jZOTGjA03Zm5ULdZ0xtBlYCKML5zAj0+0fzpy5Lzt4jrdmhvS6ytUJX6ff2klwjG34dhE+zFDm8lnqdatqdOvW6UqWVqzvVwdTIQJyuy3HjOcPG44OeAcXqUurdKurtKuErMSei5ga8B+ZOTYUUOryW+uT63enLG+QJ4X76CSHUyECcIZcjePtjSPHTX5rNXqNTdk1RenFCRyqw9BET3WvrbJM8cn2nESr9ZV1urWlioLEznm6wYTYaLzhL2tk6fbJk+fMXemCnVrNasqNeUrpNmJU4s85jacmDx1zNA26TVXaVbV69ZV6iviOEkENBNMhInG6DMf7G9qs5yxBx01usoafWW5MoFafZwhV7vx7Cljx2lTp5qvrNKurtFX5kgy4x3XwoKJcMnAyUi3tfeU8Wy7sWMq5CxXFleoSkqVhalC3eIH48cDXdbzp4zn2o0dBEVU69bU6NaWKgsJnACLtR4hNBswESag6GnX6DMfN5w8MXFqwDlcklJYqSmvUJXohdrFj8eH+bttfefM3WfMnbbAVIWqpFJTsUZTLudIFz+YuICJcEmyBx0dps6zlu5Oy3mMxAvleYWK/HxZTo40i0NnL9BOHUFnr/3iedvFblvvqNuQL1uxWl1WqanIFKdNv2cxF+aFZgMmwgR02WnXg3k7TF3tprMdpk4SUOXKogJ5XqEiP0OcukAtKyRFTXiNF6cG+uz9PbYLJp85X76iLKWoXFWSJ8telo2fVwcH1C9Jco70xsyNN2ZuBABY/Lbztgu99v7Px48Pu8aVPEWmOC1DlJYm0mkFai1fPdf5n6J8mN/oM497JkbdhiHn6CXnEBbBChS5K+W5d5f9e4Esdzl1lUNQHAmZgg1pNRvSagAAk15Tl7W3x9b3waUDFr8tU5yWI8nMEKemCnU6gUbBlV1HliIp0uy3Gr1mg3dyxDU+4h4fdA6LWMI8WU505cUV0mw6CvsyYgYmwjhQ8hRKnqIhvQ4AQFDEqNsw6hofdo81j7ZMek2TPhMTZcq4UgVHJmQJhCwBm8YSsPgAAAQgPAbXh/sBAAE86Mf9PszvDLmngg5bYIqkSDVflSrUpov0WzI33i+9KzkniYCgxaQVqLUC9U1ZmwAAATw46BwedI6OuMaPjrdOeE2ukFvBlck4EjFbLGYJ+Uwen8FDUZSO0imKIigCABCKhDECc4c9nrDXFXZb/HZXyCXlSLUClV6gzRCnrk+tzpFmCpiwkWChwEQYZzSEliVOzxKnN8zY6Al7bcEpe2DKHfZ6wt4wEfaGfQAAClBGn5nP4AEAOAy2hK0VMHkStljGkco5UiFrKZV3Q9Dyw2VwilMKilMKprfgBG4N2B0hpzPkdoc9Pszvw3wkReEkjiJotB2VRWeKWEK9UCtkCsRsUQpXLuNI4Q3fYoKJMBFFbwSzxOnxDgSCoHlh0BjRW8Z4BwJdTdJ1sUIQBEHQTDARQhAEQUkNJkIIgiAoqcFECEEQBCU1mAghCIKgpEYHAJw6deq1114DANx+++2VlZXRF959992DBw8qlcr77rtPp4vDlGAQBEEQtAjQrq6uTZs25ebm5ufnb968+ezZswCA559//qGHHmpoaAgGg9XV1YFAIN5xQhAEQdCCoIVCoZqamkceeaSystLlcjU1Ne3cufO222576qmnbr755q1bt77zzjssFqusrOwrP//aa69t2LAhPT19lvvDMIzJZCbhXKMJiyAIAACNBkfvJhAcx+HsrwkFwzAWa8kvZrucxPYbQdva2jZs2BB9sn79+tbWVovFMjQ0NL1xw4YNra2tsdofBEEQBCUUutlslsvl0ScKhcJkMpnNZhaLNT2xt0Kh6O7uvtLn/X7/ww8/PP0TorZt2/atb33rq98/MexveR/eESYOkiQBACgKy6YSBUVRBEEE6HDWpwQSiUSC8BtJDNzt30OF0mAwOMt2LDabfc3zG53FYkUX4gEAhMNhDofDZrNxHCdJMvrhUCjE4VxxMQQmk3nDDTcUFBTM3JiTk8Nmf/XSQrhETl/dABNh4ohEIgAAOjzIEwZFUeFw+EpHEBQXwWDwKqdBaDGxRBKExcZxfJbHyGyu8uk6nc5gMESfGAwGrVar0WgoijIajdFi0ejGK32ewWCsX79+9usR0nhCbmktTISJA65HmGgoiqL8fi5cjzCREF4vN3ar30Hzh6JoDNux0F27dr355psURVEU9dZbb+3evVsoFDY0NLzxxhsAAI/Hs3///t27d8dqfxAEQRCUUOj33HPPu+++W1NTg6Koy+V69dVXAQCPP/54Y2PjiRMn+vv76+rq6urq4h0nBEEQBC0IulQq7ejoaGtroyiqqqoq2kS2atWq/v7+9vb2lJSUKw2cgCAIgqBlAAUAMBiMurq6+vr6mR1FYrH4xhtvjHkWNBgMPp8vtj8Tmg+bzWaz2eIdBfQln8833W0PJQKCIAYGBuIdBfRPhoeHp8s852+xi+bvv//+48ePL/JOoat47rnnXnjhhXhHAX2ppaXlxz/+cbyjgL7kcrk2bdoU7yigf7J79+6xsbFY/TQ4egyCIAhKajARQhAEQUltvsOoI5FIV1fX7N/vdDp7enq4XO489wvFytjYGIqiR48ejXcg0BfOnz/vcDjgN5I4PB5PJBKB30hCCYVCp06dMhqN13znqlWreDze1d+DUBQ1n2jefvvtZ555Zvbzkvh8PjabDecxSRyhUAgAAOcxSRyRSCQUCvHhgPqEQVGUx+MRiUTxDgT6ksfj4fP5sxlT//LLL2dnZ1/9PfNNhBAEQRC0pME+QgiCICipwUQIQRAEJTWYCCEIgqCkBhMhBEEQlNQWtnoTx/EjR454PJ4NGzZctngvAODkyZOjo6Pl5eUrVqxY0DCgaR6Pp62tzev1FhUV5ebmznypo6Nj+rFCoUhNTV306JJRb29vtHAXACAQCC47FoaGhk6fPq3X66urq+MRXTKa+Y0AAEQi0XTN4eTkpNlsnn6ptLR0lmvDQtfHYrFMTExkZWWJxeLoFq/X29zcjKJoQ0PDZYMiCIL47LPPHA5HbW2tWq2e044WsGo0HA5v3LiRoqi0tLSmpqbm5ubi4uLpVx944IF9+/bV1dV98sknv/3tb7/zne8sUBjQtIsXL1ZWVq5evVqhUBw6dOiee+559NFHp1+l0Wi1tbUMBgMA0NjYCGf5WhwrV67k8XjR47ykpOR3v/vd9Evvvvvuvffe29jY2NraunHjxueffz5+YSaRO+64Y2JiIvr47Nmzu3fvfvHFF6NPH3744TfffDMnJyf69MMPP4SjXBZOYWHhyMgIjuN79+7dsWMHAMBkMlVVVRUVFUUikcHBwdbWVoVCEX0zSZJbt2612+0rV648ePDgvn371q1bN4edUQvm9ddfLykpwXGcoqif/vSn3/jGN6ZfGh4e5nK5k5OTFEV9+umnGo0m+jZoQTmdTqPRGH187tw5BEFsNtv0qyiKznwKLY78/Pzjx4//63aCIDIyMj788EOKoqxWq0Ag6OvrW/ToklogEBCLxceOHZve8l//9V979uyJY0hJ5cKFCziO5+bmfvTRR9EtDz300K233hp9vGPHjp///OfTb/74448zMzMDgQBFUU8++WT0Hmz2FrCPcN++fbt27YqOnb/55pv3799P/ePu88CBA1VVVRqNBgCwadOmUCh05syZhYsEihKLxdMtBnq9nqKomU1AAICOjo7W1la32x2P6JLXhQsXWlpaLlsDpKenx2KxNDY2AgAUCsX69ev3798fpwCT1N69e1NSUi5rlLbb7c3NzZcuXYpXVMkjLy/vsqlX9u3bd8stt0QfR3PK9Ev79+/fvn07h8OJvvTZZ595vd7Z72sBE+Hk5KROp4s+1uv14XDYbrf/60soimo0msnJyYWLBPpXjz/++KZNm6a/BQCAVCr9zW9+88ADD6Snp7///vtxjC2pcLnc119//Wc/+1lGRsYf/vCH6e1Go1GlUkVbqgEAOp1uNrNJQTH0yiuvfPe730UQZHoLgiBdXV1PPPFEdXX1tm3bwuFwHMNLQpfllJlZY+ZLWq0WQZA5HS8LWCwTiUSme5KjD3Acjz7FcXzm1Dh0On36JWgRvPLKK3v37r1sPSyj0Rg97b711lt33nnnli1brjlBHzR/bW1t0V/7yZMn169f39jYGK2XwXF8ZiEGnU6P4epr0DUNDw+3tra+/fbbMzf+/Oc/f+yxxwAAHo9n3bp1f/zjH/fs2ROnAJPRZTll5hEx8yUEQVAUnVNOWcA7QpVKNd3aY7VaaTRaSkpK9KlarZ7ZEGS1Wuda5ANdtzfffPORRx45fPiwXq+fuX365uOb3/xmKBSCK5Eujulf+9q1a7Oysjo7O6NPo8fIdG+C1WqNdiVAi+Oll17aunXrZeel6S9LKBTu3LlzZqE1tAhmJo7LjoiZ6cbhcEQikTnllAVMhHV1dU1NTdHHTU1N1dXVdDo9HA6HQqG6urrjx49HO6h6e3u9Xm95efnCRQJNe++99x588MFPP/10ukwfw7BgMDjzPdE+aq1WG48Ak5fD4TAYDFqtNhKJ+P3+wsJCFEWjfec4jh89erSuri7eMSYLgiDeeOONO++8c3qL1+slSXLme7q6umb2LECLoK6u7vDhw9HHTU1N9fX1AIBgMIjjeDTdRC8cm5qaCgsLZTLZ7H/yAg6fcDgcRUVFu3btysrKeuyxx954442tW7f+8Ic/9Pv9f/nLX2688UYEQXbs2PHss89+7Wtf++///u8FCgOa1tvbW1paumXLloKCguiWu++++/XXXz927Ni99977wQcflJeXu93uF1988etf//rTTz8d32iTQVdX1y9/+cuqqiqKov785z+npaUdOHDgr3/9609+8pPR0dHoUXP//fcfPHjQ7Xa3tLTEO95kceDAge9+97sGgyF6C4hhGIvF6ujo+MUvflFWViaXy48cOXLixImOjo7LmlWgGHrhhRdGRkZefPHF9evX5+Tk3HPPPW63u6amZs+ePZFI5Omnn25vb8/Nza2trd2+fft9991XXFxcW1tbVlb2+OOPP/nkk9/61rdmv6+FXX3CaDS++uqrLpdr165d0VEdR48exXE8Win60ksvDQ8PV1ZW3nLLLTN7pKEFYjKZLqs83L59+8TEhMlkqqys3Lt378jICJ/Pr62t3bRpU7yCTCper3fv3r0XLlxgMBjl5eW7d+9GUXRgYKCtrS06svb9998/ceJEWlraXXfdBVfxXDRtbW2BQKChoSH6lCTJZ5999tZbbz1z5kxbW5vP58vOzr7tttskEkl841ze9u/fbzKZpp/u3LkzJSWlt7f3rbfeQlH029/+dnRKkA8++CAjI6OsrMxms73yyitWq7WxsXHjxo1z2hdchgmCIAhKanCuUQiCICipwUQIQRAEJTWYCCEIgqCk9v8BM/NdRPUoXl0AAAAASUVORK5CYII=",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 13
}
],
"cell_type": "code",
"source": [
"ρ = real(scfres.ρ.real)[:, 1, 1] # converged density\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": 13
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.3"
},
"kernelspec": {
"name": "julia-1.5",
"display_name": "Julia 1.5.3",
"language": "julia"
}
},
"nbformat": 4
}