{ "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 1D example\n", "and we show how to define local potentials attached to atoms, which allows for\n", "instance to compute forces.\n", "The custom potential is actually already defined as `ElementGaussian` in DFTK, and could\n", "be used as is." ], "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 CustomPotential <: DFTK.Element\n", " α # Prefactor\n", " L # Width of the Gaussian nucleus\n", "end" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Some default values" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "CustomPotential() = CustomPotential(1.0, 0.5);" ], "metadata": {}, "execution_count": 3 }, { "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::CustomPotential, r::Real)\n", " -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)\n", "end\n", "function DFTK.local_potential_fourier(el::CustomPotential, q::Real)\n", " # = ∫ V(r) exp(-ix⋅q) dx\n", " -el.α * exp(- (q * el.L)^2 / 2)\n", "end" ], "metadata": {}, "execution_count": 4 }, { "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": 5 }, { "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", "positions = [[x1, 0, 0], [x2, 0, 0]]\n", "gauss = CustomPotential()\n", "atoms = [gauss, gauss];" ], "metadata": {}, "execution_count": 6 }, { "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", "model = Model(lattice, atoms, positions; n_electrons, terms,\n", " spin_polarization=:spinless); # use \"spinless electrons\"" ], "metadata": {}, "execution_count": 7 }, { "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 `gauss` 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 log10(ΔE) log10(Δρ) Diag\n", "--- --------------- --------- --------- ----\n", " 1 -0.143565650196 -0.42 8.0\n", " 2 -0.156033427558 -1.90 -1.10 1.0\n", " 3 -0.156768190389 -3.13 -1.56 2.0\n", " 4 -0.157046016501 -3.56 -2.32 2.0\n", " 5 -0.157056301334 -4.99 -3.41 2.0\n", " 6 -0.157056359749 -7.23 -3.59 1.0\n", " 7 -0.157056406037 -7.33 -4.35 1.0\n", " 8 -0.157056406911 -9.06 -5.39 2.0\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380298 \n AtomicLocal -0.3163471\n LocalNonlinearity 0.1212609 \n\n total -0.157056406911" }, "metadata": {}, "execution_count": 8 } ], "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": 8 }, { "cell_type": "markdown", "source": [ "Computing the forces can then be done as usual:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "2-element Vector{StaticArraysCore.SVector{3, Float64}}:\n [-0.055686637674769585, 0.0, 0.0]\n [0.05568538350058255, 0.0, 0.0]" }, "metadata": {}, "execution_count": 9 } ], "cell_type": "code", "source": [ "compute_forces(scfres)" ], "metadata": {}, "execution_count": 9 }, { "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": 10 }, { "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+gvaeTAAAgAElEQVR4nOzdd3wURd8A8Nly/S65S3LpvYcQQgq9JlIChKIEUSmCgPBgRX1UBBVFVFTURxQLIlh5VBAUQu8gNb1AOun9er/b8v5xvJGHEkJyuTrfD39clr3dubm9+e3MTkFomgYQBEEQ5KpQWycAgiAIgmwJBkIIgiDIpcFACEEQBLk0GAghCIIglwYDIQRBEOTSYCCEIAiCXBoMhBAEQZBLg4EQgiAIcmkwEEIQBEEuDQZCCIIgyKXZXSBcs2aNTqfr4c4EQcAp4mzIZDLZOgkuDea/bcH8ty0L5r/dBcLvv/9eKpX2cGeTyURRVL+mB+qGXq+3dRJcGsx/24L5b1sWzH+7C4QQBEEQZE0wEEIQBEEuDQZCCIIgyKXBQAhBEAS5NBgIIQiCIJcGAyEEQRDk0mAghCAIglwaDIQQBEGQS8NtnQCof5koUK6gKxT0dRUgKGCkgDsTeLNBlDsyyANhwBshCLIdigbX5HS5gm7RAoURYAhAERDMB1FuSKwQ4cLi2VpgTjunVh34bzV1pIn6u5X24yKxQiRcAJgYYKKgWgnOt4EyOVWlpAd7IrNC0axQJFSA2DrJEOQq2nVgTx31x3XqYjvty0UGihA/LnBnApIGJAUud4BKBVWtopM8kYkB6CPhSJQ7/Hn2LxgInc3RJvqzUvLvNvrBEHRpDPrzeFTEuvOeWgKcb6N311LD/iIHeSCrBmJTghD4g4Oge/rggw9ycnJIksQw7L7eqDCCejUtMQBPNvBhg3TWjVaZdgDa/38fBIBoACJoIDeA3Qb6Cx3g4iCQB3w4rv7zHD9+/MqVK/vjyIi9TVodGBh46dKlgICAnuys0+mYTOb9XovO6mADvS6PVJvAvwehc8JRXo9vckwU+K2G+riEMlHgg6FYRmBPf24qlUogEPQyuVCfwfy3ldGjR2dkZMTExNg6IS4kLy+vpKRk3759XVsseP3DGqEzKJPTqy6StWqwIRWdFYKi93nfyEDBvEh0XiT6Zx216iL5CQ98MQqLdHPxu08I6k5aWtqoUaNsnQoXwuFwSkpK+ungsLOEYzNR4O18aux+YnIgWvQQ/lDofUfBm80MQYsfwqcEoSP+IjYVU6R9NRZAEAT1CxgIHViJjB76J3Gpncp/EH9+IGqRLqA4Cp4fiF6aiWfXU+nZRJMGBkMIgpycZQLhTz/9FBUV5e3tvWTJkm6W1T1w4EBqaur+/fstclIX9205lZ5NPBuPZk/GA3gWbsYMFyDHpuKTAtHUvcThRhgLIQhyZhYIhNeuXXvqqad27NhRXl5eV1e3fv36O+6mUCheeeUVqVQqkUj6flJXpiXA/FPk5lLqTCa+OLq/6vQoAtYMRn99AF9yltxUDFc/hiDIaVmgGP3uu+9mzZo1atQokUi0du3abdu23bEn6osvvvj888+LxeK+n9GV1avpMfsJDAEXZ+Cxwn7vzzLWF7kwA/u5inriDGmE0RCCIGdkgUBYVlaWmJhofp2YmNje3i6TyW7Z58SJE9XV1U888UTfT+fKzrfRw/8i5kei34/DONbq8BvEQ85Ox2UGkHmYUJmsdFIIgiCrsUBpKpVKuwZzuLm5AQA6Ozs9PDy6dtBoNM8+++zevXuRHowH7ejoCAwM7PozMzNz586dd9vZpcYR7m1AV+Xg3wwnJvoZ1Gprn33HcPBCDj5uH7V7rEnMvlHj12g0PflOoX4C899WKMpJmkfa2tpSU1MbGhrMf37zzTdtbW2vv/763fY/ePDg3r17v/76a2sl8H+QJKm+qezr4fXPZrNx/B6RzgKB0MPDQ6lUml+bX3h6et68w7p164YMGaJQKHJzczUaTW1tbV1dXUhIyB2PJhaLez6gHsMwFwmEn5RQnxRTx6ZhiR5MW6VhWxpYl0dOPYUem4KZu+fQNM3n822VHgjmv62gqJP0tydJsrOz0/xap9OtX78+Ly+vm/0zMjL+/e9/l5aWxsfHWyWB/wPDsJsveAte/xYIhFFRUaWlpebXJSUlnp6eN1cHAQAURRUXFy9fvhwAUFtbu2PHDgzD1q5d2/dTuwIagNeukPvq6fMzsEBL9w69X+uSMT6DGpdNHp+KhfBhXQSC7EJdXV1hYWFERMTvv/8+bty4tLS0srKyffv2EQQxc+bMAQMGAAA0Gk12dnZpaSmXy50+fbp548127dqVmppq7sZx/PhxX1/f+Pj4HTt2ZGVlqVSqc+fOzZkzB0GQefPmffXVV5s3b7bB5+w3FrivWbx48R9//FFQUKDVat9///3Fixebq6vr168/ePAgAGDTpk05/y8+Pn7dunUwCvYQSYNlZ8nTLfSZTNzmUdDspQT0+YHo+GyyVgWHVUCQXSgoKFi5cuUzzzwjEokQBNm/f//UqVNRFOVwOBkZGWfPngUA5ObmnjlzJjg4mKKo9PT0nJycWw6SnZ2dnp5ufr1t27YTJ04AAFatWiWRSCorK9etW2f+r7S0tAMHDljvs1mFBWqECQkJH374YWZmplqtnjJlyptvvmneXl5eHhQUdMvOAQEBcHbEHjJSYP5JUmGkj07Fez5xqBU8PQBFAUg/QGaPR+LglwlBALyVR51rs9KDQwwBP4zDvTn/s1GlUu3du9fcSyMqKmrbtm1paWkAAB8fn/fee2/MmDFjx44dO3aseWcmk7l9+/bU1NSbj1BcXLxo0aJ7nj06OrqmpsbJ5rm1TPm6bNmyZcuW3bLxp59+un3PP/74wyJndHo6AmQdJ5go8tcknGV/z0BXDkBJGkw7yTwznbaTqioE2dAjEcgoHyv9UDEUeLJv3RgbG2uOggqFoqqq6u2333733XcBACqVqqOjAwDQ3Nz8r3/969q1awwGQ6fTDRw48JYjqFQqLpd7z7ObH8vBQAj1Ow0BZhwh/LjIjrEYbq9P5Z+JR1U6cuJB9PS0W29OIci53f5UIMYdiXG3QUq6sFisrhcYhn3zzTdeXl7mLeauPWvWrImLizP33v/8888PHTp0yxHEYnHXyDc2m63X67v+S6/Xs9k3Yq9EIsEw7JYekY7OXktZF6YhwPTDRJgA+WGc/UZBs2diiLnhyMSDhNRg66RAkLWQNGjU2DoRd8dms8ePH//LL7+I/p/RaAQAtLe3h4eHIwhiNBrvOCZtyJAhRUVF5teJiYmHDx8mCML85759+wYPHmx+XVhYmJiY2BV3nYN9F7SuR2kCEw8QsUJk6xisL+tIWM26ZOwBf2TGEUJL2DopENT/aAAWnSbtbRnXW2zdujU7Ozs5OfnBBx+Mj4/fuHEjAGDFihWrV69+8MEHU1NT7zg+bc6cOV3VxBUrVri7u0dHR6vV6gceeCA/P79r7sxDhw7NmTPHap/FOmDTqB1RGEHGISLRA/liFOYIQfCGTcOxZWfJmUeJ7Mk4E95ZQU7tpUtktdLunotPnTq1q8MnACAsLOzixYt1dXVtbW1hYWE+Pj4AgOnTp1+7dq26ujo6OlogEJirib6+vvX19eZ3paenv/TSS8XFxQkJCSwWa/fu3e3t7ZGRkT///POwYcPM+2g0mv3791+4cMHqH7F/wXLLXsgMYOJBYrg38uVoR4qCAAAEgK9GYRwMWXIGrmAIObMNBdTxJvpABm5vrTUMBuOWrisoioaFhQ0fPtwcBc28vb1HjBjh6enJZDLNfV5QFO2a/xlBkM2bN1+6dOnm/TEM8/X17dqSk5OzZs0a55syGgZCu2COgmN8kU+GO1gUNMNR8Gs6Vqum/32JtHVaIKhfbCunviunDk3BhTab3KnfjR49eunSpTdvee6559zd/+kFNG7cuMWLF1s9Xf0OBkLbkxrAhINEmh+yaZj9jZPoMQ4O/pqIH2qk4ZpNkPPZV0+9nkseysB8XayD9Lp164RCoa1T0e9gILQxiQFMOEBM8Ec+dOQoaCZigUMZ2Gel1C/VMBZCzuNCO730LPnnRDzK3U7ba7KzsydNmtRPB9+3b99TTz3VzQ7Xrl2bPHlyP53dOmAgtCWZAWQcJCYEIBuHOnwUNAvkIQczsBcukkea4ONCyBlUKemsY+T34/AhYjuNggCAhISE5557rj+OTNP0K6+88vTTT3ezT1xcHIPB2L9/f38kwDpgILSZdh0Yn01MDkQ+cJYoaDZAiPz2AL7gFFEohbEQcmwtWjDpILkhFc0ItN8oCAAgCEKn0wEADAbDxo0bm5qa1q1bt3r16vr6eq1W+8knn7z88sv5+fnmnTs7O7du3bpq1ar169dXVFR0HaSuru7tt99es2ZNdXX1p59+qlAoAAAnT57k8XhxcXEAgNOnT588eRIAsH379oaGBqVS+fHHH5vfu3Dhwi1btlj5U1sQDIS20aoDadnEg6HIO6lOFQXNxvoiX4zEMg+T9WoYCyFHpTKBaYeJJTHoomh7LyeLi4s/++wzAIDBYHj11VeXLFkSGBgok8kmTpxoXgVBKBSmpaWZ51o7f/58Y2Pj0KFDmUzm6NGjq6urAQAdHR3Dhw83GAyxsbHLly9fvXq1eZaZAwcOjB8/3nyWw4cPm6fb/uyzz2pqauRy+RtvvGH+r/Hjx588edIcjB0RHEdoA00a+oED5MIo9LXB9v4D67WsMLRJA6YcIs9Nx0VONQcF5BJMFJhznBgqRtb07Eeq2L/dUN7dSn6WhDM9n1iLCUR3+/8NGzakpKQsXrzYy8tr4sSJ5o6gR44cOX36dFZW1owZM2bMmEHTtFwub2ho2Llz59q1a7du3Tp+/PgNGzYAAJKSkhISEsyHKi0tnTVr1j1T5O3tzWQyq6urb5/C1CHAQGhttSp6wkFyRRz6UoLTRkGz5waiDRp61lHi8BSc7YT1Xshp0QAsO0syUeSLUT29cPmjMzmJo/s1VV0QBMX43fXkNEcjDMO8vLy6IpO3t7dEIgEAlJSULF26VC6X8/n81tbW6dOnAwAqKiqSkpLMe8bHx3fNoKbVajmcHvWU5fF4KpWqt5/JxmAgtKpKBT3xIPnvQehTA5w8Cpp9MBSbd4pccIr8Nd0xZoyDIADA6zlkmYI+MRXv+aheTCjGhPYyzBzD/onf5hm3AQAIgphnhnvuueeeeOKJJ598EgDw6quvmptAhUJh14zbarXaPO8MAEAsFkulUvNrDoejVCq7jqzT6bpiJEVRUqn05sH7jsUlimM7USyl0w6QbyS7ShQEAKAI2DEW69TTL8CB9pCD+Ooa9dt1et8knOuk1YSuiCWXy3/77TfzxmnTpv38888tLS0AgE2bNnXtPHTo0MLCQvPr5OTko0ePdlX7du3alZycbH5dXl7u7u4eFhZmtU9hWa5SItvcpXZ64kHi42HoE3b/4N2yWBjYMxE/0Ux/BAfaQ3bvrzpqfT51KAMT37bgn9N46aWXlixZMmXKlDFjxqSkpJg3Tpw4cfny5YMHDw4NDaVpmsPhmNcmzMrKOnr0KEmSAIBp06ZNnz49JiamoqJi3rx5f/311+bNm81v379//5w5cxDEYZt9aDsTEBDQ2NjYw521Wi1BEP2aHos40kh5/2Q8UE/ZOiEWplQqe7hno5oK2Wn6uYrs1/S4mp7nP9QTf7dS3j8Zczru/TsdNWrUuXPnrJCkHjIPn6Bp2txE2bVdLpebTCbza7Vardfrza8lEklhYaFOp9PpdBqNpmt/iqIoiiopKREKhRR1Ix+ysrL++OOPrn00Gk1sbOzevXtvPvvAgQPLy8v77fPRNE3v27cvMzPz5i0WvP6dtPJvT/bUUiv/Jnc/gI/2ddjbpT4L4CEHMrD0bMKThUy27yFZkGuqUtJzjpM/jMNTvBzv+sQwzPxcEEEQkeif3qQ3TxPK4/G6Xnt4eHh4eNxykKVLlyYkJJhMpq+++mrt2rVd1buNGzfePFiey+Wy2eyb512rqqpauXJldHS0RT+TVcFA2L/+U0JtKqaOTsUHihzv12VZA4TI7gn4Q8eI7Ml4qgOWNZATa9TQEw6QG4eirnyXNm/evLy8PARBfvzxxxEjRnRtDw8Pf/bZZ2/ec9GiRcHBwV1/xsTExMTEWC+h/QAGwv5CA/BWHvl7DX12OhbCd91f181G+SDbxuDTDxOnMvEYe522EXI1CiPIPEw+E4/Oj3St5/e3SEtLS0tL68me/TSdmw259Bfff3QEmHucPNlMn5uOwyh4s8xgZMMQbMohslkLJ52BbE9LgGmHiYkByIvOPq4X6gb87i2vVQfGZxMsDByZAidVuYMnotHlsejkg6TMYOukQK7NPH1MpBvygeOv/QL1BQyEFpYvoYf/SUwLRn8Yj7Hgj+suXklEMwKRaYcJDWHrpECuiqLBotMkjiLfjnHI1bBvQVHUV199tXz58g0bNpjny4Z6DgZCS/qthso4RHw0DH0jCc6jcg8fDMPihMhDRwkDHGoP2cIzF8gmLf1rOoY7RSn45JNPZmdnT5kypaamJi0tzTzyD+oh2FnGMggKrMkhf79OH52CD/KAQfDeEAC+GYM9coJ87CTpNIUR5ChWXyGvdNDHplpsFtwv87bntxVb5lj3giDIe+PWenD+GSbR0NCwe/fuhoYGPp8/c+bMuLi4Q4cOTZs2zTrpcQIwEFpAmw48eoJgYeDKLNwTPhTsMQwBP6dhs44ST5whd4yDk5FCVvJuAbW/nj6VibsxLHbMOXEz00PHWOxw3WJizJujIACgrKwsMjKSz+cDABAESU5OLi8vh4Gw52Ag7KvjzfTjp8llMejrSSgsyu8XEwW7H8CnHiZW/E1+PdoZHtVAdu6TEur7SurUNAvfs3pxPLw4tw5RtxqtVnvzdNhKpdLNzc1WiXFEsEGq90wUWJNDPn6a/GEc9mYyjIK9xMHBvkn4VRn93AX4VAPqX1uuUp+XUsenYn5cWyfF0ioqKs6cOQMAqK6uPnPmTA9HBEJmsEbYS2VyesEp0ocD8h/EnXh+XuvgM0D2ZHzSQeKFi+THw2FfW6hffFNGfVBEnZqGBfKc8KY1KSnprbfe0ul0FRUV7777bkREhK1T5EhgILxvJA0+K6XeKyDfScWejIVVastwZ4LDU/CJB4gXL5Gb4KAuyNK2llEbCqgTU7FQgRNGQQAAn88/fvx4bW2tWCy+eVpRqCdgOX5/SmX06H3EvjrqwgwcRkHLEjLB0an4mRb6xUsknHUGsqCtZdQ7BdSJqViEm3NGwS6hoaEwCvYCLMp7SkuA166QadnE4mj0+DTc6X9RNiFkgiNT8HOt9LPnYSyELOPzq9S7hU4eBQcMGLB8+XJbp8KBwUDYI7uuU/G7iXoNKJrNeDIWdovpRyIWODoVz5fQy8+RFAyGUN98VEx9WkKdmubMURAAEBUVNW/evF6/naZpvV5vkZSYlzy0yKGsCQbCe7jSQY/bT2wooLaPxX4aj/lybJ0gF+DGAIcy8CoFveAUaYLL2kO99WYuua2cOj3NJZZ/KS4uLioq6n6fK1euVFRU3L49Pz9/wIABFklGbGxsQUFBr99uMpl+//1364dSGAjv6qqcnn2MfOgYuTAKzZ2Fj/dz/t+S/eAzwIEMXGUCs4+RejiqArpPNADPXyT31dNnMvEAZ+wjerudO3f+9NNP3e/zxRdfZGdnWyc9vaPT6R5++GHrzw8He43eQZGUfreAOtVC/XsQ9tN4jAMzyRbYGNg9AVt0msw4RPw5EXdn2jpBkIMwUWDJGfK6mj4xDRe6xmVz9erVkydPUhT16quvRkZGLl26lCCIr776Kjc3NyAg4Omnn/b19b148WJubm5DQ0NbW1tSUtLcuXPveCij0fjll1/m5+eHhIQ888wzXl5e5u2nTp3au3evVCpNTk5+7rnnSJLcsWPH5cuXSZIcM2bMwoULUfSu1aqPPvooPT39t99+6+joyMrKmjJlinl7eXn5tm3bJBLJ+PHj58+fjyDIli1bAACvvfYaiqJPP/10YGCgpbPqzmCN8H8ca6KnHiamHCKHiJHquYwXE1AYBW2IgYIfx2OJHsi4/USL1tapgRyBhgAzjxJyIziS4SpREAAgFAp9fX19fX1TUlKioqIAAI899tj+/fvnzJlDkmRKSopCoRCLxZ6enkFBQSkpKWFhYXc71EMPPXTixImHH35YpVINGTJEo9EAALZt2/b4448nJSU9+uijnZ2d5meKlZWVM2bMmDVr1tatWzds2NBN8r755pt58+YlJCRMnjx56dKlf/31FwCgqqpq5MiRPj4+M2fO/PTTT1977TUAQGxsLAAgOTk5JSWFy7XerAewmAcAALkR/FhJfV1GIQCsSkD3TEDhCkp2AkXAf0Zg7xVSo/YRByZjsUKXaOaCeqddB2YcJeKFyNejrT2Ne/UfzbJrauucC2UgA1eEMd3+Kb39/f3j4uIIgpgzZw4AoKqqKjs7u6mpSSgUTp069cqVKzt27HjuuedCQ0MTExPN+9xRYWHhuXPnmpqaeDze1KlTL1y4sHPnzqVLl65du/bHH3+cMGECAMBcn+Pz+Rs3btTpdG1tbS+//PKaNWtef/31btK8YsUKc3cetVq9adOmGTNmbN68OSsr68UXXwQAREREpKSkvPnmm+np6QCArKwsHLdqbHLpQGiiwNEm+qcq6mADlRGEfj4SG+eHwILWDq1ORAN5IC2b+PUBfKwv/IqgO6hQ0FMPk/MjkTeTbTBpbcgUn4BxXlY6GQJujoK3q6qqCgsLEwqF5j+Tk5Pv2Efmjm+Mjo7uGolofqNMJmttbR0+fPjNe2o0mkceeaS8vDwkJESn07W0tHR/5MTERPOLwYMHr1mzBgBgrlCaN8bHxyMIUldX5+fn15N0WpwrBkINAY41UXvr6H11VIwQmR+Jfj6S4QFXjbBvCyLRAC7y8HHiw2HYgkjYpA/9j9Mt9CMniPeGYIuibXNt4BwM59i4Hamrs6VQKJTL5V3bZTKZj49PT44gFAplMtnNbwwMDBQIBAwGQyKRmFe3MNu2bRuDwTDH10uXLnU99rubrvTI5XKRSHRLIjUajcFgMG+3CVcpUAwkONdKbyigJh4k/H42fX6VSvJECh7C/56O/ysOhVHQIaT7Iyen4W/lUWty4BBD6B/byqlHThA/p+G2ioL2QCwW19fXm18nJibSNL1r1y4AQENDw549e8yB6uZ97mjIkCFyufzAgQMAgOrq6gMHDmRkZOA4PmPGjHfeecfcmbOjo4OmaZ1OhyAIAIAkyf/85z/3TN7WrVsJgqAo6ssvv8zIyAAATJky5fvvvzfHws2bN6ekpHh7ewsEAjab3X0i+4Nz1ghJGjRp6ColuCqnS6R0bid9VU7Hi5Cxvshz8VjaRITnnJ/b+cUJkYsz8NnHiNnHyB/GYwLLrScHOSKCAv++TB5ooM9k4lHuLt1mPnfu3D179oSGho4ZM+bHH3/cuXPnokWL1q1b197evnr16lGjRgEAlixZsmDBgtDQ0NmzZ2/atKnrvSiKstlsAICbm9vOnTuXLFkiEAg6Ojo2bNiQlJQEAPjiiy8WL14cGhoqEon0en1ZWdmiRYt++OGHuLg4iqKysrK6mmHZbPYdu49GREQMHDjQYDCEh4d/8803AIB58+bl5ubGxsYKhUIWi7Vz504AAIIg69evT09PxzBsz549gwYN6v+cAwAAxN5mAQgMDLx06VJAQEBPdj5wXX+6A0cQRGMCShNQGEG7jm7VgWYtLWYjkW4gTogMFCFJXkiSJ2KppaihLiqVSiAQWP+8Rgo8e54810bvnYhFOvWMId2zVf7bCYkBzD1OMFCwM93aHURHjx69ceNGc3SxZx0dHR4eHhh232Vfe3u7l5fXLSHNvOqhr69v15bm5mZPT08W6x5NatHR0T/88MPgwYPVanXXeAwzk8mkVCo9PT3vmaT9+/c/+/5XI9/5040BnhqAxosQC17/jl0zYqFAyAAoigRwgRsTCJlAzEb9uMCfi8Bun06MiYKvRmNfl1Gj9xHfjsEzg103FrqsvE56znEyKwx5dwgGF3S+G7FY3Ls3ent7376Ry+XeMqTB39+/58dks9nmeufNGAxGT6KgmR8XZAQiahPgW7opyLED4UhvenwggmGu+2DAlS2PRRM9kLknyEsdyLpkWBq6kG3l1Gs55JaR2Oww+Nt3AOvXrw8PD+/7cTxYyPz+6SgHLyPIgQ33RnJm4Rfb6YkH4Ih7l6A2gYWnyE9LqDOZOIyCjmLu3Ll3rGLaD3glQY5NzAaHM/A0fzRlr+lwo3098IYsq1BKD/mTYGLg0kw8xrW7xkCWBQMh5PBQBLyehP43HX/yHLnqIpyk2wnRAHxSQk06SLyRhH47BuM69iMdyO7AQAg5ibG+SMGDeLMWDP2TKJTCqqHzaNDQkw8Su69Tl2bgj0bAIguyPHhVQc5DxAK/pmP/HoROOkhsKKAIuJah4/u+kkrdS4z3Q09n4qEC2BwK9QvYxAA5mwWRaJofsvQsubeW2jYWG+QBS0+H1KihV5wjm7Tg6BTc3r5EFov18MMPs1gsBE5ObC1arbb/Bm7CQAg5oUAecjAD31FBTTxIPBmLrhmMwekUHAhFg6/LqDdzyWfjsVcSUYb9tVvt2bOns7NTo9F0zU8NWcHNY/ktCwZCyDkhACyORjMC0WcvkIP+ILaMxCYEwJt3B1AkpZefI3EUnMrEB9jrqltubm5ubm4uPrOPM4GBEHJmflzw+wNYdgP95DlyqBj5aBgayLPTshVSGMG6PPKXauqdFGxpLAq/J8hq7K/RAYIsbVoQUjobjxWCpD3EuwWUjrB1gqD/RdFgWzkVt8ukJcDVLMYyGAUh64KBEHIJHBysS8Yuz8QLJPSA3cTOariOk7040Uyn7iV+qKT2TcK/Ho15wjXRIKuDTaOQCwkTIL89gJ1rpV+8RG4qpt4fAh8c2lKBhF59haxSgveGoFlwvjTIdmAghFzOaF/k4kx813XqqfNkABe8k4qN9IHh0KrK5PS6POp0C7U2CXsy1h77hUIuBV6AkCtCACNzjqkAACAASURBVJgThpbOxudHofNOkVMOEefbYFupNZTJ6YWnyHHZxGBPpGou46kBMApCtgevQch14Sh4Ihotn4M/FIrOP0U+cIA42gTDYX/J66QfOUGOzyZihUjlw4xXE1EebJCC7AO8EiFXx0TBslh0cTT6SzX1wkWSiYGXEtCsMFhTsQwagCON9KZiskwOViWg345hWHxVVQjqIxgIIQgAAHAULIxCF0Sh2fX0JyXky5eppwagS2JQ8a1LakM9pSHAT1XUZyUUEwOrBqKPRsB7C8hOwUAIQf9AAMgMRjKD8UIpvbmUivndlBmELo9DR8HeNPejVEZ/U0b9XEWN80M/H4Wl+cHcg+waDIQQdAeJHsi3Y7APh2LbK6hlZ0kEgCdi0HmRqC/H1imzY0oT+K2G2l5B1anBE9FI/kN4EJzHB3IEMBBC0F2JWOCFBPSFBPRsK729ghqwyzTSG3ksEp0ZAjt6/MNEgSNN9C9V1IEGKt0ffTURnRqEYjACQo4D/poh6N7G+CJjfLHNBLanlvqpinrqbzIjCJ0ThmQEoi67WjpBgZMt9O/Xqb21VLQ78lgE+tlIBpwXBnJErvojhqD7x8PB/Eh0fiTaqQd7aqmvr1FPnCHT/NAZIcjUINTHNVpNFUZwpIn6q44+2EBFuSNzwtCcWXgwH1YAIQcGAyEE3TcvNlgWiy6LRWUGcKCB+quefvGSKUKAZAQhkwLQEd4I7lzdIykaFEjoo0304UYqp5Me7YvMCEY3DsX9uTD+Qc7AkoHQYDCwWN21jJhMJgYDjiGCnIeIBeZFovMiAUFhf7fRR5qoFy6SFQp6pA8yzg8d44ukeiEsx1wTmKRBkZQ+20qfbqFPt1DeHGRSIPLSIGy8H+KyrcGQs7LMFd3Y2Pjoo48WFRUxmcxPPvlk/vz5t+zw5ptvbt++vaOjg8fjPf300+vWrbPIeSHITuAoGOeHjPPDNqQCmQGcaaVOtdDPX6DKFHSiBzJUjAwVIyleSKQ7Ys91qAYNndtJ53TQF9vpKx10IA8Z7YvMDkM+H8nw49o6cRDUbywTCJ9//vmBAweePn06Nzc3PT09LS0tICDg5h3EYvHRo0djYmLKysrGjh0bHx8/Z84ci5waguyNiAVmhqAzQwAAQG0COZ30pXb69+v06hxKqqcHeSIJIiTBA4kVInFCxIbjMaQGUK6gr8roq3K6SEoXSmgUASleyBAx8mICNswb8YA9XyDXgNB0XydXlMvlYrG4srIyNDQUAJCZmTl27NiXX375bvtnZWXFx8e/9dZbd/zfwMDAS5cu3RJH70an0zGZTAxzzLYnx6dSqQQCga1T4UikBlAopYul9FUZfU1OX5PTehJEuSPhAiSUD4L5SBAf+HMRfy4Qs+/dptqT/Cco0K6n23SgSQMaNXS9mq5Vg+squlJBEzSIcUfiRUisEEn0QBI8AHzmd1/g9W9bFsx/C9QI6+rqmEymOQoCAAYMGFBdXX23naVS6blz51asWHG3HWiaVigUXO6Nhhg2m83huEZvPMgFeLBAmh9y80wrciOoVNDXVXStGlQq6WPNoFVLNWtBu45mYcCLjXiwgJAJ3JkIBwc8HHBwwP7/AGk04kwmaX5tooDaBPQk0BFAZqQVRiA3gk49rTIBMRv4cJBAHgjgIoE8ZFoQCBOgkW6IN/xhQRAAwCKBUC6X83i8rj8FAsH169fvuKfJZFq4cOHkyZMnTJhwt6NJJJIRI0ag6I1edxMnTvzuu+/utjOsEdqWRqNB7PqZlwPAAYjjgDgOAN63/pfShEgMQG4EciOiNAEdCXQkoiWAkbqR5zhtZIEbrzGMDmUDJkpzceDOAO5M4M6gPVjAg3mXJh8SqNX99aFcBLz+bauH+c9ms3H8HpHOAoHQy8tLoVB0/SmTyXx8fG7fjSTJBQsWAAC2bt3a/dF63jSKYRgMhDZE0zSfz7d1KpwWHwD/bndQqUwCAZwU3Gbg9W9bFsx/Cwx3CgkJwXG8tLTU/Gd+fn5sbOwt+1AUtXjxYqlUumvXLiaT2feTQhAEQZBFWCAQ8vn8xx57bM2aNS0tLTt37szPz3/ssccAAEVFRZmZmeZ9Vq5cefz48WeeeebcuXPHjh27du1a388LQRAEQX1nmeETmzZtWrVq1ciRI/38/P78808PDw8AAE3TBEGYd1Cr1QMGDPjss8/Mf06dOjUuLs4ip4YgCIKgvrDA8AnLgsMnHAjsPm5bMP9tC+a/bVkw/51rSkQIgiAIuk9w0kCnJdPLSzqulXSU1SkbWtRtEp3MRBF6Qu/GFLiz3by5XhGisChRWIrvYE+OyNaJhSDXojZq8tqKyiWV1fLaZnWbQq9UGpU4ymBhTHeWmx/fJ8jNf4BXTIJ4gC/vtoE1kKXBQOhsOrWSo7WnzzZcrFc2DhTHxnvFzojK8OP5eHI9mCiDjbMVBqXCoGrVtFXLav9uvLw591tvrtfYoJEZ4ek+PLGtkw9BzkxpUB2tPX2y7myNvG6Q94A4z5jMyMmBAn93lsCd5WaiCCNplOkVreq2WmXDuYZLW/K2u7PcxgYNfyBkbIh7kK2T77TgM0LnkdtauKciu7CtdGzwiPHBo5J8BuHovXOGoqlrkoqj10+fqDsb4xn5aNxDyb6DenhG+IzEtmD+29Z95X+V7PrOq39cbM4Z4T9kUvj4JO8EBnbvpXgomi6TVJyuP3+09lSQIGBm9JTxwaNQBD7SAsCi1z8MhM7gfNPlH0t+05p0c2JnPBA6joP3ZpC1iTQdqzvz36t/sHH2ssQFqX6D7/kWWBDbFsx/2+ph/lfKarYW/Fgtr50TO2N65GQeozcLeRAUea7x4u7y/VKdbH581qTwNAxx9XIPBsIbYCC82lm+Je87HaF/POGR0YHD0T5P+EQD+kz9hW8Lf/TmiZ9KXhIuDOlmZ1gQ2xbMf9u6Z/53aDu/zv8hr63w8YRHpkZMZKAWeBRV2F76ffF/O7SS5UmPjw4c1vcDOi4YCG9w5UAo1cm25G0vbC9Zkjh/Ulha30PgzUia3F91ZHvRL5PC0hYPeuxuVUxYENsWzH/b6ib/SZrcVbbv59Jds6KnPjrgod410nTjcnPel/nbhWz3VUNWBLsFWvbgjgIGwhtcMxBSNL234sD3xf/NjJq0IP5hNt5fq8bJDYqv8nbktRW/MvyZFN/E23eABbFtwfy3rbvlf7W89r3zn4rYwueHLA8Q+PXT2Sma2lNx4IfiXzMjJz6e8AgTc7mpK2EgvMEFA2G9sumDi5+hCPrSsKeD3XqUS32U21q48eLm4f4pK5MXs//3xhYWxLYF89+2bs9/iqZ+Lt21u3zfiqTFGeHpVkiDVCfbnPttlazm5eHPJohda7ouGAhvsFQgJA2UUWEyaUjSQJEGEgCAMlCMibJEDJaIgaB2sdIKDehdZft+Kvl90aBHZkZNtWxbaPc0Ju3mnK2lnWVvjP53lCi8azssiG0L5r9t3ZL/bZqOd85/zMQYq4c/58X1tGZKzjZc/DTn6wdCxixNnG8/VUOD3KSXGCkTbS5UMSaKMlEGH2fycZxngdoLDIQ39CIQmtSEts2g6zDqOw26DqNeYjTITDRJM91xBh9HmSjOxgAApJGijJReajSpCK4/WxTNF8UJ3MJ6093LItq1ne9d+NREEq+NfN6f72uTNByvO7s555sFAx+eHTPdvAUWxLYF89+2bs7/sw0XN13+Ym7cg3PjHrTmTWoXpUG16fKWOmXj2pGrIm+6W7UybYteelUlK1er6nU4C2V7MTEWirFuFKqkgTKpCZOKoEiaLWKwPZlsLyZHzOKImRwxiyVkgPvJORgIb+gmEFIEbVSYDDKTXmbUS0z6TqOu06BvNwIMcMQsrjeLI2ayvVhsLyZbxMC5dw2lFEGr6rTyCrWkSIlgiP8YT+9UIYJZ9UI/Vf/3p1e+mhM789EBD9l2CFGLuu3Nsxv9+D4vD3+Gx+DCgti2YP7bljn/CYr8uuD7M/Xn1415Oc4z2rZJOnL91Ja8bY/FZ82JnYHcV1TpIxpISpTNZyX6TqNXopswmu8WzsNYdy2sSANlkBp1EqO+06jrMOo6DLp2A6GnOF5MjpjJ9rxRMrNEDKaQgTHvfBwYCG9oLZCoKvQIgtAUTeop0kgROpLQkiYVQZkophvD3LZ5477Di8nxYvW+Sk4DeaW66VSnQWaKfDjAOrVDPWHYnLu1oK3k9VEvxnpGWeGM92QiTZtzv81rK1o/drUXKoIFsQ3BQGhbKpXKhBNvnN3IZ3JfG7nKjWkX30WLum3935t4DO5rI58XsYVWOKOu3VD1ezNppALTvTwT3Hr9LIk0ULoOg77TqOu80VZnkJuMchMAgCHAGTwM52IYC8NYqP8YT14AGwbCGyTVckMrgaIoQBCcg6IMFOdiDC7GEOA4p7960EiKlDV7WzwHuYXN8O3Xx4fX5XXrzn0Q4xm1asgKi3e/7iPzjefygY9PiZlg67S4LhgIbetyXd4HeZ/PiJq8YODDVq1+3QtJkzuK/nug5thrI56/Y39vC2o+I2k41hE8Sew3yrOf8oA0UEYlQWhJQkeSBpI0UKIYPtOdAQPhDbbqNUroyfIfGgCCxC4M6qb63xf7qg5/W/DTyuTFk63S96wXqmTX15zakB46ZtngBXDOJ5uAgdCGzL/Q10Y+P8w/xdZpubP8tuIN5z+ZHJ72xKDH+mMaGpqia/a0KK9rBywJYYnuPV2cxcFAeIMNh0/QFF39R4uqTpvwr7BuHjH2gsak3XTpizpl45ujX7bOAIlea5a0fFSwBUfxN0a9xGfybJ0clwMDoU0QFPlZzjeF7SVrUldF+0baOjndkRsUG85/oicMb4x6Ucz1suCRaZK+tqOeJunYx4P7qTJwT3A9QttDUCQyy18Yzb/6XT1FWOxmokJavezgKj6T/+XkD+08CgIABEz+h+nrgtz8/3X4pXplk62TA0H9Tm5QvHB8rUQn/XLyR348H1sn5x6ELPcP0t4c7p/y5KEXLzTlWOy4NKj8rQkAMGBpiK2ioGU5w2ewobBMX7YHo+z7epqyQCz8ozz75ZPrliUueGHov+xnMFD3MAR7JmXZowNmP3t09eXmPFsnB4L6UbW8dsWhlwZ5x68f+xqXwbF1cnoEAci8+Ky3x7z6yZUvv8zbTlBk349Zd7BN22aIWRBkJ2Os+w4Gwr5BQOTDAaSeajjW0ZfDqIzqtWfeO3T9+JbJH6aFjLZU6qxmasSEd8a+tvHiZ79d22vrtEBQvzjbcOHF468/OXjh0sT5Nhkp2BcJ4rhvp3xar2x65ugrLeq2vhyqs1DZUaCIXxpyt1ENjsh5PomtoDgS+3hQ63mpqk7XuyMUd1xbcuB5P573F5M+sNVg+b4bKI79MuOjo7Wn37vwqYk02To5EGQxNKC/L/51c+63H6StSw8ZY+vk9JIbS/Du+DUPhIz71+GXTtad691BjEqi+o/mmPlBDL5TLeoOA6EFMPh4xGz/8p8aSAN1X2+kaOr74v++cfb9VUNWPJWyxCKrtNiQN9dr88T3DaTxuWNrJDqZrZMDQRagJ/Trzn5wuSX3q8kfRXtE2Do5fYIAJCt2+gdp67YV/bzx4md6Qn9/76dB5a9N/qM8BMGO0SzcczAQWoZngptbGLc2+z7aHFo17c8efa2wvfTbKZ+OCEjtv7RZExtnvTn63yMCUlccerFMUtl/J6IIWttmkFeozf+0rQbSeH93IZAjogha125QVGnM37umSd+v33urpv2pI69wGZxPH9jgwRH134msKdojYuuUT2iaXnZwVbmkqudvbLkgJbRk4ARx/6XNVhy7CmJXwh/0y32/0neEiOd378Hvh2tObMnb/mj8Qw/HznK45w3dQwCyYODD4cLQV0+9/a8kC4+D1LUbOguVkiKFts3AEjFvTE5IA6PSpJeaGDxMEMIRhHJF0XxuD74FyCEYZCZZmUp5Xauq0xlkRpaIwXS/MRW+SUXoOgwsEdMzQeA1yJ1v0ZpKflvx239/ND8+q2tyXafBwdmvjnjuZN25V069lRUz47H42fccCkxoyfpD7Qkrw5ymg8zN4DhCS2o+K5FeVQ1cHtrNPnKDYtOlLQ2q5rUjX4gUhVkraf2i+3E8dYqGNWfeHeqX/FTKE30fz6uXGusOtisq1F5J7ub69+0/SL3EqKrTKq9rZdfUNE17JriJB7sLQrj2NOmHJTn3OEJtq6GjQCEpUpg0pCiG7x7BE4RwuT6sW79NGqib9ZIiRUeegu3JDM304QdaIBz+XvbXL1d3vzHqpSSfhLvt4wT536HtfO/Cf/SE4bWRzwcK/LvZs2ZvC0XQkVnd7WNlcED9DfYWCGmSzvuwKnymryjuzl/P6frz/8n5enJY+hODHmNgNpiLwbLueSGqjZp3zm/SmnRvjXmlL9MeNp3ubDjW4T/GM2CcVw/HLWnbDJJCRUeBgjLS3kOE3qlCtqdjjEjpOScoiG9nUhPtufL2K3JCR3oNdvdKdBcEcXpyK0OTdNslWf2Rdq9E97AZvr2eGV9PGD689Hm9svGdsat9eN7d7Okc+U8Dek/5gR3FOx9PmPtgdOYdG6h0Hcaiz6qTX4myqz4yMBDeYG+BEAAgLVXVZrcmvRR5S31Fppd/euXr64r6V4Y/E+8Va6vkWVZPLkSKpr8v/m929dF1o18eKL7vD04aqapfm3QSY9yiYJawN7cOmiZ92xVZR56C58/2HSHyHOhm5cVD+o9zFMQ30EBeqW69IJVXaDwGCnyGitzDeb2oypN6quKXRpOWjH08iCm471K7SdXy+pn3ojzCXxi6knWvsbzOlP+NquaNFzfTNPXy8Gdvn8rj2vZ6QQgnMN2+ng7CQHiDHQZCAEDR5hr/MZ5eg93Nf9KAPlB9bGvBD9MiJj6e8IijjJTviZ5fiBeacjZe/GxefFZW7PSeT09MGqiSr65zfdgRWf4o3qfoRRG0tETZekGqbTP4DBX5jvCwyeyIluUcBTGhIduuyFrOS3EW6jvCQ5wsxNh968RHg4ZjHW2XZAlPh93XzdPZhgubLm9ZPOixmVFTerK/c+R/F4qm/6w8sL1o5+yY6Y/Fz+7qxK5p1pdurUtdE93H36DFwUB4g30GQulVVf2h9sEvRAAArivqP778JUGZXhz6lKM/EbzdfV2IrZr2N86+78vzfmX4szzGvRexogj66rd1bA9G5JwACz7k07UbWs5L23PlbqFcv1Geohi+4z5BdPSCWFWva/lbIi1ReQwU+I30FIRYsqtL0+nOtouyhKfDGT1Yec1EEd/kf3+m4cJbY17p+Xpnjp7/d9Sm6fhPzteNqpYXhqwY7JMAAKj4uZHrxw5Mt+RUpRYBA+EN9hkIAQ3yPqz0z/TYrdt7tPbU4kGPTY/McLKuoWb3eyGaSNOW/O0Xm3LeGP3SPZYwpUHZDw0A0P00jRNlpDryFS1/Swkd6TvSw2eoqCfFpb1x0IKYNFIdeYrW81JCT/qN8PAZKur9KqHdqjvQJqtQD1oZhnY7B0qLuu2tcx96ckSvjHj2vtYUdND874mzDRc/z/023itmWfjCui+lQ9bE9LWa3g9gILzBPgMhRVNHs//uzFO0TK5bNniBO8vN1inqL727EM82XPj4ylcPx86cG/fg3e4Pms9IOgsUCU+F9ffzPFW9rvW8VFKsFMXxfUd49O65lK04XEGsada3XpB25Cvcw3m+Iz2sUB2v+LkRZaKRc+7a1/F43dnNOVvnx8+ZHZt5v2sKOlz+3xcDafy5dJfqsDHCK2TCghFsO1sSFcBA2MUOA+G5xkvfFvzoxfKc/ffDCcsi+AF2d/VYUK8vxHZt54a/P0ZR9LURz9++OoymWV/yVW3i8+FsDys9TyV0ZHuOvPWClCaBzzCR9xBhL/pZWJ+jFMSEnuzMV7RdkhmVhM8wkc8wUe/6PfUCaaDyP6oKm+HrmXDr/ajGpP30ytfl0qrXR70YJQrvxcEdJf97jdCQl98tP5d+8rIyd1HC3CkRE/pjXcNeg4HwBrsKhLmthdsKfzaQhmWDFwz3T2082alrN0TNtfellPqiLxciRdO/XN29q+zPlclLJoWN/2c7QRd8XB2Y7uWd2vvhFr2mqtO2XpRJipRuYVzvISKPeIG9dRC4mZ0XxDRFKyo1bTly2VWVMJrnM8waVcDbqeq0176rH/xiJNPtn5ub/Lbi9y/8Z5h/ysrkJ9g4q5dHtu/877umU53aFkPUowFlkspvCn5o13YuHvRYWvBoO3nQAwPhDXYSCHNaCr4v+VVhUD4+cG5ayBjzVWJSE7nvVQ55I8Y51uu6o75fiJWymg3nPwkS+L8w9F/mgYZ1h9p17YbYhUEWSmNvkEZKUqRsz5Grm3SeA93ESe7ukTw7nFDDbgtiVZ2uI1/eWaBgCRniFKF3srCfngL2UP3hdm2LPnZRMABATxi2Fv5wuv7Cy8OeHuqf3JfD2m3+W0rexsrIuQFuoTe6tuW2Fm4v+kVlVC8cODctZPQ9J6PpbzAQ3mDbQEjS5Jn6C79c3W2iiPnxWekhY2+5USrbUS+MFfgOd5IpCm9nkQvRRJp2FO/Mrj62MvmJcaJRBZuqBr8YabWms+4ZlURHvrwzX6GXmTwT3LwGublH8OxnGKJ9FcQ0UNVrJUXKzkIlgiPiJHdxsjtH3MvKlmVRBJ33fmXUowE1/JoPLm4eKI57JnXpffWLuSP7yn9LU17XVv3WlPzKrX1oc1oKfij5tUMrmRs3KyP8gV7Xp/sOBsIbbBUIlQZVdvXRPRXZfjyfuQMeHBGQesfH7LJrqvoj7YnPOfaM9d2w4IVYIa3eePGziWWTEyJj42f05oFNv9JLjZ2FSkmRUtdpEMUKPOMFwhg+zrFxU4Q9FMSUkVJUayQlKulVFc5BPQe5ew1y4/nb3aPxhiutpYeub0/YtmroiuH+lpnj3h7yv/9U7mzi+rMCxt151ERpZ9nOq3uK269OjZgwK3pK91Pw9BMYCG+wciCkAV3UfnVf5eELzVfGBI2YHZPZ/TN2mqJz3qmIXxbirBNAW7YgkNdqCrdVfpbwnwcHTJ0bN8s+p6AzKgnpVZW0VKmo1vD82KIYvjCazw/m2KTh1GYFMQ00rXp5hVperlZe1/KDOB7xAs94N7aXPU4WQQP6UM2Jrfk/Li9fHv9AROAwixXZThwIST11ZX15yup7zKnWom77oyL7UM3xeK+Y6ZEZwwNSrNmbBgbCG6wWCBtVzcdqTx+uOcnCmNMiJ00OT+thu0rdoXZST4bP8uvvFNqEZQuC4i3XvYcI6QGmz3O/vS6vfyZ1qaXu3PsDRdDKao2sXC2vUBtkJrcwrnsEzy2Myw/iWK3t1JoFMU3R2laDskarqNEoqjUYCxVG80XRfGE03w5HmHUpl1R9mvM1AOD5Icv9VQFl39enrom21BfkxIGw9YJUXqGOfTy4JzsbSOOJurPZVUcaVS0TQ8dNCBsX4xHZ3ykEMBB26e9AWK9sOtd48WTdOYlOOj549OTwtPv9gvVSY+F/aoa+GWOHXS36zoIXoqpOV/5jQ8prUeaMutyStzlnqy/PZ2Xy4jBhiEVO0X9MakJRo1VWaxQ1Gn2nkefPFoRwBcEcfhCH7cHsv36S/V0QGxUmdaNeVa9V1evUdTqGG26O9+6RPDt5iNuNDm3n1sKfclrylw1emBGebn54UfpNrddgd5+hlnls78SBsGhzTdAE8d0WD7ibRlXzkesnj9WeQQAyPmTU2MAR0Z4R9zs6s+dgILyhPwKhkTQWd1y71Jx7oSlHR+hHBw4bFzwy0Tu+112kCj+tDpnmI4ziWzCRdsKCF+K1HfXCSJ7faM+uLQRF/ll54MeS30YFDluc8KgX17Obt9sP0kCp6nVqc/Bo0JEGiufP5vmzuX5srg+L68uy4JNFyxbEpIHStRu0bQZti17drNc06wEN+IFsfjBXEMwRhHAdZfIdjUm78+off1YenBGVMW9AFpfxz8xt8kpNzR/NyS9HWaRwdtZAaFSY8j+qGrouttdV53Jp1en682cbLuoJ/TD/lGH+KSm+iTd/ERYBA+ENlgqEBtJ4TVJR1F5a0FZyTVIRLgwxf3nRHha4nWk62amTGO1qHS9LsdSFqOswFG2+PmRt9O1TYamNml+u7t5XdXhaxMRHBzzkcNP0mDSkplmnadZrWw3aVr2uzYhggCNmscVMtieT7cFkiRgsEYPpxujFgMXe5T9N0SYVoZeaDHKTQWrUS4y6TqOuw0hoSa4Pi+PD4vmyzJHb/qt9t9AThr2VB/579Y8RAUMWD3rM+7a5GgAABZ9UB0/y9oi3wHXrrIGw+YxE06K3yBjoBmXTheacS825VzvLw9xDknwTEsRxCeIBPZlt+J5gILyh14FQTxiuK+qqZNcrpNVlksp6ZWO4MDRBHDfYJyHRO94iX9I/53Le1lFLXYhVvzcz3fDgyXftxdCpk/5Y8tuJurPTIyfPiZ0pYrv3/aS2YlQRunaDXmLUS4wGqUkvNRrkJqOSwDkYU4AzzP+4GM7DcA6GczCMhWIsFGdjAAXm2iSCIebBqWq1ms/nAwAoE02ZKAAAaaBoiiZ0JGWkST1J6ClCR5rUBKElTRrSpDQZVYRJQzJ4GMuDyRIy2B4MtieT7cnkiFksIcOBppe7hZ7Q/1V5aOe1PYne8YsTHg1xv+s41M4CRfNZyaBnLNAz2VkDYdHmmqCJ3qJYSzZiGUljaWd5fltxUXtpubTKm+sV6xkV7RERIQqLFIbxmbxeHNOC+e8A80j1kcKg7NBK2rUdLeq2RlVLk6qlXtko1clC3IMiRWFRoogp4RMiRaH9tzoS24PJFjEU1VphVG++bKdnUhOdhYqU1d3Nwe3F8Vg1ZMW8+KyfS3ct3LdyYtj4uXGzfHj2tTpaDzEFOFOAu0fcejGYVIRRTZhUhFFFEBrSpCV1HUZST5J6ijRShI4EsgzA/gAAIABJREFUNCB0JACAJmnSQAEAaJpGEAQAgDIQlIECADAWiqAIzkFRJoqxMJyN4lyMJWLwAzgMPsZwYzAFOIOPOdM9mdKo2lN+YE/F/sE+CR+nv33PJ8qeg9xqs9tUdTrLLnbhNIwKk67dYPHCiokxk3wSknwSAAAkTdYqGsokleWSqhN1Z2vkdVycE+IeFCDwCxD4+fN9vble3lwvIVtotSlsHLtGmNNYUCavRFFUa9KRFKkn9XrCoDXpVEa1wqBSGJRyvZzD4Ig5nj48sQ/PO0DgFyjwD3EP9OV5W3NahMaTnXpnbB21yB1Z0+lObYsh6pGetsNIdbLfyv7Mrj46xC9pbtws6/RPs0/OWiPpoWZ16+9lfx67fmZ00LBHB8y+fTnZu2k82anvMEQ+3NemP6fMfwu2i/Zcu7azXtHYqGpuUrW0aNraNZ3t2k6VUeXOcnNnubmx3NyYfC6Dw8JYPAYXQZCpERMCBf6wRniDgTCojRoEQbgMDsbAvDBPNs7iMbh8Bs+d7SZkuQnZwq7lJW3IK9Gt8NOaiIf8nOlO3FLaLski59zHr86DI1qRtGjBwIf3Vx15/cz7XhyP2TGZY4NH2sMXDVkBRdM5Lfl/VOy/1lmZGTlxR+bnnpz76wXqkyrMfb8ybKafE09/2GudhYqgidYeHW+uAqb6Db55I0mTcr1SblAoDSqFQakj9AbCoCV0/ZEAxy47kn0GDQ9Ktflco/dk7hOhvK69vUHMxanqtDQFuiYz7Dkegzs3btac2Bl/N17eU5G9OffbjPD0zMhJgQJnq3ZDXSQ62cGa49lVR/hM3oPR094a8yqrV080GALcPZLXWaDwGea00x/2jlFFaNss3y7aOxiCeXJE93uX0zuOHQgdiMcAgeyaCgbCW7RdkvkME/W6jwaKoGOCho8JGt6oat5fdeSZo6sD+L4Z4Q+MDx7Vu8fvkB0yksbzTVcO1Rwv6SgbFzxy3eiXYzz72h7uM0zUeLwDBsJbyK6qhNF8+5lN12pgILQSUayg8tem0Exbp8OekEaqs0iZ/PKts/r2QqDAf0XSomWDF1xqzj1Uc2JL3nfJvokTQscO90+14aTAUF8QFJnXVnii9uy5xkvRHhEZ4enrRr9iqW9TFMuv+q1J22bg+sDL4x+yMpXHAAcboWQRMBBaiSCYY1ITBpmJJXKwsVn9p7NQ4R7Ou3mVuD7CEGxkwNCRAUM1Ju2Z+vPZ1Uc/uLh5iF/S2KARwwNSLTsqBuonRtKY01pwtuHi342Xg9z8xwePXjZ4ocXbxxAU8R4iar8sC53ua9kjOy6apOUVmoiHXPHhAgyE1oIAUSxfdk3lO9LD1kmxF535/fWQhsfgTomYMCVigtKg+rvx0rHaM5sub4n1jBoekDrcP7Xn3Qshq+nQdpqnc8pvK472iBgTNPxuI+ItxTvZvXRrXWimr+OOnrQs5XUtW8xkCFwxKLjiZ7YVUaygI18OA6EZoSNVdTrzWqn9x40lMEdEPaHPbS260HRlV9lfCECG+CWl+CYm+w5yuKlqnImO0Be0leS1Fl5pLZDqZEP8ktJCRr8y4tm+rxTYE1w/NsZCVQ06QTAcUAiAuV30PicXdRowEFqPKI5fvauZIuheTKblfCRFSmE0D7ttTrV+wsbZowKHjgocCgCoUzRcaSk4cv3Uh5c+9+GJB/sMHCSOTxDHOcp0pg5NaVCVdF4rar9a2F56XV4X5xWd4pP46vBnoz0irTZ6uotnorukUAEDoZn0qrrnw3mdDAyE1oNzMK4fS1mtEcY44QTc96uzSOmdKrTJqUPcg0Lcg7Jip1M0VSGtLmwvPVp76tMrX7Nw1kCvmDivmDjP6EhRGOxlYxEERdbIa8sklVc7y69KKjq1kjiv6EHiASuSFsV5RvXfjE494TXI7dp39bB1FABgkJtMakIQ5KL3BDAQWpUoTiC9poKBkNCRyuua2IV3nRDSOlAEjfWMivWMmhs3CwDQoGy6Jqm42llxvPbMdUW9H98nShQWJQqPEIWFC0MdeoJTa9KYtDXy2mpZbaWsplJaU6dsDOD7xnhGxotjs2JnhAtDrDmpU/d4/mwEQ9SNOr6rBoAusmsqUSzfZW8IYCC0KmEUv+r3JlunwvakJSphFN/e5vUIcgsIcguYFJYGACAo8rqirlJaUymr+bvpSo2sFkXQcFFIsFtgiFtQsFtAoJu/N1ds/dY8e9OpkzaqmhuUTfXKplp5fa2iXm3ShLoHRwhDoz0ipkZMiBDadd3aK9Gts1AJA6G8UuOyDwgBDIRWxg9iG2Qmk5pg8F065zuLFOLBdl3BwlEsShQeJfpnjQKJTlarqK9VNNQq6s82XmhUNssNSn++jz/fz1/g48vz8eV5+/K9vbleTtkBR2PStms6WjTtbZr2VnV7s7q1Wd3WpGpm46xAQUCIe2CQW0Cyz6AwYbAPT9x/a7FanGeie9mO+tBMH1snxKZooKjShLnwSBKXLo6tD0ERtzCuolrjlWjXYaBfUSZKUa2Jnhdo64TcH/NsTym+iV1b9IShWd1q/teibstvK2rTdLZrOgykQcz18uJ4iLleIra7F9dTyHITst2FLHch292dJWDjbBt+kDsykkalQSU3KGV6udygkOuVEp1UopNJdFKJTtqu7QQA+HDFPjxvH57Yj+8T5xXtz/cNEPg5+uhMfgCbpmgXH1mvbdVjbNSVhzjDQGht7lE8RZVLB0JFlYYfwMHZ9j5D7D2xcVa4MCT8tnV/9IShXdsp1Uk7dBKZTt6pk1bLa2V6uUwnVxiUSqOKomk3Jl/A5POZfD6Ty2PwODibz+RxcDYTY/KZPCbKYOEsJsZkYUwMQTn/v7S3gPk/T5c1Wo0K0dy8RWvSkjRlTgNBEQRF6Ai9iTLpCYPOpDOQRq1JpzFpdYRObdRqTVqVUa00qlVGFUlR7iyBO8tNyHb3YIuEbDcPjijMPdiT4+HF9RRzPR094HVDFCeQlalcORDKqzTCKJfuuAADobUJI/llFxpsnQpbkl5Ti+Kc+VfHxlnBbgHdDNs3kEaVQaUyqlVGtdqk1Zq0OkKvMqr1hF6uVzSpWgyk0UgaDaTBSJpImtKZdAAAGtBq4/+EPYqiUPR/nrNyGVwMQQEALJzFQHEMxbg4h4Ey2DiLg7OZOFPA4vvyvTk4m8vg8hlcPpMvYPLdWAKO/VVSrcYjTtB8VhIwrh8H79s5RZXGy74fVfQ3GAit7f/au/PoNsr7XeCjdTZJXmLLkuXYDolDSIAkJARCuAQSlnAgbaCsaaAltAQCPaw9pxRuoKU9QOGy9gIlTeDQhSWUnBJyf4VyIJA2e1iSkMWJ43iVbHnTMjPaRnP/EHWNHduyNdJsz+cvWR5r3igaPTPzft/3ZSupFJdKhJL2IoPeiOg9HDntlvyOo1c50mInmQm5D1vU5Xp4hVdUxx75S4sYT6uteqtAJCJ8nJusu9VSx8SQ//HKMhGuyWyogRt9Sz0SOuNSSmI9xr3+ALWx2M3OaqbvaFTphigj2ibYHFa7IWdW64cgVEDxFDZ0zKBB2Hs4WnKaUztFhWAIJac5eg8ZNAhDR7kiY3cQEghCRRTVOfqOGjQIezLjdgHUpPQ0Z+/hiNKtUEaogSuaYvTFOxGECmDcpBhPx/uSSjek0MREOtLEF09FEIK60G7SZDbxgbjSDSk4iQif4ItO0W1JcJYQhEowEc4aOtLEK92OQgs3cI4q2qAlCaBuJac6eo8Y7qKQD8RsrMXg83sQCEKluCax4UbDBWHoGFds+JswoE6ZAb5Kt6LQwo28axIOSQShQlyTGAMGYd8xdMuDShVNZsPHeSktKd2Qggo38s5ao98XJRCESnFMpIXOuBhPK92QwknFRKEzjtmNQZ1sDqu9yMa1xZRuSEGFT/CuSQhCBKFCzFYT66WiLYLSDSmccAPvrGGwKDGoVnEda6hy7mQkleJFxm3cueX6yRaEPM83Nzen08Ne4oii2NTUFIsZ64RrBE6D3R3tOxotrkNvBKhX0RQ2dMxAowlDjbxrEoNBvYRcQfjCCy/4fL5FixZNmzbt0KFDQzfYs2fPKaeccskll1RWVr7xxhuy7FTrXLVMuNFAp5+hYxiuBKpWNJkNn+Al0SjdhJFGDvdFM2QIwqampl/+8pfbtm07duzY8uXL77777qHbrFq16v7776+vr//oo49Wr17d3d2d+361znUKGzkhGKRzPsWLsZ6EowodhKBeVsZCTbAbp8MClTL9ZAjCt95666KLLjrttNMIgli9evUnn3zS0dExcINDhw4dPHjwJz/5CUEQc+fOnTlz5saNG3Pfr9bZWIvNaeE7DDGGN3SMc01iTRbchQFVK57C9hljEEU6meYDcSeK1wiCkGX1icbGxqlTp2Yeu91up9PZ1NRUUVExcAOfz8cw35561NXVnThxYrhXE0Vx3759gUAg86PX662sHHZadCnOJzubRbNWS37Y8lTvl022lFbbL/J8gsnqjLLnq5RjginRcjTfTTKU7N9/yBJbnA58KVac2pfNxpp+/yPtElUqpToalG7ImJlMZqu3xmSRcxIAGV4rEomUlf13KS+WZUOh0KANaJoeYYOBQqHQz3/+c5vt2yWK5syZ8+yzzw63cfTf/0N8s238TVccf3pPY5nlwBal2zFO6XRayO4spK/7unLn593HAvlukqFk//5DltISGe28qevN/2siRu+z0PT738edaRGLut/cqnRDxoO69IfWyWdGo1lVNlEU1R8ow5EhCN1ud29vb/+Pvb29Ay8HMxv09fUN3GD69OnDvVppaemHH37o8w27qOlA1kU/sC+50WLR6lrnzmbh6IY27/3XKd2QccpyPTwxkT6x5nD1L/43xk7IC+sR5kPHk0eLfvh/WN/oK4Vp+v0P/7nFM83pnnuj0g3JiVzvvwynMzNnzty1a1fm8YEDB8xm8+TJkwduMH369M7Ozra2tsyPu3btmjlzZu771QHWR8WCCTGh82H1kSae9VFIQdAEZ60hxjVFmgRnNToIvyVDEF577bXNzc1PP/30119/fc8999xyyy0syxIE8Ytf/OKJJ54gCKKiouKaa66566679u3b98gjj0iSdPnll+e+Xx0wWUyMh9T9ZBaRRsxeAZrhmsSET+g8CFO8mOJFuhxD6b8lQxCyLPvxxx9v27Zt1apVc+fO/d3vfpd5vqqqyuv1Zh6//PLL1dXVP/3pT+vr6z/88EOr1eiTnfdzVDORZp0fdeETqNIGzXDVMhG9B2GkWXBU0RhK30+eQDrjjDPee++9QU/edddd/Y9dLtfzzz8vy750xllN63xFUImINAtTlyMIQRvoclJMpON9SbJ4lAoL7Yo28w7cFx1AqyVPuuGYSEea9TyAlwvEbKwVC56BZhhgudBIi4ARhAMhCBXGuMkUJ6Y4UemG5EvkBO/CfVHQFN2vkhZtEXBFOBCCUGkmwlFFRfQ7q1O4kXeiUgY0Rd/dhPHeJCEROr7xOw4IQuU5JjI6nt4wjCtC0BrHRJoPxNM6HdcUweXgEAhC5Tmqab0WjqY4MRUVmQpUaYOWmG1mxkNGdTquKYoOwiEQhMpzVFF6HUoYaREcE1GlDdrjqKb1ep8mmjkqYQAEofKoEruYSCcjKaUbIr8IqrRBm5wTGb2Wc3NtAuvDUfkdCEIVMBEOH6XL+zDRZkzjBJqk1w6LeG+SMJvsLgxn+g4EoSqwPppr0+HpZ7QVN2FAk/Q6rinaKmB97KEQhKrg8FHRVr1dEaJKGzTMRLA+Ktqqt9NTri3mqBp9YQ2jQRCqAltFR3V3RRhpFpw1GDgBWuWs1mE3YbRNcKCDcAgEoSowbjIZSaViuroPE21BpQxomC67CaOtMRZXhEMgCNXBRLBevQ2iiKBSBrTMWU1H9XVFmIym0ok0VWJXuiGqgyBUC1ZnowklgmuN4SYMaBdZbCNMRLwvqXRDZBNtjTmqKIzrHQpBqBYOH62nnnm+M25zWqysRemGAIyfY6KuLgq5NoFFyejJIAjVgq2i9VQ4Gm1BlTZonmMiracJ8aNtMUclOghPAkGoFqyHjPUk0kmdzPMbbcW5J2ieo0pXA3xxRTgcBKFamCwmuszOB+JKN0QeGK4EOuDwUbqZcVRMpBOhFF2OSpmTQBCqCFtJce26uDsqEVwbKmVA8+xFNpPZpI96Gb49RntIkxmlMieBIFQR3QSh0BW3MhYrg0oZ0Dy2itZHOTfXjg7CYSEIVYStpDi/Lg65NgzaBZ1w6GWiNa49xnhxVJ4cglBFdHNFGMUIQtALVi/jmrj2GIsrwmEgCFXE5rCaLXrokIi2CqiUAX3QybrZEsEFEITDQhCqC+vTw0UhVv4E3aBK9bBudqwnYaUtVhrd9ieHIFQX1qv5IMTKn6Arulg3m2uPseggHB6CUF100E2IlT9BZ3SwbjY6CEeGIFQXRvtBiKH0oDM6WDcbQTgyBKG6MG4y3pfU9ERrUXQQgr7oYN1szo8gHAmCUF1MFhNdru2J1rj2mMOHQw70gy63J8MpMabV01Mxnk6GU1QZJlcbFoJQdTRdL5MSxBQvUqU45EA/TGYT7SH5gFaPSt4foyswudpIEISqw3op7R5y3xan4YgDfdH06Snui44KQag6jJfi/Fq9NYo+edAltpKKajYI+UCc8ZBKt0LVEISqw3hIXrMzjnLtMRYdhKA7mh7XxPkxiHAUCELVIYttaVFKRjU5kwXXhitC0CG2kuL9MUJSuh3jwgfimG57ZAhCNWI8lBYLR6W0JHTGGQ8OOdAbK22xMpZYd0LphoxZIpIiJMnuxExPI0EQqhHrIbW4HpMQTNiLrBYSHyrQIY3eHeX9WH1pdPjOUiPGq8krQlTKgI6xPlqTQRiIo4NwVAhCNWK9mqyX4dpjbCXmlAF90ugVIeePoWR0VAhCNWK8FKfBnnlcEYKOObQZhLg1mg0EoRpZaYuFsmhuhV6uTUAQgl5RE+xJLpUSRKUbMhYSwXfEmQpcEY4CQahSrFdj9TIpTkwnJbLYpnRDAPLDRDBeSlt9FrGehJXBeryjQxCqlOYOOS5zBwaTq4F+sV6K01QVG4bSZwlBqFKsR4uHHO7AgJ6xXpLXVDchhtJnCUGoUozWCkf5AIbSg84xHkpbHRZ8ACWjWUEQqhTjJoVgQkprpnKUa0elDOgcW0nx/riGyrl5f5zF6WkWEIQqZbab7S6rZqZ0yhSn4dwTdM3KWMx2UzykjXJuKS0J3QnajcVBR4cgVC/GQ2plfplYT8JKozgN9O/b2be1INaVIIusZhu+5EeH90i9NDT1NorTwCAY7azQy/nRbZ8tBKF6sR5SK0vV834Up4EhsNpZNxuVMtlDEKqX1q4IcciB/mlogC8KubOHIFQvuoIUurVROIr5DMEgmMxRKWrhqMQVYdYQhOpltprIIpsQVHvhaDolxXqTjBuHHOif2Woii21Cp9pv1UiiFOtJ0uU4KrOCIFQ1TRSOCh1xqtRusmB2NTAETXQTCsEEWWIzW3FUZgVBqGqMFupluAA6CMFAGK8Gjko+EGNxXzRrCEJV00S9DEpGwVBYLUy0xmGW0bFAEKqaJq4I0ScPhsJ4NdBhgaNyTBCEqsa4yVhPUuUlahhND4ZCTbAnIykxnla6ISPB2IkxQRCqmsliokpVXaImxtMpXqRKMZ8hGIXJbKLdJN+h3qMynZLivUm6DEdlthCEasdUqPqQ4/0xuoLEerxgKIxH1aukCZ1xaoINhdzZQxCqHe2h1ByEXADrvIDhqLyKje+IMxU4KscAQah2TIWqe+bRJw8GxHpJvkO9V4RYE22sEIRqp/LCUfTJgwExHlWPqcfp6VghCNVO5YWjPKbbBuMhi23pRFoUVFo4ygdwa3RsEIRqZ7KYyBKb0KXGGUdTfDqdkuxFNqUbAlBYJoKpIGOdajwqJVGK9ybpcpSMjgGCUANUe3c01pHA7BVgTIyXEjqTSrfiJPjOOFWKktGxQRBqgGrrZeLBJOYzBGNiPGQ8qMYgFFApM3YIQg1gKlQ6gkLoTKJSBoyJ8VCxDjUGIerXxgFBqAGqvjWKc08wJNZLCqrsI+Q7YkwFjsqxkS0IBUHYs2dPc3PzcBt0dXXt3bu3tbVVrj0aB+0mY91qLByNBZMIQjAmm8NqMpkSkZTSDRmMD8QRhGMlTxDu3bt38uTJ99xzz9lnn/3ggw8O3eCGG26YOnXq6tWrZ8+efdVVVyUSajyTUq1vF8VWWeFoIpQ0W0w2h1XphgAog3Lb1DbRWmZhegoL04+RPEH4wAMP3H333f/617+++OKLl19++eDBg4M2WLlyZSAQ2Llz5/Hjxw8cOPDaa6/Jsl/jUOFS9Zw/TlWgRBuMiyy3qe2oFDrjVCkWph8zGYKwo6Pjs88+u/XWWwmC8Pl8S5YseeeddwZtc+mll9rtdoIgnE7njBkzOjo6ct+voaiwm5DviJFluBwE46IrbGpboZfvwH3R8ZDhi6ylpcXpdJaVlWV+nDRp0gg9hQ0NDZ999tmjjz463AbJZHLr1q39r1ZZWTl9+vTcG6l1TAXVfSCsdCu+g/fHKQ+uCMG4KLc9tF9lRyVKRscl2yB8/PHH9+/fP+jJuXPn3nfffTzPZ672MiiK4jjupC/S09Nz9dVXP/DAA7NmzRpuRxzHvfjiiyT57UnNzJkzH3vsseE2FgTBbrdbLJYs/xUa5hKj7UI0GlW6Hf8VaeNK62hVNcloOI4zmXATTDFpR5JrF6KRqHqWIQu3RYtnsAY5KrP8/FMUZbWOknTZBuGCBQsmT5486Emfz0cQREVFRV9fXzqdNpvNBEF0d3d7PJ6hrxAKhZYsWXLxxRc/9NBDI+youLj4nXfeybzyqCwWi0GCkKmVjvYGWJpVy4QREhHvShXXuBwOh9JNMS5JkvD+K0iSJCsdtqVIskQtswwmuvylNUWMwxAXhTJ+/rMNwgsuuGC4X9XW1hYXF+/cuXP+/PkEQfz73/++9957B23Dcdz3vve9WbNmPf300+Nuq5H1F46qpAMg1pOwMhYLiXGoYGisl+QDMZUEIUpGx02GLzKSJFevXn3nnXf+4x//ePDBB4PB4DXXXEMQxLZt26qrqzPbXHPNNc3NzXPmzFm7du2rr776+eef575fo1FV4Si6IgAIgmAqVLQeE0pGx02eqr81a9ZMmDDhpZde8vl8n332GUVRBEFUVFRcf/31mQ0WLlw4c+bMxsbGzI/9tTCQPaYisxaoS+mGEARB8AGsvgRAMB4y1HDykojCQ8nouJkkSV3zlVRVVe3cuTPLPkIDFcsQRPCLUPeB8LSbJyrdEIIgiCN/aS051UGfanE6nUq3xbgikQjefwVFIhGi19rwbvus+wbXTyii+R+dBEFUL3Er3ZACkfHzjz4ezVDXrVE/lsAGIBgPyXfGpbQqLif4DhyV44Qg1AzaTca6E2qYcVRKS0IwQbtxyIHRWexmu8Ma71HFMhSYZXTcEISaYbaqZan6WFfCXmS12PHhASAYL6mG+WUyJaM4PR0ffJdpiUpW6OX8cRYL0wMQBEEQjIdSw1EpYGH6HCAItUQlM47ygRiDIAQgCIIgWI8qrghRMpoLBKGWqGSpet4fY9EnD0AQBEEwXkodp6cY2jt+CEItUckVIYdDDuA/VLJuNkpGc4Eg1BI1HHLplBTvTdLlWHcCgCAIwmw1UaU2vlPhWzWcH7dGxw9BqCXfFo4qesjxgRhdbkefPEA/xkspu1T9t6enKBkdLwShxrBeUtluQnRFAAyi+GQXQmecmoDT0/FDEGoM46E4RQ85zDIKMAjrpZQtHMVMTzlCEGoM41H4JgznxxUhwHcoPpSQC8RRyJ0LBKHGKF44yvtjDK4IAQagy+zJaEqMp5VqAIb25ghBqDF0uT0eSqWTyhxyqZiYEkSqBCWjAAOYCLpcyW5C9NznCEGoMSaziS6zK1Uvw/vjjIck0CUP8F2MV7FbNWIinYikqAk2RfauDwhC7VGwQ4IPxDDLKMBQrEex+WX4QJxxkyYzzk/HD0GoPayXVKpeBpUyACfFeCmuXbHTU5SM5ghBqD0KjqDg2mNsJYIQYDC2kuLaBUV2jQ7C3CEItUfBwlEMVwI4KbvLSphMiUiq8LvG0N7cIQi1hyq1p3gxFRMLvN94b9JsM9kc1gLvF0ATWC/JtytwhooOi9whCDXIRNAVpFDwu6OcH/dFAYbFeCmu4EGYEkQxLpLFKBnNCYJQkxSZ0olrR8kowLCUPCpRMZobBKEmKXLI8X7MXgEwLJyeaheCUJPYSgVqtXFrFGAEjJcUgokCLxfK+2MMjsqcIQg1ia2k+PYYUcAjLp2SYt1Y8AxgWGabmSy2CcFEIXfK+XFFKAMEoSZZGYvZbor3JQu2R74jTpfZzVb0RQAMiy1wvYyUGUSI09NcIQi1iq0s6CHHt6ODEGAUrJcsZDdhrDthY61W2lKwPeoVglCr2MqC9sxzfgzaBRgF4y3ocqFcO9ZEkweCUKsKPGgJlTIAo1Lg9BRHpRwQhFrFFvbcE7dGAUZFldpTgpgSCjTrE49KGZkgCLWKqSBjPcnCrNCbCKfSaQmzVwCMwkSwXoprK9AZahSnpzJBEGqVyWKiywu0Qi/XJjh8dAF2BKB1jio62laIZSjERDoZTtHl9gLsS/cQhBpWsJksom0xhw8nngCjYysLdEXI++M01uOVCYJQwxgvVZjZ7rk29MkDZIX1FSgIMZReRghCDXP4qGhBDrlom8BW4dYowOhYLyV0J9KpvE/7xLUJLO7TyARBqGFsFc215X2iNTGGrgiAbJksJrrMXoCK7mhbDEEoFwShhtlYi4U0x3rzO7ch1y4wXgpdEQBZcvjovN+qkQjejyCUDYJQ21gfxbXm95DDiSfAmLA+istz4SjfGbc5rVYKk6vJA0GobQ5f3mu1uXaUjAKMQQHqZTCiSV4IQm0rxCHXKrDHo4Y2AAATHElEQVQ45ACy5vDRnD+/nffRtpijCqenskEQapujio625vGKUBIlPphgsc4LQNYslNnmsArBPE52gdNTeSEItY0ssaVFKRFJ5en1+Y44VWIz2/E5ARiDfA9tiraj515O+ILTPEclzeXtojDaIjgm4sQTYGwcE/N4qybWkzBbTHanNU+vb0AIQs1zVOXx3DPSLDirEYQAY+OYSEeb8xWEXFvMgQkuZIUg1DzWR+evVjvaIjgmMnl6cQC9clTT0TZBSuelYCbaFsNMT/JCEGqeo4qK5mcoYTolCZ1xthKVMgBjY6UsdqdVCOZlsguuVcCIJnkhCDWPLieTXCofa4Fy7TG63G624UMCMGaOiUye7o5GWwXcGpUXvuO0z0Q4qvLSIRFt5h3VuC8KMB7OajrSwsv+svG+JCERZAlWyZYTglAPnDVMpEn+Qw6VMgDjlqd6mUiT4KzB6anMEIR64KymI/m4IsTYCYDxcvgoLhCXfT2maDPvwOmp3BCEeuCsYWQPQjGejvclGcwpAzAuZruZLrPzAZkL2SLNgrMGQSgzBKEe2F1Ws9UU65GzRC3aIrCVWH0JYPxkvzsqpSVUyuQDglAnHNV0pEnOQy6C+6IAuXFW05EWOY9KPhAni21WGqsvyQxBqBPOaibaLGe9TLSZR6UMQC4c1TJXsUWaBXQQ5gOCUCeccl8Rhht51yRWxhcEMBrWS8b7kilOtjG+kSYeJaP5gCDUCUc1zfljkihPiZoQTBAmE8YqAeTCZDY5q5nwCdkuCqMY0ZQfCEKdsNjNVKmd88tTohZu5Iom48QTIFeuU5hwozxBKMbTse4E68XkavJDEOqHs4aOyHTuGTnBu2pxXxQgV65JbPg4J8tLRZoF1keZLCjklh+CUD9ck9lQgzxBGDrOu07BFSFArlw1NOePyTKsPtzAFU3G6WleIAj1o2gyG2rgiJyPuCQnJsMpDKUHyJ3ZbmYqSFkqukMIwrxBEOoHWWyz2M18ZzzH1wkf55yTGAylB5CFaxKbezdhOiVFWwVnLe7T5AWCUFeKJrPhhlw7JMKNvGsSjjcAebgmMaHjuQZhpIlnPKSFxDd2XuBt1RXXFCYkQxByCEIAubhOYSKNfI6r1eO+aF4hCHXl227CHIjxNB+IY6wSgFxsDqvNZeX9OfVZhBt4BGH+IAh1hSq1mywmITj+2bdDxzhnDYNV6QFkVHKqo/dIdNx/LolSpJl34j5N3uD7Tm9yvCjsPRIpmeaQsT0AUDLN0XckMu4/jzQLtJu0UphrO18QhHqTaxAejhafiiAEkJNrMhtpEcR4enx/Hmrgik7BfdE8QhDqTfGpjt7DkfH1zMd6EulEmvVgDicAOVnsZufE8Rey9R6M4PQ0r+QJwnQ6vX379s2bN/f19Y2wWTwe37t3b3d3tyw7hZMii21kkW18K1H0HoqWTHMSGEAIILfiaY7ew+PpJkxxIh+IF03BFWEeyRCEoiheeeWVt91226uvvjpt2rT9+/cPt+WaNWvmzZv3wQcf5L5TGEHpdGfPwfF0SPQeRgchQF6UnDrObsKeQ5GiqazZivPTPJIhCD/44INjx47t2rXr73//+2233fbII4+cdLNdu3Zt3bp19uzZue8RRlY6w9XzTXisfyWJUvg4jxNPgHxgvVRm+Yix/mHPN+HSGc58NAn6yRCEGzduvOqqq2iaJghi+fLlmzZtSqVSg7ZJJBK33377Sy+9ZLGg8CnvnNV0khNjPWM75MKNPF1utzmseWoVgKGZiJJTnWO9OyqJUl89VzoNQZhfMnzrtba2zp07N/O4uro6lUoFAoGqqqqB2/zmN79ZunTprFmzRn21eDz+/vvvl5aWZn6sqqo699xzh9tYFEVRlG31Zz0pmebo2h/ynl+a/Z8Ev+4rnu4Y0/uJ919ZeP+VNdb3v3g66/9Xj/vcouz/JHSUo8rtZsaE/+ihsnz/zWazyTTKjeWsgvDrr7++7777hj6/fv36mpqaRCJhs327lHnmQTz+nTkU9u3bt3Hjxt27d2ezr1gstmnTpsz1JUEQZ5111gh3U+PxuCRJuMocylFHBneGS8/O9j6nlJa6vg6fuqpy0P/dyBKJxJi2B3nh/VfWWN9/5hQbvyEW6eTsRdlegXTtD7lOpfG/fFJZvv92u91qHeUNz+r/o7a2ds2aNUOfLy8vJwjC4/F0dXVlngkGgwRBeL3egZv99re/9fl8jz76KEEQLS0tGzZscDqdV1999Un3VVRUtHbtWp/Pl03DTCaT3W5HEA5FnkE1v9dlN5FWOqs3p/dQhC4jS3xjOFclCEIURYbBbBeKwfuvrHG8/xPOKOIOJ4ovcmW1tUSEj7Sedms1w2BE00nI+PnPKgiLiooWLlw43G8XLFjw/vvvP/TQQwRBfPrpp7Nnzx7UuJUrVzY1NWUekyRZXl5eVlaWQ5thdBbSXDzNEfwi5F2Q1d3R4Bch91ljS0EAGKvys4pPbAr4LsrqCzDUwFlIM+tFCuadDH2EN99885NPPnn//fdPnz794Ycffu655zLPL1iw4KqrrnrggQcuu+yy/o3XrVt34YUXXnDBBbnvF0bmObek8f1ANkGYTqR7DkYmfc9TgFYBGFnxFDYRSfEdcaZi9IWvAzt6K84tKUCrQIaq0ZKSkh07dtjt9j179qxfv/7666/PPH/HHXcMDbw777xzzpw5ue8URlVc5xAT6WjL6CPruw9EnDW0zYl6UYA8MxHls1xdX4ZG3TDFib2HIu65xQVoFMjz3VddXf34448PenLFihVDt7z55ptl2SOMzkR4zikJ7OidMnGUNZU69/SWn4XjDaAQyucUH369ZeKl5SbzSKWMnXv7Smc4s+zjhxxhrlE9c88r6foqNPJUv9FWge+Il81CByFAITiqaLLUFvxilIvCwI4eD+6LFgqCUM/sTmvRFLZzz0gTwDb/o7NqcTkmcAIomOrL3C3/7BxhZvzQMU5KE65JmOapQBCEOle9xN38UWeSO/mw02irwPljFfNw4glQOEWTWXvRsBeFUlo6/nd/zRI3pr8vGAShzrFeqmxmUdP/dJz0t7gcBFBEzeUVzR+d/KLQ/68eG2tFb0UhIQj1r+Zyd8+BcLR5cPlox65evjOOy0GAwnNNYugye9PmwWeoiUiq5ePg5Ku9J/0ryBMEof5ZaUvtlZ76N1sHznzfczDS9P86Zvy0BpeDAIo4dcXEnoORts+6+p9JcWL9n1sqzimh3aOPMgQZYeiYIbjnFqcE8evnj5+yzEOW2KMtQsvHwek/qaHLcbwBKMPKWGasqt334vGUkC6azBAScfTttrKZRTVL3Eo3zXAQhEZR+b8muCYxR99qM1tNrI8+bWW1s3qU8YUAkFdkse30VbX+f/e0fBSMh5KTf1BZOh0rLikAQWggjip69gNTlG4FAPwX7SZPuQo9ggpDHyEAABgaghAAAAwNQQgAAIam7SB8//33Dx48qHQrjOu1117r7OxUuhXG9eyzzyaTSaVbYVCpVOqZZ55RuhXG1dXVtW7dOrleTdtBuHnz5t27dyvdCuN68803jx49qnQrjOvFF1/s6xtpIlnIn0gk0r/2KhReQ0PDX/7yF7leTdtBCAAAkCMEIQAAGBqCEAAADM0kScOuiaUIh8NRVlZmsWS1LnMwGKRp2uFw5LtVcFJ+v7+0tJQkMU+bMpqbm6uqqsxmnM4qIJ1Ot7S01NTUKN0Qg0okEl1dXZWVlaNuuXz58scee2zkbVQXhO3t7bFYLMuNk8mkxWLBF4FS4vE4UlBBeP+VhfdfWVm+/16vl6ZHmU5SdUEIAABQSLiWAgAAQ0MQAgCAoSEIAQDA0BCEAABgaJpZjzAYDK5bty4YDF5xxRWLFi0auoEoim+88cb+/funTZt2yy232Gy2wjdSr4LB4KZNmw4dOlRcXHzdddfV1dUN2iCdTv/xj3/s//H0008/77zzCttGPduyZUt9fX3msdVqXbly5dBt2tvb169f39fXt2zZsvPPP7+wDdS5tWvXDiwqHPrxrq+v37JlS/+Py5Ytc7uxynxOOjs79+7d29LScuGFF06dOrX/+dbW1tdeey0cDv/gBz8499xzh/5hIpFYt27d0aNHZ82atWLFiizHFGjjipDn+fnz5x85cqSmpmb58uVvv/320G3uvPPOl156qa6u7s9//vOPf/zjgrdRz+66664PP/zQ6/V2dnbOmjVr27ZtgzZIpVKrVq06fPjw8ePHjx8/3t3drUg79er1119/++23M+9tY2Pj0A36+vrmzZvX1tZWVVX1/e9/f/PmzYVvpI41NjYe/4977713//79gzbYvn37U0891b9NPB5XpJ16snjx4l//+tcPPvjg9u3b+5/s6uo6++yzg8Gg1+u9/PLL//nPfw79wxtvvPHtt9+uq6t7/vnn77vvvmz3J2nB+vXrzz777HQ6LUnSX//61zPPPHPQBn6/nyTJlpYWSZJ6enooimpoaFCgoTolCEL/41WrVq1cuXLQBpkjPxKJFLZdRvGjH/3omWeeGWGD5557btGiRZnHr7zyyoIFCwrSLsPZs2cPTdO9vb2Dnn/99deXLl2qSJP0ShRFSZLOOeec119/vf/JJ5988vLLL888fv755/s/8/0OHTrEMEwoFJIkqbGxkabprq6ubHanjSvCzz///JJLLjGZTARBXHLJJfv27evt7R24wfbt26dMmVJVVUUQRElJyVlnnbV161Zl2qpHFEX1P47FYsNN5fOHP/zhhRde+OKLLwrVLgPZsWPHU089tWHDhpOuu/T5559feumlmceXXHLJ9u3bsTxTPqxbt+7aa68tLi4e+qvW1tannnpq/fr1wWCw8A3Tn5Pe0hz0Od+6dWs6nR60wbx581wuF0EQtbW1EydO3LlzZ1a7y7nBheD3+8vLyzOPJ0yYYLVa/X7/cBsQBFFRUdHe3l7QJhrD9u3bN27c+LOf/WzorxYvXtzd3X3o0KFFixY99dRThW+bjlVXV5eXl/f09Dz++OPz5s3jeX7QBgM//263O51OBwKBgjdT5wRBeOutt07aQVtUVHT66aeHQqGNGzdOmzZt6L1TkMWgz3kymezq6hq4QSAQGBgEbrc7yyDQRrGM1WpNpVKZx6IoiqJot9sHbSCKYv+PyWRy0AaQuyNHjlx77bVr166dMmXKoF/Z7faPP/448/jGG2+8+OKLV69ezbJswduoT7/+9a/7H5x11lnr16+/6667Bm4w8ADJPMDnX3Z/+9vfSkpKLrjggqG/WrZs2bJlyzKP77jjjl/96lfvvvtuYVtnCKN+zscdBNq4IvT5fP3Bnnng9XoHbdDW1tb/Y1tbWzaTsUL26uvrFy9e/MQTT1x33XUjbzl//vxUKtXa2lqYhhmKzWabN2/e0HqZgQdIW1ubzWYrKysreOt0bv369bfeemumg2YE55133vHjxwvTJKMZ9DlnGGbQbepxB4E2gnDp0qWbNm3KTMb97rvvLlq0KHO18dVXXzU1NREEsXDhwq6urr179xIEcfTo0cOHD/ffSobcNTU1XXbZZWvWrFmxYsXA53fv3p352AmC0P/kBx98wDBMbW1tgRupY/1vbyQS2bJly4wZMwiCSCQSn3zySaZMaenSpRs3bsz0C27YsOGKK67Icv0WyFJjY+PWrVtvuumm/md4nv/kk08y1yX9/0GSJG3evPn0009XppV6t3Tp0vfeey9zzbdhw4alS5dmnt+9e3cmIC+77LKvv/46c6a4Y8cOnucXLFiQ1UvLUuGTb6lU6tJLL50zZ85NN900YcKEbdu2ZZ5ftGjRY489lnn83HPPeb3elStXTpw4sf9JkMWSJUtYlp3zH7fffnvm+dmzZ//+97+XJGnt2rUzZsz44Q9/uGTJEpfL9ac//UnR9upNaWnp0qVLly9f7vV6r7zyymQyKUlS5sg/ceKEJEmxWOz888+fP3/+8uXLy8rKvvzyS6WbrDcPPfTQFVdcMfCZw4cPEwTR3d0tSdKSJUsWL168YsWKM844o66urqmpSaFm6se99947Z84clmVra2vnzJmT+c7nef6cc845//zzb7jhBrfbfeDAgczGZ5555iuvvJJ5/PDDD9fU1KxcudLj8bz88stZ7k4zq0+Iorhly5aurq6FCxd6PJ7Mk0eOHHE6nf0XvwcOHMgMqJ89e7ZyLdWh+vr6SCTS/6PT6cwMcf3mm2/Ky8szvdZ79+5tbGwsKiqaO3cuRhPL68SJE1999VUikZg6deqsWbMyTyaTyS+//HLWrFmZXpBkMvnpp5/29fVddNFFA+sFQBb19fUul6v/m4cgiFgstm/fvjlz5lgslp6enl27dvX29vp8vvnz52M2j9w1NDT09fX1/1hXV5epBc3cCAmHw4sXL54wYULmtwcOHKioqOj/2O/du7e+vn7mzJnTp0/PcneaCUIAAIB80EYfIQAAQJ4gCAEAwNAQhAAAYGgIQgAAMDQEIQAAGBqCEAAADA1BCAAAhoYgBAAAQ0MQAgCAoSEIAQDA0BCEAABgaP8fiClnQd7uknsAAAAASUVORK5CYII=", "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": 11 } ], "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", "ψ = ifft(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": 11 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.3" }, "kernelspec": { "name": "julia-1.8", "display_name": "Julia 1.8.3", "language": "julia" } }, "nbformat": 4 }