{ "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", " LocalNonlinearity(ρ -> 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.143657376032 NaN 3.81e-01 0.80 7.0\n", " 2 -0.156049075383 -1.24e-02 7.89e-02 0.80 1.0\n", " 3 -0.156758202565 -7.09e-04 2.75e-02 0.80 2.0\n", " 4 -0.156930937634 -1.73e-04 1.27e-02 0.80 2.0\n", " 5 -0.156865730718 6.52e-05 1.53e-02 0.80 2.0\n", " 6 -0.157056387519 -1.91e-04 1.99e-04 0.80 1.0\n", " 7 -0.157056406771 -1.93e-08 2.33e-05 0.80 2.0\n", " 8 -0.157056406859 -8.83e-11 1.17e-05 0.80 1.0\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380281 \n AtomicLocal -0.3163440\n LocalNonlinearity 0.1212595 \n\n total -0.157056406859" }, "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.055663643422245127, 0.0, 0.0]\n [0.055663860923571756, 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+gvaeTAAAgAElEQVR4nOzdd3xTVd8A8HNHdtIkTfege7eUtpRdoOxCQTaCgILwyKM4cIHgQMHHifKI4gAEFHEgiELZG2SX0tLSTfdu9k7ueP8Ib+VhlNKmSZqc74c/0nByc3Jyc373nHsGQtM0gCAIgiBXhdo7AxAEQRBkTzAQQhAEQS4NBkIIgiDIpcFACEEQBLk0GAghCIIglwYDIQRBEOTSYCCEIAiCXBoMhBAEQZBLg4EQgiAIcmkwEEIQBEEuzeECYV5entls7mBiiqLgEnF2RJKkvbPg0mD52xcsf/uyYvk7XCAcP358c3NzBxMbjUaKoro1P1A7dDqdvbPg0mD52xcsf/uyYvk7XCCEIAiCIFuCgRCCIAhyaTAQQhAEQS4NBkIIgiDIpcFACEEQBLk0GAghCIIglwYDIQRBEOTSYCCEIAiCXBpu7wxA3ctMgWIlXaKkK9XATAGSBnwG8OaAcDektzvCgBdCEGQ/JA2KFHSxkq7XAbkRMFGAIKAXD0QKkWgRwoXVs63AknZOjXrwSzl1pI76u5H25SJRIiRMAJgYwBDQpAfnm0CRgipT0X0kyORgdHowEixA7J1lCHIVLQbwRyW1u4K62Ez7cpE4MeLHBWIW0BKApMDVFlCipMrVdJIEGe2PPh6KRAjhz7N7wUDobI7W0V8UkH830VOC0EVR6E/DUTHr/il1BDjfRO+upPr/RfZ2R5bFYxmBCPzBQdBdPvnkkytXrtz7PEmSGIY90qGUJlCtoaVG4MEGXmwwkoXgKAAANAHQdEeySADCaKAwgt1G+is94OIggAe8Oa7+8xw+fPizzz7bHUdGHG3R6oCAgEuXLvn7+3cksV6vZzKZj3ouOquDNfTqa6TGDF7rjc4IRXkdvsgxU+C3W9Rn+ZSZAh/3w8YFdPTnplarBQJBJ7MLdRksf9tIS0sbM2ZMdHS0vTPi0q5du5afn79v3762Z6x4/sMWoTMoUtDLLpKVGvB+X3RyEIo+4nUjAwVPhKNPhKN/VlHLLpKf88BXg7FwNxe/+oSgf6Snpw8ZMsTeuXBpHA4nPz+/mw4OB0v0bGYKvJdDDd1PjA1A86biU4MfOQre6bEg9MZUPCMQHfgXse4GRTpWZwEEQVC3gIGwB8uX0/3+JC41UzlT8JfiUasMAcVR8FI8eukxPKuaGpFF1GlhMIQgyMlZJxDu2LEjIiLCy8vr6aef1uv1D0p24MCBvn377t+/3ypv6uI2F1MjsogX4tCssbg/z8rdmKEC5Nh4fEwA2ncvcbgWxkIIgpyZFQJhYWHhc889t23btuLi4qqqqjVr1tw3mVKpXL58uVwul0qlXX9TV6YjwNxT5IYC6kwmviCyu9r0KAJW9UF/HYk/fZZcdwPufgxBkNOyQjX6/fffT548efDgwWKx+M0339yyZct9R6K+8sory5Yt8/Dw6Po7urJqDZ22n8AQcHESHi3q9vEsQ32QC5Own8qohWdIE4yGEAQ5IysEwqKiosTERMvjxMTE5uZmuVx+V5oTJ07cunVrwYIFXX87V3ahmR7wF/FEGLp9GMax1YDfQB5ydiIuN4LMw4TabKM3hSAIshkr1KYymaxtMoebmxsAoLW11d3dvS2BVqt94YUX9u7di3RgPqhWqw0ICGj788cff5w8efKDErvUPMK9NejL2fg3/YgxfkaNxtbvvm0AePkqPmwftXuo2ZN9u8Wv1Wo78p1C3QSWv21QlJN0hjQ3NycnJ9fW1lr+3LRpU0NDw9tvv/2g9IcOHdqzZ893331nqww+BEmSmjvqvg6e/1wuF0Uf0uSzQiB0d3dXqVSWx5YHEonkzgSrV69OTU1VKpXZ2dlarbaysrKqqiooKOi+R+PxePn5+R2cUI9hmIsEws/zqc9vUEcysD4Spr3ysCUdrL5GZpxEj4/HLMNzaJrm8/n2yg8Ey982HlqN9hQURbW2tloe6/X69957Lzs7u530Y8eOfe211woKCuLi4mySwYfAMOzOE96K578VAmFERERBQYHlcX5+vkQiubM5CACgKOrGjRvPPPMMAKCysnLbtm0Yhr355ptdf2tXQAOw8gq5r5o+PwkLsPbo0Ee1OhnjM6hhWeSxDAwuTwpBdlFVVZWbmxsWFrZr166hQ4eOGDGiqKho//79ZrN50qRJlqCl1WqzsrIKCgo4HM7EiRPvjWS7d+9OSUnx8vICAJw4ccLLyys+Pn779u1Tp07VarVnzpyZOXMmgiBPPPHE119//eWXX9rhc9qQFa50FixYsGfPnpycHJ1O98EHHyxYsMDSXF2zZs3BgwcBAOvWrbv6/+Li4lavXg2jYAeRNFh8ljzdQJ/JxO0eBS1eTUCXxaPDs8gKNZxWAUF2kJeX9+yzzy5dulQsFqMompWVlZGRgSAIl8vNyMg4c+YMAODatWtnzpzp1asXTdMjR468evXqXQfJyspKT0+3PN6yZcuJEycAAMuWLZNKpaWlpatXr7b8V3p6elZWlu0+m51YoUWYkJDwySefTJo0Sa1Wjx8//p133rE8X1JSEhgYeFfi6OjouzpOoQcxUWDuSVJpoo+Oxzu+cKgNPBeLoggYeYDcPxyJhUtdQq7nWiu9/Appm/ca4IWsSbn77o9Kpdq7d69QKAQAREZGbtmyZcSIEQAAb2/vDz74YOjQoWlpaWlpaZbELBbr+++/79u3751HuHHjxvz58x/67pGRkZWVlSqVyjL+w1lZp35dvHjx4sWL73ryxx9/vDfl9u3brfKOTk9PgOnHCSaK/DUGZznePdB/x6AEBTJPMs9MpB2kqQpBNhMsQJb3ttHP0oN9nyejoqIsUVClUpWWlq5Zs+aDDz4AAKjV6ubmZgBAQ0PDkiVLCgsLGQyGXq+/t2tUrVZzudyHvrvlJpxarYaBELI1LQEmHSF8uci2oRjuqPfpn49D1Xpy9EH09ATci2Pv3ECQDbmzwCh/e17/sVi3N1ezjBb89ttvPT09Lc9Yhva8+eab0dHRlrH6X3311YEDB+46gqenZ9s8NzabbTAY2v7LYDCw2bfDr1QqRVHU6bvxHLWWdWFaAkw8TIQIkB+GOW4UtHg+ipgViow+SMiM9s4KBLkkNpudnp6+c+dO8f8zmUwAgKamptDQUARBTCbTzz//fO8LU1NT8/LyLI8TExOPHDlCEITlz3379vXp08fyODc3NzExsS0uOivHrmhdj8oMRh8gokXIpjSsK/tI2MzqZGykHzLpCKEj7J0VCHJJmzZtOnjwYHJy8pQpU+Lj4z/88EMAwJIlS1auXDllypS+ffvedzbajBkzLIMZAQDPPPOMUCiMjIzUaDSjRo3Kzs5uWynz0KFD06dPt9lnsRfYNepAlCYw7hCR6I58NRjrCUHwtnUDsMVnyceOElljcSa8soKgbjZu3Ljhw4e3/RkcHHzx4sXKysqmpqaQkBBvb28AQGZmZmFhYXl5eWRkpEAgsDQTvby8ampqLK9KT083Go15eXm9e/dmsVi7d+9ubm6OiIjYsWPHgAEDLGl0Ot2+ffvOnz9v609oc7DechQKExh9kOjvhXw9pCdFQQAAAsA3gzEOhjx9Bu5gCEHdjsFg3LUzO4IgISEhAwYMsERBCy8vr4EDB0okEiaTaRnzgqJo261EBEE2bNhw6dKlO9NjGObj49P2zJUrV1auXGmZa+jcYCB0CHIjGH2ASPNB1g/oYVHQAkfBryOwSg392iUbjSmHIKiLBg8efNdo/xdffNEyGNVi2LBhCxcutHm+7AAGQvuTGcGog8RwX2Rdf8ebJ9FhHBz8NRo/VEvDPZsgqId65513xGKxvXNhBzAQ2pnUCEYdIEb5IZ/05ChoIWaBQ+OwLwqoneUwFkLOw9EGgh08eHD06NHddPD9+/c/++yz7SQoKioaM2bMfffa67lgILQnuRGMO0iM8kc+6tfjo6BFAA85OA57+SJ5pM6pfieQyypT0dUaxzqZ4+PjX3rppe44Mk3Ty5cvX7p0aTtpoqOj2Wz2vn37uiMD9gIDod0068HwLGJsAPKxs0RBi1gR8ttIfN4pIlfmWNUHBD2qBh0Yc5D04TjWjXuCIHQ6HQDAZDJ99NFHdXV1q1evXrFiRVVVlU6n+/zzz1977bVr165ZEre2tm7atGnZsmVr1qwpKSlpO0h1dfV77723atWqsrKy9evXK5VKAMDp06c5HE5sbCwA4MyZM5YFSLdt21ZdXa1Wq9etW2d57bx5877++msbf+puBQOhfTTqQXoWMSUYWdvXqaKgxVAf5KtBWOZh0tEupSGo49RmMOEw8XQUKmLZOyv/Kz8//4svvgAAGI3GFStWPP300wEBAUqlcvTo0ZbNz93d3dPT01taWgAAFy5cqK2t7devH5PJHDJkSFlZGQCgpaWlf//+BoMhOjp6yZIlb7zxhkwmAwAcOHCgbWLG4cOHLevRfPHFF7du3VIoFG+99Zblv4YPH37y5Em9Xm+HD9894DxCO6jT0iMPkPMj0JV9nPZCZHoIWqcFGYfIcxNxsYPVIxD0UGYKTD9G9PdCVvVBD93zv6aqIsXvX9kmJ8zQeNGUZ9pJsHbt2r59+y5cuNDDw2PUqFGWgaBHjhw5derUjBkzJk6cOHHiRJqmFQpFbW3tzz///NZbb23evHnYsGH/+c9/AADJycnx8fGWQ+Xn50+aNOmhWfL09ORwOGVlZQkJCdb4iPYHA6GtVarpUQfJJTHoqwlOGwUtXoxHa7T05KPE4Qyc7YTtXshp0QAsPkuyceTLQfc/cXHvXqKZL9gmMyiH134CSxhDUdTDw6MtpHl5eUmlUgBAQUHBokWLZDKZQCBobGzMzMwEAJSUlCQlJVlSxsbGtq1cqtPpOJwOLRzM5XLVanWnPpAjgoHQpkqV9OiD5Gu90edinTwKWnzcD3viFDnvFPnriJ6xYhwEAQDevEoWK+nj4/EHzepF2VxmYIRtM/VAGPZPtLasuA0AQBDEMrDzxRdffOqppyz7or/xxhuW6CgSidpW3NZqtZZ1ZwAAnp6elj5SAACHw1EoFG1H1uv1bTGSoiiZTHbn5P2eziWqYwdxQ0anHyDfTnaVKAgAQBGwbSgmNdDLLsKJ9lDP8E0h9XsFvW8MznWKZkJbxFIoFL/99pvlyQkTJuzcubO+vh4A0DYEBgDQv3//3Nxcy+Pk5ORjx461Nft+//33lJQUy+OSkhKBQBAaGmqzT9HdXKVGtrtLzfTog8Rn/dGFka5V5iwM/DEaP9lAfwon2kMO788qak0OdXAcdt9dAHuiV155ZdGiRRkZGWlpacnJyZYnR40a9cwzzyQlJQUHB1MUxeFwLHsTzpw5s20bivHjx0+aNCkqKqq4uHju3Ll79+7dsGGD5eVZWVkzZsxAECfq5KEdjL+/f21tbQcT63Q6giC6NT9WcaSW8tphOlBN2TsjVqZSqTqYslZDBf1s3lFKdmt+XE3Hyx/qiL8bKa8dpuyWu3+nQ4YMOXv2rF2ydF8EQej1epqmLV2Ubc8rFAqz2Wx5rNFoDAaD5bFUKs3NzdXr9Xq9XqvVtqWnKIqiqIKCAqFQSFG3P/WMGTN2797dlkan08XExPzxxx9tz5AkmZCQUFRU1G2f7/727duXmZl55zNWPP+dovHv2P6opJ79m9w9Eh/i40QXUI/In4cczcCGZ5EebGRsgOuWA+SwylT0jOPkD8PwZA9HPz8xDLPcF0QQ5M4V0e5cJpTH+2eIjbu7u7u7+10HWbx4cXx8vNls/uabb1atWtXWvPvwww/vnCzP4XDYbLZIJGp7pqysbMmSJVFRUVb9THYGA2H3+m8+te4GdXQ8Hi929F9Xd4sQIr+NxKYeI7LG4n0dvq6BXEqtlh51gPyoH+o6V2lPPPFEdnY2giA//PDDoEGD2p4PDQ198cUX70z51FNP9erVq+3PyMjIyMhI22XUJmAg7C40AO9eI3fdos9OxIL4rvLrat9gb2RLGj7xMHEqE48SwjKBHILSBDIPk8/HoXPDXej+/fDhw+/c1LAdL7xgo4kiduRCX7wt6Qkw6zh5sp4+NxGHUfBOmb2Q91OxjENkvQ4uOgPZn44AEw4To/2RV5x9Xi/UDvjdW1+jHgzPIlgYOJIBF1W5j4WR6DPR6NiDpMxo76xArs1MgRnHiXA35OOev/cL1BWwa9TKcqT0lKPkwij0rSQ4g/yBlieiUiOdeZg4Oh7nwXMQsgeKBk+dJhkosjmt5+2GTVHUd999l5OTExgYuHTp0jsHs0CdAFuE1vTbLWrcIeLT/ujbMAo+zEf9sBgRMvUoYYRT7SGbowFYep6s19G/jMDwHlgLPvPMM/v378/IyKiqqkpPTydJ+CvqEng1bh0EBVZdJXdV0Ecz8N7uMAg+HALAd2nY7BPk7JPkbz2zMoJ6rjeukNmt9LHxnVwF92Zryfor31g7U/fX2yt2acqiO5+pra3dtWtXTU2NQCB47LHHYmNjDx48aFlEFOocGAitoEkPZp8gWBi4MhmXwJuCHYYhYEc6NvkosfAMuW0YXIwUspH3r1MHauhTE3ABo5NHCBIGvNK/vW3crYjPuHvR7eLi4rCwMIFAAABAECQ5Obm4uBgGwq6AgbCrjtfTT54mF0ehbyWhsCp/VEwU7B6Jjz9MLPmb/HZIz7tVA/U4n+dTP5RSpzNx9y5cs/IY3Cj3cOtl6tHodDqVStX2p0qlcnNzs1dmnAPskOo8MwVWXSWfPE3+MAx7JxlGwU7i4GDfGPymnH7xArzPAXWvjTepLwuo4+Mxnw7tNeS4ysrKTp06BQCoqKg4ffp0enq6vXPUs8FA2ElFCnrQX0SulM6Zgo/wgzGwS/gMkDUWv9RMvww3qYC6zXdF1Md51PHxWACvx/9g+/Tps2bNmkGDBqWmpr7//vvh4XZrnjoH2DX6yEgafFFAfXCdXNsX+1c0vJKwDiETHM7AxxwkXrlEroOTuiBr21RE/ec6dWICFizo8VEQAMDj8Y4fP15ZWenp6XnnsqJQ58B6/NEUyOkh+4h9VdSFSTiMgtYlYoIjGfiZBvrliyRcdQayou+KqLXXqePjsVCniIJtgoODYRS0CliVd5SOACuvkOlZxIJI9PgEPMzNqX5RDkLEBEfH4+eb6efPw1gIWceGAuqDXOrkBMxpfrMxMTFLliyxdy6cCgyEHfJ7BRW3m6jWgrxpjH9Fw2Ex3cjSLsyV0f86S1IwGEJd8+kN6osC6tQEp2oLhoeHz507t9Mvp2naYDBYJSeWLQ+tcij7goHwIa600MP2E+9fp7YOxXYM7/GDzXoENwY4NA6/pabnniLNcFt7qLPezia/L6ZOTXDC7V/y8/Nzc3PbT3PlypWSkpJ7n79+/Xp0dLRVshEbG5uTk9PplxMEsWvXLoqy/48cBsIHuqmgpx0jpx4j50eg2ZPx4b7O9ltyZDwcZI3FtQSYeowwwJGk0COiAXjxAplVQ5/OxP17/hjRe/3yyy8//vhj+2k2btx45xa7Dkiv18+cOdMR1oeDo0bvI09G/+c6daqBeq03tmM4xoGFZA9sDPw+Eltwhhx7kPhzDC5i2jtDUA9hpsDCM2Slhj4xHhc642lTWFh44sQJgiBWrFgRFha2ePFigiC+/fbbq1ev+vv7L1261MfH59KlS9nZ2dXV1S0tLUlJSbNmzbrvoUwm09dff52TkxMUFLR06VJPT0/L86dOndq7d69MJktOTn7xxRdJkty2bdvly5cJgkhLS3vyySdR9IGNqHXr1qWnp+/atau5uXn69OkZGRmW50tKSjZv3iyVSocPHz537lwEQTZu3AgAWLVqFYqizz33XGBgoLWLqqNgi/B/HKujxx8mMg6RqZ5I+SzGKwkojIJ2xEDBD8OwJA9k2H4C7l8IdYSWAJOOECozODLOOaMgAEAoFPr6+np7e6ekpFg2i587d+6+fftmzJhBkmRKSopCofDw8JBIJAEBASkpKSEhIQ861LRp006cODFz5ky1Wp2amqrVagEA33///fz585OSkubMmdPa2kpRlMFgKCkpmTRp0tSpU7ds2bJ27dp2srdp06Y5c+bExcWNHTt20aJFf/31FwCgvLx84MCBXl5ejz322Pr169944w0AgKWHNjk5OSUlxb7DX2E1DwAAChP4sZT6tohCAFiWgP4xCmXBmWyOAUXA+gHYh7nUkH1k1lgsRuSE3VyQtTTrwcQjRII78s3g7l3GXVWhK9lZ241vcAdRBC98pv+dz/j5+cXExBgMhhkzZgAAysvL//rrr7q6OrFYPH78+KtXr27btu2ll14KDg6Oj4+3pLmvGzdunDlzpq6ujs/njx8//uLFizt37ly8ePGbb765ffv20aNHAwDGjRsHAODz+R9//LFer29qalq+fPmKFSvefvvtdvK8ZMkSy3AerVb76aefTpo0acOGDdOmTXv11VcBAOHh4UlJSe+8886IESMAANOmTWMwOrvqq5W4dCA0U+BoHb2jjDpYQ40LRL8chA3zRWBF64BWJKIBPJCeRfw6Ah8Gb9ZC91OspCccJudH2GIrUH4AO35JcDe/yW0o8yEhvaysLDg4WCwWW/5MTk6+7xiZe5WWlkZGRvL5/DtfqFAoGhoaBgwYcGdKrVY7e/bsoqKioKAgvV7f0NDQ/pETExMtD/r06WPJTGlpaduy4LGxsRiGVVVV+fv7P/AQtuWKgVBLgGN11N4qel8VFSVC5oajXw5idGUFXsgG5oaj/jxk1gnik/7YvHDYpQ/9j1MN9OwTxIf9sCcjbHFuoAyULbFzx2vbvAWRSKRQKNqel8vlbbf62icSieRy+Z0v9PPz4/P5TCZTKpVadrew2Lp1K4ZhlpB25cqVMWPGtH/ktvwoFApLhBaLxW1P6nQ6g8HQFrkdgatUKEYSnGuk379OjT5I+P5k/vImlSRBrk/F/56I/zsGhVGwR0j3RU5OwN+9Rq28AqcYQv/YUkzNPkHsHIHbJgo6Ak9Pz+rqasvj3r17oyi6a9cuAEBNTc2ePXss41PuTHNfqampKpUqKysLAHDr1q2srKxx48bhOD5p0qS1a9daBnO2tLTQNK3X6xEEAQCQJLl+/fqHZu+7774jCIKiqG+++cbSuZqRkfHDDz9YYuGGDRuSk5O9vb35fD6Hw2k/k7bhnC1CkgZ1WrpMBW4q6HwZnd1K31TQcWJkqA/yYhyWPhrhOefndn4xIuTSY/jUo8S0Y+QPw7FO7ycHOQeCAq9dJg/U0Gcy8QihC/WZz5o1a/fu3cHBwUOGDNmxY8fPP//81FNPvfvuu83NzStWrBgyZAgAYOHChfPnzw8ODp46depnn33W9loURdlsNgBAIBDs3Lnz6aefFggELS0ta9euTU5OBgB89dVXCxYssHS3GgyGwsLCJ598cvv27TExMRRFzZgxQyQSWQ7FZrPvO3w0IiIiPj7eZDIFBwd/8803AIDZs2dfvXo1KipKLBazWKydO3cCABAEWbNmzahRo1AU3bNnT1uHqu0hjrYuQEBAwKVLlzrYd1yvNPy3GEcQRGsGKjNQmUCzgW7QgXod7clGwt1AjAiJFyNJHkiSBOncVtRQO9Rq9Z39JzZjosCLF8gzDfTe0ZhLVX93sVf5OwipEcw6TjAxsDO9e2fXpKWlffDBB5bo4shaWlrc3d0x7JFruubmZolEctcLLbseent7I/8/cKK+vt7d3d0SRNsRHR29devWpKQkjUbj4eFx53+ZzWaVSiWRSB41hwCA/fv3v/jhN4PW/unNAZ/2x4BVz/+e3TLCEFrEACiK+HOBGxMImcCLjfpygR8XgcM+nRgTBV8Pxr4rotL2E5vSsIm9XKVDDGqT3UrPOE7ODEXe74vBDZ0tOnhr8F5eXl73Psnlcrlc7p3P+Pn5dfyYbDb73pDJYDA6FwUtfLhgTADSHbdFenYgdGOA13sjGAbrQVf0r2g00R2ZeYK82Ey/lwJrQxeyqYh6M5v8ejA2NRj+9h3Oe++9Fxoa2h1Hdmch3TRQDp5GUA/W3wu5Ohm/0kKPOgBn3LsEjRnMO0VuuEmdzcRhFHRMM2fO9Pb2tncuHg08k6CezZMNDo3DR/qhffcSh2phLHRm16V0370EGwOXJuGRLnxvGLI6GAihHg9FwJtJ6K8j8CXnyGUXSbhIt/OhAfjsBjX2ELE6Gd2UBpf/hawMBkLISaT5IDlT8Hod6PcnkSuDTUPnUaOlxx4k/qiiLk3CHw+DVRZkffCsgpyHmAV+HYG91hsdc5B4/zpF2H+bM6irtpdSffcSw33RUxPwYCfaXBdyKLCLAXI288LRdF9k0VlybyW1ZSjW2x3Wnj1SrZZeco6s04GjGbh9v8SUlJRZs2bdOxmApmkELk5sKzqdbvDgwd108J49oV6v1zOZzE5MIIWswpEndNMAbCuhVlwh/xWNruqDOeVyCo5c/l1B0eDbIuqdbPKFOGx5Isqwd7/Vg5aZ1mq19t08yNX4+PjcObURTqiHoIdAAFgQiY4LQF+4QCbsJr4ejI3yhxfvPUCujF5yjsRRcCoTj3WMXbc4HM59J8Y564WIC4KBEHJmvlywaySWVUP/6xyZ6oms648G8ByiboXupTCB1dfIX8qptX2xp6O6fSslCGpj704HCOp+EwKRgml4rAhJ+oN4/zqlJ+ydIeh/UTTYXEzF/m42EKBgOmMRjIKQbcFACLkEDg7eSUavPIbnyeiY34md5RTcyMlBHK+nU/YSO8qo/WPwb4ZgErgnGmRzsGsUciHBAuTXEdi5RvrVy+S6G9SHqdhoeOPQfq5L6RVXyHIV+DAVnRYCL8ohu4GBEHI5Q3yQC5Pw3RXU0vOkHxes7YsN9obh0KYKFfTqa9TZRurNPtjiaPuPC4VcHDwBIVeEADA9BC2Yhs+PQOeeIscdIv5ugl2ltlCooOedIodnEckSpHQm49lYGAUh+4PnIOS6cBQsiESLZ+DTQ9D5p8gRWcSROowaU4wAACAASURBVBgOu0t2Kz3zOJmeRcSKkbKZjOWJKA92SEGOAZ6JkKtjomBRFPpUBPrzLerVSySOgFd7ozNCYEvFOmgADtfS626QJUqwLB7dOowB4x/kaOApCUEAAICjYF44OjccPVBDf36DfP0y9VwsuigK9bx7XS2oozRmsKOM+qKAYmNgWQL6eCi8toAcFAyEEPQPBIAJgciEQDxPRm8ooKJ2mScEos9Eo0N84GiaR5Avp78ronaWUcP90I2DseG+sPQghwYDIQTdR293ZFMa9nE/bFsp9a9zJABgYRQ6Nxz14dg7Zw5MaQK/3aK2llA1WrAwErk+FYfr+EA9AgyEEPRAYhZYFo8ui0fPNdJbS6jY380DvZA5YejkYDjQ4x9mChyupXeWUwdrqJH+6Mo+WEYggsEICPUc8NcMQQ83xAcZ4oNtILA/Kqmd5dTS8+TYAHRGCJIRiHJd9TdkpsCpBnpXBbW3kooUInPC0A2DGHBdGKgnctUfMQQ9Oi4OnghHnwhHpUbwRyW1qZhaeIZM90MnBSHjA12l11RpAkfqqD+r6IM1VJQQmRGKvjkZ78WHDUCoB4OBEIIemYQFFkWhi6JQhQkcqKH+rKJfvWQOEyBjA5DR/uggb8TJhkdSNMiR0kfq6CO1VHYrneaDTApCP+nH8OU+/LUQ5PisGQiNRiOL1V7PiMlkYjKZVnxHCLIvERPMCUPnhAGCwv5uog/XUq9eIkuU9EBvZLgvOsQbSfVEWD1zT2CSBtel9LlG+nQjfbqB8uYgo/2RV3tj6b6Iy/YGQ87KOmd0bW3t7Nmz8/LymEzmZ599Nm/evLsSvPPOO1u3bm1paeFyuc8///zq1aut8r4Q5CBwFAzzRYb5Yv9JBXIjONNInWqgl12kChV0ogTp74mkeiIpHki4G+LIOwxVa+jsVvpqK32xmb7aQgfykSHeyIwQZONghot0/EKuyTqB8KWXXoqPjz99+nR2dvaIESNGjBjh7+9/ZwJPT8+jR49GRUUVFRUNHTo0Li5uxowZVnlrCHI0YhZ4LAh9LAgAALQEuNpCX2qhd1fSq65SUgOd4I4kuCMJYiRahMSIEDv2LkqNoFhB31TQBXL6hoy+LqUZKEjxQFI90VcT0AFeiBiOfIFcA0LTXV1cUaFQeHp6lpaWBgcHAwAyMzOHDh36+uuvPyj99OnT4+Li3n333fv+b0BAwKVLl+6Kow+i1+uZTCaG9cy+p55PrVYLBAJ756InkRtBnozOk9EFcrpISd+U03oSRLghIQIkRAB68ZFAHvDlIv484MV+eJ9qR8rfTIEWA92oA3U6ukYDarR0pQZUqOlSJU3RIEqExImQaBGSKEESxPaMyj0RPP/ty4rlb4UWYVVVFZPJtERBAEBsbGx5efmDEstksnPnzi1ZsuRBCWiaViqVXO7tXySPx4O3FSGnIWZZelD/6R5VmkCpiq5Q05VqUKaiT9aDBj1VrwXNBpqBAgkL8WADIRMImQgHA3wGYGGg7RadyYQzmaTlsZEEOgLoCaAngcJEK4xAYQItBlpLAC824sMFflwQwEMCeMjEXiBEgIa7IXD1OAiysEIgVCgUPB6v7U+BQFBRUXHflGazef78+WPHjh01atSDjqbT6QYOHIiit0fdffvtt+PHj39QYtgitC+tVosgDnzLqyfAAIhmg2g2AJ53/5eaQKRGIDcCpRlRmYGeADoSMZJAT94uc5w2scDtx2IGzeUANkazMSBiADcGLWQCCQuImQ/o8iGARtNdH8pFwPPfvjpY/lwuty2gPIgVAqGHh4dSqWz7Uy6Xe3t735uMJEnLIJpNmza1czQej9fxrlEMw2AgtCOapvl8vr1z4bT4APi2m0CtNgsEsFlnN/D8ty8rlr8VpjsFBwfjOF5QUGD5MycnJyYm5q40FEUtWLBAJpP9/vvvsKsTgiAIchxWCIQ8Hm/OnDmrVq1qaGj4+eefc3JyZs+eDQDIy8vLzMy0pPn3v/99/Pjx559//ty5c8eOHSssLOz6+0IQBEFQ11ln+sS6deuWLVs2aNAgX1/fP//8093dHQBA0zRBEJYEWq02Njb2iy++sPw5fvz4e1uNEARBEGR7Vpg+YV1w+kQPAoeP2xcsf/uC5W9fVix/51oSEYIgCIIeEVw00GnJDYr8lsL8lqIqVU2Dpkmql5spwkAY3JgCIdvNi+sRJg6JEIek+PSRcMT2ziwEuRaNSXutKa9YWlquqKzXNCkNKpVJhaMMFsYUstx8+d6Bbn6xHlEJnrE+PC97Z9b5wUDobFp10qOVp8/WXKxW1cZ7Rsd5RE+KGOfL85Zw3Zkog42zlUaV0qhu1DaVyyv/rr28IXuzF9djaOCgcaEjvHn3zGWDIMh6VEb10crTJ6vO3lJU9faKjZFEZYaPDRD4CVkCIcvNTBEm0iQ3KBs1TZWqmnM1lzZe2ypkuQ0NHDAyaGiQMNDe2Xda8B6h88huzP2jJCu3qWBor4HDew1O8u6Now8vGYqmCqUlRytOn6g6GyUJnx0zNdmndwffEd4jsS9Y/vb1SOVfJq/4+eaei/VXB/qljgkdnuSVwMAYD30VRdNF0pLT1eePVp4KFPg/FpkxvNdgFIG3tACw6vkPA6EzOF93+cf833Rm/YzoSSODh3HwzkyyNpPmY1Vnfrm5h42zFyXOTfVNeuhLYEVsX7D87auD5V8iK9+cu6NcUTkjetLE8LE8RmdWdCUo8lztxd3F+2V6+dy46WNC0zHE1es9GAhvg4HwZmvxxmvf6wnDkwmPDwkYgHZ5wSca0GeqL2zO/dGT6/FcytNhouB2EsOK2L5g+dvXQ8u/Wdf6bc72nKa8JxMeHx82moFa4VZUbnPB9hu/tOikzyQ9OSSgf9cP2HPBQHibKwdCmV6+8drW3Ob8pxPnjglJ73oIvBNJk/vLjmzN+3l0yPCFvec8qIkJK2L7guVvX+2UP0mTu4r+2lmwe0rkhMdjp3Suk6Ydl+uvfZ2zVcQWLktd0sstwLoH7ylgILzNNQMhRdN7Sw5sv/FLZsSYeXEz2Xh37RqnNKq+vrb1WtON5QOeT/FJvDcBrIjtC5a/fT2o/MsVlR+cXy/miJalLvHj+3TTu1M09UfJgR9u/JoZPvrJhMeZmMstXQkD4W0uGAirVXUfX/wCRdBX+y/t5dahUuqi7Mbcjy5uGOCX8mzyAvb/XtjCiti+YPnb173lT9HUTwW/7y7etyRpwbjQETbIg0wv35C9uUx+6/UBLyR4utZyXTAQ3matQEgaKZPSbNaSpJEijSQAAGWgGBNliRksMQNBHWKnFRrQvxft25G/66nejz8WMd66faHt05p1G65uKmgtemvwq5HuYW3Pw4rYvmD529dd5d+kbVnz9zoWznxjwIseXIktc3K25uL6q9+ODEpblDjXcZqGRoXZIDVRZtpSqWJMFGWiDD7O5OM4zwqtFxgIb+tEIDRrCF2TUd9iMrQa9S0mg9RklJtpkmYKcQYfR5kozsYAAKSJokyUQWYyqwmuH1scyRfHCNxC7LaBd7Ou9YML680ksXLQS93X2dK+41VnN1z9bm7czGnRmQhAAKyI7Q2Wv33dWf5nay6su7xxVsyUWTFTbHmR2kZlVK+7vLFKVfvmoGXh4lDbZ8BC12CQ3VTLizXqaj3OQtkeTIyFYqzblSpppMwawqwmKJJmixlsCZPtweR4sjieTI4niyVigEcpORgIb2snEFIEbVKajXKzQW4ySM2GVpO+1WhoNgEMcDxZXC8Wx5PJ9mCxPZhsMQPnPjCUUgStrtIpSjTSPBWCIX5pEq++IgSz6Yl+qvrv9Ve+mRH92OzYqfadQtSgaXrn3Ec+PK/lA17gMbiwIrYvWP72ZSl/giK/zdl2pubCu2nLoyUR9s3SkYpTG69tmRM3fUb0JOSRokoX0UCar6o/KzW0mjwS3USRfLdQHsZ6YGVFGimjzKSXmgytJn2LSd9i1DcbCQPF8WByPJlsye2amSVmMEUMjHn/48BAeJtaqm08rkAQhKZo0kCRJorQk4SONKsJykwx3RiWvs3b1x0eTI4Hq/NNchooSjV1p1qNcnP4TH/btA4NhHFD9qbrTflvDX7F7r8xCzNp/vLalqsN19cOW+mBimFFbEcwENqXWq024cQ7Zz8SMHlvDHrJjekQ30WDpmnN3+t4DO7KQS+J2SIbvKO+2Vi2q540UQEjPCQJbp2+l0QaKX2L0dBq0rfe7qszKswmhRkAwBDgDB6GczGMhTEEeNhUXwADYRuNQqsq0KMoChAE56AoA8W5GIOLMQQ4zumuETTSPNWtvQ2S3m4hk3y69fZhhaJq9bmPoyQRy1KXWH34dRcdrTz1VfaWZ+KfzIgaZe+8uC4YCO3rctW1j699+VjEuLnxM2za/HoYkia35f1y4NaxlQNfuu94byuqPyOtOdbSa4yn72BJN5UBaaRMKoLQkYSeJI0kTQHPJCGAgbCNvUaNEgay+IcagCDR8wPbaf53xb6yw5uv73g2ecFYm4w964QyecWqU++PCE5b3GceXPPJLmAgtKO/Sg9tyf1p1aBl/fyS7Z2X+8tpuvH++c/HhqYv7D2nO5ahoSn61h8Nqgpd7NNBLPHDl4uzOhgIb7Pj9Amaosv3NKirdAn/DmnnFmMnaM26dZe+qlLVvjPkddtMkOi0emnDp9c34ij+9uBX+UyevbPjcmAgtAuCIr+4+l1ec8HKvi9F+oTbOzvtURiV75//3EAY3x78iifXw4pHpkm6cFs1TdLRT/bqpsbAQ8H9CO0PQZHw6X6iSP7N76spwmoXEyWy8sUHl/GZ/K/HfuLgURAAIGDyPxmxupeb/5LDr1ar6uydHQjqdgqj8uXjb0r1so1jP/Hleds7Ow8hYgk/Tn9ngF/Kvw69cqHuqtWOS4PS3+oAALGLguwVBa3LGT6DHYVk+rDdGUXbq2nKCrFwT3HW6ydXL06c93K/fzvOZKD2YQi2NGXRE7HTXjj6xqX6bHtnB4K6UZm84pmDryR6xa8ZupLL4Ng7Ox2CAOSJuOnvpa34/MrXX1/bSlBk149ZdbBJ12SMmhfoIHOsuw4Gwq5BQPhMf9JA1Rxr6cph1CbNm2c+OFRxfOPYT9KDhlgrdzaTETZq7dCVH1/68tfCvfbOCwR1izM1F1498faSpKeeTnzCLjMFuyLBM2ZzxvpqVd3zR5c3aJq6cqjWXFXLdWXcoqAHzWroiZznk9gLiiPRTwY2npepq/SdO8KNlsKnD7zky/P6aszH9pos33XxntFfj/3kWOXpDy6sN5Nme2cHgqyGBvT2G798mb354/TVPfE61cKNJfjP8FUjg4b9+/CrJ6vOde4gJhVRvqc+am4gg+9Um7rDQGgFDD4eNs2veEcNaaQe6YUUTW2/8cvbZz9clrrkuZSnrbJLix15cT02jP7QSJpeOLZSqpfbOzsQZAUGwrD67MeXG659M/bTO9cX7IkQgEyPnvhx+uoteT99dPELA2F4tNfToPTXOr/B7oJePaNbuONgILQOSYKbWwi3MusR+hwatc0vHF2Z21ywOWP9QP++3Zc3W2LjrHeGvDbIP3XJoVeKpKXd90YUQeuajIoSjeWfrsHwqFchUE9EEbS+2ago1Vq+d21d937vjdrm544s5zI460e+784Rd98b2VKke9imjM9pml58cFmxtKzjL2y4ICN0ZMAoz+7Lm7307CaIQwmd4pv9YanPQDHP9+GT3w/fOrHx2tbZcVNnRk/ucfcb2ocAZF78zDBx8IpT71l9DX59s7H1urL1hkrfZGS5M1nC24sTmlRmg8yMczBBEMctmCuK4nfkW4B6BIPMpCjSqCp06mqdUW5miZlMIW4ZpmFWE/pWE0uES+LdJIlC67ZUrjXmrTm/bm7cjGlRmVY8rCPg4OwVA188WXVu+al3p0dNmhM37aFTgQkdWX2oOeHZEKcZIHMnOI/QmurPSmU31fHPBLeTRmFUrru0sUZd/+agl8PFIbbKWrdofx5PlbJm1Zn/9PNNfi5lYdfn8xqkpqqDzcpSjUeS0KO3myCYe+8P0iAzqSv1qkqtvFBDk7Qkwc0jSegWxHWkRT+sybnnEeoaDC3XldI8lVlHiqP5wjCeIIjL9WLd/W3SQFtvaL2haslWsCXM4ExvfoAVwuGuor923tz9zuBX+3gnPCiNE5R/i671gwv/NRDGlYNeChD4tZPy1t4GiqDDp7eXxsbghPrbHC0Q0iR97ZOy0Md8xDH3/3pOV5//79Vvx4aMWNh7DgOzw1oM1vXQE1Fj0q49v05n1r+btrwryx7WnW6tOdbilybxH+bRwXlL+mZja66qJUdJGknvVLFXqogt6RkzUjrOCSrie5nVRPM1RfNlBWEgPfsIPRKF/EBORy5laIpuuiivPtLskSgMmeTT6ZXxDYTxk0tfVqtq1w59w5vn1U5K5yh/GtB/FB/YduPnJxNmTYnMvG8Hlb7FlPdFefLyCIcaIwMD4W2OFggBALICdWVWY9Kr4Xe1V+QGxfor31Yoq5cPeD7OI9pe2bOujpyIFE1vv/FLVvnR1UNej/d85A9OmqiyX+v0UlPMU71Yos5cOmjrDU2X5S3XlDxfts9AsSTBzcabh3Qf56iIb6OBokTTeEGmKNVK4t28+omEobxONOVJA1Wys9asJaKf7MV0e+Rau07d8NaZDyLcQ1/u9yzrYXN5nan8a9X1H13cQNPU6wNeuHcpj8Kt1YIgTsAIx7o7CAPhbQ4YCAEAeRtu+aVJPPoILX/SgD5QfmzT9R8mhI1+MuHxnjJTviM6fiJeqLv60cUvnoibPj16YseXJyaNVP43FVxvdth0PxTvUvSiSVp6Q9V4Qa5rNHj1E/sOcrfL6ojW5RwVMaElmy7LG87LcA7mM1DsmSzq6mIlNKg53tJ0UZ6wNOSRLp4sewou6D3nsYiMjqR3jvJvQ9H0n6UHtub9PC1q4py4aW2D2LX1hoJNVX1XRXbxN2h1MBDe5piBUHZTXX2ouc/LYQCACmX1Z5e/JijzK/2e6+l3BO/1SCdio7b57bMftu1l+ND0FEHf3FzFdmeEz/C34k0+fYux8bys+apCEMz1HSwRR/F77h3Enl4Rq6v0DX9LZQVqSbybj7UH5dedbm28IO/9fCijAzuvmSniu5ztj7qnYE8v//tq0rb89+q3teqGl1OXWO6PlvxUy/VlB4yw5lKlVgED4W2OGQgBDa59UuqX6b5bv/do5akFvedMDB/nZENDLR71RDST5o05Wy/WXX17yKsxksj2ktKg6IdqAJCoeQHdMUqNMlMtOcqGv2WEjvQZ6O7dT+RQNz86qIdWxKSJarmmbDwvIwyk7yB3735i6y5b36bqQJO8WJPwXEj7a6A0aJrePfeJhCNePvCFR9pTsIeWf0ecrbn4ZfbmOI+oxaHzq76Wpa6KwtgON9cOBsLbHDMQUjR1NOvv1mvKhrFVi/vME7Lc7J2j7tK5E/FszYXPrnwzM/qxWTFTHnR9UH9G2npdmfBcSHffz9NU6xsuyKR5KnEM32eAuzCsM/el7KXHVcTaekPjeVnLdaUwjOc7yF0U2e3N8ZKdtSgDDZ/xwLGOx6vObri6aW7cjGnRmY+6p2CPK/9HYiRNPxX8rj5sCvMIGjVvINvBtkQFMBC2ccBAeK720ubrP3qwJNP+npmwOIzv73BnjxV1+kRs1rW+//dnKIquHPjSvbvDaOsN+d9UJr4Uyna30f1UQk82X1U0XpDTJOXdX+zVV9yJcRa211MqYsJAtuYoGy/JzWrCu7/Yp7+YKbTRDVrSSOWsKwuZ6CNJuPt6VGvWrb/ybbGs7K3Br0SIQztx8J5S/p1GaMnL/yk+N+LkZVX2UwmzMsJGdce+hp0GA+FtDhUIsxtzt+T+ZCSNi/vMG+DXt/Zkq77ZGDHL0bdS6oqunIgUTe+8ufv3oj+fTX56TMjwf54n6OuflQeM9PBK6fx0i05TV+mbLslac1VuIVyvVJF7nJujDRC4k4NXxDRFK0u1TVcU8kK1KJLv3V9slzuy6ip94fdVfV4Jv/PiJqfpxocX/tvfL+XZ5IVsnNXJIzt2+Xdd3alWXYMxYrZ/kbT0u+s/NOtaFyTMTg9Kc5AbPTAQ3uYggfBKQ872G78ojeqnEh5vO0vMGiL7g9LUt6OcY7+u++r6iVgqv/X++c8DBX4v9/u3ZaJh1aFmfbMxen6glfLYGZSJas1TNWcrNDV6SYKbZ5JQGM5zwAU1HLQipoG6Wt96XdGSo2SJGV4pIs9kUTfdBeyg6sPNugZD9FO9AAAGwrgp94fT1Rde77+0i5vLO2j5W8+1j0rDZ/m7Bd8e2pbdmPt93k9qk3Z+/Mz0oCF2bx3CQHibfQMhSZOnq8//fHOPmSLmxk0fEZR21zJFRduqRdECnwFOskThvaxyIppJ87YbP2eVH3s2eeEw8eDr68r6vBLeuSmDVmdSES05itbrSoPULElw80h0E4bxHGcaomNVxDRQV+ta81TSXBXKQDz6CD2TRRxPh5gsRBH0tQ9LI2b7l/PKP7n0ZYJX7NKUpx9pXMx9OVb5W5uqQlf2W13y8rvH0F5tuP5j/m/NutZZMZPHhY7sdHu662AgvM1egVBlVGeVH/2jJMuX5z0rdspA/773vc0uL1RXH2lOfLFnr1jfDiueiCWy8o8ubhhVOKZ3ZHTcxM7csOlWBplJmqtqzVPpW43iKIEkXiCK4uMc57ki7jTKRCnKtLICtaxAhXMxj95CSW83np/D3RqvvtJ481DF1oQtL/f7d3+/FKsc0xHKv/uU/lzH9WP5D7v/rImbrcU7b+650XxzfNioyZEZ7S/B001gILzNxoGQBnRuc8H+0iMX6q+kBQ6cFpXZ/j12mqKvri2JWxzEddIFoK1bESirtNc3l36R8N/JsRmPx0xxzCXoTCpCdlMtK1Apy7U8X7Yoii+K5At6cezScWq3ipgG2gaDolSjKNaoKnT8QI57nEAS58b2cIj2311oQB8qP/5dzo9LSpbEjQgLGGC1KtuJAyFpoK6sKU554yFrqjVomvaUZB26dTxWEjUpYlx/vxQctd3VIQyEt9ksENaq649Vnj586yQLZ2WGjRkbmi5g8jvywqpDzaSBDJ3s2905tAvrVgQ3NlZ4pYroWPOX2ZsrFNXP9100wM9xd6eiCFp1Sysv1ihKNEap2S2U6xbGcwvhCgI5Nus7tWVFTFO0rtGouqVVluuU5Vqcg4oi+JbrAEe+C14sLVt/9VsEIC+lPuOr9ivaXt13VaS1viAnDoSNF2SKEk30k706kthImk5WndtfdqRWXT86eNiokGFR7uHdnUMAA2Gb7g6E1aq6c7UXT1adk+plw3sNGRua/qhfsEFmyv3vrX7vRDngUIuus+KJqK7SF/9Yk7IywlJQlxuubbi6yYfn/WzyghBRkFXeovuYtaSqXKss16oqdPoWI8+PLejFFQRx+AEctoTZfeMku7siNirM2jqDulqnrtKrq3UsIUMQzBWG8YThPAe5iduOFl3rptwdVxty/tXnybGh6ZabFwXfVXr0EXr3s85teycOhHkbbgWO8nzQ5gEPUqduOFxx4ljlGQQgw4MGDw0YGCkJe9TZmR0HA+Ft3REITaTpRkvhpfrsC3VX9YRhSED/Yb0GJXrFPXS/rgfJXV8eNMFbFNGhFmTPYsUTsXBbtSic5ztE0vYMQZF/lh74Mf+3wQH9FyTM9uBK2nm54yCNlKZGr67Sqav1mlo9qad4/myuL4vny+b6sLneLCuOn7RuRUwaKX2TUdtk0DUYtfUGbb0BAMAP5AgCOYIgjqAXF+/AWmWOQGvW/Xxzz5+lBydFjHsidjqX8c/KbYpS7a099cmvR1ilcnbWQGhSmnM+Leu3OrrTTediWdnp6vNnay4aCEN/v5T+fikpPol3fhFWAQPhbdYKhAbCWCQtyW2+mducXygtCRUFWb68SHcrXM7UnWzVS00OtY+XtVjrRNS3GPM2VKS+GYnesxSWxqTdeXP3vrLDE8JGPx47RcQSdv3tbInQkpp6va7BqGs0aBuN+iYjQAHHk8XxYLIlTLaEyRIzWCIGU8joxITFzpU/TdEmFWGUm41ys0FmMkhNBqlJ32wkDBTXi8X1ZlnCNs+PbbNp79ZiIAx/lBz4tfCPgf6pC3rP8bpnrQYAwPXPy3uN8XKPs8J566yBsP6MVNtgsMoc6BpV3YX6q5fqs2+2FgcLeyV7JyR4xcZ7xPCZvK4fHAbC2zodCA2E4ZaiulxRUSIrL5KWVqtqw0QhCZ4xid7xfbzirXvl4sS9o9Y6Ect+q2MKGb3GPnAUQ6tetiN/1/GqMxPDx86IfkzM7mHh8E5mNaFvMeqlJoPUbJSZDDKzSWE2qcwYG2PwcaYbzuDjDB6GczGci+FsDGOhGBvF2BiCgtvjVJHbDzQaDZ/PBwCQBoqmaAAAYaAARRN6kjRRpJEi9RShJ81agtCRZi1pUpnNasKsJRl8nC1msMQMlpjJljDYHkyOB4slYvSg5eXuYiAMf5Ye+qXwj0SvuAUJs4OED5yH2pqrrD8j7f28FUYmO2sgzNtwK3C0lzjamp1YJtJU0Fqc03TjRsvNImmpJ9cjRhIR6R4WJg4JF4V0Li5asfx7wDpSXaQ0qlp0rc261np1U52moVZdX62sVRiVQW6BYeLgCHFYRuioCHFI941RZLsz2WKGslwnirDCRZDzMWuI1jxVyhvtrcHtwXF/KfWZOXHTfir4ff6+Z0eHDJ8VM9mb51i7o3UQQ4AzBLhb6N0ng1lDmNWESU2Y1YRZRxJa0tBqIvQkaaRII0UaSJoChJ4EAAD69gOaphEEAQBgLNTSi4WzUYAgOBdDGQjGxnA2inMwlpjBD+AweBjTjcEQ4EwB3nMD3r1UJvUfxQf+KNnfxzvhsxHvPfSOsiTBrXJ/k7pKJwh6+BYoLsikNOubjVaviciQxgAAIABJREFUrJgYM8k7Ick7AQBA0VSFsrpIWlosLTtRdbZCUc3B2UHCQH+Brx/fx4/v483z9OJ6iNgimy1h07NbhFK17GDVcRRFdWY9SZEG0mAgjDqzXm3SKI1qpVGlMCg4DI4nR+LN8/The/nxfQMEfkHCAB+ety1XCao92Wpwxt5Rq1yR1Z1q1TUaIx7vaD+MTC//rejPrPKjqb5Js6InR0lsMT7NMTlri6SD6tQNvxf/dazizJDA/rNjp927neyDWGv5Q6csfyv2i3Zcs661Wllbq66v0zQ2aJqatS3Nula1SS1kuQlZbm4sNzcmn8vgsDAWj8Fl46wnEx4HsEXYhqIpjUmLIAiXwcEYmAcmYeMsHoPLZ/CEbDcRy03EFrVtL2lHHoluuetvhU31db7e0a5ruiQPn/kIvzp3jnhJ0lPz4mdmlR15++yHEo77tKjMob0GOcIXDdkARdNXG3P2FO8vkpZOCB+zPfNLd86jjQL17ivK/rA0dLKvI0/8sJfWXGXgaFvPjvfienhxPfr69rnzSZImFQaVwqhUGdVKo0pPGIyEUUfo2Zj1p2X37LqDz+AtTpxn97VGH4rtzmSJGaoKnTAM9o7+D3WVjqZB22KGHcdjcGfGTJ4ePenv2st/lGRtyN48LnTEhLDRgR1uFkA9jlQvP3jreFbZEQGTPzly/HtpK5hYZ6bwMwS4MJzXel3p3d9plz/sHJOa0DVZv1+0czAEk3DEkke8yumcnh0IexD3WIG8UA0D4V2aLsm9+4s7fcsKRdC0wAFpgQPq1A37yg6/cGylP99nXOjI4b0GW2VYGuQITKTpfN2VQ7eO57cUDe81ePWQ17veH+7dX1x7vAUGwrvIb6pFkXzHWU3XZmAgtBFxtKD017rgTHvnw5GQJqo1T5X8+t2r+naCv8B3SdJTi/vMu1SffejWiY3Xvk/2SRwVPHSAX187LgoMdQVBkdeack9Unj1XeynSPWxc6IjVQ5Zb69sUR/PLdtXrmoxcb3h6/ENepHaPddqNxNsBA6GNCHpxzBrCKDezxD1sblb3ac1VCkN5VtwCF0OwQf79Bvn305p1Z6rPZ5Uf/fjihlTfpLTAAQP8+sI2Yo9gJE3ZjdfP1lz8u/ZyoJvf8F5DFveZb/X+MQRFvPqKmi/Lgyf6WPfIPRdN0ooSbdhUZxvT1xEwENoKAsTRfHmh2meQu72z4iharim7aY8qHoObETYqI2yUyqj+u/bSiaqzn13+OkoSPtA/dYBf344PL4RspkXXalnOKafpRpQkfEhA/wfNiLcWr2Rhwaaq4EwfZ5pM0hWqCh3bk8kQuGJQcMXPbC/iaEFLjgIGQgtCT2qq9OIFHVrVt9PcWAJLRDQQxuzG3At1V34v+gsBSKpvUopPYpJPQo9bqsaZ6Mz63OaCa425Vxqvy/TyVN+kEUFpywe+0PWdAjuC68vGWKi6Wi8IsvLSXz2UvEjt/oiLizoNGAhtRxzDL/+9niLoTiym5XykeSpRFA+7Z021bsLGWYMD+g0O6AcAqFLWXGm4fqTi1KeXv/LkeiR5x/f2jEvwjOkpy5n2aCqjOr+1MK/5Zm5zfoWiOtYjKtm794oBL0S6h9tyaq+FJFHYmqeEgdBCdlNj4+mDjgMGQtvBORjXl6Uq14qinHAB7kfVmqfyThXZ5a2DhIFBwsDp0RMpmiqV3brenH+08tT6K9+ycFa8R1SMR1S0JCJCHApH2VgFQZHliooiaWlha8lNaUmrThrjEdnbM3ZJ0oJYSaR9d5306O1W+H11COwdBcCoMJs1hKCXi14TwEBoU+IYgaxQDQMhoSfVFbro+Q9cENI2UASNkoRHScJnxUwGANSo6gqlJTdbS05UnbmlqPble0eIQ8LFoeHikFBRcI9e4NSWNCbtLUVluaKyTF5RKrtVpar15/tEScLjPKOnR08KFQV1eiMXq+P5sREM0dTq+YEuGgDayAvV4ii+y14QwEBoU6IIftmuOnvnwv6k+SphBM/R1vUIdPMPdPMfE5IOACAoskJZVSa7VSqvuFB3pVxeiSJoiKhXkDAwyC2wl5t/gJufF9fT9r15jqZVL6tV19eo6qqVtRXK6ipljcasDREGhYmDI93DJoSNDhUFO3Lb2iPRrTVXBQOholTrsjcIAQyENsYPZBvlZrOGYPBduuSleSrPPg7dwMJRLEIcGiEOzfj/Z6R6eaWyukpZW6WqOVd7sUZVpzCqfPnefnxvP76vD9/Ll+flzfPy4nk45QAcjUnbrGtt1DY3apobtU31mqYGTWOduoGNswIE/kHCgEA3/xTfxGBhL2+eZ/ftxWp1kkRh0bbq4Exve2fErmigLNOGuPBMEpeujm0PQRG3EK6yXOuR6IR1ZQdRZkpZro18IsDeGXk0ltWeUnwS254xkqZ6dUO9pqle09iobcptym/StjTpWgyE0ZMr8eC4e3I93Nkid47YnS0SsYUillDEFgpZAjZu/cUSu8hEmlRGtcKokhsUCqNSblDK9HKpXi7Vy1r1smZtC4qgXlwPH76XN8/Lh+cV6xFl2SuAx+jZezjw/dmABi4+s17XaMBYqCtPcYaB0NaEETxlmUsHQmWZlh/AwdmOvkLsQ7EwZogo6N59f4ykqUXXKtXLW3StMoNCqpNVKKvlBoVcr1AaVSqTmqJpAZMnYPIFTD6PweMxOFwGl8/ksXEWC2PxmTwcxTk4m4kxWBgLRdC2DTK5OAdD/yk3rU6rRrR3vrXapLE80BMGkiLNlNlAGE2kyUiatGadiTTpzQaNWas3G7Rmrcas05g0apNWbVKTFCVkCYQsNxFb6M4Wi9lCd444VBQk4bhLOO5eXA+rby/uOETRfHmh2pUDoaJM6yDri9oLDIS2JgznN16osXcu7ElWqBFHO/PdCBbGDBD4BQgeuEKHkTSpTRq1SaMxaTQmnc6s0xF6jUmrJwxKo6pO3XBnAKNoSmfWW16oI/QkRbYdh6IoFP2f+6wC5u1xWGychaM4A2WwcRYTY7IwJpfBYWJMAYvvw/fi4Bw+k8tj8Pj/H485jtdItRn3GEH9Gan/8G6cvO/gXPzSHMBAaHt8PzahJUxKM1Pooh0R8iJ1TDfPo3dwLIzJ4rh7cLq6tIJT7odne8IIXvFPNaSRcrTRWzZCA9UtbZjT7Zb6SFzyi7cvBLiF8ZTl2oendEb6ZiNN0Dwf121/QI4GY6KCXlxFqcbeGbEPTZ2ewceZLrmyWhsYCO1AFM5TlrloIJQXacQxgp4zqBByCeIYvrzQRQOhslQrjHD1mc0wENqBMOL/2rvzKCmqO17g1Xst3T0LMz3d07Mw7ILKsCriEwURPErEuOMWMU8MkqMSc46GhCRmUYMvB02eGgngMXmJipETwfeiMYqSAKIIAhEYlmHW7tlnurqqequu90ebyThrz3R11/b9/NXT1HRf2675Vt37u/c6u08bNAg7T7AF04x+1oHaFF7g6jrJKt0KZfSc5fImGbpShkAQKoL2OMRoMtodV7ohuSbGkmwdnz8FQQjqQnkcJrOJD0aVbkjOSUToPJ83QdtzYDKHIFSCiXBVUmwdr3Q7ci10lnOWUQYtSQB1K5jq7DpluJtCPhixMRaDr+9BIAiV4q5iQrWGC8KeM1y+4TthQJ1SE3yVbkWuhWp5dxVOSQShQtxVtAGDsPsMhuVBpfImMqFzvJSUlG5IToVqedd4o/eLEghCpTjLKaE1KkaTSjckdxIRUWiNYnVjUCeb02rPs3GNEaUbklOh87y7CkGIIFSI2WpifGS4QVC6IbkTOsu7KmlsSgyqlT+Z6TZS72icTSR4kfYYd225XrIFIc/z9fX1yeSQtziiKNbV1UUixrrgGobLYL2j3afDBl/PEFQubxLTc8ZAswl7anl3FY1JvYRcQfj888/7/f7FixdPmzbtxIkTAw/47LPPJkyYsHTp0tLS0ldffVWWN9U693g6VGugy8+eM5iuBKqWN5EJnecl0SjDhGwth37RFBmCsK6u7gc/+MG+ffvOnDmzatWqhx9+eOAxa9as+d73vldTU/Pee++tXbu2o6Mj8/fVOvcEhj0vGGRwPsGL0c64swwDhKBeVtpCjbMbZ8AClTK9ZAjC11577aqrrrrgggsIgli7du0//vGPlpaWvgecOHHiyy+//Pa3v00QxNy5c2fOnLlz587M31frbIzF5rLwLYaYw9tzhnNV0SYLemFA1fImOQ0yTJiMJ/lg1IXiNYIgZNl9ora2dsqUKanHHo/H7XbX1dWVlJT0PcDv99P0V5cekydPPn/+/FCvlkwmjx49GgwGUz9OmDChoKBgqIMlMRFvrBfNWi35YYoTXYfrbAmttl/k+Rid1hVl55GEs8gUazid7SYZSvqfP6SJyU8GD4slU7vTOVjTnz/bLJGFUqLlrNINGTWT2WLzT5D3NWUIQpZli4r+u5UXwzA9PT39DqAoapgD+opGo9///vdttq+2KHryySevvPLKoQ7mW5uTf3luzC1XHn9hZ22R5fgepdsxRslkUkjvKqS749Zi18cdp4PZbpKhpP/5Q5qSkiPcenf7n/+3iRh5zELTn383d7FFzOv4816lGzJqJsrJfOtHBEGEw2lVNtE0bbGMsA24DEHo8Xi6urp6f+zq6up7O5g6oLu7u+8B06dPH+rVKIp69913/X5/Om9ttVbYv/ebEf8jVctVL5ze0eT73q1KN2SM0twPT4wlz288WfH4jzB3Ql7YjzAbWp45nXfn/2L8I+8UpunPP/THBu80l2fuHUo3JCNyff4yXM7MnDnz4MGDqcfHjx83m80TJ07se8D06dNbW1ubmppSPx48eHDmzJmZv68OMH4y0hYTYzqfVs/W8YyfRAqCJhhk1Se2TnBVYIDwKzIE4S233FJfX//ss88eOXLkkUceue+++xiGIQji8ccff/rppwmCKCkpufnmm9etW/fFF19s3LhRkqRrr7028/fVAZPFRHsdXJPO51aytVi9AjTDNZ4Ondd5ECZ4McGLVDGm0n9FhiBkGOb9998/cODAgw8+OG/evF/96lep58vLy30+X+rxiy++WFlZ+cADD5w5c+bdd9+1Wo2+2HkvZwXN1uv8rAudR5U2aIZ7PM3qPQjZesFZRmEqfS95Aumiiy568803+z350EMP9T52u92bN2+W5b10xlVB6XxHUIlg64UpqxCEoA1UsUOMJaPdcUe+Tem2ZEu4nneiX7QPrZY86YaznGLr9TyBlwtGbIwVG56BZhhgu1C2QcAMwr4QhAqjPY4EJyY4UemGZAt7nnejXxQ0Rff1MuEGAXeEfSEIlWYinGUkq99VnUK1vAuVMqApbl3Xy0S74oRE6LjjdwwQhMpzltM6Xt4whDtC0BpnOSUEo0mdzmticTs4AIJQec4KSq+FowlOTIRFugRV2qAlZpuZ9jrCOp3XFMYA4QAIQuU5y0i9TiVkGwRnOaq0QXucFZRe+2nCqbMS+kAQKo8ssIuxZJxNKN0Q+bGo0gZtcpXTei3n5poExo+z8msQhCpgIpx+Upf9MOF6LOMEmqTXAYtoV5wwm+xuTGf6GgShKjB+imvS4eVnuBGdMKBJep3XFG4UsD/2QAhCVXD6yXCj3u4IUaUNGmYiGD8ZbtTb5SnXFHGWjbyxhtEgCFWBKaPCursjZOsFVyUmToBWuSp0OEwYbhKcGCAcAEGoCrTHEWcTiYiu+mHCDaiUAQ3T5TBhuDHC4I5wAAShOpgIxqe3SRQsKmVAy1wVVFhfd4TxcCIZS5IFdqUbojoIQrVgdDabUCK4xgg6YUC7HPk2wkREu+NKN0Q24caIs4zEvN6BEIRq4fRTehqZ51ujNpfFyliUbgjA2DnLdXVTyDUJDEpGB4MgVAumjNJT4Wi4AVXaoHnOckpPC+KHmyLOUgwQDgJBqBaM1xHpjCXjOlnnN9yIa0/QPGeZrib44o5wKAhCtTBZTFSRnQ9GlW6IPDBdCXTA6Sd1s+KoGEvGehJUMSplBoEgVBGmlOSaddE7KhFcEyplQPPseTaT2aSPehm+OUJ5HSYzSmUGgSBUEd0EodAetdIWK41KGdA8pozSRzk314wBwiEhCFWEKSW5gC5OuSZM2gWdcOploTWuOUL7cFYODkGoIrq5IwxjBiHoBaOXeU1cc4TBHeEQEIQqYnNazRY9DEiEGwVUyoA+6GTfbIngggjCISEI1YXx6+GmEDt/gm6QhXrYNzvSGbNSFiuFYfvBIQjVhfFpPgix8yfoii72zeaaIwwGCIeGIFQXHQwTYudP0Bkd7JuNAcLhIQjVhdZ+EGIqPeiMDvbNRhAOD0GoLrTHEe2Oa3qhtTAGCEFfdLBvNhdAEA4HQaguJouJKtb2Qmtcc8TpxykH+kEV2+OhhBjR6uWpGE3GQwmyCIurDQlBqDqarpdJCGKCF8lCnHKgHyazifI6+KBWz0o+EKFKsLjacBCEqsP4SO2ecl8Vp+GMA33R9OUp+kVHhCBUHdpHcgGtdo1iTB50iSklw5oNQj4Ypb0OpVuhaghC1aG9Dl6zK45yzREGA4SgO5qe18QFMIlwBAhC1XHk25KiFA9rciULrgl3hKBDTCnJByKEpHQ7xoQPRrHc9vAQhGpEe0ktFo5KSUlojdJenHKgN1bKYqUtkY6Y0g0ZtRibICTJ7sJKT8NBEKoR43VocT8moS1mz7NaHPhSgQ5pdJiQD2D3pZHhb5Ya0T5N3hGiUgZ0jPFTvBaDMBjFAOGIEIRqxPg0WS/DNUeYUqwpA/qk0XoZLhBByeiIEIRqRPtIToMj87gjBB1zajMI0TWaDgShGlkpi4W0aG6HXq5JQBCCXpHj7HEukRBEpRsyGhLBt0TpEtwRjgBBqFKMT2P1MglOTMYlR75N6YYAZIeJoH2ktsYsIp0xK439eEeGIFQpzZ1yXKoHBourgX4xPpLTVBUbptKnCUGoUozWphJygQjjQw8M6Bnjc2ircBRT6dOEIFQpWmtdo3wQU+lB52gvqbWzEiWjaUEQqhTtcQhtMSmpmcpRrhmVMqBzTCnJB6IaKufmA1EGl6dpQBCqlNlutrutmlnSKVWchmtP0DUrbTHbTVop55aSktARozzYHHRkCEL1or0OrQwTRjpjVgrFaaB/TKlmekcj7TFHntVswx/5keEzUi8NLb2N4jQwCA2Vc3MBDNunC0GoXozXoZWt6vkAitPAEBjt7JuNSpn0IQjVS2t3hDjlQP80dEeIQu70IQjViypxCB3aKBzFeoZgEHTqrBS1cFbijjBtCEL1MltNjjyr0Kb2wtFkQop0xWkPTjnQP7PV5Mi3Ca1q76qRRCnSGaeKcVamBUGoaproHRVaomSh3WTB6mpgCJoYJhTaYo4Cm9mKszItCEJVo7VQL8MFMUAIBkL7NHBW8sEIg37RtCEIVU0Td4QoGQVDYbSw0BqHVUZHA0Goapq4I8SYPBgK7dPAShc4K0cFQahqtMcR6YyrvEQNs+nBUMhx9jibEKNJpRsyHMydGBUEoaqZLCayUNUlamI0meBFshDrGYJRmMwmyuPgW9R7ViYTUrQrThXhrEwXglDt6BJVn3J8IEKVOLAfLxgK7XWoeVq90Bolx9lQyJ0+BKHaUV5SzUHIBbHPCxiOyqvY+JYoXYKzchQQhGpHl6h6ZB5j8mBAjM/Bt6j3jhB7oo0WglDtVF44ijF5MCDaq+o59bg8HS0EodqpvHCUx3LbYDyOfFsylhQFlRaO8kF0jY4OglDtTBaTo8AmtKtxxdEEn0wmJHueTemGAOSWiaBLHJFWNZ6VkihFu+JUMUpGRwFBqAGq7R2NtMSwegUYE+0jhda40q0YBN8aJQtRMjo6CEINUG29TLQtjvUMwZhoryPapsYgFFApM3oIQg2gS1Q6g0JojaNSBoyJ9pKRFjUGIerXxgBBqAGq7hrFtScYEuNzCKocI+RbInQJzsrRkS0IBUH47LPP6uvrhzqgvb390KFDDQ0Ncr2jcVAeR6RDjYWjkbY4ghCMyea0mkymWCihdEP644NRBOFoyROEhw4dmjhx4iOPPDJv3rwnnnhi4AG33377lClT1q5dO3v27BtvvDEWU+OVlGp9tSm2ygpHYz1xs8Vkc1qVbgiAMkiPTW1dNamN6UlsTD9K8gThY4899vDDD//zn/88fPjwSy+99OWXX/Y7YPXq1cFg8JNPPjl37tzx48e3b98uy/saB+1VXb0MF4iSJSjRBuNyFNvUdlYKrVGyEBvTj5oMQdjS0vLRRx/df//9BEGUlpYuW7bsjTfe6HfMNddcY7fbCYJwuVwzZsxoaWnJ/H0NRYXDhHxLxFGE20EwLqrEprYdevkW9IuOhQx/yBoaGlwuV1FRUerHqqqqYUYKz5w589FHH/3kJz8Z6gBRFPfu3dv7ahdffLHH48m8kVpHl5Adx0NKt+Jr+ECU9OKOEIyL9Nh7jqnsrETJ6JikG4RPPfXUsWPH+j05d+7c9evX8zyfuttLIUmS47hBX6Szs/Omm2567LHHqqurh3qjWCz2m9/8xuH46qLm8ccfv+yyy4Y6WBAEu91usVjS/K/QMLcYbhbC4bDS7fgvtpkvnEyqqklGw3GcyYROMMUknXGuWQizYfVsQxZqCufPYAxyVqb5/adp2mweoe8z3SBcuHDhxIkT+z3p9/sJgigpKenu7k4mk6k36+jo8Hq9A1+hp6dn+fLlV1999YYNG4Z5I4qi3njjjdQrj8hisRgkCOnx0umuIEMxalkwQiKibfH8So/T6VS6KcYlSRI+fwVJkmSlQraEw1GgllUGY+2Bwso82mmIm0IZv//pBuEVV1wx1D9VVVXl5+cfOHAgdev2r3/969FHH+13DMdx3/jGN6qrq5999tkxt9XIegtHVTIAEOmMWWmLxYF5qGBojM/BByMqCUKUjI6ZDH/I7Hb72rVr161b97e//e2JJ55oa2u7+eabCYLYt29fRUVF6pibb765vr5+zpw5W7Zsefnllz/++OPM39doVFU4iqEIAIIg6BIV7ceEktExk6fqb+PGjePGjXvhhRf8fv9HH31EkiRBECUlJbfddlvqgEWLFs2cObO2tjb1Y28tDKSPLkntBepWuiEEQRB8ELsvARC019FzdvCSiNxDyeiYmSRJXeuVlJWVffLJJ2mOERqoWIYg2j7v6TgemnZPudINIQiCOPV/GgumOqmpFpfLpXRbjItlWXz+CmJZluiynn2zuXp9//oJRdT/rZUgiIrlRimzl/H7jzEezVBX12gAW2ADELTXwbdGpaQqbif4FpyVY4Qg1AzK44h0xNSw4qiUlIS2GOXBKQdGZ7Gb7U5rtFMV21BgldExQxBqhtmqlq3qI+0xe57VYseXB4CgfQ41rC+TKhnF5enY4G+Zlqhkh14uEGWwMT0AQRAEQXtJNZyVAjamzwCCUEv+UziqMD4YoRGEAARBEATjVcUdIUpGM4Eg1BKVXHvygQiDMXkAgiAIgvaRalgQH1N7M4Eg1BKV7EHB4ZQD+I/UvtnJhMJVbCgZzQSCUEvUsFV9MiFFu+JUMfadACAIgjBbTWShTWhTuKuGC6BrdOwQhFryVeFom5KFo3wwQhXbMSYP0Iv2kbyiw4RfXZ6iZHSsEIQak1rkV8EGYCgCoB9G6cUuhNYoOQ6Xp2OHINQY2ktyip5yWGUUoB/aRypbOIqVnjKEINQY2utQthOGC+COEOBraC/JK7oHBReMopA7EwhCjaG9Ctdq84EIjTtCgD6oInucS4jRpFINwNTeDCEINYYqtkd7Esm4MqdcIiImBJEsQMkoQB8mgipWcpgQI/cZQhBqjMlsoorsfIsypxwfiNJeB4EheYCvo5WrYhNjyRibIMfZFHl3fUAQao+C68vwwQhWGQUYiPEqVi/DB6O0x2Ey4/p07BCE2sP4FKuXQaUMwKBon2L1MnwQJaOZQhBqD+0lleoa5ZojTCmCEKA/ppTkmgVF3hoDhJlDEGoPrdRq9xKmKwEMzu62EiZTLJTI/Vtjam/mEITaQxbaE7yYiIg5ft9od9xsM9mc1hy/L4AmKDVmgQGLzCEINchEUCUOIef1MlwA/aIAQ6J9JNec6yBMCKIYFR35KBnNCIJQkxgllnTCACHAMBQ7K30kZjRlCEGoSYqccnwAq1cADIkpVS4IITMIQk1iSkmuOeddozjlAIZGex1CWyzH24XygQiNfpqMIQg1iSkl+eYIkcMzLpmQIl1x7PwJMBSzzezIz/V2oVwAl6cyQBBqkpW2mO2maHc8Z+/It0QpbHgGMCwmx/UyUmoSIS5PM4Ug1CqmNKenHN+MHhiAETC+nM7xjXTEbIzVSlly9o56hSDUqhyPzHMBTNoFGAFdSuZyKiHXjD3R5IEg1KocT1pCpQzAiHLcNYqpvXJBEGoV48vptSfmTgCMiCy0JyJigs/Rqk88KmVkgiDUKrrEEemM52aH3lgokUxKWL0CYASmnN4UhptxeSoPBKFWmSwmqjhHO/RyTYLTT+XgjQC0zllGhZtysQ2FGEvGQwmq2J6D99I9BKGG5Wx9mXBTxOnHhSfAyJhSkmvKxVnJB6IU9uOVCYJQw2gfyeekE4ZrijAIQoA0MH4y3JiTsxIDhPJBEGqY00+Gc3LtGW4SGHSNAqSB8ZGRzlgykfVln7gmAZenckEQahhTRnFNWV9oTYxgKAIgXSaLiSqy56CiO4x+GvkgCDXMxlgsDnOkK7trG3LNAu0jMRQBkCann8p6V41E8AEEoWwQhNrG+EkuywMSuPAEGBXGT3JZLhzlW6M2l9VKYnE1eSAItc3pz3qtNteMklGAUWD8WS8cxYwmeSEItS0Xp1wjKmUARsHpp7hAdgfvw00RZxkuT2WDINQ2ZxkVbsziHaEkSnxbjME+LwBps5Bmm9MqtGVxsQtcnsoLQahtjgJbUpRibCJLr8+3RMkCm9mO7wnAKGR7alO4GSP3csIfOM0mvq1iAAASiElEQVRzllJc1m4Kww2CsxwXngCj4yzPYldNpDNmtpjsLmuWXt+AEISa5yzL4rUnWy+4KhCEAKPjLKfC9dkKQq4p4izDWSknBKHmMX4qe7Xa4QbBWU5n6cUB9MpZQYWbBCmZlYKZcFOEQRDKCkGoec6ybK1tmExIQmuUKUWlDMDoWEmL3WUV2rKy2AXXKGBGk7wQhJpHFTviXCIhyL8XKNccoYrtZhu+JACj5iyns9Q7Gm4U0DUqL/yN0z4T4SzLyoBEuJ53VqBfFGAsXBUU28DL/rLR7jghEY4C7JItJwShHrgqabZO/lMOlTIAY5alehm2TnBV4vJUZghCPXBVUGw27ggxdwJgrJx+kgtGZd+PKVzPO3F5KjcEoR64KmnZg1CMJqPdcRprygCMidluporsfFDmQja2XnBVIghlhiDUA7vbaraaIp1ylqiFGwSmFLsvAYyd7L2jUlJCpUw2IAh1wllBsXVynnIs+kUBMuOqoNgGOc9KPhh15NusFHZfkhmCUCdcFXS4Xs56mXA9j0oZgEw4K2j2vJxnJVsvYIAwGxCEOuGS+44wVMu7qxgZXxDAaBifI9aTSHCyzfFl63iUjGYDglAnnBUUF4hIojwlakJbjDCZMFcJIBMms8lZSYXkuykMY0ZTdiAIdcJiN5OFdi4gT4laqJbLm4gLT4BMuavoUK08QShGk5GOGOPD4mryQxDqh6uSkmtAgj3Pu8ejXxQgU+4qJnSOk+Wl2HqB8ZMmCwq55Ycg1A/3RKbnrDxB2HOOd0/AHSFAptyVFBeIyDKtPnSWy5uIy9OsQBDqR95EpucsR2R8xsU5MR5KYCo9QObMdjNd4pClorsHQZg1CEL9cOTbLHYz3xrN8HVC5zhXFY2p9ACycFcxPecyDcJkQgo3Cq7x6KfJCgShruRNZEJnMx2QCNXy7iqcbwDykKVehq3jaa/D4sBf7KzAx6or7kl0jwxByCEIAeTinkCztXyGu9WjXzSrEIS68tUwYQbEaJIPRjFXCUAuNqfV5rbygYzGLEJneQRh9iAIdYUstJssJqFt7Ktv95zhXJU0dqUHkFHBVGfXqfCYf10SJbaed6GfJmvw905vMrwp7DrFFkxzytgeACiY5uw+xY7519l6gfI4rCTW2s4WBKHeZBqEJ8P5UxGEAHJyT2TYBkGMJsf26z1nubwJ6BfNIgSh3uRPdXadZMc2Mh/pjCVjScaLNZwA5GSxm13lYy9k6/qSxeVpVskThMlkcv/+/e+88053d/cwh0Wj0UOHDnV0dMjypjAoR77NkWcb204UXSfCBdNcBCYQAsgtf5qz6+RYhgkTnMgHo3mTcEeYRTIEoSiK119//QMPPPDyyy9Pmzbt2LFjQx25cePG+fPn7969O/M3hWEUTnd1fjmWAYmukxggBMiKgqljHCbsPMHmTWHMVlyfZpEMQbh79+4zZ84cPHjwr3/96wMPPPDjH/940MMOHjy4d+/eWbNmZf6OMLzCGe7Of4dG+1uSKIXO8bjwBMgGxkemto8Y7S92/jtUOMOVjSZBLxmCcOfOnTfeeCNFUQRBrFq1ateuXYlEot8xsVjswQcffOGFFywWFD5lnauCinNipHN0p1yolqeK7TanNUutAjA0E1Ew1TXa3lFJlLpruMJpCMLskuGvXmNj49y5c1OPKyoqEolEMBgsKyvre8zPf/7zFStWVFdXj/hqiUTi7bffLiwsTP146aWX9nupvkRRFEXZdn/Wk4JpzvZjPb7LC9P/lbYvuvOnO0f1eeLzVxY+f2WN9vPPn84E/tnpuTQv/V/pOc2RxXYzbcL/6IHS/PzNZrPJNELHclpBeOTIkfXr1w98fvv27ZWVlbFYzGb7aivz1INo9GtrKBw9enTnzp2ffvppOu8liuKuXbtS95cEQZSUlBQXFw91cDQalSQJd5kDOSc72j4JFc5Lt59TSkrtX4Smrint9/9ueLFYbFTHg7zw+StrtJ8/PcHG74iwrZw9L907kPZjPe6pFP4vDyrNz58kSXmCsKqq6ic/+cnA51MR5fV629vbU8+0tbURBOHz+foe9otf/MLv96deoaGhYceOHS6X65vf/Oag7+VwOLZs2eL3+9NpmMlkstvtCMKBHBeR9W+1200OK5XWh9N1gqWKHAX+UVyrEgQhiiJNY7ULxeDzV9YYPv9xF+VxJ2P5V7nTOloiQqcaL7i/gqYxo2kQMn7/0wrCvLy8K664Yqh/Xbhw4dtvv71hwwaCID788MNZs2b1a9zq1avr6upSjx0OR3FxcVFRUQZthpFZHOb8ac62z3t8C9PqHW37vMcze3QpCACjVTw7//yuoP+qtP4A9pzlLA4z40MKZp0MY4T33HPPM888s379+unTp//oRz/avHlz6vmFCxfeeOONjz322LJly3oP3rp165VXXjlMrIJcvJcW1L4dTCcIk7Fk55ds1Te8OWgVgJHlT2JibIJvidIlI298HTzQVXJpQQ5aBTJUjRYUFBw4cIAkycOHD2/fvv22225LPb927dpFixb1O/ihhx6aM2dO5m8KI8qf7BRjyXDDyDPrO46zrkrK5kK9KECWmYjianf74Z4RD0xwYtcJ1jM3PweNAnn+9lVUVPzyl7/s9+Sdd9458Mh77rlHlneEkZkI7yUFwQNdk8pH2FOp9bOu4tk43wByoXhO/slXGsqvKTaZh6vgaD3UXTjDleYYP2QIa43qmWd+QfuRnuGX+g03CnxLtKgaA4QAueAsoxyFtrbPR7gpDB7o9KJfNFcQhHpmd1nzJjGtnw23AGz931rLlhRjASeAnKlY5mn4e+swK+P3nOGkJOGuwjJPOYIg1LmK5Z7691rj3ODTTsONAheIlMzHhSdA7uRNZOx5Q94USknp3F8Dlcs9WP4+ZxCEOsf4yKKZeXX/r2XQf8XtIIAiKq8tqX9v8JvCwD87bYwVoxW5hCDUv8prPZ3HQ+H6/uWjLQe7+NYobgcBcs9dRVNF9rp3+l+hxthEw/ttE7/pG/S3IEsQhPpnpSzjr/fW/Lmx78r3nV+ydf+3Zcb/rMTtIIAipt5V3vkl2/RRe+8zCU6s+WNDySUFlGfkWYYgI0wdMwTP3PyEIH7x3LkJK72OAnu4QWh4v236tyupYpxvAMqw0pYZa8Yf/c25hJDMm0gTEnH69aaimXmVyz1KN81wEIRGUfo/xrmr6NOvNZmtJsZPXbC6wlUxwvxCAMgqR77twjXjA//qbHivLdoTn3hTaeF07LikAAShgTjLqFmPTVK6FQDwX5THMeFGjAgqDGOEAABgaAhCAAAwNAQhAAAYmraDsL29nWVZpVthXIFAIBKJKN0K42poaBDFwdcMgmxLJpP19fVKt8K4otFoc3OzXK+m7SDcsGHDW2+9pXQrjOvee+89dOiQ0q0wrkWLFnV2dirdCoPq6em5/PLLlW6FcR05cuSuu+6S69W0HYQAAAAZQhACAIChIQgBAMDQTJI05J5YirjhhhuOHj1qNqeV0CzL2mw2kiSz3SoYVHd3N8MwNptN6YYYVEdHR0FBQZonC8hLkqSOjo6ioiKlG2JQiUSCZdmCgpH3DNi9e/cFF1ww/DGqC0KWZdva2pRuBQAA6EFZWZndbh/+GNUFIQAAQC6hUwUAAAwNQQgAAIaGIAQAAENDEAIAgKFpZj/Ctra2rVu3trW1XXfddYsXLx54gCiKr7766rFjx6ZNm3bfffehpl9GbW1tu3btOnHiRH5+/q233jp58uR+BySTyd///ve9P1544YWXXXZZbtuoZ3v27KmpqUk9tlqtq1evHnhMc3Pztm3buru7V65ciaW/5LVly5a+RYUDv941NTV79uzp/XHlypUeD3aZz0hra+uhQ4caGhquvPLKKVOm9D7f2Ni4ffv2UCh00003XXrppQN/MRaLbd269fTp09XV1XfddVeak4u0cUfI8/yCBQtOnTpVWVm5atWq119/feAxDz300AsvvDB58uQ//vGP3/rWt3LeRj1bt27du+++6/P5Wltbq6ur9+3b1++ARCKxZs2akydPnjt37ty5cx0dHYq0U69eeeWV119/PfXZ1tbWDjygu7t7/vz5TU1NZWVlN9xwwzvvvJP7RupYbW3tuf949NFHjx071u+A/fv3b9q0qfeYaDSqSDv1ZMmSJU8++eQTTzyxf//+3ifb29vnzZvX1tbm8/muvfbav//97wN/8Y477nj99dcnT5783HPPrV+/Pt33k7Rg27Zt8+bNSyaTkiT96U9/uvjii/sdEAgEHA5HQ0ODJEmdnZ0kSZ49e1aBhuqUIAi9j9esWbN69ep+B6TOfJZlc9suo7j33nt//etfD3PA5s2bFy9enHr80ksvLVy4MCftMpzPPvuMoqiurq5+z7/yyisrVqxQpEl6JYqiJEmXXHLJK6+80vvkM888c+2116YeP/fcc73f+V4nTpygabqnp0eSpNraWoqi2tvb03k7bdwRfvzxx0uXLjWZTARBLF269OjRo11dXX0P2L9//6RJk8rKygiCKCgomD179t69e5Vpqx71XbsnEok4nc5BD/vd7373/PPPf/7557lql4EcOHBg06ZNO3bsiMfjA//1448/vuaaa1KPly5dun///kEPgwxt3br1lltuyc/PH/hPjY2NmzZt2rZtG9YDkcWgXZr9vud79+5NJpP9Dpg/f77b7SYIYvz48eXl5Z988klab5dxg3MhEAgUFxenHo8bN85qtQYCgaEOIAiipKRExq2qoNf+/ft37tz53e9+d+A/LVmypKOj48SJE4sXL960aVPu26ZjFRUVxcXFnZ2dTz311Pz583me73dA3++/x+NJJpPBYDDnzdQ5QRBee+21QQdo8/LyLrzwwp6enp07d06bNm1g3ynIot/3PB6Pt7e39z0gGAz2DQKPx5NmEGijWMZqtSYSidRjURRFUey3ZI7Vau27Q2k8Hh9xTR0YrVOnTt1yyy1btmyZNGlSv3+y2+3vv/9+6vEdd9xx9dVXr127lmGYnLdRn5588sneB7Nnz962bdu6dev6HtD3BEk9wPdfdn/5y18KCgquuOKKgf+0cuXKlStXph5/5zvf+elPf/rmm2/mtnWGMOL3fMxBoI07Qr/f3xvsqQc+n6/fAU1NTb0/NjU1lZaW5rKFuldTU7NkyZKnn3761ltvHf7IBQsWJBKJxsbG3DTMUGw22/z58wfWy/Q9QZqammw2GxaDlt22bdvuv//+1ADNMC677LJz587lpklG0+97TtN0v27qMQeBNoJwxYoVu3btikQiBEG8+eabixcvTt1tHDlypK6ujiCIRYsWtbe3p3ZLP3369MmTJ3u7kiFzdXV1y5Yt27hxY78toT/99NPU104QhN4nd+/eTdP0+PHjc9xIHev9eFmW3bNnz4wZMwiCiMViH3zwQapMacWKFTt37kyNC+7YseO6666zWCwKNlh/amtr9+7de/fdd/c+w/P8Bx98kLov6f0fJEnSO++8c+GFFyrTSr1bsWLFW2+9lbrn27Fjx4oVK1LPf/rpp6mAXLZs2RdffJG6Ujxw4ADP8wsXLkzrpWWp8Mm2RCJxzTXXzJkz5+677x43bty+fftSzy9evPhnP/tZ6vHmzZt9Pt/q1avLy8t7nwRZLF++nGGYOf/x4IMPpp6fNWvWb3/7W0mStmzZMmPGjDvvvHP58uVut/sPf/iDou3Vm8LCwhUrVqxatcrn811//fXxeFySpNSZf/78eUmSIpHI5ZdfvmDBglWrVhUVFR0+fFjpJuvNhg0brrvuur7PnDx5kiCIjo4OSZKWL1++ZMmSu+6666KLLpo8eXJdXZ1CzdSPRx99dM6cOQzDjB8/fs6cOam/+TzPX3LJJZdffvntt9/u8XiOHz+eOvjiiy9+6aWXUo9/+MMfVlZWrl692uv1vvjii2m+nWZ2nxBFcc+ePe3t7YsWLfJ6vaknT5065XK5em9+jx8/nppQP2vWLOVaqkM1NTUsy/b+6HK5UlNc//3vfxcXF6dGrQ8dOlRbW5uXlzd37lzMJpbX+fPnjxw5EovFpkyZUl1dnXoyHo8fPny4uro6NQoSj8c//PDD7u7uq666qm+9AMiipqbG7Xb3/uUhCCISiRw9enTOnDkWi6Wzs/PgwYNdXV1+v3/BggVYzSNzZ8+e7e7u7v1x8uTJqVrQVEdIKBRasmTJuHHjUv96/PjxkpKS3q/9oUOHampqZs6cOX369DTfTjNBCAAAkA3aGCMEAADIEgQhAAAYGoIQAAAMDUEIAACGhiAEAABDQxACAIChIQgBAMDQEIQAAGBoCEIAADA0BCEAABgaghAAAAzt/wPXqEp/lcXgcAAAAABJRU5ErkJggg==", "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "execution_count": 10 } ], "cell_type": "code", "source": [ "ρ = scfres.ρ[:, 1, 1, 1] # converged density, first spin component\n", "ψ_fourier = scfres.ψ[1][:, 1]; # first k-point, all G components, first eigenvector\n", "ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n", "ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\n", "\n", "using Plots\n", "x = a * vec(first.(DFTK.r_vectors(basis)))\n", "p = plot(x, real.(ψ), label=\"real(ψ)\")\n", "plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n", "plot!(p, x, ρ, label=\"ρ\")\n", "plot!(p, x, tot_local_pot, label=\"tot local pot\")" ], "metadata": {}, "execution_count": 10 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.1" }, "kernelspec": { "name": "julia-1.7", "display_name": "Julia 1.7.1", "language": "julia" } }, "nbformat": 4 }