{ "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.143566993365 NaN 3.81e-01 0.80 8.0\n", " 2 -0.156033237910 -1.25e-02 7.92e-02 0.80 1.0\n", " 3 -0.156768052529 -7.35e-04 2.74e-02 0.80 2.0\n", " 4 -0.157046025712 -2.78e-04 4.79e-03 0.80 2.0\n", " 5 -0.157056320778 -1.03e-05 3.58e-04 0.80 2.0\n", " 6 -0.157056353872 -3.31e-08 2.70e-04 0.80 1.0\n", " 7 -0.157056406077 -5.22e-08 4.34e-05 0.80 1.0\n", " 8 -0.157056406911 -8.34e-10 4.17e-06 0.80 2.0\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380298 \n AtomicLocal -0.3163472\n PowerNonlinearity 0.1212609 \n\n total -0.157056406911" }, "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.05568574649946116, 0.0, 0.0]\n [0.055685969057733764, 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+gvaeTAAAgAElEQVR4nOzdd3xTVd8A8HNHdtIk3XvvQWnLhgJFVrEMmYKAovLAg6ggKigOFNzyyAOKIiCoCCqCKJQNQkFmBy0tdO/dJs2ed7x/hLfyMEpHmqTJ+X74Iw0nNycnN+d3z7lnIDRNAwiCIAhyVKi1MwBBEARB1gQDIQRBEOTQYCCEIAiCHBoMhBAEQZBDg4EQgiAIcmgwEEIQBEEODQZCCIIgyKHBQAhBEAQ5NBgIIQiCIIcGAyEEQRDk0GwuEObl5RmNxk4mpigKLhFnRSRJWjsLDg2Wv3XB8rcuM5a/zQXCSZMmNTc3dzKxXq+nKKpX8wN1QKPRWDsLDg2Wv3XB8rcuM5a/zQVCCIIgCLIkGAghCIIghwYDIQRBEOTQYCCEIAiCHBoMhBAEQZBDg4EQgiAIcmgwEEIQBEEODQZCCIIgyKHh1s4A1LuMFCiS08VyulIJjBQgacBnAA8OCHFC4p0RBrwQgiDrIWlQKKOL5HS9BrTpARMFCAL8eSBciESKEC6sni0FlrR9atSCn8uok3XU3420FxeJECEhAsDEAIaAJi241ARuy6gyBd3fBZkagM4MQoIEiLWzDEGOokUHfq+kDlRQV5ppLy4SI0a8uUDMAmoCkBTIbAHFcqpMSSe4ION80CeDkTAh/Hn2LhgI7c2pOnpzAfl3E/1EAPp8BPrTaFTMenBKNQEuNdG/VVBD/iT7OSMrY7FUPwT+4CDoHp999tn169fvf54kSQzDunQouQFUq2iJHriygTsbPMZCcBQAAJoAaLorWTgAITSQ6cEBPf2VFnBx4MsDHhxH/3mOHj162bJlvXFkxNYWrfb19b169aqPj09nEmu1WiaT2dVz0V4dq6HXZZMqI3itHzorGOV1+iLHSIFfy6n/5FNGCnw6CJvo29mfm1KpFAgE3cwu1GOw/C0jOTl5/PjxkZGR1s6IQ8vOzs7Pzz98+HD7M2Y8/2GL0B4UyuiVV8hKFfhgADotAEW7eN3IQMFToehToegfVdTKK+QXPPDVcCzUycGvPiHoHykpKSNGjLB2Lhwah8PJz8/vpYPDwRJ9m5EC7+dQI48QE3zRvOn49MAuR8G7TQ1Ab07HU/3QoX8SG29SpG11FkAQBPUKGAj7sPw2etAfxNVmKucJfEUsapYhoDgKVsSiV6fi6dXUmHSiTg2DIQRBds48gXDPnj1hYWHu7u7PPfecVqt9WLKjR48OGDDgyJEjZnlTB7ejiBqTTrwUg6ZPwH14Zu7GDBYgpyfh43zQAYeIE7UwFkIQZM/MEAhv3779wgsv7N69u6ioqKqqav369Q9MJpfLV69eLZVKJRJJz9/UkWkIMP8cuaWAykjDF4X3VpseRcBbCejPY/BnM8iNN+HuxxAE2S0zVKPffffdtGnThg8fLhaL33rrrZ07dz5wJOqqVatWrFjh5ubW83d0ZNUqOvkIgSHgyhQ8UtTr41lGeSFXpmJ7SqlFGaQBRkMIguyRGQJhYWFhfHy86XF8fHxzc3NbW9s9ac6ePVtWVvbss8/2/O0c2aUmesifxPxQ9PtRGMdSA379eMjFybhMD9JOEEqjhd4UgiDIYsxQm0ql0vbJHE5OTgCA1tZWZ2fn9gRqtfqll146dOgQ0on5oGq12tfXt/3PH3/8cdq0aQ9L7FDzCA/VoCsz8W+HEOO89CqVpd999xDwSiY+6jB1YKTRjX2nxa9WqzvznUK9BJa/ZVCUnXSGNDc3JyYm1tbWmv7cvn17Q0PDO++887D0x48fP3jw4LfffmupDD4CSZKqu+q+Tp7/XC4XRR/R5DNDIHR2dlYoFKbHpgcuLi53J1i3bt3AgQPlcnlWVpZara6srKyqqgoICHjg0Xg8Xn5+ficn1GMY5iCB8It86oub1OnHsXhnprXysDMFrMsmJ51DT6dipuE5NE3z+Xxr5QeC5W8Zj6xG+wqKolpbW02PtVrt+++/n5WV1UH6CRMmvPbaawUFBTExMRbJ4CNgGHb3CW/G898MgTAsLKygoMD0OD8/38XF5e7mIACAoqibN28uWbIEAFBZWbl7924Mw956662ev7UjoAF48zp5uJq+NAXzNffo0K5al4jxGdSodPLMJCyAD9siEGQFVVVVubm5ISEh+/fvHzly5JgxYwoLC48cOWI0GqdMmWIKWmq1Oj09vaCggMPhTJ48+f5IduDAgaSkJHd3dwDA2bNn3d3dY2Njv//+++nTp6vV6oyMjNmzZyMI8tRTT3399ddffvmlFT6nBZnhSmfRokUHDx68ceOGRqP5+OOPFy1aZGqurl+//tixYwCAjRs3Zv6/mJiYdevWwSjYSSQNFl8gzzfQGWm41aOgyatx6IpYdHQ6WamE0yogyAry8vKWLVu2fPlysViMomh6enpqaiqCIFwuNzU1NSMjAwCQnZ2dkZHh7+9P0/Rjjz2WmZl5z0HS09NTUlJMj3fu3Hn27FkAwMqVKyUSSUlJybp160z/lZKSkp6ebrnPZiVmaBHGxcV99tlnaWlpKpUqNTX13XffNT1fVFTk5+d3T2IfHx+4OmInGSgw/y9SbqBPTcI7v3CoBSyPRlEAxhwl00cjUfDLhBxPdiu9+jppmfca4o6sT7r37o9CoTh06JBQKAQAhIeH79y5c8yYMQAADw+Pjz76aOTIkcnJycnJyabELBbru+++GzBgwN1HuHnz5sKFCx/57uHh4ZWVlQqFwjT+w16Zp35dvHjx4sWL73lyz54996c8ePCgWd7R7mkJMPMMwUSRP8fjLNu7B7osGiVp8PhfzIzJtI00VSHIYgIFyOp+FvpZurIf8GRERIQpCioUipKSkvXr13/00UcAAKVS2dzcDABoaGhYunTp7du3GQyGVqu9v2tUqVRyudxHvrvpJpxSqYSBELI0NQGmnCS8uMjukRhuq/fpX4xBlVpy3DH0/OO4O8fauYEgC3JmgbE+1rz+Y7HubK5mGi24bdu29inapqE9b731VmRkpGms/ldffXX06NF7juDm5tY+z43NZut0uvb/0ul0bPad8CuRSFAUvWf8o/2x1VrWgakJMPkEESRAfhhlu1HQ5MUIYk4wMu4YIdVbOysQ5JDYbHZKSsrevXvF/89gMAAAmpqagoODEQQxGAz79u27/4UDBw7My8szPY6Pjz958iRBEKY/Dx8+3L9/f9Pj3Nzc+Pj49rhor2y7onU8CiMYd5SIFCHbk7Ge7CNhMesSsce8kSknCQ1h7axAkEPavn37sWPHEhMTn3jiidjY2I8//hgAsHTp0jfffPOJJ54YMGDAA2ejzZo1yzSYEQCwZMkSoVAYHh6uUqnGjh2blZXVvlLm8ePHZ86cabHPYi2wa9SGyA1g4nEi3hn5ajjWF4LgHRuHYIsvkFNPEekTcCa8soKgXjZx4sTRo0e3/xkYGHjlypXKysqmpqagoCAPDw8AQFpa2u3bt8vKysLDwwUCgamZ6O7uXlNTY3pVSkqKXq/Py8vr168fi8U6cOBAc3NzWFjYnj17hgwZYkqj0WgOHz586dIlS39Ci4P1lq2QGcC4Y8Rgd+TrEX0pCgIAEAC+GY5xMOS5DLiDIQT1OgaDcc/YewRBgoKChgwZYoqCJu7u7kOHDnVxcWEymaYxLyiKtt9KRBBky5YtV69evTs9hmGenp7tz1y/fv3NN980zTW0bzAQ2oQ2PRh3lEj2RDYN6WNR0ARHwS9jsEoV/dpVC40phyCoh4YPH37PaP+XX37ZNBjVZNSoUQ6yQDQMhNYn1YOxx4jRXsjGwbY3T6LTODj4cxx+vJaGezZBUB/17rvvisVia+fCCmAgtDKJHow9Soz1Rj7ry1HQRMwCxydimwuovWUwFkL2w9YGgh07dmzcuHG9dPAjR44sW7asgwSFhYXjx49/4F57fRcMhNbUpgcTjxFjfZBPBvX5KGjiy0OOTcReuUKerLOr3wnksEoVdLXKtk7m2NjYFStW9MaRaZpevXr18uXLO0gTGRnJZrMPHz7cGxmwFhgIraZZC0anExN8kU/tJQqaRIuQXx/DF5wjcqW2VX1AUFc1aMD4Y6Qnx7Zu3BMEodFoAAAGg+GTTz6pq6tbt27dmjVrqqqqNBrNF1988dprr2VnZ5sSt7a2bt++feXKlevXry8uLm4/SHV19fvvv7927drS0tJNmzbJ5XIAwPnz5zkcTnR0NAAgIyPDtADp7t27q6urlUrlxo0bTa9dsGDB119/beFP3atgILSORi1ISSeeCEQ2DLCrKGgy0hP5ahiWdoK0tUtpCOo8pRE8foJ4LgIVsaydlf+Vn5+/efNmAIBer1+zZs1zzz3n6+srl8vHjRu3aNEiAICzs3NKSkpLSwsA4PLly7W1tYMGDWIymSNGjCgtLQUAtLS0DB48WKfTRUZGLl269I033pBKpQCAo0ePtk/MOHHihGk9ms2bN5eXl8tksrffftv0X6NHj/7rr7+0Wq0VPnzvgPMIraBOTT92lFwYhr7Z324vRGYGoXVqkHqcvDgZF9tYPQJBj2SkwKwzxCA3ZG1/9Ph9/2uoKpT99pVlcsIMjhU9saSDBBs2bBgwYMCzzz7r6uo6duxY00DQkydPnjt3btasWZMnT548eTJN0zKZrLa2dt++fW+//faOHTtGjRr14YcfAgASExNjY2NNh8rPz58yZcojs+Tm5sbhcEpLS+Pi4szxEa0PBkJLq1TSY4+RS6PQV+PsNgqavByL1qjpaaeIE6k42w7bvZDdogFYfIFkoshXwx984uIe/qLZL1kmMyiH13ECUxhDUdTV1bU9pLm7u0skEgBAQUHB888/L5VKBQJBY2NjWloaAKC4uDghIcGUMjo6un3lUo1Gw+F0auFgLperVCq79YFsEQyEFlUip8cdI1/rh74QbedR0OTTQdhT58gF58hfxvSNFeMgCADwdiZZKKfPTsIfNqsXZXOZfmGWzdRDYdg/0dq04jYAAEEQ08DOl19++ZlnnjHti/7GG2+YoqNIJGpfcVutVpvWnQEAuLm5mfpIAQAcDkcmk7UfWavVtsdIiqKkUundk/f7Ooeojm3ETSmdcpR8J9FRoiAAAEXA7pFYq45+BU60h/qIb25Tv1bQh8fjXLtoJrRHLJlM9uuvv5qefPzxx/fu3VtfXw8AaB8CAwAYPHhwbm6u6XFiYuLp06fbm32//fZbUlKS6XFxcbFAIAgODrbYp+htjlIjW93VZnrcMeI/g9Fnwx2rzFkY+H0cfrae/hxOtIds3p9V1Poc6vhEzM1etltYtWrV888/n5qampycnJiYaHpy7NixS5YsSUhICAwMpCiKw+GY9iacPXt2+zYUkyZNmjJlSkRERFFR0fz58w8dOrRlyxbTy9PT02fNmoUgdtTJQ9sYHx+f2traTibWaDQEQfRqfsziZC3lvsdwtJqydkbMTKFQdDJlrYoK2Gf8qZTs1fw4ms6XP9QZfzdS7nsMmS33/k5HjBhx4cIFq2TpgQiC0Gq1NE2buijbn5fJZEaj0fRYpVLpdDrTY4lEkpubq9VqtVqtWq1uT09RFEVRBQUFQqGQou586lmzZh04cKA9jUajiYqK+v3339ufIUkyLi6usLCw1z7fgx0+fDgtLe3uZ8x4/ttF49+2/V5JLfubPPAYPsLTji6gusiHhxydiI1JJ1xYyARfxy0HyGaVKuhZZ8gfRuFJrrZ+fmIYZroviCDI3Sui3b1MKI/3zxAbZ2dnZ2fnew6yePHi2NhYo9H4zTffrF27tr159/HHH989WZ7D4bDZbJFI1P5MaWnp0qVLIyIizPqZrAwGwt7133xq403q1CQ8Vmzrv67eFi1CDozFp58m0ifgA2y+roEcSq2aHnuU/GQQ6jhXaU899VRWVhaCID/88MOwYcPanw8ODn755ZfvTvnMM8/4+/u3/xkeHh4eHm65jFoEDIS9hQbgvWxyfzl9YTIWwHeUX1fHhnsgO5PxySeIc2l4hBCWCWQT5AaQdoJ8MQadH+pA9+9Hjx5996aGHXjpJQtNFLEiB/riLUlLgDlnyL/q6YuTcRgF75bmj3wwEEs9TtZr4KIzkPVpCPD4CWKcD7LK3uf1Qh2A3735NWrB6HSChYGTqXBRlQd4NhxdEolOOEa26a2dFcixmZaPCXVCPu37e79APQG7Rs0sR0I/cYp8NgJ9OwHOIH+o1fFoq45+/ARxahLOg+cgZA0UDZ45T+IosiO57+2GTVHUt99+m5OT4+fnt3z58rsHs0DdAFuE5vRrOTXxOPH5YPQdGAUf5dPBWJQImX6K0MOp9pA1vHiZrNPQv4zB8D5YCy5ZsuTIkSOpqalVVVUpKSkkCX9FPQKvxs2DoMDaTHJ/BX0qFe/nDIPgoyEAfJuMPXmWnPcX2UcrI6jveuM6eb2FPj2pm6vg3mot3nT9G3Nn6sH6uUcvT3r+7mdqa2v3799fU1MjEAimTp0aHR197Ngx0yKiUPfAQGgGTVow9yzBwsD1abgLvCnYaRgCfkrBpp0ins0gd4+Ci5FCFvLhDepINX0uDXdidPMIAULfVYM72sbdjPiMexfdLioqCgkJEQgEAAAEQRITE4uKimAg7AkYCHvqTD399HlycQT6dgIKq/KuYqLgwGP4pBPE0r/JbSP63q0aqM/5Ip/6voQ693iPrll5DG6Ec6j5MtU1Go1GoVC0/6lQKJycnKyVGfsAO6S6z0iBtZnk0+fJH0Zh7ybCKNhNHBwcHo/faqNfvgzvc0C9a+st6ssC6swkzItr7az0TGlp6blz5wAAFRUV58+fT0lJsXaO+jYYCLupUEYP+5PIldA5T+BjvGEM7BE+A6RPwK82069cgbEQ6i3fFlKf5lFnJmG+vD7/g+3fv//69euHDRs2cODADz74IDTUas1T+wC7RruMpMHmAuqjG+SGAdi/IuGVhHkImeBEKj7uKLHqKrkRTuqCzG17IfXBDersJCxQ0OejIACAx+OdOXOmsrLSzc3t7mVFoe6B9XjXFLTRIw4Th6uoy1NwGAXNS8QEpybhGQ30qqskXHUGMqPthdSGG9TZSViIkz1EwXaBgYEwCpoFrMo7S0OAN6+TKenEonD0zOO4nf2ibISICU6m4hcb6ZcuwVgImceXt6gPc+0qCkZFRS1dutTaubArMBB2ym8VVMwBoloN8mYw/hUJh8X0IjELnJqE50joJRdJCgZDqGc+v0ltyqfOPW4/URAAEBoaOn/+/G6/nKZpnU5nlpyYtjw0y6GsCwbCR7jeQo86Qnxwg9o1EtszGvPkWDtDDsCJAY5PxEvl9IJzpBFuaw9117tZ5M4i6vzjdrj9S35+fm5ubsdprl+/XlxcfP/zN27ciIyMNEs2oqOjc3Jyuv1ygiD2799PUdb/kcNA+FC3ZPSM0+T00+TCMDRrGj7ay95+S7aMzwBHJ+JKI5hxmtTBkaRQF9EArLhCHq6mM9Jwn74/RvR+P//8848//thxmq1bt969xa4N0mq1s2fPtoX14eCo0QfIk9If3qDONVCv9cP2jMY4sJCsgY2BA2OxZ86TE48Tf4zDhUxrZwjqI4wUeC6DrFDRZx/HRfZ42ty+ffvs2bMEQaxZsyYkJGTx4sUEQWzbti0zM9PHx2f58uWenp5Xr17Nysqqrq5uaWlJSEiYM2fOAw9lMBi+/vrrnJycgICA5cuXu7m5mZ4/d+7coUOHpFJpYmLiyy+/TJLk7t27r127RhBEcnLy008/jaIPbURt3LgxJSVl//79zc3NM2fOTE1NNT1fXFy8Y8cOiUQyevTo+fPnIwiydetWAMDatWtRFH3hhRf8/PzMXVSdBVuE/+N0HT3pBJF6nBzohpTNYayKQ2EUtCIGCn4cjcU7I6OOEA0aa+cG6gvUBJh6ipAZwMmJ9hkFAQBCodDLy8vDwyMpKcm0Wfz8+fMPHz48a9YskiSTkpJkMpmrq6uLi4uvr29SUlJQUNDDDjVjxoyzZ8/Onj1bqVQOHDhQrVYDAL777ruFCxcmJCTMmzevtbWVoiidTldcXDxlypTp06fv3Llzw4YNHWRv+/bt8+bNi4mJmTBhwvPPP//nn38CAMrKyoYOHeru7j516tRNmza98cYbAABTD21iYmJSUpJ1h7/Cah4AAGQG8GMJta2QQgBYGYf+PhZlwZlstgFFwH+HYh/lUsMPE0cnYJEiO+zmgsylWQumnCJiRMi2Eb27jLuiQlO8t7YX3+AuojBe6Gyfu5/x9vaOiorS6XSzZs0CAJSVlf355591dXVisXjSpEmZmZm7d+9esWJFYGBgbGysKc0D3bx5MyMjo66ujs/nT5o06cqVK3v37l28ePFbb731/fffjxs3DgAwceJEAACfz//000+1Wm1TU9Pq1avXrFnzzjvvdJDnpUuXmobzqNXqzz//fMqUKVu2bJkxY8arr74KAAgNDU1ISHj33XfHjBkDAJgxYwaD0d1VX83EoQOhkQKn6ug9pdSxGmqiH/rlMGyUFwIrWhv0RjzqywMp6cQvj+EjPeFXBD1AsZyedIKcH4q8m9jri9byfdmxSwN7+U3uQJmPCOmlpaWBgYFisdj0Z2Ji4gPHyNyvpKQkPDycz+ff/UKZTNbQ0DBkyJC7U6rV6rlz5xYWFgYEBGi12oaGho6PHB8fb3rQv39/U2ZKSkralwWPjo7GMKyqqsrHx+ehh7AsRwyEagKcrqMOVdGHq6gIETI/FP1yGMMZ7hph2xaEoj5cZPYZ4rPB2IJQ2KUP/Y/zDfSTZ4mPBmLPhFvi3EAZKNvFyh2v7fMWRCKRTCZrf76tra39Vl/HRCJRW1vb3S/09vbm8/lMJlMikZh2tzDZtWsXhmGmkHb9+vXx48d3fOT2/MhkMlOEFovF7U9qNBqdTtceuW2Bo1QoehJcbKQ/uEGNO0Z4/WT88haV4ILcmI7/PRn/dxQKo2CfMMYb+etx/L1sam0mnGII/WNnEfXkWeKnFNwyUdAWuLm5VVdXmx7369cPRdH9+/cDAGpqag4ePGgan3J3mgcaOHCgQqFIT08HAJSXl6enp0+cOBHH8SlTpmzYsME0mLOlpYWmaa1WiyAIAIAkyU2bNj0ye99++y1BEBRFffPNN6bO1dTU1B9++MEUC7ds2ZKYmOjh4cHn8zkcTseZtAz7bBGSNKhT06UKcEtG50vprFb6loyOESMjPZGXY7CUcQjPPj+3/YsSIVem4DNOEzNOkz+MxgRWvrMAWRlBgdeukUdr6Iw0PEzoQH3mc+bMOXDgQGBg4IgRI/bs2bNv375nnnnmvffea25uXrNmzYgRIwAAzz777MKFCwMDA6dPn/6f//yn/bUoirLZbACAQCDYu3fvc889JxAIWlpaNmzYkJiYCAD46quvFi1aZOpu1el0t2/ffvrpp7///vuoqCiKombNmiUSiUyHYrPZDxw+GhYWFhsbazAYAgMDv/nmGwDA3LlzMzMzIyIixGIxi8Xau3cvAABBkPXr148dOxZF0YMHD7Z3qFoeYmvrAvj6+l69erWTfcf1ct1/i3AEQdRGoDAChQE06+gGDajX0G5sJNQJRImQWDGS4IokuCDd24oa6oBSqby7/8RiDBR46RJ5sYk+NA4LtaMVQ7rKWuVvIyR6MOcMwUDBvjG9O0A0OTn5o48+MkUXW9bS0uLs7IxhXa7pmpubXVxc7nmhaddDDw8P5P8HTtTX1zs7O5uCaAciIyN37dqVkJCgUqlcXV3v/i+j0ahQKFxcXLqaQwDAkSNHXv74m2Eb/vDggM8HY8Cs53/fbhlhCC1iABRFfLjAiQmETODORr24wJuLwGGfdoyJgm9GYNsKqRGHiR3JeJq/48ZCh5XdSs86Q84MQj4ciMENnU06eWvwfu7u7vc/yeVyudz/2bbR29u788dks9n3h0wGg9G9KGjiyQXjfZHeuC3StwOhEwO83g/BMEe5MQDdbUkkGu+MzDlLXm1B1iXC2tCB7Cyi3swktw7DZgTB377Nef/994ODg3vjyM4spJcGysHTCOrDhrgjmdPwK830uKNwxr1DUBnBwnPkpnwqIw2HUdA2zZ4928PDw9q56Bp4JkF9mxsbnJiIp3ijSYeMJ2pt64Y3ZF65UnrgHwQTA1en4hGONDQG6m0wEEJ9HoqAtxPQn8fg/7pIrrwCF+m2QzQAX+RT448R7ySgO5Ixbt++pQPZHBgIITsx0hO58QRerwGD/iBypbBpaD9q1PSEY8SBCurqFHxuCKyyIPODZxVkP8Qs8MsY7LV+6PhjxAc3KML625xBPfV9CTXgEDHaCz2fhgcKYHco1CtgFwNkbxaEoileyPMXyEOV1M6RWD9nWHv2SbVqeulFsk4DTqXi1v0Sk5KS5syZc/9kAJqmEbg4saVoNJrhw4f30sH79oR6rVbLZDK7MYEUMgtbntBNA7C7mFpznfxXJLq2P2aXyynYcvn3BEWDbYXUu1nkSzHY6niUYe1+q4ctM61Wq627eZCj8fT0vHtqI5xQD0GPgACwKByd6Iu+dJnsd5DYOgwb6wMv3vuAPCm95CKJo+BcGh5tG7tucTicB06Ms9cLEQcEAyFkz7y4YP9jWHoN/a+L5CA35PPBqC/PJupW6H5yA1iXTe4tozYkYc9HovB7gizG2p0OENT7HvdDCmbgkSKQ8Dvx4Q1KS1g7Q9D/omiws4iK+s2oIcCtmYzFMApClgUDIeQQODhYl4hdm4rfkNDRB4h9ZXAfJ1txtp4ecIj4oYQ6PB7fNgJzgXuiQRYHu0YhBxIkQH59DLvYSK+6Sm68SX08EN44tKYbEvqN62SpAnw0EJ0J10uDrAcGQsjhjPBErkzFf6ugXrhE+nDBhgHYMA8YDi2qUEavy6bON1BvJWD/irT+uFDIwcETEHJECACzgtCCGfj8MPSpc2TqceJSE+wrtYRCGb3wHDkqnejvgpTOYbwQDaMgZH3wHIQcF46CZ8PRoln49EB0/jnysaPEqToYDgtLjy0AACAASURBVHtLdiv95FlydDoRKUJKZjPWxKM82CEF2QZ4JkKOjomCxZHoonB0bxn1yhWSiYFX49CZQbClYh40ACdr6Y03yUIZWBmH7khm8BnWzhME/S8YCCEIAABwFCwMQxeEoenV9Bf55OvXqBei0eciULd719WCOktNgD2l1OZ8iomBlbHo3BB4bQHZKBgIIegfCABp/kiaP54rpbcUUBH7jWl+6JIodDgcTdMVBW30t4XUT6XUKC/0y+FYihcsPcimwUAIQQ8Q74zsSMY+G4TtKqYWXyARAJ6NQJ8KRT051s6ZDVMYwa/l1K5iqkoFng1HcqbjfnAdH6gvgIEQgh5KzAKvxKGvxKEXGuldxVT0b8Zh7si8UHRqABzo8Q8jBU7W0XtLqaM11BhvdE08OskPxWAEhPoO+GuGoEdL9kSSPbEtBPZ7JfVTKfXC3+QEX3RWEJLqhzrsbukEBc410PsrqN8rqTAhMi8E3TyMAdeFgfoiR/0RQ1DX8XAwPxSdH4pK9OD3SurbQuq5C2SKFzolAJnkh3o4Rq+p3ABO1lF/VtHHaqgwITIrCL0+DQ/gwwYg1IfBQAhBXebCAs9HoM9HoG16cLSG+rOaXnXVGCJAJvoh433Qoe4Ibl/DIyka3JDQp+roE7VUZis9whOZ4o9+Mgj35sL4B9kDcwZCvV7PYnXUM2I0GhkMOIcIsh9iFngqFH0qFBAU9ncTfbKOeuUKWSynh3kgo7zQZE9kgCvC6pt7ApM0yJPSFxrp8w30+QbKnYOM90Ve7YeN9kIctjcYslfmOaNra2vnzp2bl5fHZDK/+OKL+fPn35Pg3Xff3bVrV0tLC4/HW758+bp168zyvhBkI3AUjPJCRnlhHwwAbXqQ0Uida6BXXKYK5XS8MzLIDRnkhiS5IqFCxJbbUDVqOquVzmyhrzTT11toXx4ywhOZEYR8OYzhxX30yyGojzJPIFyxYkVsbOz58+ezsrLGjBmTkpLi4+NzdwI3N7dTp05FREQUFhaOHDkyJiZm1qxZZnlrCLI1YhaYGoBODQAAAJURZLbSV5vp/RX0G5mUVEf3c0HixEicMxIpQqJEiBXnY0j1oEhO32qjb8noPCmdK6FRBCS5IgPdkFVx2GB3xBmOfIEcA0LTPV1cUSaTubm5lZSUBAYGAgDS0tJGjhz5+uuvPyz9zJkzY2Ji3nvvvQf+r6+v79WrV++Jow+j1WqZTCaG9c2+p75PqVQKBAJr56IvkepBrpS+KaVvtdG3ZfRtGa0jQZgQCRYggXzgz0f8+MCbi3hzgRv70X2qnSl/ggLNOrpJC+rUoFZNV6voShWoUNIlcpqgQYQQiREjkSIk3hmJcwbwnl+XwPPfusxY/mZoEVZVVTGZTFMUBABER0eXlZU9LLFUKr148eLSpUsfloCmablczuXe6Yjh8XhMJrPnmYQgW+DMAileyN0rrcgMoEROVyjpShUoUdCn60GjhqrXgGYtzcKAKxtxZgEREwiZCAcHPBxwcMD+/wBpMOBMJml6bKSAygh0JNASoM1Ayw1AZgCtOlppBG5s4MFBfHnAh4v48pDH/UCQAA11QtwdY5grBD2SGQKhTCbj8XjtfwoEgoqKigemNBqNCxcunDBhwtixYx92NI1GM3ToUBS9M+pu27ZtkyZNelhi2CK0LrVajdj0Pa8+AAcgigOiOAC43/tfCiMi0QOZAcgMiMIItCTQkoiGAAbqTpnjtIEF7jzGMDqQDZgozcWBkAGETCBk0M4s4Mx8SJcPCVSq3vpQDgKe/9bVyfLncrntAeVhzBAIXV1d5XJ5+59tbW0eHh73JyNJcsGCBQCA7du3d3A0Ho/X+a5RDMNgILQimqb5fL61c2G3+AB4d5hAqTQKBHBRcKuB5791mbH8zTDdKSAgAMfxgoIC0585OTmRkZH3pKEoatGiRVKp9LfffoNdnRAEQZDtMEMg5PP58+bNW7t2bUNDw759+3JycubNmwcAyMvLS0tLM6VZtmzZmTNnXnzxxYsXL54+ffr27ds9f18IgiAI6jnzTJ/YuHHjypUrhw0b5uXl9ccffzg7OwMAaJomCMKUQKVSRUdHb9682fTnpEmToqKizPLWEARBENQTZpg+YV5w+kQfAoePWxcsf+uC5W9dZix/+1oSEYIgCIK6CC4aaLfadLL8ltv5LYVVipoGVZNE22akCB2hc2IKhGwnd65riDgoTByU5NnfhSO2dmYhyLGoDOrsprwiSUmZrLJe1STXKRQGBY4yWBhTyHLy4nv4OXlHu0bEuUV78u6bWAOZGwyE9qZVIzlVef5CzZVqRW2sW2SMa+SUsIlePA8XrjMTZbBxtlyvkOuVjeqmsrbKv2uvbcna4c51Hek3bGLwGA+em7WzD0H2TKFXnqo8/1fVhXJZVT/36CiXiLTQCb4CbyFLIGQ5GSnCQBradPJGVVOlouZizdWt2buELKeRfkMeCxgZIPSzdvbtFrxHaD+yGnN/L07PbSoY6T90tP/wBI9+OProkqFo6rak+FTF+bNVFyJcQudGTU/07NfJd4T3SKwLlr91dan8S9sq9t06eKU+c6j3wPHBoxPc4xjYo7fioWi6UFJ8vvrSqcpzfgKfqeGpo/2Howi8pQWAWc9/GAjtwaW6az/m/6oxamdFTnkscBQH784kayNpPF2V8fOtg2ycvTh+wQCv/o98CayIrQuWv3V1svxL2sq33/ixTFY5K3LK5NAJPEZ3NvIgKPJi7ZUDRUek2rb5MTPHB6dgiKPXezAQ3gED4a3Woq3Z32kJ3dNxT47wHYL2eMEnGtAZ1Zd35P7oznN7IfG5YFFAB4lhRWxdsPyt65Hl36Jp3ZbzQ3ZT7tNxT04KGcdAzXArKre54PubP7doJEsSnh7hO7jnB+y7YCC8w5EDoVTbtjV7V25z/nPx88cHpfQ8BN6NpMkjpSd35e0dH5SyqN+8hzUxYUVsXbD8rauD8idp8rfCwz8V/DYtfNLc6Ond66TpwLX67K9zdonYwpUDl/o7+Zr34H0FDIR3OGYgpGj6UPHR72/+nBY2fkHMbDbeW7vGyfTyb7J3ZzfdXD3kxSTP+PsTwIrYumD5W9fDyr9MVvnRpU1itmjFwCU+Aq9eeneKpn4vPvrDzV/SQsc9HfckE3O4pSthILzDAQNhtaLu0yubUQR9dfByf6dOlVIPZTXmfnJlyxDvpGWJi9j/e2ELK2LrguVvXfeXP0VTPxX8dqDo8NKERRODx1ggD1Jt25asHaVt5a8PeSnOzbGW64KB8A5zBUJSTxnkRqOaJPUUqScBACgDxZgoS8xgiRkIahM7rdCA/q3w8J78/c/0e3Jq2CTz9oV2TG3UbMncXtBa+M6I18LEwe3Pw4rYumD5W9c95d+kbtlw6T9MjPHGkJdduS6WzMmFmiubMrc9FpD8fPx822ka6mVGncRAGWlTpYoxUZSJMvg4k4/jPDO0XmAgvKMbgdCoIjRNem2LQdeq17YYdBKDvs1IkzRTiDP4OMpEcTYGACANFGWgdFKDUUlwvdnicL44SuAU1J3hXmbRrGn96PImI0m8OWyFN9/TKnk4U3VhS+a3C2Jnz4iYbHoGVsTWBcvfuu4u/ws1VzZe+2pO1BNzop6w5EVqO4VeufHa1ipF7VvDVobedbVqYZoGnfSWsq1IpazW4iyU7crEWCjGulOpknrKqCKMSoIiabaYwXZhsl2ZHDcWx43JcWOxRAzQlZKDgfCODgIhRdAGuVHfZtS1GXQSo67VoG3V65oNAAMcNxbXncVxY7JdWWxXJlvMwLkPDaUUQSurNLJilSRPgWCId7KL+wARgln0RD9X/fem69/Mipw6N3q6dacQNaia3r3wiRff4/UhL/IYXFgRWxcsf+sylT9BkdtufJ9RfWld8utRLuHWzdLJinNbs3fOi5k5K3IK0qWo0kM0kOQr6i9IdK0G13gnUTjfKZiHsR5aWZF6Si81aCUGXatB22LQtui1zXpCR3FcmRw3JtvlTs3MEjOYIgbGfPBxYCC8QylRN56RIQhCUzSpo0gDRWhJQkMalQRlpJhODFPf5p3rDlcmx5XV/SY5DWQlqrpzrfo2Y+hsH8u0DnWEfkvW9htN+W8PXxXpEmaBd3wkI2nckrUjuylv/cg3XFExrIitCAZC61IqlUaceOfCJ3wm981hK52YNvFdNKia1v+9kcfgvjlshZgtssA7apv1pfvrSQPlO8bVJc6p2/eSSD2lbdHrWg3a1jt9dXqZ0SAzAgAYApzBw3AuhrEwhgAPme4FYCBsp5KpFQVaFEUBguAcFGWgOBdjcDGGAMc5vTWCRpKnKD/U4NLPKWiKZ6/ePqyQVa27+GmES9jKgUvNPvy6h0wXnktin06NGGvtvDguGAit61pV9qfZX04Jm7AgdrZFm1+PQtLk7ryfj5affnPoigeO9zaj+gxJzekW//FuXsNdeqkMSD1lUBCEhiS0JKknaQq4JQgBDITtrDVqlNCRRT/UAASJXOjXQfO/Jw6XnthxY8+yxEUTLDL2rBtK2yrWnvtgTGDy4v4L4JpPVgEDoRWZfqFvDlsx2DvJ2nl5sJymmx9c+mJCcMqz/eb1xjI0NEWX/96gqNBEPxfAEj96uTizg4HwDitOn6Apuuxgg7JKE/fvoA5uMXaD2qjZePWrKkXtuyNet8wEiW6rlzR8fmMrjuLvDH+Vz+RZOzsOBwZCqyAocnPmt7nN+WsHrAz3DLV2djoi08s/uPSFjtC/M3yVG9fVjEemSfr27mqapCOf9u+lxsAjwf0IrQ9BkdCZ3qJw/q3vqinCbBcTxdKyxcdW8pn8ryd8ZuNREAAgYPI/G7POz8n73yderVbUWTs7ENTrZHr5K2fekmilX0/43IvnYe3sPIKIJfw05d0h3kn/Or7qcl2m2Y5Lg5Jf6wAA0c8HWCsKmpc9fAYrCkrzZDszCr+vpikzxMKDRemv/7VucfyCVwb923YmA3UMQ7AXkxbPjZ7x0qk3rtVnWzs7ENSLymSVS4+/2s89Zv3IN7kMjrWz0ykIQJ6Kmfl+8povrn/9dfYugiJ7fsyqY02aJn3EAj8bmWPdczAQ9gwCQmf7kDqq5nRLTw6jNKjeyvjoeMWZrRM+SwkYYa7cWcykkLEbRr75yZXNv94+ZO28QFCvuFBzedWZt//Vf+Hz8fOtMlOwJ+LconakbqpW1L14anWDqqknh2rNVbTckMc8H/CwWQ19kf18EmtBcSTyab/GS1JllbZ7R7jZcvu5oyu8eO5fjf/UWpPley7WLfLriZ+fqjz/0eVNRtJo7exAkNnQgP7+5i9bsnZ8mrJuTECytbPTTU4swYej1z4WMOrfJ179q+pi9w5iUBBlB+sj5vsx+Ha1qTsMhGbA4OMhM7yL9tSQeqpLL6Ro6vubP79z4eOVA5e+kPScWXZpsSJ3ruuWcR/rScPLp9dKtG3Wzg4EmYGO0K278Om1hqxvJnwe7hxi7ez0CAKQmZGTP01ZtzPvp0+ubNYRuq69ngYlv9R5D3cW+PeNbuHOg4HQPFzinJyCuJXpXehzaFQ3v3Tqzdzmgh2pm4b6DOi9vFkSG2e9O+K1oT4Dlh5fVSgp6b03ogha06SXFatM/zSNetLQtasQqC+iCFrbrJeXqk3fu7pO16vfe6O6+YWTq7kMzqbHPnDmiHvvjSwp3Dlke+oXNE0vPraySFLa+Rc2XJYSGtJ3rFvv5c1a+nYTxKYEP+GV9XGJ51Axz+vRk99PlJ/dmr1rbsz02ZHT+tz9ho4hAFkQOztYFLjm3Pv/TjDzPEhts741VyHJk2ua9Cwx887ihDQwKIw6qZHBwwQBHEEgVxzO53biW4D6BH2bsa1QqajQKKu0+jYDS8xgCu8shW9UEtoWPUvMdIkTuPYT8s3aUslpuvn+35/Pj5nZvriu3eDg7DVDX/6r6uLqc+/NjJgyL2bGI6cCExqy+nhz3LIguxkgczc4j9Cc6i9IpLeUsUsCO0gj08s3Xt1ao6x/a9groeIgS2WtV3Q8j6dKXrM248NBXokvJD3b8/m8Oqmh6lizvFjlmiA0tb/v/0HqJAZllUZRoWm7raJp2iXOya2/UBDAtaVFP8zJvucRahr1LTfkkjy5UU2KI/jCEJ4ggMv1YN37bdJAVa+T5MlbsuVsF2Zgmgff1wzhcH/hn3tvHXhn+KsJHnEPS2MH5d+iaf3o8n91hP7NYSt8Bd4dpCw/1EARdOjMjtJYGJxQf4etBUKapLM/Kw2e6imOevDXc7760n8zt00IGvNsv3kMzAprMZjXI09ElUG94dJGjVH7XvLqnix7WHe+teZ0i3eyi88o107OW9I06SW58pYbcspAuw8UuQ8QsV36xoyUzrODivh+RhXRnCVrvi4jtKRrf6FrvFDgx+nMpQxN0k1X26pPNrvGC4OmeHZ7ZXwdof/s6pfVitoNI9/w4Ll3kNI+yp8G9O9FR3ff3Pd03JwnwtMe2EGlbTHkbS5LXB1mU2NkYCC8w9YCIQBAWqCsTG9MeDX0nvZKm0626fq2Cnn16iEvxrhGWit75tWZE5Gi6e9v/pxedmrdiNdj3br8wUkDVfpLnVZiiHrGnyXqzqWDuk7XdL2tJVvO82Z7DhW7xDpZePOQ3mMfFfEdNJCVqBovS2XFaudYgccgsTCY142mPKmjivfWGjVk5NN+TEGXa+06ZcPbGR+FOQe/MmgZ61Fzee2p/GuV9Z9c2ULT1OtDXrp/KY/bu6oFARzfMbZ1dxAGwjtsMBACAPK2lHsnu7j2F5r+pAF9tOz09hs/PB4y7um4J/vKTPnO6PyJeLku85Mrm5+KmTkzcnLnlycm9VT+NxVcD3bITG8U71H0oghamq9ovCzVNOk9Bok9hzpbZXVE87KPiphQk03X2xouSXEW6jnU2S1RhLF7NoiPBjWnW5qutsUtD+rSxdOFmssbr21d1G/e1LDUzqS3j/JvR9H0HyVHd+XtmxExeV7MjPZB7Op6XcH2qgFrw3v4GzQ7GAjvsM1AKL2lrD7e3P+VEABAhbz6P9e+JijjqkEv9PU7gvfr0onYqG5+58LHnjz31UNe4jEevYkVRdC3dlSxnRmhs3zMeJNP26xvuCRtzpI5BXK9hruII/h99w5iX6+IldXahr8l0nylc6zAa5iLIMCcQ13qzrc2XWmLWx7M6MTOa0aK+Dbn+4yay+8lr+78fmd9vfwfqEnd8t/MbbXKhlcGLu3vEQcAKP6pluvF9h1jzqVKzQIGwjtsMxACGmR/VuKd5nxAe+hU5blF/eZNDp1oZ0NDTbp6IhpJ49acXVfqMt8Z8eojtjClQeEPNQDQvbSME2WgWnLkDX9LCS3pOczZY5C4M9WlremjFTFpoFqy5Y2XpISO9Brq7DFI3P1dQjtUdbSprVjVb1kQ2uEaKA2qpvcufubCEa8e+lKX9hTso+XfGRdqrnyZtSPGNWJx8MKqr6UD10b0tJneC2AgvMM2AyFFU6fS/27NljdMqFrcf4GQ5WTtHPWW7p2IF2ou/+f6N7Mjp86JeuJh1wf1GZLWG/K4F4J6+36eslrbeEkquakQR/E9hzp3776UtfS5ilhdr2u8LG3JkQuDeZ7DnC3QHC/+qRZloqGzHjrW8UzVhS2Z2+fHzJoRmdbVPQX7XPl3iZ40/FTwm/KEIcQ1YOyCoWwb2xIVwEDYzgYD4cXaqztu/OjKcpnx9+y4xSF8H5s7e8yo2ydis6b1g7//g6Lom0NX3L87jLpel/9NZfyKYLazhe6nElqyOVPWeFlKk8BjsNh9oKgb4ywsr69UxISObM2RN11tMygIj8Fij8Hi7o176gZST+V8Xho0xdMl7t7rUbVRs+n6tiJp6dvDV4WJg7tx8L5S/t1GqMlrHxZdHPPXNUXWM3FzUkPG9sa+ht0GA+EdNhUIsxpzd+b+pCf1i/svGOI9oPavVm2zPmyOrW+l1BM9OREpmt5768BvhX8sS3xufNDof54n6Bv/KfMd4+o+oPvTLbpNWaVpvNImyVM4BXHdB4qcY5xsbYDA3Wy8IqYpWl6qbroua7ulFIXzPAZbogl4P2WV5vZ31f1XhTKd/rm4yWm6+fHl/w72TlqW+CwbZ3XzyLZd/j1Xd65V06APm+tTKCn59sYPzZrWRf3mpfiPsJEbPTAQ3mEjgTCz4cb3+b/I9YqnY+ekBCSbzhKjisj6qGTgOxH2sV/XA/X8RCxpK//g0hd+Au9XBv3bNNGw6niztlkfudDPTHnsDtJASfIUzZkyVa3WJc7JLUEoDOXZ4IIaNlsRK6u1rTnylhwZS8RwSxK5J4p66S5gJ1WfaNY06CKf8QcA6Aj99twfzldffn3w8kHeiT05rM2Wv7lkf1ISOsfHKfDO0LasxtxdeXuVBtXC2DkpASMeuRhNb4OB8A7rBkKSJjOqL++9dcBIEfNjZo4JGHnPhVLh7mpRpMBziJ0sUXg/s5yIRtK4++a+9LLTyxKfHSUefmNjaf9VoRbrOuuYQUG05Mhac+S6NqNLnJNrPydhCM92piHaVkVMA2W1RpKnaM1VIDjiliB0SxRy3LrZ2DIviqCzPy4Jm+tTzi//9MqWWLeoFwc836VxMQ9kW+VvbooKTemvdYmr7x1Dm9lw44f8X1o0kjlR0yYGP9bt9nTPwUB4h7UCoUKvTC879XtxuhfPY070E0N9BjzwNnvbbWX1yeb4l/v2ivUdMOOJWCwt++TK5nGFE+JCI2OmdOeGTa/SSQyteQpJnkLbqhdHClxiBKIIPs6xcleELVTElIGSl6kl+UrpLSXOQV36CV37OfG8be7WeM31xoLjFbvidq4ctHSIt3nWuLeF8u89JfvquN4sn1EPnjVR0Fq479bvN5tvTQoZOy08teMleHoJDIR3WDgQ0oDOa751uOTE5frryX5DZ0SkdXyPnabozA3FMYsD7HUBaPNWBLJKde7Oks1x/30ietKcqGm2uQSdQUFIbymlBQp5mZrnxRZH8EXhfL4/xyodp1ariGmgbtTJilWyIpWiQsP34zjHCFxinNiutrhYBA3o4+Vnt+f8uKRoScxjIb6DzVZl23EgJHXU9fVFSW88Yk21BlXTweL04+VnYlwjJodOHOKTZMnRNDAQ3mGxQFirrD9def5E+V8sjPl46PgJwSmd7FepOt5M6sjgaV69nUOrMG9FcHNrhftAER1t/DJrR4Ws+sUBz5vryr03UAStKFO3FalkxSp9m9EpiCsM4TkFcfl+HIv1nVqyIqYpWtOoV5Rr5OVqeZkaY6GicL44nC8K59vgDLN2RZLSTZnbAAArBi7xVvoUfl89YG24ub4gOw6EjZelsmJV5NP+nUmsJw1nqy6kl56sVTaMCxw1NmhUhHNob+cQwEDYrrcDYbWi7mLtlb+qLkq00tH+IyYEp3T1C9ZJDbn/LR/0boQNDrXoOTOeiMoqbdGPNUlvhpkK6lpD9pbM7Z48j2WJi4JEAWZ5i95jVBHyco2iTC0vV+taDTxvtiCAK/Dn8P04bGdm742T7O2K2CA3qmp1ymqNslqrqtIynHBTvBeG8mzkJm4HWjSt23P3ZDbkLO6/cGLwGNPNi4JvK137Cz0Gmee2vR0Hwrwt5X5j3R62ecDD1CrrT1b8dboyAwHI6IDhI32HhruEdHV2ZufBQHhHbwRCA2m42XL7an3W5bpMLaEb4Tt4lP+wePeYbg+Ryt1UFvC4hyiMb8ZM2ggznoi3d1eLQnleI1zanyEo8o+Soz/m/zrcd/CiuLmuXJcOXm47SD2lrNaqTMGjRkvqKZ43m+fN5nqxuR4srifLjHcWzVsRk3pK26zXNOk1DTpVvU5drwM04Puy+f5cgT9HEMDtK4vvqI2afbcO/lFybErYxKeiZ3IZ/6zcJitRlx+sT3w9zCyVs70GQoPcmPN56aB1kd1uOhdJS89XX7pQc0VH6AZ7Jw32TkryjL/7izALGAjvMFcg1JOG25LivOaCG035tyXFwaIA05cX7myGy5m6v1q1EoNN7eNlLuY6EbUt+rwtFQPfCr9/KSyVQb331oHDpSceDxk3N3p6n1umx6gm1fVadb1O06jXNOq0TQYEAxw3FtuNyXZhsp2ZLDGDJWYwnRjdmLDYvfKnKdqoJHRSo15m1EsNOolB22rQthgIDcn1YHE8WDxPlily236z7x46Qn+o5OjPtw4O9Rm4qN889/vWagAA3PiizH+8u3OMGc5bew2E9RkSdYPOLHOgaxR1l+szr9Zn3WotChIGJHjGxblFxblFd2a14UeCgfCObgdCHaGvkFeVtlUUS8sKJSXVitpgUWCcW1R/j7h49xizfEn/vJf99o6a60Qs3V/PdML9Jzx0FEOrVvpj/q9nqy5MDp0wK3KqmC3s+Ztai0FJaJv1OolBJzHopUad1KCXGQ0KAudgTAHOMP3jYjgPwzkYzsEwFoqxUJyNARSYWpMIhpgmp6pUKj6fDwCgjDRlpAAApJ6iKZrQkpSBJnUkoaMILWlUEYSGNKpJo8JoUBJGNcngYSxnJkvEYDsz2C5MtguT48ZiiRh9aHm5e+gI3Z8lx/fd/j3ePWZR3NwA4UPnobbekNdfkPR70Qwjk+01EOZtKfcb5y6ONGcnloE0FLQW5TTdzGsuKJKWunNdI13Cwp1DQsRBoaIgPpPXjWOasfz7wDpSPSTXK1o0kmZNS4OqqVbZUKdsqFbUSrVtAUK/UHFQmDgkNXhsqDiw93ZHYjsz2WKGvEwjCuvOl233jCqiNVee9EZHa3C7cpxXDlz6VMzMnwp+W3h42big0XOipnnwbGt3tE5iCnCmABeG3HsyGJWEQUUYlYRBSRBq0qghtS0GUkeSOoo0UISWBDQgtCQAgCZpUk8BAGiaRhAEAIAyEJSBAgAwFoqgCM5BUSaKsTCcjeJcjCVm8H04DD7GcGIwBTiDj9nTNZnCoPy96OjvxUf6e8T9Z8z7Avh9FQAAIABJREFUj7yj7NLPqTK9SVmlNe9mF3bDIDdqm/Vmr6yYGDPBIy7BIw4AQNJkpbymUFJSJCk9W3WhXFbFxTkBQj8fgZePwMub7+nOdXXnuorYIostYdO3W4QSpfRY1RkURTVGLUmROlKnI/Qao1ZpUMn1SrleIdPJOAyOG8fFg+fmwXP3EXj5CrwDhL6ePHdLLotQ+1erzh57R81yRVZ3vlXToA97srP9MFJt26+Ff6SXnRrolTAnapplxqfZJnttkXRSvapxf+EfpysyRvgNnhs94/7tZB+m9q9WXYs+dHZPu/7ssvzN2C/aec2a1mp5ba2yvk7Z0KBuala3NmtalQalkOUkZDk5sZycmHwug8PCWDwGl42zno57EsAWYTuKplQGNYIgXAYHY2CumAsbZ/EYXD6DJ2Q7iVhOIraofXtJK3KNd8rdVB4y3cuersTNpelqW+isLvzqnDnipQnPLIidfaT05NsZH7tynGdEpI30H2YLXzRkARRNZzbkHCw+cru1JC103O60L104XRsF6jFAlPVxSdBULzte/rDbWnPlfuMsPTve1AQc4NX/7idJmpTpFDK9XKFXyvUKLaHTE3oNoWVj5p+W3bfrDj6Dtzh+gdXXGn0k05gIRYXm/g4xB6es0tAUaF/MsPN4DO6cqGmzIqf8XXvt9+L0LVk7JgaPSQsd7yuwt2Y31E6ibTtWfia99CSfyXsi/PH3ktewunVHgyHAhaG81htyj8F2u/xh9xiUhKbJ/P2i3YMhmAtH3NWrnO7p24GwD3GOFrTdVsJAeI+mq20eg8XdHqOBImiy35BkvyG1yvojpSdfPPWGD99zYvBjo/2Hd+/2O2SDDKThUt314+Vn8lsKR/kPWzfi9QiXnvaHewwW155pgYHwHm23lKJwvu2spmsxMBBaiDhSUPJLXWCatfNhS0gD1ZqnSHz93lV9u8FX4L004ZnF/Rdcrc86Xn52a/Z3iZ7xYwNHDvEeYMVFgaGeICgyuyn3bOWFi7VXw51DJgaPWTditbm+TXEkv/TXOk2TnusBT49/tBUqnaP72Awls4CB0EIE/hyjitC3GVniPjY3q/e05sqFwby7d4nrIQzBhvkMGuYzSG3UZFRfSi879emVLQO9Ekb6DR3iM8C8s2KgXmIgDZmNNy7UXPm79pqfk/do/xGL+y80e/8YgiLuA8XN19oCJ3ua98h9F03SsmJ1yHRHvLkAA6GlIEAcyW+7rfQc5mztrNiK1pzeuknDY3BTQ8amhoxV6JV/1149XZmx8drWSJewIT4DhngP6PzwQshiWjStpuWccppuhjuHJPsNediMeHNxTxQWbK8KTPPsu7MnzUtRoWG7MRkCRwwKjviZrUUcKWjJkcFAaEJoSWWV1rRXau9xYglMEVFH6LIa8y7XXf+t8E8EIAO9EpI84xM9+/W5pWrsiZbQ3WjKz27Mvd54Q6ptG+iVkBIwYvXQl3q+U2BncL3YGAtV1mgF/nBCIQCmftEuLi5qN2AgtBxxFL/st3qKoLuxmJb9keQpROE87L411XoJG2cP9x003HcQAKBKXnO94cbJinOfXf3Sg+fW3yO2n1tMnFtUX1nOtE9T6JX5rbfzmm/lNhdUyKqiXMOTPOLXDHkp3DnUYrOn27nECyW5chgITaS3VJ2fzmtnYCC0HJyDcb1YijK1KMIOF+DuqtY8hfsAkVXeOkDoFyD0mxk5maKpYmlZbnPBqcpzm65vY+GsWNeIKNeIKJfwUHEQHGVjFgRFlssqCyUlt1qLbkmKWzWSKNfwfm7RSxOeiXIJ670VnTrDtZ/T7e+qYe8oAEAvMxpVhMDPQa8JYCC0KHGUQHpbCQMhoSUVFerIhQ9dENIyUASNdAmLdAmbEzUNAFCjqLstKb7VWnymMqNCXu3F9wgTB4WJg0PEQcGiwD69wKklqY2aclllWVtlSVt5ibS8SlHrw/eMcAmNcYucGTklWBRgyUWdOsbzZiMYoqrV8h01ALRru60UR/Id9oIABkKLEoXxS/fXWTsX1ifNV4rC+La2roefk4+fk8/4oBQAAEGRFfKqEml5SVv533XXy9sqUQQNFgf4O/kGOPn5O/n4Onm7c90s35tna1q10lplfY2irlpRVymrrpRXq4zqQKF/iCgw3DlkUsjYEJFNt61d451acxUwEMpK1A57gxDAQGhhfD+2vs1oVBEMvkOXfGue3K2/TTewcBQLEweHif/Zo0CibauUV1fKayrl1RdqL9cq6mV6hTffw5vv5S3w8OR5ePLcPfnu7lxXuxyAozZqmtUtDermJnVzo6q5XtVYr2qqU9azcZavwCdA6Ovn5JPo0S9I5O/Bc+u9vVjNziVeWLi7OjDNw9oZsSoayEvVQQ48k8Shq2PLQ1DEKYgrL1O7xtt0GOhVlJGSl6nDn/K1dka6xrTaU5JnfPszOkJfr2o0/WtQNeU05TWpW5vVLXpS78Z1deU4u3FdxWyhK9dFxHISsYUillDEFgpZAjZu/sUSe8hAGhR6pUyvaNPJZHq5TKeQaKUSbZtEK5Vopc2aVgCAB9fNg+fuwXPz4ntEuYZ78z19BF59fXYm34dNU7SDz6zXNOowNurIU5xhILQ0YRhPXurQgVBequb7cHC2ra8Q+0hsnBUsCgi+b98fHaFv1rRKtdIWraRNK2vVSstklW06WZtWJtcrFAYlRdNOTL6Ayecz+Xwml8fgcXA2n8nj4GwmxuQzeUyUwcJZTIzJwpgYgnL+f2tvAfN/7i6rNWolor77GY1RQ9KUKQ8ERRAUoSV0RsqoI/Rao1ZPGjRGrdqo0RJalUGjMWqUBpXCoFIalCRFCVkCIctJxBY6s8UitpMzRxwk9HfhOLtyXdy4Ln094HVAHCVoK1Q6ciCUlapFYQ49cAEGQksThfILL9dYOxfWJL2tEkfZ86+OjbP8nXw6mLavJw1KvVJpUCkNKpVRozFqtIROaVDpCJ1MJ69TNuhJg4E06Em9gTSSNKU1agEANKBVhv8JexRFoej/3GflMrgYggIAWDiLgeIYinFxDgNlsHEWB2czcaaAxffku3NwNpfB5TO4fCZfwOQ7sQQc22ukWoxzlKD+gsRnVC9O3rdx8lK1q23fquhtMBBaGs+bTagJg9zIFDpoR0RboTJqUe/Oo7dxLIzJ4rr0fNqiXe6HZ3nCMF7RTzWknrK10VsWQgNFuTrE7nZL7RKH/OKtCwFOITx5mfrRKe2RtkVPEzTP03HbH5CtwZiowJ8rK1H9X3t3HiVFefcLvHqtrbtnYaaX6dlYBhBUBkEU8YqCKB4lonElLhETUSTHNedo9MW8Jjlq8M1Rk6tGAnpMbqJi5ETw3mg8ipKwCSpLBIZlmLV79pnurqrequv+0WYymbVnurpr+37+6ump6XpouvpbVc/veR6lG6KMSItgc1jthpxZrR+CUAGF09i+kwYNwp6jkaKznNopKgRDKDrL0XPUoEHYd4IrMHYHIYEgVERBjaP3hEGDsDs9bhdATYrPcvYcCyvdCmX0neIKphl98U4EoQIYNynGUrHehNINyTcxngo38IXTEYSgLrSbNJlNfDCmdEPyTiJCZ/iCKbotCc4QglAJJsJZRYcbeKXbkW+hU5yjnDZoSQKoW9EMR89xw10U8sGojbUYfH4PAkGoFNdkNlRvuCDsO8kVGv4mDKhTeoCv0q3It1A975qMQxJBqBDXZMaAQdh7Et3yoFIFU9nQaV5KSUo3JK9C9byz2uj3RQkEoVIcFbTQHhNjKaUbkj/JqCi0xzC7MaiTzWG1F9i4lqjSDcmr0BneNRlBiCBUiNlqYn1UpElQuiH5EzrFO6sYLEoMqlVYwxqqnDsRTiZ5kXEbd265frIFIc/zjY2NqdSIlziiKDY0NESjxjrhGoXTYHdH+05GCmvQGwHqVTCN7TtpoNGEffW8azKDQb2EXEH40ksv+f3+JUuWzJw58+jRo0M32L9//5QpU5YtW1ZWVvbmm2/KslOtc1UzoXoDnX72nsRwJVC1gqls6AwviUbpJgzXc7gvmiZDEDY0NPzkJz/ZtWvXyZMnV61a9cADDwzdZs2aNY888khdXd1HH320du3arq6u7Perda4pbPiMYJDO+SQvRrvijnJ0EIJ6WRkLNclunA4LVMr0kyEI33rrrcsuu+yss84iCGLt2rWffPJJW1vbwA2OHj36zTff/OAHPyAIYv78+XPmzNm6dWv2+9U6G2uxOS18myHG8Pad5FyTWZMFd2FA1Qqnsb3GGESRSqT4YMyJ4jWCIGRZfaK+vn769Onpx2632+l0NjQ0eDyegRv4/X6G+fbUo6am5syZMyO9WiqVOnToUDAYTP84ZcqUoqKikTaWxGSiuVE0a7Xkhy1N9nzVYEtqtf0iz8eZjM4ou79OOiaZ4k0nct0kQ8n8/YcMsYWp4FeiZ0ZvJhtr+v0Pt0pUsZRsO6V0Q8bNZLbY/FPkfU0ZgjAcDpeU/HspL5Zl+/r6Bm1A0/QoGwwUi8V+/OMf22zfLlH09NNPX3rppSNtzLe3pv784oRbrjz+7O76EsuRHUq3Y4JSqZSQ2VlIb9dNpc7Pu04Gc90kQ8n8/YcMpSQy0n5755/+t4kYu89C0+9/L3euRSzo+tNOpRsybibawX7/vwiCiEQyqmxiGMZiGWMZcBmC0O129/T09P/Y09Mz8HIwvUFvb+/ADWbNmjXSq9E0/eGHH/r9Iy5qOpDVWml/5Ndj/iNVy9konNjS4nvkJqUbMkEZrocnxlNn1h+rfOy/MHZCXliPMBfanjtR8L3/Yf1jrxSm6fc/9Icm70yne/6tSjckK3K9/zKczsyZM2ffvn3px0eOHDGbzVOnTh24waxZs9rb21taWtI/7tu3b86cOdnvVwdYPxXtiItxnQ+rDzfwrJ9CCoImOKsNMa4p3CA4K9FB+C0ZgvDGG29sbGx8/vnnDx48+OCDD951110syxIE8dhjjz377LMEQXg8nhtuuGHdunWHDh166qmnJEm66qqrst+vDpgsJsZL6n4yi3A9Zq8AzXBNZkJndB6ESV5M8iJdiqH035IhCFmW/fjjj3ft2rVmzZr58+f/8pe/TD9fXl7u8/nSj1955ZXKysof/vCHdXV1H374odVq9MnO+zkqmXCjzo+60BlUaYNmuKqZsN6DMNwoOMppDKXvJ08gnXPOOe+9996gJ9etW9f/2OVyvfiilqtacsZZSet8RVCJCDcK01chCEEb6FJSjKdivQmy0KZ0W3Il0sg7cF90AK2WPOmGo4ION+p5AC8XjNpYKxY8A80wwHKh4SYBIwgHQhAqjHGTSU5McqLSDcmV8BnehfuioCm6XyUt0iTginAgBKHSTISjnArrd1anUD3vRKUMaIq+uwljPQlCInR843cCEITKc1QwOp7eMIQrQtAaRwXNB2MpnY5rCuNycAgEofIclbReC0eTnJiMiIwHVdqgJWabmfGSEZ2Oa4qgg3AIBKHyHOWUXocShpsERwWqtEF7HJW0Xu/TRNJHJQyAIFQeVWQX46lEOKl0Q+QXRpU2aJOzgtFrOTfXIrB+HJX/AUGoAibC4ad0eR8m0ohpnECT9NphEetJEGaT3YXhTP8BQagKrJ/mWnR4+hlpxk0Y0CS9jmuKNAtYH3soBKEqOPxUpFlvV4So0gYNMxGsn4o06+30lGuJOsrHXljDaBCEqsCW0xHdXRGGGwVnFQZOgFY5K3XYTRhpERzoIBwCQagKjJtMhJPJqK7uw0SaUCkDGqbLbsJIc5TFFeEQCEJ1MBGsT2+DKMKolAEtc1bSEX1dESYiyVQ8RRXZlW6I6iAI1YLV2WhCieCao7gJA9pFFtoIExHrTSjdENlEmqOOcgrjeodCEKqFw0/rqWeeb4/ZnBYra1G6IQAT56jQ1UUh1yKwKBkdDoJQLdhyWk+Fo5EmVGmD5jkqaD1NiB9piTrK0EE4DAShWrBeMtodTyV0Ms9vpBnnnqB5jnJdDfDFFeFIEIRqYbKY6BI7H4wp3RB5YLgS6IDDT+lmxlExnor3JelSVMoMA0GoImwZxbXq4u6oRHAtqJQBzbMX2Exmkz7qZfjWKO0lTWaUygwDQagiuglCoTNmZSxWBpUyoHlsOa2Pcm6uFR2EI0IQqghbRnEBXRxyLRi0Czrh0MtEa1xrlPHhqBweglBFdHNFGMEIQtALVi/jmrjWKIsrwhEgCFXE5rCaLXrokIg0C6iUAX3QybrZEsEFEYQjQhCqC+vXw0UhVv4E3aCK9bBudrQ7bqUtVhrd9sNDEKoL69N8EGLlT9AVXaybzbVGWXQQjgxBqC466CbEyp+gMzpYNxsdhKNDEKoLo/0gxFB60BkdrJuNIBwdglBdGDcZ601oeqK1CDoIQV90sG42F0AQjgZBqC4mi4ku1fZEa1xr1OHHIQf6QZfaE6GkGNXq6akYSyVCSaoEk6uNCEGoOpqul0kKYpIXqWIccqAfJrOJ9pJ8UKtHJR+I0h5MrjYaBKHqsD5Ku4fct8VpOOJAXzR9eor7omNCEKoO46O4gFZvjaJPHnSJLaMimg1CPhhjvKTSrVA1BKHqMF6S1+yMo1xrlEUHIeiOpsc1cQEMIhwDglB1yEJbSpQSEU3OZMG14IoQdIgto/hAlJCUbseE8MEYptseHYJQjRgvpcXCUSklCe0xxotDDvTGSlusjCXaFVe6IeMWDycJSbI7MdPTaBCEasR6SS2uxyR0xO0FVguJDxXokEa7CfkAVl8aG76z1IjxafKKEJUyoGOsn+a1GITBGDoIx4QgVCPWp8l6Ga41ypZhThnQJ43Wy3CBKEpGx4QgVCPGR3Ea7JnHFSHomEObQYhbo5lAEKqRlbZYKIvmVujlWgQEIegVNcme4JJJQVS6IeMhEXxbjPHginAMCEKVYn0aq5dJcmIqIZGFNqUbApAbJoLxUdrqs4h2x60M1uMdG4JQpTR3yHHpOzCYXA30i/VRnKaq2DCUPkMIQpVivVo85HAHBvSM9ZHaKhzFUPoMIQhVitFa4SgfxFB60DnGS2mrw4IPomQ0IwhClWLcpNARl1KaqRxFySjoHltG8YGYhsq5+UCMxelpBhCEKmW2m+0uq2amdJIIvg3nnqBzVsZitpu0Us4tpSShK067sTjo2BCE6sV4Sa3MLxPtiVtpFKeB/rFlmrk7Gu2MkwVWsw1f8mPDe6ReGpp6+9v1eAH0TkPl3FwA3faZQhCqF+sltbJUPR9AcRoYAquddbNRKZM5BKF6aemKEGMnwBg0dEWIQu7MIQjVi/aQQpc2CkcxnyEYBJM+KkUtHJW4IswYglC9zFYTWWATOtReOJpKStGeBOPGIQf6Z7aayEKb0K72WzWSKEW7E3QpjsqMIAhVTROFo0JbjCq2myyYXQ0MQRPdhEJHnCyyma04KjOCIFQ1Rgv1MlwQHYRgIIxPA0clH4yyuC+aMQShqmmiXgYlo2AorBYmWuMwy+h4IAhVTRNXhOiTB0NhfBrosMBROS4IQlVj3GS0O6HyEjUs9QKGQk2yJ8JJMZZSuiGjwdiJcUEQqprJYqKKVV2iJsZSSV6kijGfIRiFyWyi3STfpt6jMpWUYj0JugRHZaYQhGrHeFR9yPGBKO0hsR4vGArjVfUqaUJ7jJpkQyF35hCEakd7KTUHIRfEOi9gOCqvYuPbYowHR+U4IAjVjvGoumceffJgQKyP5NvUe0XIt8VwVI4LglDtVF44ij55MCDGq+ox9Tg9HS8EodqpvHCUx3TbYDxkoS0VT4mCSgtH+SBujY4PglDtTBYTWWQTOtU442iST6WSkr3ApnRDAPLLRDAeMtquxqNSEqVYT4IuRcnoOCAINUC1d0ejbXHMXgHGxPgooT2hdCuGwbfHqGKUjI4PglADVFsvE+tIYD5DMCbGS8Y61BiEAiplxg9BqAGMR6UjKIT2BCplwJgYLxVtU2MQon5tAhCEGqDqW6M49wRDYn2koMo+Qr4tynhwVI6PbEEoCML+/fsbGxtH2qCzs/PAgQPNzc1y7dE4aDcZ7VJj4Wi0I4EgBGOyOawmkykeSirdkMH4YAxBOF7yBOGBAwemTp364IMPnn/++Y8//vjQDW655Zbp06evXbt27ty51113XTyuxjMp1fp2UWyVFY7G+xJmi8nmsCrdEABlUG6b2m7VpBemp7Aw/TjJE4SPPvroAw888Pe///3LL7985ZVXvvnmm0EbrF69OhgM7t279/Tp00eOHHn99ddl2a9xqHCpei4Qozwo0QbjIkttajsqhfYYVYyF6cdNhiBsa2v77LPP7r77boIg/H7/8uXL33nnnUHbXHHFFXa7nSAIp9M5e/bstra27PdrKCrsJuTbomQJLgfBuGiPTW0r9PJtuC86ETJ8kTU1NTmdzpKSkvSPkydPHqWn8NSpU5999tlPf/rTkTYQRXHnzp39r3buuee63e7sG6l1jIfqOhJSuhX/gQ/EKC+uCMG4KLe977DKjkqUjE5IpkH4zDPPHD58eNCT8+fPf/jhh3meT1/tpVEUxXHcsC/S3d19/fXXP/roo7W1tSPtKB6P//rXvybJb09qHnvssYsuumikjQVBsNvtFoslw3+FhrnESKsQiUSUbse/hVu44hpaVU0yGo7jTCbcBFNMypHgWoVIOKKeZchCLZHC2axBjsoMP/8Mw5jNY9z7zDQIFy1aNHXq1EFP+v1+giA8Hk9vb28qlUrvrKury+v1Dn2Fvr6+5cuXX3755U888cQoO6Jp+p133km/8pgsFotBgpCplk70BFmaVcuEERIR60wWVrkcDofSTTEuSZLw/itIkiQrHbIlSbJILbMMxjsDxVUFjMMQF4Uyfv4zDcJLLrlkpF9VV1cXFhbu3bt34cKFBEH84x//eOihhwZtw3Hcd77zndra2ueff37CbTWy/sJRlXQARLvjVsZiITEOFQyN9ZF8MKqSIETJ6ITJ8EVGkuTatWvvv//+v/71r48//nhHR8cNN9xAEMSuXbsqKyvT29xwww2NjY3z5s3buHHja6+99vnnn2e/X6NRVeEouiIACIJgPCpajwkloxMmT9Xf+vXrJ02a9PLLL/v9/s8++4yiKIIgPB7PzTffnN5g8eLFc+bMqa+vT//YXwsDmWM86bVAXUo3hCAIgg9i9SUAgvGSfaeGL4nIP5SMTphJktQ1X0l5efnevXsz7CM0ULEMQXR82dd1JDTzjgqlG0IQBHH8/zQXzXDQMyxOp1PpthhXOBzG+6+gcDhM9FhPvdta+/Dg+glFNP61nSCIyuVGKbOX8fOPPh7NUNet0QCWwAYgGC/Jt8eklCouJ/g2HJUThCDUDNpNRrviaphxVEpJQkecduOQA6Oz2M12hzXWrYplKDDL6IQhCDXDbFXLUvXRzri9wGqx48MDQDA+Ug3zy6RLRnF6OjH4LtMSlazQywViLBamByAIgiAYL6WGo1LAwvRZQBBqyb8KRxXGB6MMghCAIAiCYL2quCLksTB9FhCEWqKSc08+EGVxyAEQBEEQjI9Sw4T4fDDGeHB6OkEIQi1RyRoUHEbTA/xLet3sVFLhKjaUjGYDQaglaliqPpWUYj0JuhTrTgAQBEGYrSaq2CZ0KHyrhgvg9HTiEIRa8m3haIeShaN8MEqX2tEnD9CP8VG8ot2EOD3NEoJQY9KT/CrYAMwyCjCI4pNdCO0xahJOTycOQagxjJfiFD3kMMsowCCsj1K2cJQP4KjMCoJQYxgvqexNGHRFAAyieDk36teyhCDUGMarcK02H4gyOPcEGIAusSciSTGWUqoBfBAlo1lBEGoMXWqP9SVTCWUOuWRUTAoiVYQ+eYABTARdqmQ3IY/7NNlBEGqMyWyiS+x8mzKHHB+IMV6SQJc8wH9ilKtiE+OpRCRJTbIpsnd9QBBqj4IdEnwwillGAYZileuz4IMx2k2azDg/nTgEofawPsXqZVApAzAsxkdxrYqdnqKDMEsIQu1hvJRSt0a51ihbhiAEGIwto7hWQZFdo4MwewhC7WGUm+0eC9MDDMvushImUzyczP+u+TYMIswWglB7qGJ7kheTUTHP+431JMw2k81hzfN+ATSB9ZF8qwJnqOiwyB6CUINMBO0hhbzXy3AB3BcFGBHjo7i8B2FSEMWYSBaiZDQrCEJNUmRKJ64VJaMAI1LyqETFaHYQhJqkyCHHB7AwPcCI2DKcnmoVglCT2DIFarVxaxRgFIyXFDrieV4ulA9EGRyVWUMQahJbRvGtUSKPR1wqKUW7ErQbxWkAwzPbzGRhvpcL5QK4IpQBglCTrIzFbDfFehN52yPfFqNL7GYr+iIARsTmuV5GSq8PitPTbCEItYoty+shx7eigxBgDKwvr2N8o11xG2u10pa87VGvEIRaleeeeQ4rfwKMhSmj8jn9IdeKNdHkgSDUqjwPWkKlDMCY8nxrFEelXBCEWsX68nruiVujAGOiiu3JqJjk8zTrE49KGZkgCLWK8ZDR7kR+VuiNh5KplITZKwDGYMrrRWEEp6cyQRBqlcliokvztEIv1yI4/HQedgSgdY5yOtKSj2UoxHgqEUrSpfY87Ev3EIQalrf5ZSItUYcfJ54AY2PLKK4lH0clH8B6vLJBEGoY46PyM9s914I+eYCMsH4q0pyXoxIdhPJBEGqYw09F8nLuGWkR2HLcGgUYG+ujot3xVDLn0z5xLQKL+zQyQRBqGFtOcy05n2hNjKIrAiBTJouJLrHnoaI70hJFEMoFQahhNtZiIc3RntzObci1CoyPQlcEQIYcfjrnt2okgg8gCGWDINQ21k9xOe6QwIknwLiwforLceEo3x6zOa1WCpOryQNBqG0Of85rtblWlIwCjAPrz3nhKEY0yQtBqG35OOSaBRaHHEDGHH6aC+S28z7SEnWU4/RUNghCbXOU05HmHF4RSqLEd8RZrPMCkDELZbY5rEJHDie7wOmpvBCE2kYW2VKiFA8l3WZEAAASrElEQVQnc/T6fFuMKrKZ7ficAIxDroc2RVrRcy8nfMFpnqOM5nJ2URhpEhwVOPEEGB9HRQ5v1US742aLye605uj1DQhBqHmO8hyee4YbBWclghBgfBwVdKQxV0HItUQdmOBCVghCzWP9dO5qtSNNgqOCydGLA+iVo5KOtAhSKicFM5GWKGZ6kheCUPMc5bma2zCVlIT2GFuGShmA8bFSFrvTKnTkZLILrlnAiCZ5IQg1jy4lE1wyKci/FijXGqVL7WYbPiQA4+aoYHJ0dzTSLODWqLzwHad9JsJRnpMOiUgj76jEfVGAiXBW0uEmXvaXjfUmCIkgi7BKtpwQhHrgrGLCDfIfcqiUAZiwHNXLhBsEZxVOT2WGINQDZyUdzsUVIcZOAEyUw09xwZjs6zFFGnkHTk/lhiDUA2cVI3sQirFUrDfBYE4ZgAkx2810iZ0PylzIFm4UnFUIQpkhCPXA7rKaraZot5wlapEmgS3D6ksAEyf73VEpJaFSJhcQhDrhqKTDDXIecmHcFwXIjrOSDjfJeVTywRhZaLPSWH1JZghCnXBWMpFGOetlIo08KmUAsuGolLmKLdwooIMwFxCEOuGU+4owVM+7JrMyviCA0bA+MtabSHKyjfENN/AoGc0FBKFOOCppLhCVRHlK1ISOOGEyYawSQDZMZpOzkgmdke2iMIIRTbmBINQJi91MFdu5gDwlaqF6rmAqTjwBsuWawoTq5QlCMZaKdsVZHyZXkx+CUD+cVXRYpnPP8BneVY37ogDZck1mQ6c5WV4q3CiwfspkQSG3/BCE+uGayvadkicI+07zrim4IgTIlquK5gJRWYbVh05xBVNxepoTCEL9KJjK9p3iiKyPuAQnJkJJDKUHyJ7ZbmY8pCwV3X0IwpxBEOoHWWiz2M18eyzL1wmd5pyTGQylB5CFazKbfTdhKilFmgVnNe7T5ASCUFcKprKhU9l2SITqeddkHG8A8nBNZvpOZxuE4Qae8ZIWEt/YOYG3VVdc05g+GYKQQxACyMU1hQnX81muVo/7ojmFINSVb7sJsyDGUnwwhrFKAHKxOaw2l5UPZNVnETrFIwhzB0GoK1Sx3WQxCR0Tn3277yTnrGKwKj2AjIpmOHqORyb855IohRt5J+7T5Ay+7/Qmy4vCnuPhopkOGdsDAEUzHb3HwxP+83CjQLtJK4W5tnMFQag32QbhsUjhDAQhgJxcU9lwkyDGUhP7875TXMEU3BfNIQSh3hTOcPQcC0+sZz7aHU/FU6wXczgByMliNzsrJl7I1vNNGKenOSVPEKZSqd27d3/wwQe9vb2jbBaLxQ4cONDV1SXLTmFYZKGNLLBNbCWKnqORoplOAgMIAeRWONPRc2wi3YRJTuSDsYJpuCLMIRmCUBTFa6655p577nnttddmzpx5+PDhkbZcv379ggULtm/fnv1OYRTFs5zd30ykQ6LnGDoIAXKiaMYEuwm7j4YLprNmK85Pc0iGINy+ffvJkyf37dv3l7/85Z577nnqqaeG3Wzfvn07d+6cO3du9nuE0RXPdnX/MzTev5JEKXSax4knQC6wPiq9fMR4/7D7n6Hi2c5cNAn6yRCEW7duve6662iaJghi1apV27ZtSyaTg7aJx+P33nvvyy+/bLGg8CnnnJV0ghOj3eM75EL1PF1qtzmsOWoVgKGZiKIZzvHeHZVEqbeOK56JIMwtGb71mpub58+fn35cWVmZTCaDwWB5efnAbX7+85+vWLGitrZ2zFdLJpPvv/9+cXFx+scLL7xw0EsNJIqiKMq2+rOeFM10dB7u811cnPmfdBzsLZzlGNf7ifdfWXj/lTXe979wFhv4e7f7woLM/6TvBEeV2s2MCf/RQ2X4/pvNZpNpjBvLGQXhwYMHH3744aHPb968uaqqKh6P22zfLmWefhCL/cccCocOHdq6desXX3yRyb5EUdy2bVv6+pIgCI/HU1paOtLGsVhMkiRcZQ7lqCE79oaKz8/0PqeUkjoPhmasKRv0fze6eDw+ru1BXnj/lTXe95+ZYuO3RMPtnL0g0yuQzsN9rhk0/peHleH7T1GUPEFYXV29fv36oc+nI8rr9XZ2dqaf6ejoIAjC5/MN3OwXv/iF3+//6U9/ShBEU1PTli1bnE7n9ddfP+y+SJLcuHGj3+/PpGEmk8lutyMIhyLPoRrf67SbSCud0ZvTczRMl5BF/nGcqxIEIYoiw2C2C8Xg/VfWBN7/SecUcMfihZe5MtpaIkLHm8+6u5JhMKJpGDJ+/jMKwoKCgsWLF4/020WLFr3//vtPPPEEQRCffvrp3LlzBzVu9erVDQ0N6cckSZaWlpaUlGTRZhibhTQXznR0fNnnW5TR3dGOL/vc540vBQFgvErPKzyzLei/LKMvwL5TnIU0sz6kYM7J0Ed4xx13PPfcc4888sisWbOefPLJF154If38okWLrrvuukcfffTKK6/s33jTpk2XXnrpJZdckv1+YXTeC4vq3w9mEoSpeKr7m/Dk73jz0CoAIyucxsbDSb4txnjGXvg6uKfHc2FRHloFMlSNFhUV7dmzx26379+/f/PmzTfffHP6+fvuu29o4N1///3z5s3LfqcwpsIahxhPRZrGHlnfdSTsrKJtTtSLAuSYiSitdXV+1TfmhklO7Dkads8vzEOjQJ7vvsrKymeeeWbQk7fddtvQLe+44w5Z9ghjMxHeC4qCe3qmVYyxplL7/p7S83C8AeRD6bzCY280VVxRajKPVsHRfqC3eLYzwz5+yBLmGtUz94Kizq/7Rp/qN9Is8G2xklp0EALkg6OcJottHV+OcVEY3NPtxX3RfEEQ6pndaS2YxrbvH20C2Ma/tpcvLcUETgB5U3mlu+lv7aPMjN93kpNShGsypnnKEwShzlUudzd+1J7ghh92GmkWuEDUswAnngD5UzCVtReMeFEopaTTfwlULXdj+vu8QRDqHOujSuYUNPy/tmF/i8tBAEVUXeVp/Gj4i8LA37ttrBW9FfmEINS/qqvc3UdCkcbB5aNt+3r49hguBwHyzzWZoUvsDR8MPkONh5NNH3dMvd437F9BjiAI9c9KW6qv8db9qXngzPfd34Qb/m/b7B9W4XIQQBEzbqvo/ibc8lln/zNJTqz7Q5PngiLaPfYoQ5ARho4Zgnt+YVIQD754espKL1lkjzQJTR93zPpBFV2K4w1AGVbGMntN9aFfn04KqYKpDCERJ95uKZlTULXcrXTTDAdBaBRl/2uSazJz4q0Ws9XE+umzVlc6K8cYXwgAOUUW2s5eUx34R3fTRx2xvsTU75YVz8KKSwpAEBqIo5ye++g0pVsBAP9Gu8kp16FHUGHoIwQAAENDEAIAgKEhCAEAwNC0HYSdnZ3hcFjpVhhXIBCIRqNKt8K4mpqaRHH4OYMg11KpVGNjo9KtMK5YLNba2irXq2k7CJ944on33ntP6VYY15133nngwAGlW2Fcixcv7u7uVroVBtXX13fxxRcr3Qrj+vrrr4dd4GhitB2EAAAAWUIQAgCAoSEIAQDA0EySNOKaWIq49tprDx06ZDZnlNDhcNhms1EUletWwbB6e3tZlrXZbEo3xKC6urqKiooyPFhAXpIkdXV1lZSUKN0Qg0omk+FwuKho7DUDtm/fftZZZ42+jeqCMBwOd3R0KN0KAADQg/LycrvdPvo2qgtCAACAfMJNFQAAMDQEIQAAGBqCEAAADA1BCAAAhqaZ9Qg7Ojo2bdrU0dFx9dVXL1myZOgGoii++eabhw8fnjlz5l133YWafhl1dHRs27bt6NGjhYWFN910U01NzaANUqnU7373u/4fzz777Isuuii/bdSzHTt21NXVpR9brdbVq1cP3aa1tXXz5s29vb0rV67E1F/y2rhx48CiwqEf77q6uh07dvT/uHLlSrcbq8xnpb29/cCBA01NTZdeeun06dP7n29ubn799ddDodB3v/vdCy+8cOgfxuPxTZs2nThxora29rbbbstwcJE2rgh5nl+4cOHx48erqqpWrVr19ttvD93m/vvvf/nll2tqav7whz98//vfz3sb9WzdunUffvihz+drb2+vra3dtWvXoA2SyeSaNWuOHTt2+vTp06dPd3V1KdJOvXrjjTfefvvt9HtbX18/dIPe3t4FCxa0tLSUl5dfe+21H3zwQf4bqWP19fWn/+Whhx46fPjwoA127969YcOG/m1isZgi7dSTpUuXPv30048//vju3bv7n+zs7Dz//PM7Ojp8Pt9VV131t7/9begf3nrrrW+//XZNTc2LL7748MMPZ7o/SQs2b958/vnnp1IpSZL++Mc/nnvuuYM2CAQCJEk2NTVJktTd3U1R1KlTpxRoqE4JgtD/eM2aNatXrx60QfrID4fD+W2XUdx5552/+tWvRtnghRdeWLJkSfrxq6++umjRory0y3D2799P03RPT8+g5994440VK1Yo0iS9EkVRkqQLLrjgjTfe6H/yueeeu+qqq9KPX3zxxf7PfL+jR48yDNPX1ydJUn19PU3TnZ2dmexOG1eEn3/++bJly0wmE0EQy5YtO3ToUE9Pz8ANdu/ePW3atPLycoIgioqKzjvvvJ07dyrTVj0aOHdPNBp1OBzDbvbb3/72pZde+vLLL/PVLgPZs2fPhg0btmzZkkgkhv72888/v+KKK9KPly1btnv37mE3gyxt2rTpxhtvLCwsHPqr5ubmDRs2bN68GfOByGLYW5qDPuc7d+5MpVKDNliwYIHL5SIIorq6uqKiYu/evRntLusG50MgECgtLU0/njRpktVqDQQCI21AEITH45FxqSrot3v37q1bt/7oRz8a+qulS5d2dXUdPXp0yZIlGzZsyH/bdKyysrK0tLS7u/uZZ55ZsGABz/ODNhj4+Xe73alUKhgM5r2ZOicIwltvvTVsB21BQcHZZ5/d19e3devWmTNnDr13CrIY9DlPJBKdnZ0DNwgGgwODwO12ZxgE2iiWsVqtyWQy/VgURVEUB02ZY7VaB65QmkgkxpxTB8br+PHjN95448aNG6dNmzboV3a7/eOPP04/vvXWWy+//PK1a9eyLJv3NurT008/3f/gvPPO27x587p16wZuMPAAST/A5192f/7zn4uKii655JKhv1q5cuXKlSvTj++7777//u//fvfdd/PbOkMY83M+4SDQxhWh3+/vD/b0A5/PN2iDlpaW/h9bWlrKysry2ULdq6urW7p06bPPPnvTTTeNvuXChQuTyWRzc3N+GmYoNpttwYIFQ+tlBh4gLS0tNpsNk0HLbvPmzXfffXe6g2YUF1100enTp/PTJKMZ9DlnGGbQbeoJB4E2gnDFihXbtm2LRqMEQbz77rtLlixJX218/fXXDQ0NBEEsXry4s7MzvVr6iRMnjh071n8rGbLX0NBw5ZVXrl+/ftCS0F988UX6YycIQv+T27dvZximuro6z43Usf63NxwO79ixY/bs2QRBxOPxTz75JF2mtGLFiq1bt6b7Bbds2XL11VdbLBYFG6w/9fX1O3fuvP322/uf4Xn+k08+SV+X9P8HSZL0wQcfnH322cq0Uu9WrFjx3nvvpa/5tmzZsmLFivTzX3zxRTogr7zyyoMHD6bPFPfs2cPz/KJFizJ6aVkqfHItmUxeccUV8+bNu/322ydNmrRr167080uWLPnZz36WfvzCCy/4fL7Vq1dXVFT0PwmyWL58Ocuy8/7l3nvvTT8/d+7c3/zmN5Ikbdy4cfbs2d/73veWL1/ucrl+//vfK9pevSkuLl6xYsWqVat8Pt8111yTSCQkSUof+WfOnJEkKRqNXnzxxQsXLly1alVJSclXX32ldJP15oknnrj66qsHPnPs2DGCILq6uiRJWr58+dKlS2+77bZzzjmnpqamoaFBoWbqx0MPPTRv3jyWZaurq+fNm5f+zud5/oILLrj44otvueUWt9t95MiR9Mbnnnvuq6++mn785JNPVlVVrV692uv1vvLKKxnuTjOrT4iiuGPHjs7OzsWLF3u93vSTx48fdzqd/Re/R44cSQ+onzt3rnIt1aG6urpwONz/o9PpTA9x/ec//1laWprutT5w4EB9fX1BQcH8+fMxmlheZ86c+frrr+Px+PTp02tra9NPJhKJr776qra2Nt0LkkgkPv30097e3ssuu2xgvQDIoq6uzuVy9X/zEAQRjUYPHTo0b948i8XS3d29b9++np4ev9+/cOFCzOaRvVOnTvX29vb/WFNTk64FTd8ICYVCS5cunTRpUvq3R44c8Xg8/R/7AwcO1NXVzZkzZ9asWRnuTjNBCAAAkAva6CMEAADIEQQhAAAYGoIQAAAMDUEIAACGhiAEAABDQxACAIChIQgBAMDQEIQAAGBoCEIAADA0BCEAABgaghAAAAzt/wPn/Es426/WpwAAAABJRU5ErkJggg==", "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "execution_count": 10 } ], "cell_type": "code", "source": [ "ρ = scfres.ρ[:, 1, 1, 1] # converged density, first spin component\n", "ψ_fourier = scfres.ψ[1][:, 1]; # first k-point, all G components, first eigenvector\n", "ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n", "ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\n", "\n", "using Plots\n", "x = a * vec(first.(DFTK.r_vectors(basis)))\n", "p = plot(x, real.(ψ), label=\"real(ψ)\")\n", "plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n", "plot!(p, x, ρ, label=\"ρ\")\n", "plot!(p, x, tot_local_pot, label=\"tot local pot\")" ], "metadata": {}, "execution_count": 10 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.3" }, "kernelspec": { "name": "julia-1.6", "display_name": "Julia 1.6.3", "language": "julia" } }, "nbformat": 4 }