{ "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": [ { "output_type": "execute_result", "data": { "text/plain": "2-element Vector{Main.##339.ElementGaussian}:\n Main.##339.ElementGaussian(X)\n Main.##339.ElementGaussian(X)" }, "metadata": {}, "execution_count": 5 } ], "cell_type": "code", "source": [ "x1 = 0.2\n", "x2 = 0.8\n", "positions = [[x1, 0, 0], [x2, 0, 0]]\n", "gauss = ElementGaussian(1.0, 0.5)\n", "atoms = [gauss, gauss]" ], "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", "model = Model(lattice, atoms, positions; n_electrons, 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 `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.143574048788 -0.42 8.0\n", " 2 -0.156034088097 -1.90 -1.10 1.0\n", " 3 -0.156768131103 -3.13 -1.56 2.0\n", " 4 -0.157046022272 -3.56 -2.32 2.0\n", " 5 -0.157056275417 -4.99 -3.37 2.0\n", " 6 -0.157056363555 -7.05 -3.60 1.0\n", " 7 -0.157056406034 -7.37 -4.35 1.0\n", " 8 -0.157056406912 -9.06 -5.42 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.157056406912" }, "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-element Vector{StaticArrays.SVector{3, Float64}}:\n [-0.05568508591796781, 0.0, 0.0]\n [0.055685498201180614, 0.0, 0.0]" }, "metadata": {}, "execution_count": 8 } ], "cell_type": "code", "source": [ "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+gvaeTAAAgAElEQVR4nOzdd3xTVd8A8HNHdtKM7r13S2nLhgJFVrEMmYKAovLAg6ggKioOFNzyyCOKIiCoCCqCKJQNQkFmC7S0dO89kmbPO94/wlt5GKUjTdLkfD/8kYaTm5OTm/O759wzEJqmAQRBEAQ5K9TWGYAgCIIgW4KBEIIgCHJqMBBCEARBTg0GQgiCIMipwUAIQRAEOTUYCCEIgiCnBgMhBEEQ5NRgIIQgCIKcGgyEEARBkFODgRCCIAhyanYXCHNzc00mUycTUxQFl4izIZIkbZ0FpwbL37Zg+duWBcvf7gLhpEmTmpubO5nYYDBQFNWr+YE6oNVqbZ0FpwbL37Zg+duWBcvf7gIhBEEQBFkTDIQQBEGQU4OBEIIgCHJqMBBCEARBTg0GQgiCIMipwUAIQRAEOTUYCCEIgiCnBgMhBEEQ5NRwW2cA6l0mChQp6GIFXakCJgqQNOAzgCcHhLogCRKEAS+EIMh2SBoUyukiBV2vBW0GwEQBgoAAHogQIlEihAurZ2uBJe2YGnXg5zLqeB31dyPtzUUiRUioADAxgCGgSQcuNIECOVWmpPu7IlMD0ZnBSLAAsXWWIchZtOjB75XUvgrqUjPtzUVixYgPF4hZQEMAkgJZLaBYQZWp6ERXZJwv+ngIEi6EP8/eBQOhozlRR3+RT/7dRD8WiD4bif40GhWz7p9SQ4ALTfRvFdSQP8l+EmRlHJbmj8AfHATd5dNPP7169eq9z5MkiWFYlw6lMIJqNS01ADc28GCDR1gIjgIAQBMATXckiwAglAZyA9hnoL/SAS4O/HjAk+PsP8/Ro0cvW7asN46M2Nui1X5+fpcvX/b19e1MYp1Ox2Qyu3ouOqojNfTaa6TaBF7ph84KQXmdvsgxUeDXcuo/eZSJAp8Mwib6dfbnplKpBAJBN7ML9Rgsf+tISUkZP358VFSUrTPi1K5du5aXl3fw4MH2Zyx4/sMWoSMolNMrL5GVavD+AHRaIIp28bqRgYInwtAnwtA/qqiVl8jPeeCr4ViYi5NffULQP1JTU0eMGGHrXDg1DoeTl5fXSweHgyX6NhMF3rtOjTxETPBDc6fj04O6HAXvNDUQvTkdT/NHh/5JbLhJkfbVWQBBENQrYCDsw/La6EF/EJebqeuP4SviUIsMAcVRsCIOvTwVz6imxmQQdRoYDCEIcnCWCYS7du0KDw/38PB45plndDrdg5IdPnx4wIABhw4dssibOrltRdSYDOKFWDRjAu7Ls3A3ZogAOTkJH+eLDjhAHKuFsRCCIEdmgUBYUFDw3HPP7dy5s6ioqKqqat26dfdNplAoVq9eLZPJpFJpz9/UmWkJMP8MuSmfykzHF0X0VpseRcCbiejPY/CnM8kNN+HuxxAEOSwLVKPffffdtGnThg8fLhaL33zzze3bt993JOqqVatWrFjh7u7e83d0ZtVqOuUQgSHg0hQ8StTr41lGeSOXpmK7SqlFmaQRRkMIghyRBQJhYWFhQkKC+XFCQkJzc3NbW9tdaU6fPl1WVvb000/3/O2c2YUmesifxPww9PtRGMdaA379ecj5ybjcANKPESqTld4UgiDIaixQm8pksvbJHC4uLgCA1tZWiUTSnkCj0bzwwgsHDhxAOjEfVKPR+Pn5tf/5448/Tps27UGJnWoe4YEadGUW/u0QYpy3Qa229rvvHAJeysJHHaT2jTS5s2+3+DUaTWe+U6iXwPK3DopykM6Q5ubmpKSk2tpa859bt25taGh4++23H5T+6NGj+/fv//bbb62VwYcgSVJ9R93XyfOfy+Wi6EOafBYIhBKJRKlUmh+bH7i6ut6ZYO3atQMHDlQoFNnZ2RqNprKysqqqKjAw8L5H4/F4eXl5nZxQj2GYkwTCz/Ooz29SJx/FEiRMW+VheypYe42cdAY9mYaZh+fQNM3n822VHwiWv3U8tBrtKyiKam1tNT/W6XTvvfdednZ2B+knTJjwyiuv5Ofnx8bGWiWDD4Fh2J0nvAXPfwsEwvDw8Pz8fPPjvLw8V1fXO5uDAACKom7evLlkyRIAQGVl5c6dOzEMe/PNN3v+1s6ABuCNq+TBavrCFMzP0qNDu2ptEsZnUKMyyFOTsEA+bItAkA1UVVXl5OSEhobu3bt35MiRY8aMKSwsPHTokMlkmjJlijloaTSajIyM/Px8DoczefLkeyPZvn37kpOTPTw8AACnT5/28PCIi4v7/vvvp0+frtFoMjMzZ8+ejSDIE0888fXXX3/55Zc2+JxWZIErnUWLFu3fv//GjRtarfajjz5atGiRubm6bt26I0eOAAA2bNiQ9f9iY2PXrl0Lo2AnkTRYfI4820BnpuM2j4JmL8ejK+LQ0RlkpQpOq4AgG8jNzV22bNny5cvFYjGKohkZGWlpaQiCcLnctLS0zMxMAMC1a9cyMzMDAgJomn7kkUeysrLuOkhGRkZqaqr58fbt20+fPg0AWLlypVQqLSkpWbt2rfm/UlNTMzIyrPfZbMQCLcL4+PhPP/00PT1drVanpaW988475ueLior8/f3vSuzr6wtXR+wkIwXm/0UqjPSJSXjnFw61guUxKArAmMNkxmgkGn6ZkPO51kqvvkpa572GeCDrku+++6NUKg8cOCAUCgEAERER27dvHzNmDADA09Pzww8/HDlyZEpKSkpKijkxi8X67rvvBgwYcOcRbt68uXDhwoe+e0RERGVlpVKpNI//cFSWqV8XL168ePHiu57ctWvXvSn3799vkXd0eDoCzDxFMFHkz/E4y/7ugS6LQUkaPPoXM3MybSdNVQiymiABsrqflX6Wbuz7PBkZGWmOgkqlsqSkZN26dR9++CEAQKVSNTc3AwAaGhqWLl1aUFDAYDB0Ot29XaMqlYrL5T703c034VQqFQyEkLVpCDDlOOHNRXaOxHB7vU//fCyq0pHjjqBnH8U9OLbODQRZkYQFxvra8vqPxbq9uZp5tOCWLVvap2ibh/a8+eabUVFR5rH6X3311eHDh+86gru7e/s8Nzabrdfr2/9Lr9ez2bfDr1QqRVH0rvGPjsdea1knpiHA5GNEsAD5YZT9RkGz5yOJOSHIuCOEzGDrrECQU2Kz2ampqbt37xb/P6PRCABoamoKCQlBEMRoNO7Zs+feFw4cODA3N9f8OCEh4fjx4wRBmP88ePBg//79zY9zcnISEhLa46Kjsu+K1vkoTWDcYSJKhGxNwXqyj4TVrE3CHvFBphwntIStswJBTmnr1q1HjhxJSkp67LHH4uLiPvroIwDA0qVL33jjjccee2zAgAH3nY02a9Ys82BGAMCSJUuEQmFERIRarR47dmx2dnb7SplHjx6dOXOm1T6LrcCuUTuiMIKJR4kECfLVcKwvBMHbNgzBFp8jp54gMibgTHhlBUG9bOLEiaNHj27/Mygo6NKlS5WVlU1NTcHBwZ6engCA9PT0goKCsrKyiIgIgUBgbiZ6eHjU1NSYX5WammowGHJzc/v168disfbt29fc3BweHr5r164hQ4aY02i12oMHD164cMHan9DqYL1lL+RGMO4IMdgD+XpEX4qCAAAEgG+GYxwMeSYT7mAIQb2OwWDcNfYeQZDg4OAhQ4aYo6CZh4fH0KFDXV1dmUymecwLiqLttxIRBNm0adPly5fvTI9hmJeXV/szV69efeONN8xzDR0bDIR2oc0Axh0mUryQjUP6WBQ0w1HwyxisUk2/ctlKY8ohCOqh4cOH3zXa/8UXXzQPRjUbNWqUkywQDQOh7ckMYOwRYrQ3smGw/c2T6DQODv4chx+tpeGeTRDUR73zzjtisdjWubABGAhtTGoAYw8TY32QT/tyFDQTs8DRidgX+dTuMhgLIcdhbwPBjhw5Mm7cuF46+KFDh5YtW9ZBgsLCwvHjx993r72+CwZCW2ozgIlHiLG+yMeD+nwUNPPjIUcmYi9dIo/XOdTvBHJapUq6Wm1fJ3NcXNyKFSt648g0Ta9evXr58uUdpImKimKz2QcPHuyNDNgKDIQ206wDozOICX7IJ44SBc1iRMivj+ALzhA5MvuqPiCoqxq0YPwR0otjXzfuCYLQarUAAKPR+PHHH9fV1a1du/a1116rqqrSarWff/75K6+8cu3aNXPi1tbWrVu3rly5ct26dcXFxe0Hqa6ufu+999asWVNaWrpx40aFQgEAOHv2LIfDiYmJAQBkZmaaFyDduXNndXW1SqXasGGD+bULFiz4+uuvrfypexUMhLbRqAOpGcRjQcj6AQ4VBc1GeiFfDcPSj5H2dikNQZ2nMoFHjxHPRKIilq2z8r/y8vK++OILAIDBYHjttdeeeeYZPz8/hUIxbty4RYsWAQAkEklqampLSwsA4OLFi7W1tYMGDWIymSNGjCgtLQUAtLS0DB48WK/XR0VFLV269PXXX5fJZACAw4cPt0/MOHbsmHk9mi+++KK8vFwul7/11lvm/xo9evRff/2l0+ls8OF7B5xHaAN1GvqRw+TCcPSN/g57ITIzGK3TgLSj5PnJuNjO6hEIeigTBWadIga5I2v6o0fv+V9jVaH8t6+skxNmSJzosSUdJFi/fv2AAQOefvppNze3sWPHmgeCHj9+/MyZM7NmzZo8efLkyZNpmpbL5bW1tXv27Hnrrbe2bds2atSoDz74AACQlJQUFxdnPlReXt6UKVMemiV3d3cOh1NaWhofH2+Jj2h7MBBaW6WKHnuEXBqNvhzvsFHQ7MU4tEZDTztBHEvD2Q7Y7oUcFg3A4nMkE0W+Gn7/Exf3DBDNfsE6mUE5vI4TmMMYiqJubm7tIc3Dw0MqlQIA8vPzn332WZlMJhAIGhsb09PTAQDFxcWJiYnmlDExMe0rl2q1Wg6nUwsHc7lclUrVrQ9kj2AgtKoSBT3uCPlKP/S5GAePgmafDMKeOEMuOEP+MqZvrBgHQQCAt7LIQgV9ehL+oFm9KJvL9A+3bqYeCMP+idbmFbcBAAiCmAd2vvjii0899ZR5X/TXX3/dHB1FIlH7itsajca87gwAwN3d3dxHCgDgcDhyubz9yDqdrj1GUhQlk8nunLzf1zlFdWwnbsro1MPk20nOEgUBACgCdo7EWvX0S3CiPdRHfFNA/VpBHxyPcx2imdAeseRy+a+//mp+8tFHH929e3d9fT0AoH0IDABg8ODBOTk55sdJSUknT55sb/b99ttvycnJ5sfFxcUCgSAkJMRqn6K3OUuNbHOXm+lxR4j/DEafjnCuMmdh4Pdx+Ol6+jM40R6ye39WUeuuU0cnYu6Ost3CqlWrnn322bS0tJSUlKSkJPOTY8eOXbJkSWJiYlBQEEVRHA7HvDfh7Nmz27ehmDRp0pQpUyIjI4uKiubPn3/gwIFNmzaZX56RkTFr1iwEcaBOHtrO+Pr61tbWdjKxVqslCKJX82MRx2spj13Gw9WUrTNiYUqlspMpa9VU4B7TT6Vkr+bH2XS+/KHO+LuR8thlzGq5+3c6YsSIc+fO2SRL90UQhE6no2na3EXZ/rxcLjeZTObHarVar9ebH0ul0pycHJ1Op9PpNBpNe3qKoiiKys/PFwqFFHX7U8+aNWvfvn3tabRabXR09O+//97+DEmS8fHxhYWFvfb57u/gwYPp6el3PmPB898hGv/27fdKatnf5L5H8BFeDnQB1UW+POTwRGxMBuHKQib4OW85QHarVEnPOkX+MApPdrP38xPDMPN9QQRB7lwR7c5lQnm8f4bYSCQSiURy10EWL14cFxdnMpm++eabNWvWtDfvPvroozsny3M4HDabLRKJ2p8pLS1dunRpZGSkRT+TjcFA2Lv+m0dtuEmdmITHie3919XbYkTIvrH49JNExgR8gN3XNZBTqdXQYw+THw9Cnecq7YknnsjOzkYQ5Icffhg2bFj78yEhIS+++OKdKZ966qmAgID2PyMiIiIiIqyXUauAgbC30AC8e43cW06fm4wF8p3l19Wx4Z7I9hR88jHiTDoeKYRlAtkFhRGkHyOfj0XnhznR/fvRo0ffualhB154wUoTRWzIib54a9IRYM4p8q96+vxkHEbBO6UHIO8PxNKOkvVauOgMZHtaAjx6jBjni6xy9Hm9UAfgd295jTowOoNgYeB4GlxU5T6ejkCXRKETjpBtBltnBXJu5uVjwlyQT/r+3i9QT8CuUQu7LqUfO0E+HYm+lQhnkD/Q6gS0VU8/eow4MQnnwXMQsgWKBk+dJXEU2ZbS93bDpijq22+/vX79ur+///Lly+8czAJ1A2wRWtKv5dTEo8Rng9G3YRR8mE8GY9EiZPoJwgCn2kO28PxFsk5L/zIGw/tgLbhkyZJDhw6lpaVVVVWlpqaSJPwV9Qi8GrcMggJrssi9FfSJNLyfBAbBh0MA+DYFe/w0Oe8vso9WRlDf9fpV8moLfXJSN1fBvdVavPHqN5bO1P3184hZnvzsnc/U1tbu3bu3pqZGIBBMnTo1JibmyJEj5kVEoe6BgdACmnRg7mmChYGr03BXeFOw0zAE/JSKTTtBPJ1J7hwFFyOFrOSDG9ShavpMOu7C6OYRAoV+qwZ3tI27BfEZdy+6XVRUFBoaKhAIAAAIgiQlJRUVFcFA2BMwEPbUqXr6ybPk4kj0rUQUVuVdxUTBvkfwSceIpX+TW0b0vVs1UJ/zeR71fQl15tEeXbPyGNxISZjlMtU1Wq1WqVS2/6lUKl1cXGyVGccAO6S6z0SBNVnkk2fJH0Zh7yTBKNhNHBwcHI/faqNfvAjvc0C9a/Mt6st86tQkzJtr66z0TGlp6ZkzZwAAFRUVZ8+eTU1NtXWO+jYYCLupUE4P+5PIkdLXH8PH+MAY2CN8BsiYgF9upl+6BGMh1Fu+LaQ+yaVOTcL8eH3+B9u/f/9169YNGzZs4MCB77//fliYzZqnjgF2jXYZSYMv8qkPb5DrB2D/ioJXEpYhZIJjafi4w8Sqy+QGOKkLsrSthdT7N6jTk7AgQZ+PggAAHo936tSpyspKd3f3O5cVhboH1uNdk99GjzhIHKyiLk7BYRS0LBETnJiEZzbQqy6TcNUZyIK2FlLrb1CnJ2GhLo4QBdsFBQXBKGgRsCrvLC0B3rhKpmYQiyLQU4/iDvaLshMiJjiehp9vpF+4AGMhZBlf3qI+yHGoKBgdHb106VJb58KhwEDYKb9VULH7iGoNyJ3B+FcUHBbTi8QscGISfl1KLzlPUjAYQj3z2U1qYx515lHHiYIAgLCwsPnz53f75TRN6/V6i+TEvOWhRQ5lWzAQPsTVFnrUIeL9G9SOkdiu0ZgXx9YZcgIuDHB0Il6qoBecIU1wW3uou97JJrcXUWcfdcDtX/Ly8nJycjpOc/Xq1eLi4nufv3HjRlRUlEWyERMTc/369W6/nCCIvXv3UpTtf+QwED7QLTk94yQ5/SS5MBzNnoaP9na035I94zPA4Ym4ygRmnCT1cCQp1EU0ACsukQer6cx03LfvjxG9188///zjjz92nGbz5s13brFrh3Q63ezZs+1hfTg4avQ+cmX0BzeoMw3UK/2wXaMxDiwkW2BjYN9Y7Kmz5MSjxB/jcCHT1hmC+ggTBZ7JJCvU9OlHcZEjnjYFBQWnT58mCOK1114LDQ1dvHgxQRBbtmzJysry9fVdvny5l5fX5cuXs7Ozq6urW1paEhMT58yZc99DGY3Gr7/++vr164GBgcuXL3d3dzc/f+bMmQMHDshksqSkpBdffJEkyZ07d165coUgiJSUlCeffBJFH9iI2rBhQ2pq6t69e5ubm2fOnJmWlmZ+vri4eNu2bVKpdPTo0fPnz0cQZPPmzQCANWvWoCj63HPP+fv7W7qoOgu2CP/HyTp60jEi7Sg50B0pm8NYFY/CKGhDDBT8OBpLkCCjDhENWlvnBuoLNASYeoKQG8HxiY4ZBQEAQqHQ29vb09MzOTnZvFn8/PnzDx48OGvWLJIkk5OT5XK5m5ubq6urn59fcnJycHDwgw41Y8aM06dPz549W6VSDRw4UKPRAAC+++67hQsXJiYmzps3r7W1laIovV5fXFw8ZcqU6dOnb9++ff369R1kb+vWrfPmzYuNjZ0wYcKzzz77559/AgDKysqGDh3q4eExderUjRs3vv766wAAcw9tUlJScnKybYe/wmoeAADkRvBjCbWlkEIAWBmP/j4WZcGZbPYBRcB/h2If5lDDDxKHJ2BRIgfs5oIspVkHppwgYkXIlhG9u4y7skJbvLu2F9/gDqJwXths3zuf8fHxiY6O1uv1s2bNAgCUlZX9+eefdXV1YrF40qRJWVlZO3fuXLFiRVBQUFxcnDnNfd28eTMzM7Ouro7P50+aNOnSpUu7d+9evHjxm2+++f33348bNw4AMHHiRAAAn8//5JNPdDpdU1PT6tWrX3vttbfffruDPC9dutQ8nEej0Xz22WdTpkzZtGnTjBkzXn75ZQBAWFhYYmLiO++8M2bMGADAjBkzGIzurvpqIU4dCE0UOFFH7yqljtRQE/3RL4dho7wRWNHaodcTUD8eSM0gfnkEH+kFvyLoPooV9KRj5Pww5J2kXl+0lu/Hjlsa1MtvchvKfEhILy0tDQoKEovF5j+TkpLuO0bmXiUlJREREXw+/84XyuXyhoaGIUOG3JlSo9HMnTu3sLAwMDBQp9M1NDR0fOSEhATzg/79+5szU1JS0r4seExMDIZhVVVVvr6+DzyEdTljINQQ4GQddaCKPlhFRYqQ+WHol8MYErhrhH1bEIb6cpHZp4hPB2MLwmCXPvQ/zjbQj58mPhyIPRVhjXMDZaBsVxt3vLbPWxCJRHK5vP35tra29lt9HROJRG1tbXe+0MfHh8/nM5lMqVRq3t3CbMeOHRiGmUPa1atXx48f3/GR2/Mjl8vNEVosFrc/qdVq9Xp9e+S2B85SoRhIcL6Rfv8GNe4I4f2T6ctbVKIrcmM6/vdk/N/RKIyCfcIYH+SvR/F3r1FrsuAUQ+gf24uox08TP6Xi1omC9sDd3b26utr8uF+/fiiK7t27FwBQU1Ozf/9+8/iUO9Pc18CBA5VKZUZGBgCgvLw8IyNj4sSJOI5PmTJl/fr15sGcLS0tNE3rdDoEQQAAJElu3Ljxodn79ttvCYKgKOqbb74xd66mpaX98MMP5li4adOmpKQkT09PPp/P4XA6zqR1OGaLkKRBnYYuVYJbcjpPRme30rfkdKwYGemFvBiLpY5DeI75uR1ftAi5NAWfcZKYcZL8YTQmsPGdBcjGCAq8coU8XENnpuPhQifqM58zZ86+ffuCgoJGjBixa9euPXv2PPXUU++++25zc/Nrr702YsQIAMDTTz+9cOHCoKCg6dOn/+c//2l/LYqibDYbACAQCHbv3v3MM88IBIKWlpb169cnJSUBAL766qtFixaZu1v1en1BQcGTTz75/fffR0dHUxQ1a9YskUhkPhSbzb7v8NHw8PC4uDij0RgUFPTNN98AAObOnZuVlRUZGSkWi1ks1u7duwEACIKsW7du7NixKIru37+/vUPV+hB7WxfAz8/v8uXLnew7rlfo/1uEIwiiMQGlCSiNoFlPN2hBvZZ2ZyNhLiBahMSJkUQ3JNEV6d5W1FAHVCrVnf0nVmOkwAsXyPNN9IFxWJgDrRjSVbYqfzshNYA5pwgGCvaM6d0BoikpKR9++KE5utizlpYWiUSCYV2u6Zqbm11dXe96oXnXQ09PT+T/B07U19dLJBJzEO1AVFTUjh07EhMT1Wq1m5vbnf9lMpmUSqWrq2tXcwgAOHTo0IsffTNs/R+eHPDZYAxY9Pzv2y0jDKFFDICiiC8XuDCBkAk82Kg3F/hwETjs04ExUfDNCGxLITXiILEtBU8PcN5Y6LSutdKzTpEzg5EPBmJwQ2ezTt4avJeHh8e9T3K5XC73f7Zt9PHx6fwx2Wz2vSGTwWB0LwqaeXHBeD+kN26L9O1A6MIAr/ZDMMxZbgxAd1oShSZIkDmnycstyNokWBs6ke1F1BtZ5OZh2Ixg+Nu3O++9915ISEhvHFnCQnppoBw8jaA+bIgHkjUNv9RMjzsMZ9w7BbUJLDxDbsyjMtNxGAXt0+zZsz09PW2di66BZxLUt7mzwbGJeKoPmnzAdKzWvm54Q5aVI6MH/kEwMXB5Kh7pTENjoN4GAyHU56EIeCsR/XkM/q/z5MpLcJFuB0QD8HkeNf4I8XYiui0F4/btWzqQ3YGBEHIQI72QG4/h9Vow6A8iRwabho6jRkNPOELsq6AuT8HnhsIqC7I8eFZBjkPMAr+MwV7ph44/Qrx/gyJsv80Z1FPfl1ADDhCjvdGz6XiQAHaHQr0CdjFAjmZBGJrqjTx7jjxQSW0fifWTwNqzT6rV0EvPk3VacCINt+2XmJycPGfOnHsnA9A0jcDFia1Fq9UOHz68lw7etyfU63Q6JpPZjQmkkEXY84RuGoCdxdRrV8l/RaFr+mMOuZyCPZd/T1A02FJIvZNNvhCLrU5AGbbut3rQMtMajca2mwc5Gy8vrzunNsIJ9RD0EAgAiyLQiX7oCxfJfvuJzcOwsb7w4r0PyJXRS86TOArOpOMx9rHrFofDue/EOEe9EHFCMBBCjsybC/Y+gmXU0P86Tw5yRz4bjPrx7KJuhe6lMIK118jdZdT6ZOzZKBR+T5DV2LrTAYJ636P+SP4MPEoEEn8nPrhB6QhbZwj6XxQNthdR0b+ZtAS4NZOxGEZByLpgIIScAgcHa5OwK1PxG1I6Zh+xpwzu42QvTtfTAw4QP5RQB8fjW0ZgrnBPNMjqYNco5ESCBcivj2DnG+lVl8kNN6mPBsIbh7Z0Q0q/fpUsVYIPB6Iz4XppkO3AQAg5nRFeyKWp+G8V1HMXSF8uWD8AG+YJw6FVFcrptdeosw3Um4nYv6JsPy4UcnLwBIScEQLArGA0fwY+Pxx94gyZdpS40AT7Sq2hUE4vPEOOyiD6uyKlcxjPxcAoCNkePAch54Wj4OkItGgWPj0InX+GfOQwcaIOhr9UgNcAACAASURBVMPecq2Vfvw0OTqDiBIhJbMZryWgPNghBdkHeCZCzo6JgsVR6KIIdHcZ9dIlkomBl+PRmcGwpWIZNADHa+kNN8lCOVgZj25LYfAZts4TBP0vGAghCAAAcBQsDEcXhKMZ1fTneeSrV6jnYtBnIlH3u9fVgjpLQ4BdpdQXeRQTAyvj0Lmh8NoCslMwEELQPxAA0gOQ9AA8R0Zvyqci95rS/dEl0ehwOJqmK/Lb6G8LqZ9KqVHe6JfDsVRvWHqQXYOBEILuI0GCbEvBPh2E7SimFp8jEQCejkSfCEO9OLbOmR1TmsCv5dSOYqpKDZ6OQK5Px/3hOj5QXwADIQQ9kJgFXopHX4pHzzXSO4qpmN9MwzyQeWHo1EA40OMfJgocr6N3l1KHa6gxPuhrCegkfxSDERDqO+CvGYIeLsULSfHCNhHY75XUT6XUc3+TE/zQWcFImj/qtLulExT4q4HeW0EdqKQihMi8UPSLYQy4LgzUFznrjxiCuo6Hg/lh6PwwVGoAv1dS3xZSz5wjU73RKYHIJH/U0zl6TRVGcLyO+rOKPlJDhQuRWcFo1jQ8gA8bgFAfBgMhBHWZKws8G4k+G4m2GcDhGurPanrVZVOoAJnoj4z3RYd6ILhjDY+kaHBDSp+oo4/VUlmt9AgvZEoA+vEg3IcL4x/kCCwZCA0GA4vVUc+IyWRiMOAcIshxiFngiTD0iTBAUNjfTfTxOuqlS2Sxgh7miYzyRlO8kAFuCKtv7glM0iBXRp9rpM820GcbKA8OMt4PebkfNtobcdreYMhRWeaMrq2tnTt3bm5uLpPJ/Pzzz+fPn39XgnfeeWfHjh0tLS08Hm/58uVr1661yPtCkJ3AUTDKGxnljb0/ALQZQGYjdaaBXnGRKlTQCRJkkDsyyB1JdkPChIg9t6FqNHR2K53VQl9qpq+20H48ZIQXMiMY+XIYw5v78JdDUB9lmUC4YsWKuLi4s2fPZmdnjxkzJjU11dfX984E7u7uJ06ciIyMLCwsHDlyZGxs7KxZsyzy1hBkb8QsMDUQnRoIAABqE8hqpS8303sr6NezKJme7ueKxIuReAkSJUKiRYgN52PIDKBIQd9qo2/J6VwZnSOlUQQkuyED3ZFV8dhgD0QCR75AzgGh6Z4uriiXy93d3UtKSoKCggAA6enpI0eOfPXVVx+UfubMmbGxse++++59/9fPz+/y5ct3xdEH0el0TCYTw/pm31Pfp1KpBAKBrXPRl8gMIEdG35TRt9roAjldIKf1JAgXIiECJIgPAviIPx/4cBEfLnBnP7xPtTPlT1CgWU836UCdBtRq6Go1XakGFSq6REETNIgUIrFiJEqEJEiQeAmA9/y6BJ7/tmXB8rdAi7CqqorJZJqjIAAgJiamrKzsQYllMtn58+eXLl36oAQ0TSsUCi73dkcMj8djMpk9zyQE2QMJC6R6I3eutCI3ghIFXaGiK9WgREmfrAeNWqpeC5p1NAsDbmxEwgIiJhAyEQ4OeDjg4ID9/wHSaMSZTNL82EQBtQnoSaAjQJuRVhiB3Aha9bTKBNzZwJOD+PGALxfx4yGP+oNgARrmgng4xzBXCHooCwRCuVzO4/Ha/xQIBBUVFfdNaTKZFi5cOGHChLFjxz7oaFqtdujQoSh6e9Tdli1bJk2a9KDEsEVoWxqNBrHre159AA5ANAdEcwDwuPu/lCZEagByI5AbEaUJ6EigIxEtAYzU7TLHaSML3H6MYXQQGzBRmosDIQMImUDIoCUsIGE+oMuHBGp1b30oJwHPf9vqZPlzudz2gPIgFgiEbm5uCoWi/c+2tjZPT897k5EkuWDBAgDA1q1bOzgaj8frfNcohmEwENoQTdN8Pt/WuXBYfAB8OkygUpkEArgouM3A89+2LFj+FpjuFBgYiON4fn6++c/r169HRUXdlYaiqEWLFslkst9++w12dUIQBEH2wwKBkM/nz5s3b82aNQ0NDXv27Ll+/fq8efMAALm5uenp6eY0y5YtO3Xq1PPPP3/+/PmTJ08WFBT0/H0hCIIgqOcsM31iw4YNK1euHDZsmLe39x9//CGRSAAANE0TBGFOoFarY2JivvjiC/OfkyZNio6OtshbQxAEQVBPWGD6hGXB6RN9CBw+bluw/G0Llr9tWbD8HWtJRAiCIAjqIrhooMNq08vzWgryWgqrlDUN6iaprs1EEXpC78IUCNkuHly3UHFwuDg42au/K0ds68xCkHNRGzXXmnKLpCVl8sp6dZNCr1QalTjKYGFMIcvFm+/p7+IT4xYZ7x7jxbtnYg1kaTAQOppWrfRE5dlzNZeqlbVx7lGxblFTwid68zxduRImymDjbIVBqTCoGjVNZW2Vf9de2ZS9zYPrNtJ/2MSQMZ48d1tnH4IcmdKgOlF59q+qc+Xyqn4eMdGukelhE/wEPkKWQMhyMVGEkTS26RWN6qZKZc35msubr+0QslxG+g95JHBkoNDf1tl3WPAeoePIbsz5vTgjpyl/ZMDQ0QHDEz374ejDS4aiqQJp8YmKs6erzkW6hs2Nnp7k1a+T7wjvkdgWLH/b6lL5l7ZV7Lm1/1J91lCfgeNDRid6xDOwh2/FQ9F0obT4bPWFE5Vn/AW+UyPSRgcMRxF4SwsAi57/MBA6ggt1V37M+1Vr0s2KmvJI0CgO3p1J1ibSdLIq8+db+9k4e3HCggHe/R/6ElgR2xYsf9vqZPmXtJVvvfFjmbxyVtSUyWETeIzubORBUOT52kv7ig7JdG3zY2eOD0nFEGev92AgvA0GwlutRZuvfacj9E/GPz7Cbwja4wWfaEBnVl/clvOjB8/9uaRnQkSBHSSGFbFtwfK3rYeWf4u2dcv1H6415TwZ//ik0HEM1AK3onKa87+/+XOLVrok8ckRfoN7fsC+CwbC25w5EMp0bZuv7chpznsmYf744NSeh8A7kTR5qPT4jtzd44NTF/Wb96AmJqyIbQuWv211UP4kTf5WePCn/N+mRUyaGzO9e500HbhSf+3r6ztEbOHKgUsDXPwse/C+AgbC25wzEFI0faD48Pc3f04PH78gdjYb761d4+QGxTfXdl5rurl6yPPJXgn3JoAVsW3B8retB5V/mbzywwsbxWzRioFLfAXevfTuFE39Xnz4h5u/pIeNezL+cSbmdEtXwkB4mxMGwmpl3SeXvkAR9OXBywNcOlVKPZTdmPPxpU1DfJKXJS1i/++FLayIbQuWv23dW/4UTf2U/9u+ooNLExdNDBljhTzIdG2bsreVtpW/OuSFeHfnWq4LBsLbLBUISQNlVJhMGpI0UKSBBACgDBRjoiwxgyVmIKhd7LRCA/q3woO78vY+1e/xqeGTLNsX2jGNSbspa2t+a+HbI14JF4e0Pw8rYtuC5W9bd5V/k6Zl/YX/MDHG60NedOO6WjMn52oubcza8khgyrMJ8+2naWiQm/RSI2WizZUqxkRRJsrg40w+jvMs0HqBgfC2bgRCk5rQNhl0LUZ9q0HXYtRLjYY2E03STCHO4OMoE8XZGACANFKUkdLLjCYVwfVhiyP44miBS3B3hntZRLO29cOLG00k8cawFT58L5vk4VTVuU1Z3y6Imz0jcrL5GVgR2xYsf9u6s/zP1VzacOWrOdGPzYl+zJoXqe2UBtWGK5urlLVvDlsZdsfVqpVpG/SyW6q2IrWqWoezULYbE2OhGOt2pUoaKJOaMKkIiqTZYgbblcl2Y3LcWRx3JsedxRIxQFdKDgbC2zoIhBRBGxUmQ5tJ32bUS036VqOu1aBvNgIMcNxZXA8Wx53JdmOx3ZhsMQPnPjCUUgStqtLKi9XSXCWCIT4prh4DRAhm1RP9TPXfG69+Mytq6tyY6badQtSgbnrn3MfefM9XhzzPY3BhRWxbsPxty1z+BEVuufF9ZvWFtSmvRrtG2DZLxyvObL62fV7szFlRU5AuRZUeooE0T1l/TqpvNboluIgi+C4hPIz1wMqKNFAGmVEnNepbjboWo67FoGs2EHqK48bkuDPZrrdrZpaYwRQxMOb9jwMD4W0qqabxlBxBEJqiST1FGilCRxJa0qQiKBPFdGGY+zZvX3e4MTlurO43yWkgL1HXnWk1tJnCZvtap3WoJwybsrfeaMp7a/iqKNdwK7zjQ5lI06bsbdeacteNfN0NFcOK2IZgILQtlUplwom3z33MZ3LfGLbShWkX30WDumnd3xt4DO4bw1aI2SIrvKOu2VC6t540Un5j3FzjXbp9L4k0ULoWg77VqGu93VdnkJuMchMAgCHAGTwM52IYC2MI8NDp3gAGwnZquUaZr0NRFCAIzkFRBopzMQYXYwhwnNNbI2ikucryAw2u/VyCp3j16u3DCnnV2vOfRLqGrxy41OLDr3vIfOG5JO7JtMixts6L84KB0LauVF375NqXU8InLIibbdXm18OQNLkz9+fD5SffGLrivuO9Lag+U1pzsiVgvLv3cNdeKgPSQBmVBKElCR1JGkiaAu6JQgADYTtbjRol9GTRDzUAQaIW+nfQ/O+Jg6XHtt3YtSxp0QSrjD3rhtK2ijVn3h8TlLK4/wK45pNNwEBoQ+Zf6BvDVgz2SbZ1Xu7vetPN9y98PiEk9el+83pjGRqaost/b1BWaGOeCWSJH75cnMXBQHibDadP0BRdtr9BVaWN/3dwB7cYu0Fj0m64/FWVsvadEa9aZ4JEt9VLGz67sRlH8beHv8xn8mydHacDA6FNEBT5Rda3Oc15awasjPAKs3V2OiI3KN6/8LmeMLw9fJU7182CR6ZJumBnNU3SUU8G9FJj4KHgfoS2h6BI2EwfUQT/1nfVFGGxi4liWdniIyv5TP7XEz618ygIABAw+Z+OWevv4vPvYy9XK+tsnR0I6nVyg+KlU29KdbKvJ3zmzfO0dXYeQsQSfpL6zhCf5H8dXXWxLstix6VBya91AICYZwNtFQUtyxE+gw0Fp3uxJYzC76tpygKxcH9Rxqt/rV2csOClQf+2n8lAHcMQ7PnkxXNjZrxw4vUr9ddsnR0I6kVl8sqlR1/u5xG7buQbXAbH1tnpFAQgT8TOfC/ltc+vfv31tR0ERfb8mFVHmrRNhsgF/nYyx7rnYCDsGQSEzfYl9VTNyZaeHEZlVL+Z+eHRilObJ3yaGjjCUrmzmkmhY9ePfOPjS1/8WnDA1nmBoF5xrubiqlNv/av/wmcT5ttkpmBPxLtHb0vbWK2se/7E6gZ1U08O1ZqjbLmhiH028EGzGvoix/kktoLiSNST/o0XZKoqXfeOcLOl4JnDK7x5Hl+N/8RWk+V7Ls496uuJn52oPPvhxY0m0mTr7ECQxdCA/v7mL5uyt32SunZMYIqts9NNLizBB6PXPBI46t/HXv6r6nz3DmJUEmX76yPn+zP4DrWpOwyEFsDg46EzfIp21ZAGqksvpGjq+5s/v33uo5UDlz6X/IxFdmmxIQ+u26ZxHxlI44sn10h1bbbODgRZgJ7Qrz33yZWG7G8mfBYhCbV1dnoEAcjMqMmfpK7dnvvTx5e+0BP6rr2eBiW/1PkMlwgC+ka3cOfBQGgZrvEuLsHcyowu9Dk0appfOPFGTnP+trSNQ30H9F7erImNs94Z8cpQ3wFLj64qlJb03htRBK1tMsiL1eZ/2kYDaezaVQjUF1EErWs2KEo15u9dU6fv1e+9UdP83PHVXAZn4yPvSzji3nsja4qQhG5N+5ym6cVHVhZJSzv/woaLMkJL+o1177282UrfboLYlZDHvLM/KvEaKuZ5P3zy+7Hy05uv7ZgbO3121LQ+d7+hYwhAFsTNDhEFvXbmvX8nWngepK7Z0JqjlOYqtE0Glph5e3FCGhiVJr3MxOBhggCOIIgrjuRzO/EtQH2Coc3UVqhSVmhVVTpDm5ElZjCFt5fCN6kIXYuBJWa6xgvc+gn5Fm2pXG+6+d7fn82Pndm+uK7D4ODs14a++FfV+dVn3p0ZOWVe7IyHTgUmtGT10eb4ZcEOM0DmTnAeoSXVn5PKbqnilgR1kEZuUGy4vLlGVf/msJfCxMHWylqv6HgeT5WiZk3mB4O8k55Lfrrn83n1MmPVkWZFsdotUWhuf9/7g9RLjaoqnbJC01agpmnaNd7Fvb9QEMi1p0U/LMmx5xFqGw0tNxTSXIVJQ4oj+cJQniCQy/Vk3f1t0kBdr5fmKlquKdiuzKB0T76fBcLh3sI/d9/a9/bwlxM94x+UxgHKv0Xb+uHF/+oJwxvDVvgJfDpIWX6ggSLosJkdpbEyOKH+NnsLhDRJX/u0NGSqlzj6/l/P2eoL/83aMiF4zNP95jEwG6zFYFkPPRHVRs36Cxu0Jt27Kat7suxh3dnWmpMtPimuvqPcOjlvSdtkkOYoWm4oKCPtMVDkMUDEdu0bM1I6zwEq4nuZ1ERztrz5qpzQkW79hW4JQoE/pzOXMjRJN11uqz7e7JYgDJ7i1e2V8fWE4dPLX1Yra9ePfN2T59FBSscofxrQvxcd3nlzz5Pxcx6LSL9vB5WuxZj7RVnS6nC7GiMDA+Ft9hYIAQCyfFVlRmPiy2F3tVfa9PKNV7dUKKpXD3k+1i3KVtmzrM6ciBRNf3/z54yyE2tHvBrn3uUPThqp0l/qdFJj9FMBLFF3Lh00dfqmq20t1xQ8H7bXULFrnIuVNw/pPY5REd9GA3mJuvGiTF6skcQJPAeJhSG8bjTlST1VvLvWpCWjnvRnCrpca9epGt7K/DBcEvLSoGWsh83ldaTyr1XVf3xpE01Trw554d6lPAp2VAsCOX5j7OvuIAyEt9lhIAQA5G4q90lxdesvNP9JA/pw2cmtN354NHTck/GP95WZ8p3R+RPxYl3Wx5e+eCJ25syoyZ1fnpg0UHnfVHA92aEzfVC8R9GLImhZnrLxokzbZPAcJPYaKrHJ6oiW5RgVMaEhm662NVyQ4SzUa6jEPUmEsXs2iI8GNSdbmi63xS8P7tLF07maixuubF7Ub97U8LTOpHeM8m9H0fQfJYd35O6ZETl5XuyM9kHsmnp9/taqAWsievgbtDgYCG+zz0Aou6WqPtrc/6VQAECFovo/V74mKNOqQc/19TuC9+rSidioaX773EdePI/VQ17gMR6+iRVF0Le2VbEljLBZvha8yadrNjRckDVny12CuN7DXcWR/L57B7GvV8Sqal3D31JZnkoSJ/Ae5ioItORQl7qzrU2X2uKXhzA6sfOaiSK+vf59Zs3Fd1NWd36/s75e/vfVpGn5b9aWWlXDSwOX9veMBwAU/1TL9Wb7jbHkUqUWAQPhbfYZCAENrn1a4pMu2ac7cKLyzKJ+8yaHTXSwoaFmXT0RTaRp8/Udl+qy3h7x8kO2MKVB4Q81ANC9tIwTZaRarisa/pYROtJrmMRzkLgz1aW96aMVMWmkWq4pGi/ICD3pPVTiOUjc/V1CO1R1uKmtWN1vWTDa4RooDeqmd89/6soRrx76Qpf2FOyj5d8Z52oufZm9LdYtcnHIwqqvZQPXRPa0md4LYCC8zT4DIUVTJzL+br2maJhQtbj/AiHLxdY56i3dOxHP1Vz8z9VvZkdNnRP92IOuD+ozpa03FPHPBff2/TxVta7xgkx6UymO5nsNlXTvvpSt9LmKWFOvb7woa7muEIbwvIZJrNAcL/6pFmWiYbMeONbxVNW5TVlb58fOmhGV3tU9Bftc+XeJgTT+lP+b6pgx1C1w7IKhbDvbEhXAQNjODgPh+drL22786MZynfH37PjFoXxfuzt7LKjbJ2KztvX9v/+DougbQ1fcuzuMpl6f901lwooQtsRK91MJHdmcJW+8KKNJ4DlY7DFQ1I1xFtbXVypiQk+2Xlc0XW4zKgnPwWLPweLujXvqBtJAXf+sNHiKl2v83dejGpN249UtRbLSt4avCheHdOPgfaX8u43QkFc+KDo/5q8ryuyn4uekhY7tjX0Nuw0GwtvsKhBmN+Zsz/nJQBoW918wxGdA7V+tumZD+Bx730qpJ3pyIlI0vfvWvt8K/1iW9Mz44NH/PE/QN/5T5jfGzWNA96dbdJuqStt4qU2aq3QJ5noMFEliXextgMCd7LwipilaUappuipvu6USRfA8B1ujCXgvVZW24Lvq/qvCmC7/XNxcb7r50cX/DvZJXpb0NBtndfPI9l3+PVd3plXbYAif61soLfn2xg/N2tZF/ealBoywkxs9MBDeZieBMKvhxvd5vygMyifj5qQGppjPEpOayP6wZODbkY6xX9d99fxELGkrf//C5/4Cn5cG/ds80bDqaLOu2RC10N9CeewO0khJc5XNWXJ1rc413sU9USgM49nhghp2WxGrqnWt1xUt1+UsEcM9WeSRJOqlu4CdVH2sWdugj3oqAACgJwxbc344W33x1cHLB/kk9eSwdlv+lnLt45KwOb4uQbeHtmU35uzI3a0yqhfGzUkNHPHQxWh6GwyEt9k2EJI0mVl9cfetfSaKmB87c0zgyLsulAp3VouiBF5DHGSJwntZ5EQ0kaadN/dklJ1clvT0KPHwGxtK+68Ks1rXWceMSqLlurz1hkIvNbnGu7gluAhDefYzDdG+KmIaqKq10lxla44SwRH3RKF7kpDj3s3GlmVRBH3to5Lwub7l/PJPLm2Kc49+fsCzXRoXc1/2Vf6WpqzQlv5al7T67jG0WQ03fsj7pUUrnRM9bWLII91uT/ccDIS32SoQKg2qjLITvxdnePM858Q8NtR3wH1vs7cVqKqPNye82LdXrO+ABU/EYlnZx5e+GFc4IT4sKnZKd27Y9Cq9zCjNUbbmKnWtBnGUwDVWIIrk4xwbd0XYQ0VMGSlFmUaap5LdUuEc1LWf0K2fC8/H7m6N11xtzD9asSN++8pBS4f4WGaNe3so/95TsqeO68PyHXX/WRP5rYV7bv1+s/nWpNCx0yLSOl6Cp5fAQHiblQMhDejc5lsHS45drL+a4j90RmR6x/fYaYrOWl8cuzjQUReAtmxFIK/U5Gwv+SL+v4/FTJoTPc0+l6AzKgnZLZUsX6ko0/C82eJIviiCzw/g2KTj1GYVMQ00jXp5sVpepFZWaPn+HEmswDXWhe1mj4tF0IA+Wn566/UflxQtiX0k1G+wxapsBw6EpJ66uq4o+fWHrKnWoG7aX5xxtPxUrFvk5LCJQ3yTrTmaBgbC26wWCGtV9Scrzx4r/4uFMR8NGz8hJLWT/SpVR5tJPRkyzbu3c2gTlq0Ibm6u8BgoomNMX2Zvq5BXPz/gWUtdufcGiqCVZZq2IrW8WG1oM7kEc4WhPJdgLt+fY7W+U2tWxDRFaxsNynKtolyjKNNgLFQUwRdH8EURfDucYdauSFq6MWsLAGDFwCU+Kt/C76sHrImw1BfkwIGw8aJMXqyOejKgM4kNpPF01bmM0uO1qoZxQaPGBo+KlIT1dg4BDITtejsQVivrztde+qvqvFQnGx0wYkJIale/YL3MmPPf8kHvRNrhUIues+CJqKrSFf1Yk/xGuLmgrjRc25S11YvnuSxpUbAo0CJv0XtMakJRrlWWaRTlGn2rkefDFgRyBQEcvj+HLWH23jjJ3q6IjQqTulavqtaqqnXqKh3DBTfHe2EYz05u4nagRdu6NWdXVsP1xf0XTgwZY755kf9tpVt/oecgy9y2d+BAmLup3H+s+4M2D3iQWlX98Yq/TlZmIgAZHTh8pN/QCNfQrs7O7DwYCG/rjUBoJI03Wwou12dfrMvSEfoRfoNHBQxL8Ijt9hCpnI1lgY96isL5FsyknbDgiViws1oUxvMe4dr+DEGRf5Qc/jHv1+F+gxfFz3XjunbwcvtBGihVtU5tDh41OtJA8XzYPB8215vN9WRxvVgWvLNo2YqYNFC6ZoO2yaBt0Kvr9Zp6PaAB34/ND+AKAjiCQG5fWXxHY9LuubX/j5IjU8InPhEzk8v4Z+U2eYmmfH990qvhFqmcHTUQGhWm65+VDlob1e2mc5Gs9Gz1hXM1l/SEfrBP8mCf5GSvhDu/CIuAgfA2SwVCA2kskBbnNuffaMorkBaHiALNX16ExAKXM3V/teqkRrvax8tSLHUi6loMuZsqBr4Zce9SWGqjZvetfQdLjz0aOm5uzPQ+t0yPSUNq6nWaer220aBt1OuajAgGOO4stjuT7cpkS5gsMYMlZjBdGN2YsNi98qcp2qQi9DKTQW4yyIx6qVHXatS1GAktyfVkcTxZPC+WOXLbf7PvLnrCcKDk8M+39g/1Hbio3zyPe9ZqAADc+LwsYLyHJNYC562jBsL6TKmmQW+ROdA1yrqL9VmX67NvtRYFCwMTveLj3aPj3WM6s9rwQ8FAeFu3A6GeMFQoqkrbKoplZYXSkmplbYgoKN49ur9nfIJHrEW+pH/ey3F7Ry11IpburWe64AETHjiKoVUn+zHv19NV5yaHTZgVNVXMFvb8TW3FqCJ0zQa91KiXGg0yk15mNMhNRiWBczCmAGeY/3ExnIfhHAznYBgLxVgozsYACsytSQRDzJNT1Wo1n88HAFAmmjJRAADSQNEUTehIykiTepLQU4SONKkJQkuaNKRJaTKqCJOGZPAwloTJEjHYEgbblcl2ZXLcWSwRow8tL3cXPaH/s+TonoLfEzxiF8XPDRQ+cB5q6w1F/Tlpv+ctMDLZUQNh7qZy/3Ee4ihLdmIZSWN+a9H1ppu5zflFslIPrluUa3iEJDRUHBwmCuYzed04pgXLvw+sI9VDCoOyRStt1rY0qJtqVQ11qoZqZa1M1xYo9A8TB4eLQ9NCxoaJg3pvdyS2hMkWMxRlWlF4d75sh2dSE605iuTXO1qD240jWTlw6ROxM3/K/23hwWXjgkfPiZ7mybOv3dE6iSnAmQJcGHr3yWBSEUY1YVIRRhVBaEiTltS1GEk9Seop0kgROhLQgNCRAACapEkDBQCgaRpBEAAAykBQBgoAwFgogiI4B0WZKMbCcDaKczGWmMH35TD4GMOFwRTgDD7mSNdkSqPq96LDvxcf6u8Z/58xMNgCawAAIABJREFU7z30jrJrP5fKjCZVlc6ym104DKPCpGs2WLyyYmLMRM/4RM94AABJk5WKmkJpSZG09HTVuXJ5FRfnBAr9fQXevgJvH76XB9fNg+smYoustoRN324RSlWyI1WnUBTVmnQkRepJvZ4waE06lVGtMKgUBqVcL+cwOO4cV0+euyfPw1fg7SfwCRT6efE8rLksQu1frXpH7B21yBVZ3dlWbYMh/PHO9sPIdG2/Fv6RUXZioHfinOhp1hmfZp8ctUXSSfXqxr2Ff5ysyBzhP3huzIx7t5N9kNq/WvUthrDZPe36c8jyt2C/aOc1a1urFbW1qvo6VUODpqlZ09qsbVUZVUKWi5Dl4sJycWHyuQwOC2PxGFw2znoy/nEAW4TtKJpSGzUIgnAZHIyBuWGubJzFY3D5DJ6Q7SJiuYjYovbtJW3ILcElZ2N56HRvR7oSt5Smy21hs7rwq5NwxEsTn1oQN/tQ6fG3Mj9y40hmRKaPDBhmD180ZAUUTWc1XN9ffKigtSQ9bNzO9C9dOV0bBeo5QJT9UUnwVG8HXv6w21pzFP7jrD073twEHODd/84nSZqU65Vyg0JpUCkMSh2hNxAGLaFjY5aflt236w4+g7c4YYHN1xp9KPOYCGWF9t4OMSenqtLSFGhfzLDzeAzunOhps6Km/F175ffijE3Z2yaGjEkPG+8ncLRmN9ROqms7Un4qo/Q4n8l7LOLRd1NeY3XrjgZDgAvDeK03FJ6DHXb5w+4xqghtk+X7RbsHQzBXjrirVznd07cDYR8iiRG0FahgILxL0+U2z8Hibo/RQBE0xX9Iiv+QWlX9odLjz5943ZfvNTHkkdEBw7t3+x2yQ0bSeKHu6tHyU3kthaMChq0d8Wqka0/7wz0Hi2tPtcBAeJe2WypRBN9+VtO1GhgIrUQcJSj5pS4o3db5sCekkWrNVSa9eveqvt3gJ/BZmvjU4v4LLtdnHy0/vfnad0leCWODRg7xGWDDRYGhniAo8lpTzunKc+drL0dIQieGjFk7YrWlvk1xFL/01zptk4HrCU+Pf7QVqiQxfWyGkkXAQGglggCOSU0Y2kwscR+bm9V7WnMUwhDenbvE9RCGYMN8Bw3zHaQxaTOrL2SUnfjk0qaB3okj/YcO8R1g2VkxUC8xksasxhvnai79XXvF38VndMCIxf0XWrx/DEERj4Hi5ittQZO9LHvkvosmaXmxJnS6M95cgIHQWhAgjuK3Fai8hklsnRV70Xq9t27S8BjctNCxaaFjlQbV37WXT1ZmbriyOco1fIjvgCE+Azo/vBCymhZtq3k5p+tNNyMkoSn+Qx40I95SPJKE+VurgtK9+u7sSctSVmjZ7kyGwBmDgjN+ZlsRRwlarsthIDQjdKSqSmfeK7X3uLAE5oioJ/TZjbkX667+VvgnApCB3onJXglJXv363FI1jkRH6G805V1rzLnaeEOmaxvonZgaOGL10Bd6vlNgZ3C92RgLVVXDCYW3tRWqJF1cXNRhwEBoPeJoftlv9RRBd2MxLccjzVWKInjYPWuq9RI2zh7uN2i43yAAQJWi5mrDjeMVZz69/KUnz72/Z1w/99h49+i+spxpn6Y0qPJaC3Kbb+U051fIq6LdIpI9E14b8kKEJMxqs6fbuSYIW3MVMBCayW6pOz+d18HAQGg9OAfjerOUZRpRpAMuwN1VrblKjwEim7x1oNA/UOg/M2oyRVPFsrKc5vwTlWc2Xt3CwllxbpHRbpHRrhFh4mA4ysYiCIosl1cWSktutRbdkha3aqXRbhH93GOWJj4V7Rreeys6dYZbP5eC76qDYe8oAAa5yaQmBP5Oek0AA6FViaMFsgIVDISEjlRWaKIWPnBBSOtAETTKNTzKNXxO9DQAQI2yrkBafKu1+FRlZoWi2pvvGS4ODheHhIqDQ0RBfXqBU2vSmLTl8sqytsqStvISWXmVstaX7xXpGhbrHjUzakqIKNCaizp1jOfDRjBEXavjO2sAaNdWoBJH8Z32ggAGQqsShfNL99bZOhe2J8tTicL59rauh7+Lr7+L7/jgVAAAQZEViqoSWXlJW/nfdVfL2ypRBA0RBwa4+AW6+Ae4+Pq5+Hhw3a3fm2dvWnWyWlV9jbKuWllXKa+uVFSrTZogYUCoKChCEjopdGyoyK7b1m4JLq05ShgI5SUap71BCGAgtDK+P9vQZjKpCQbfqUu+NVfh3t+uG1g4ioWLQ8LF/+xRINW1VSqqKxU1lYrqc7UXa5X1coPSh+/pw/f2EXh68Ty9eB5efA8PrptDDsDRmLTNmpYGTXOTprlR3VyvbqxXN9Wp6tk4y0/gGyj083fxTfLsFywK8OS5995erBbnmiAs3FkdlO5p64zYFA0UpZpgJ55J4tTVsfUhKOISzFWUadwS7DoM9CrKRCnKNBFP+Nk6I11jXu0p2Suh/Rk9YahXN5r/NaibrjflNmlamzUtBtLgznVz40jcuW5ittCN6ypiuYjYQhFLKGILhSwBG7f8Yok9ZCSNSoNKblC26eVyg0KuV0p1MqmuTaqTSXWyZm0rAMCT6+7J8/DkuXvzPaPdInz4Xr4C774+O5PvywY0cPKZ9dpGPcZGnXmKMwyE1iYM5ylKnToQKko1fF8Ozrb3FWIfio2zQkSBIffs+6MnDM3aVplO1qKTtunkrTpZmbyyTS9v08kVBqXSqKJo2oXJFzD5fCafz+TyGDwOzuYzeRyczcSYfCaPiTJYOIuJMVkYE0NQzv9v7S1g/s/dZY1Wo0I0dz6jNWlJmjLngaAIgiJ0hN5EmfSEQWfSGUij1qTTmLQ6Qqc2arUmrcqoVhrVKqOKpCghSyBkuYjYQglbLGK7SDjiYGGAK0fixnV157r29YDXAVEUv61A5cyBUF6qEYU79cAFGAitTRTGL7xYY+tc2JKsQC2OduRfHRtnBbj4djBt30AaVQaVyqhWGdVqk1Zr0uoIvcqo1hN6uV5Rp2owkEYjaTSQBiNpImlKZ9IBAGhAq43/E/YoikLR/7nPymVwMQQFALBwFgPFMRTj4hwGymDjLA7OZuJMAYvvxffg4Gwug8tncPlMvoDJd2EJOPbXSLUaSbSgPlPqO7oXJ+/bOUWpxs2+b1X0NhgIrY3nwyY0hFFhYgqdtCOirVAVvah359HbORbGZHFdez5t0SH3w7M+YTiv6Kca0kDZ2+gtK6GBslwT6nC7pXaJU37xtoUAl1Ceokzz8JSOSNdsoAma5+W87Q/I3mBMVBDAlZeobZ2R/2vvzqOkqO6+gVevtXX3LMz0dE/PxjKAoDIIooivKIjiUSLGnbhETESRHNeco8EH85jkqMEnR01eNRLQY/ImKkZOBN83Go+iJIAIKksEhmWYtXv26aWqequu9482k8msPdPVXdv381dPU9N9bbvmW3Xv796rjEirYHNY7YZcWa0fglABhdPY4EmDBmHvsUjRWU7tFBWCIRSd5eg9atAgDJ7gCow9QEggCBVRUOvoO2HQIOxJz9sFUJPis5y9x8JKt0IZwVNcwTSjb96JIFQA4ybFWCrWl1C6IfkmxlPhRr5wOoIQ1IV2kyaziQ/ElG5I3klE6AxfMEW3JcEZQhAqwUQ4q+lwI690O/ItdIpzVNAGLUkAdSua4eg9bribQj4QtbEWg6/vQSAIleKazIYaDBeEwZNcoeE7YUCd0hN8lW5FvoUaeNdknJIIQoW4JjMGDMK+kxiWB5UqmMqGTvNSSlK6IXkVauCdNUbvFyUQhEpxVNJCR0yMpZRuSP4ko6LQEcPqxqBONofVXmDjWqJKNySvQmd412QEIYJQIWarifVSkWZB6YbkT+gU76xmsCkxqFZhLdtnpN7RRDiZ5EXGbdy15frJFoQ8zzc1NaVSI97iiKLY2NgYjRrrgmsUToP1jvadiBTWYjQC1KtgGhs8aaDZhMEG3jWZwaReQq4gfPHFF30+35IlS2bOnHn06NGhB+zfv3/KlCnLli0rLy9/4403ZHlTrXPVMKEGA11+Bk9iuhKoWsFUNnSGl0SjDBOGGzj0i6bJEISNjY0/+clPdu/effLkyVWrVj3wwANDj1mzZs0jjzxSX1//4Ycfrl27tru7O/v31TrXFDZ8RjDI4HySF6M9cUcFBghBvayMhZ5kN86ABSpl+skQhG+++eZll1121llnEQSxdu3ajz/+uL29feABR48e/eabb37wgx8QBDF//vw5c+Zs27Yt+/fVOhtrsTktfLsh5vAGT3KuyazJgl4YULWCaQ6DDBOmEik+EHOieI0gCFl2n2hoaJg+fXr6sdvtdjqdjY2NZWVlAw/w+XwM8+2lR21t7ZkzZ0Z6tVQqdejQoUAgkP5xypQpRUVFIx0siclES5No1mrJD1ua7P2q0ZbUavtFno8zGV1R9nyddEwyxZtP5LpJhpL55w8ZYgtTga/Eshl9mRys6c8/3CZRxVKy/ZTSDRk3k9li802R9zVlCMJwOFxS8u+tvFiWDQaDgw6gaXqUAwaKxWI//vGPbbZvtyh66qmnLr300pEO5jvaUn9+YcItVx5/dk9DieXITqXbMUGpVErI7Cqkr/umUudn3ScDuW6SoWT++UOGUhIZ6bi960//20SMPWah6c+/jzvXIhZ0/2mX0g0ZNxPtYL//XwRBRCIZVTYxDGOxjLENuAxB6Ha7e3t7+3/s7e0deDuYPqCvr2/gAbNmzRrp1Wia/uCDD3y+ETc1HchqrbI/8usx/yNVy9kknNja6n3kJqUbMkEZ7ocnxlNnNhyreuy/MHdCXtiPMBfanz1R8L3/YX1j7xSm6c8/9Idmz0yne/6tSjckK3J9/jJczsyZM2ffvn3px0eOHDGbzVOnTh14wKxZszo6OlpbW9M/7tu3b86cOdm/rw6wPiraGRfjOp9WH27kWR+FFARNMMiqT+FGwVmFAcJvyRCEN954Y1NT03PPPXfw4MEHH3zwrrvuYlmWIIjHHnvsmWeeIQiirKzshhtuWLdu3aFDh5588klJkq666qrs31cHTBYT4yG5Vp3PrQw3YPUK0AxnDRM6o/MgTPJikhfpUkyl/5YMQciy7EcffbR79+41a9bMnz//l7/8Zfr5iooKr9ebfvzyyy9XVVX98Ic/rK+v/+CDD6xWoy923s9RxYSbdH7Whc6gShs0w1XDhPUehOEmwVFBYyp9P3kC6Zxzznn33XcHPblu3br+xy6X64UXtFzVkjPOKlrnO4JKRLhJmL4KQQjaQJeSYjwV60uQhTal25IrkSbegX7RAbRa8qQbjko63KTnCbxcIGpjrdjwDDTDANuFhpsFzCAcCEGoMMZNJjkxyYlKNyRXwmd4F/pFQVN0Xy8TaRZwRzgQglBpJsJRQYX1u6pTqIF3olIGNMWl63qZWG+CkAgdd/xOAIJQeY5KRsfLG4ZwRwha46ikhUAspdN5TWHcDg6BIFSeo4rWa+FokhOTEZEpQ5U2aInZZmY8ZESn85oiGCAcAkGoPEcFpdephOFmwVGJKm3QHkcVrdd+mkj6rIQBEITKo4rsYjyVCCeVboj8wqjSBm1yVjJ6LefmWgXWh7PyPyAIVcBEOHyULvthIk1Yxgk0Sa8DFrHeBGE22V2YzvQfEISqwPporlWHl5+RFnTCgCbpdV5TpEXA/thDIQhVweGjIi16uyNElTZomIlgfVSkRW+Xp1xr1FEx9sYaRoMgVAW2go7o7o4w3CQ4qzFxArTKWaXDYcJIq+DAAOEQCEJVYNxkIpxMRnXVDxNpRqUMaJguhwkjLVEWd4RDIAjVwUSwXr1NogijUga0zFlFR/R1R5iIJFPxFFVkV7ohqoMgVAtWZ7MJJYJriaITBrSLLLQRJiLWl1C6IbKJtEQdFRTm9Q6FIFQLh4/W08g83xGzOS1W1qJ0QwAmzlGpq5tCrlVgUTI6HAShWrAVtJ4KRyPNqNIGzXNU0npaED/SGnWUY4BwGAhCtWA9ZLQnnkroZJ3fSAuuPUHzHBW6muCLO8KRIAjVwmQx0SV2PhBTuiHywHQl0AGHj9LNiqNiPBUPJulSVMoMA0GoImw5xbXpondUIrhWVMqA5tkLbCazSR/1MnxblPaQJjNKZYaBIFQR3QSh0BWzMhYrg0oZ0Dy2gtZHOTfXhgHCESEIVYQtpzi/Lk65VkzaBZ1w6GWhNa4tynhxVg4PQagiurkjjGAGIegFq5d5TVxblMUd4QgQhCpic1jNFj0MSERaBFTKgD7oZN9sieACCMIRIQjVhfXp4aYQO3+CblDFetg3O9oTt9IWK41h++EhCNWF9Wo+CLHzJ+iKLvbN5tqiLAYIR4YgVBcdDBNi50/QGR3sm40BwtEhCNWF0X4QYio96IwO9s1GEI4OQagujJuM9SU0vdBaBAOEoC862Deb8yMIR4MgVBeTxUSXanuhNa4t6vDhlAP9oEvtiVBSjGr18lSMpRKhJFWCxdVGhCBUHU3XyyQFMcmLVDFOOdAPk9lEe0g+oNWzkvdH6TIsrjYaBKHqsF5Ku6fct8VpOONAXzR9eYp+0TEhCFWH8VKcX6tdoxiTB11iy6mIZoOQD8QYD6l0K1QNQag6jIfkNbviKNcWZTFACLqj6XlNnB+TCMeAIFQdstCWEqVERJMrWXCtuCMEHWLLKd4fJSSl2zEhfCCG5bZHhyBUI8ZDabFwVEpJQkeM8eCUA72x0hYrY4l2x5VuyLjFw0lCkuxOrPQ0GgShGrEeUov7MQmdcXuB1ULiSwU6pNFhQt6P3ZfGhr9ZasR4NXlHiEoZ0DHWR/NaDMJADAOEY0IQqhHr1WS9DNcWZcuxpgzok0brZTh/FCWjY0IQqhHjpTgNjszjjhB0zKHNIETXaCYQhGpkpS0WyqK5HXq5VgFBCHpFTbInuGRSEJVuyHhIBN8eY8pwRzgGBKFKsV6N1cskOTGVkMhCm9INAcgNE8F4KW2NWUR74lYG+/GODUGoUpo75bh0DwwWVwP9Yr0Up6kqNkylzxCCUKVYjxZPOfTAgJ6xXlJbhaOYSp8hBKFKMVorHOUDmEoPOsd4KG0NWPABlIxmBEGoUoybFDrjUkozlaMoGQXdY8sp3h/TUDk374+xuDzNAIJQpcx2s91l1cySThLBt+PaE3TOyljMdpNWyrmllCR0x2k3NgcdG4JQvRgPqZX1ZaK9cSuN4jTQP7ZcM72j0a44WWA12/BHfmz4jNRLQ0tvf7sfL4Deaaicm/Nj2D5TCEL1Yj2kVraq5/0oTgNDYLWzbzYqZTKHIFQvLd0RYu4EGIOG7ghRyJ05BKF60WWk0K2NwlGsZwgGwaTPSlELZyXuCDOGIFQvs9VEFtiETrUXjqaSUrQ3wbhxyoH+ma0mstAmdKi9q0YSpWhPgi7FWZkRBKGqaaJwVGiPUcV2kwWrq4EhaGKYUOiMk0U2sxVnZUYQhKrGaKFehgtggBAMhPFq4KzkA1EW/aIZQxCqmibqZVAyCobCamGhNQ6rjI4HglDVNHFHiDF5MBTGq4EBC5yV44IgVDXGTUZ7EiovUcNWL2Ao1CR7IpwUYymlGzIazJ0YFwShqpksJqpY1SVqYiyV5EWqGOsZglGYzCbaTfLt6j0rU0kp1pugS3BWZgpBqHZMmapPOd4fpctI7McLhsJ4VL1LmtARoybZUMidOQSh2tEeSs1ByAWwzwsYjsqr2Pj2GFOGs3IcEIRqx5SpemQeY/JgQKyX5NvVe0fIt8dwVo4LglDtVF44ijF5MCDGo+o59bg8HS8EodqpvHCUx3LbYDxkoS0VT4mCSgtH+QC6RscHQah2JouJLLIJXWpccTTJp1JJyV5gU7ohAPllIpgyMtqhxrNSEqVYb4IuRcnoOCAINUC1vaPR9jhWrwBjYryU0JFQuhXD4DtiVDFKRscHQagBqq2XiXUmsJ4hGBPjIWOdagxCAZUy44cg1ACmTKUzKISOBCplwJgYDxVtV2MQon5tAhCEGqDqrlFce4IhsV5SUOUYId8eZcpwVo6PbEEoCML+/fubmppGOqCrq+vAgQMtLS1yvaNx0G4y2q3GwtFoZwJBCMZkc1hNJlM8lFS6IYPxgRiCcLzkCcIDBw5MnTr1wQcfPP/88x9//PGhB9xyyy3Tp09fu3bt3Llzr7vuunhcjVdSqvXtptgqKxyNBxNmi8nmsCrdEABlUG6b2rpq0hvTU9iYfpzkCcJHH330gQce+Pvf//7ll1++/PLL33zzzaADVq9eHQgEPv/889OnTx85cuS1116T5X2NQ4Vb1XP+GFWGEm0wLrLUprazUuiIUcXYmH7cZAjC9vb2Tz/99O677yYIwufzLV++/O233x50zBVXXGG32wmCcDqds2fPbm9vz/59DUWFw4R8e5Qswe0gGBddZlPbDr18O/pFJ0KGP2TNzc1Op7OkpCT94+TJk0cZKTx58uSnn37605/+dKQDRFHctWtX/6ude+65brc7+0ZqHVNGdR8JKd2K/8D7Y5QHd4RgXJTbHjyssrMSJaMTkmkQPv3004cPHx705Pz58x9++GGe59N3e2kURXEcN+yL9PT0XH/99Y8++mhdXd1IbxSPx3/961+T5LcXNY899thFF1000sGCINjtdovFkuF/hYa5xEibEIlElG7Hv4VbueJaWlVNMhqO40wmdIIpJuVIcG1CJBxRzzZkodZI4WzWIGdlht9/hmHM5jH6PjMNwkWLFk2dOnXQkz6fjyCIsrKyvr6+VCqVfrPu7m6PxzP0FYLB4PLlyy+//PL169eP8kY0Tb/99tvpVx6TxWIxSBAyNdKJ3gBLs2pZMEIiYl3JwmqXw+FQuinGJUkSPn8FSZJkpUO2JEkWqWWVwXiXv7i6gHEY4qZQxu9/pkF4ySWXjPRPNTU1hYWFe/fuTd+6/eMf/3jooYcGHcNx3He+8526urrnnntuwm01sv7CUZUMAER74lbGYiExDxUMjfWSfCCqkiBEyeiEyfCHjCTJtWvXrlu37q9//evjjz/e2dl5ww03EASxe/fuqqqq9DE33HBDU1PTvHnzNm3a9Oqrr3722WfZv6/RqKpwFEMRAARBMGUq2o8JJaMTJk/V34YNGyZNmvTSSy/5fL5PP/2UoiiCIMrKym6++eb0AYsXL54zZ05DQ0P6x/5aGMgcU5beC9SldEMIgiD4AHZfAiAYDxk8NXxJRP6hZHTCTJKkrvVKKioqPv/88wzHCA1ULEMQnV8Gu4+EZt5RqXRDCIIgjv+flqIZDnqGxel0Kt0W4wqHw/j8FRQOh4le66l32uoeHlw/oYimv3YQBFG13Chl9jJ+/zHGoxnq6hr1YwtsAILxkHxHTEqp4naCb8dZOUEIQs2g3WS0O66GFUellCR0xmk3TjkwOovdbHdYYz2q2IYCq4xOGIJQM8xWtWxVH+2K2wusFju+PAAE4yXVsL5MumQUl6cTg79lWqKSHXo5f4zFxvQABEEQBOOh1HBWCtiYPgsIQi35V+GowvhAlEEQAhAEQRCsRxV3hDw2ps8CglBLVHLtyfujLE45AIIgCILxUmpYEJ8PxJgyXJ5OEIJQS1SyBwWH2fQA/5LeNzuVVLiKDSWj2UAQaokatqpPJaVYb4Iuxb4TAARBEGariSq2CZ0Kd9VwflyeThyCUEu+LRztVLJwlA9E6VI7xuQB+jFeild0mBCXp1lCEGpMepFfBRuAVUYBBlF8sQuhI0ZNwuXpxCEINYbxUJyipxxWGQUYhPVSyhaO8n6clVlBEGoM4yGV7YTBUATAIIqXc6N+LUsIQo1hPArXavP+KINrT4AB6BJ7IpIUYymlGsAHUDKaFQShxtCl9lgwmUooc8olo2JSEKkijMkDDGAi6FIlhwl59NNkB0GoMSaziS6x8+3KnHK8P8Z4SAJD8gD/iVGuik2MpxKRJDXJpsi76wOCUHsUHJDgA1GsMgowFKvcmAUfiNFu0mTG9enEIQi1h/UqVi+DShmAYTFeimtT7PIUA4RZQhBqD+OhlOoa5dqibDmCEGAwtpzi2gRF3hoDhNlDEGoPo9xq99iYHmBYdpeVMJni4WT+35pvxyTCbCEItYcqtid5MRkV8/y+sd6E2WayOax5fl8ATWC9JN+mwBUqBiyyhyDUIBNBl5FC3utlOD/6RQFGxHgpLu9BmBREMSaShSgZzQqCUJMUWdKJa0PJKMCIlDwrUTGaHQShJilyyvF+bEwPMCK2HJenWoUg1CS2XIFabXSNAoyC8ZBCZzzP24Xy/iiDszJrCEJNYsspvi1K5PGMSyWlaHeCdqM4DWB4ZpuZLMz3dqGcH3eEMkAQapKVsZjtplhfIm/vyLfH6BK72YqxCIARsXmul5HS+4Pi8jRbCEKtYsvzesrxbRggBBgD683rHN9od9zGWq20JW/vqFcIQq3K88g8h50/AcbClFP5XP6Qa8OeaPJAEGpVnictoVIGYEx57hrFWSkXBKFWsd68XnuiaxRgTFSxPRkVk3yeVn3iUSkjEwShVjFlZLQnkZ8deuOhZColYfUKgDGY8npTGMHlqUwQhFplspjo0jzt0Mu1Cg4fnYc3AtA6RwUdac3HNhRiPJUIJelSex7eS/cQhBqWt/VlIq1Rhw8XngBjY8sprjUfZyXvx368skEQahjjpfKz2j3XijF5gIywPirSkpezEgOE8kEQapjDR0Xycu0ZaRXYCnSNAoyN9VLRnngqmfNln7hWgUU/jUwQhBrGVtBca84XWhOjGIoAyJTJYqJL7Hmo6I60RhGEckEQapiNtVhIc7Q3t2sbcm0C46UwFAGQIYePznlXjUTwfgShbBCE2sb6KC7HAxK48AQYF9ZHcTkuHOU7Yjan1UphcTV5IAi1zeHLea0214aSUYBxYH05LxzFjCZ5IQi1LR+nXIvA4pQDyJjDR3P+3A7eR1qjjgpcnsoGQahtjgo60pLDO0JJlPjOOIt9XgAyZqHMNodV6MzhYhe4PJUXglDbyCJbSpTi4WSOXp+0bceiAAASp0lEQVRvj1FFNrMd3xOAccj11KZIG0bu5YQ/cJrnKKe5nN0URpoFRyUuPAHGx1GZw66aaE/cbDHZndYcvb4BIQg1z1GRw2vPcJPgrEIQAoyPo5KONOUqCLnWqAMLXMgKQah5rI/OXa12pFlwVDI5enEAvXJU0ZFWQUrlpGAm0hrFSk/yQhBqnqMiV2sbppKS0BFjy1EpAzA+Vspid1qFzpwsdsG1CJjRJC8EoebRpWSCSyYF+fcC5dqidKndbMOXBGDcHJVMjnpHIy0Cukblhb9x2mciHBU5GZCINPGOKvSLAkyEs4oON/Oyv2ysL0FIBFmEXbLlhCDUA2c1E26U/5RDpQzAhOWoXibcKDircXkqMwShHjir6HAu7ggxdwJgohw+igvEZN+PKdLEO3B5KjcEoR44qxnZg1CMpWJ9CQZrygBMiNlupkvsfEDmQrZwk+CsRhDKDEGoB3aX1Ww1RXvkLFGLNAtsOXZfApg42XtHpZSESplcQBDqhKOKDjfKecqF0S8KkB1nFR1ulvOs5AMxstBmpbH7kswQhDrhrGIiTXLWy0SaeFTKAGTDUSVzFVu4ScAAYS4gCHXCKfcdYaiBd01mZXxBAKNhvWSsL5HkZJvjG27kUTKaCwhCnXBU0Zw/KonylKgJnXHCZMJcJYBsmMwmZxUTOiPbTWEEM5pyA0GoExa7mSq2c355StRCDVzBVFx4AmTLNYUJNcgThGIsFe2Os14sriY/BKF+OKvpsEzXnuEzvKsG/aIA2XJNZkOnOVleKtwksD7KZEEht/wQhPrhmsoGT8kThMHTvGsK7ggBsuWqpjl/VJZp9aFTXMFUXJ7mBIJQPwqmssFTHJH1GZfgxEQoian0ANkz281MGSlLRXcQQZgzCEL9IAttFruZ74hl+Tqh05xzMoOp9ACycE1msx8mTCWlSIvgrEE/TU4gCHWlYCobOpXtgESogXdNxvkGIA/XZCZ4OtsgDDfyjIe0kPiLnRP4WHXFNY0JyhCEHIIQQC6uKUy4gc9yt3r0i+YUglBXvh0mzIIYS/GBGOYqAcjF5rDaXFben9WYRegUjyDMHQShrlDFdpPFJHROfPXt4EnOWc1gV3oAGRXNcPQej0z41yVRCjfxTvTT5Az+3ulNljeFvcfDRTMdMrYHAIpmOvqOhyf86+EmgXaTVgprbecKglBvsg3CY5HCGQhCADm5prLhZkGMpSb268FTXMEU9IvmEIJQbwpnOHqPhSc2Mh/tiafiKdaDNZwA5GSxm52VEy9k6/0mjMvTnJInCFOp1J49e95///2+vr5RDovFYgcOHOju7pblTWFYZKGNLLBNbCeK3qORoplOAhMIAeRWONPRe2wiw4RJTuQDsYJpuCPMIRmCUBTFa6655p577nn11Vdnzpx5+PDhkY7csGHDggULduzYkf2bwiiKZzl7vpnIgETvMQwQAuRE0YwJDhP2HA0XTGfNVlyf5pAMQbhjx46TJ0/u27fvL3/5yz333PPkk08Oe9i+fft27do1d+7c7N8RRlc829Xzz9B4f0sSpdBpHheeALnAeqn09hHj/cWef4aKZztz0SToJ0MQbtu27brrrqNpmiCIVatWbd++PZlMDjomHo/fe++9L730ksWCwqecc1bRCU6M9ozvlAs18HSp3eaw5qhVAIZmIopmOMfbOyqJUl89VzwTQZhbMvzVa2lpmT9/fvpxVVVVMpkMBAIVFRUDj/n5z3++YsWKurq6MV8tmUy+9957xcXF6R8vvPDCQS81kCiKoijb7s96UjTT0XU46L24OPNf6TzYVzjLMa7PE5+/svD5K2u8n3/hLNb/9x73hQWZ/0rwBEeV2s2MCf+jh8rw8zebzSbTGB3LGQXhwYMHH3744aHPb9mypbq6Oh6P22zfbmWefhCL/ccaCocOHdq2bdsXX3yRyXuJorh9+/b0/SVBEGVlZaWlpSMdHIvFJEnCXeZQjlqy8/NQ8fmZ9nNKKanrYGjGmvJB/+9GF4/Hx3U8yAufv7LG+/kzU2z81mi4g7MXZHoH0nU46JpB4//ysDL8/CmKkicIa2pqNmzYMPT5dER5PJ6urq70M52dnQRBeL3egYf94he/8Pl8P/3pTwmCaG5u3rp1q9Pp/O53vzvse5EkuWnTJp/Pl0nDTCaT3W5HEA5FnkM1vdtlN5FWOqMPp/domC4hi3zjuFYlCEIURYbBaheKweevrAl8/pPOKeCOxQsvc2V0tESEjrecdXcVw2BG0zBk/P5nFIQFBQWLFy8e6V8XLVr03nvvrV+/niCITz75ZO7cuYMat3r16sbGxvRjkiRLS0tLSkqyaDOMzUKaC2c6Or8Mehdl1Dva+WXQfd74UhAAxqv0vMIz2wO+yzL6Axg8xVlIM+tFCuacDGOEd9xxx7PPPvvII4/MmjXriSeeeP7559PPL1q06Lrrrnv00UevvPLK/oM3b9586aWXXnLJJdm/L4zOc2FRw3uBTIIwFU/1fBOe/B1PHloFYGSF09h4OMm3x5iysTe+DuztLbuwKA+tAhmqRouKivbu3Wu32/fv379ly5abb745/fx99903NPDuv//+efPmZf+mMKbCWocYT0Wax55Z330k7KymbU7UiwLkmIkorXN1fRUc88AkJ/YeDbvnF+ahUSDP376qqqqnn3560JO33Xbb0CPvuOMOWd4RxmYiPBcUBfb2TqscY0+ljv29pefhfAPIh9J5hcdeb668otRkHq2Co+NAX/FsZ4Zj/JAlrDWqZ+4FRV1fB0df6jfSIvDtsZI6DBAC5IOjgiaLbZ1fjnFTGNjb40G/aL4gCPXM7rQWTGM79o+2AGzTXzsqlpZiASeAvKm60t38t45RVsYPnuSkFOGajGWe8gRBqHNVy91NH3YkuOGnnUZaBM4fLVuAC0+A/CmYytoLRrwplFLS6b/4q5e7sfx93iAIdY71UiVzChr/X/uw/4rbQQBFVF9V1vTh8DeF/r/32FgrRivyCUGof9VXuXuOhCJNg8tH2/f18h0x3A4C5J9rMkOX2BvfH3yFGg8nmz/qnPpd77C/BTmCINQ/K22pucZT/6eWgSvf93wTbvy/7bN/WI3bQQBFzLitsuebcOunXf3PJDmx/g/NZRcU0e6xZxmCjDB1zBDc8wuTgnjwhdNTVnrIInukWWj+qHPWD6rpUpxvAMqwMpbZa2oO/fp0UkgVTGUIiTjxVmvJnILq5W6lm2Y4CEKjKP9fk1yTmRNvtpqtJtZHn7W6ylk1xvxCAMgpstB29poa/z96mj/sjAUTU68vL56FHZcUgCA0EEcFPffRaUq3AgD+jXaTU67DiKDCMEYIAACGhiAEAABDQxACAIChaTsIu7q6wuGw0q0wLr/fH41GlW6FcTU3N4vi8GsGQa6lUqmmpialW2FcsVisra1NrlfTdhCuX7/+3XffVboVxnXnnXceOHBA6VYY1+LFi3t6epRuhUEFg8GLL75Y6VYY19dffz3sBkcTo+0gBAAAyBKCEAAADA1BCAAAhmaSpBH3xFLEtddee+jQIbM5o4QOh8M2m42iqFy3CobV19fHsqzNZlO6IQbV3d1dVFSU4ckC8pIkqbu7u6SkROmGGFQymQyHw0VFY+8ZsGPHjrPOOmv0Y1QXhOFwuLOzU+lWAACAHlRUVNjt9tGPUV0QAgAA5BM6VQAAwNAQhAAAYGgIQgAAMDQEIQAAGJpm9iPs7OzcvHlzZ2fn1VdfvWTJkqEHiKL4xhtvHD58eObMmXfddRdq+mXU2dm5ffv2o0ePFhYW3nTTTbW1tYMOSKVSv/vd7/p/PPvssy+66KL8tlHPdu7cWV9fn35stVpXr1499Ji2trYtW7b09fWtXLkSS3/Ja9OmTQOLCod+vevr63fu3Nn/48qVK91u7DKflY6OjgMHDjQ3N1966aXTp0/vf76lpeW1114LhULXX3/9hRdeOPQX4/H45s2bT5w4UVdXd9ttt2U4uUgbd4Q8zy9cuPD48ePV1dWrVq166623hh5z//33v/TSS7W1tX/4wx++//3v572NerZu3boPPvjA6/V2dHTU1dXt3r170AHJZHLNmjXHjh07ffr06dOnu7u7FWmnXr3++utvvfVW+rNtaGgYekBfX9+CBQtaW1srKiquvfba999/P/+N1LGGhobT//LQQw8dPnx40AF79uzZuHFj/zGxWEyRdurJ0qVLn3rqqccff3zPnj39T3Z1dZ1//vmdnZ1er/eqq67629/+NvQXb7311rfeequ2tvaFF154+OGHM30/SQu2bNly/vnnp1IpSZL++Mc/nnvuuYMO8Pv9JEk2NzdLktTT00NR1KlTpxRoqE4JgtD/eM2aNatXrx50QPrMD4fD+W2XUdx5552/+tWvRjng+eefX7JkSfrxK6+8smjRory0y3D2799P03Rvb++g519//fUVK1Yo0iS9EkVRkqQLLrjg9ddf73/y2Wefveqqq9KPX3jhhf7vfL+jR48yDBMMBiVJamhooGm6q6srk7fTxh3hZ599tmzZMpPJRBDEsmXLDh061NvbO/CAPXv2TJs2raKigiCIoqKi8847b9euXcq0VY8Grt0TjUYdDsewh/32t7998cUXv/zyy3y1y0D27t27cePGrVu3JhKJof/62WefXXHFFenHy5Yt27Nnz7CHQZY2b9584403FhYWDv2nlpaWjRs3btmyBeuByGLYLs1B3/Ndu3alUqlBByxYsMDlchEEUVNTU1lZ+fnnn2f0dlk3OB/8fn9paWn68aRJk6xWq9/vH+kAgiDKyspk3KoK+u3Zs2fbtm0/+tGPhv7T0qVLu7u7jx49umTJko0bN+a/bTpWVVVVWlra09Pz9NNPL1iwgOf5QQcM/P673e5UKhUIBPLeTJ0TBOHNN98cdoC2oKDg7LPPDgaD27Ztmzlz5tC+U5DFoO95IpHo6uoaeEAgEBgYBG63O8Mg0EaxjNVqTSaT6ceiKIqiOGjJHKvVOnCH0kQiMeaaOjBex48fv/HGGzdt2jRt2rRB/2S32z/66KP041tvvfXyyy9fu3Yty7J5b6M+PfXUU/0PzjvvvC1btqxbt27gAQNPkPQDfP9l9+c//7moqOiSSy4Z+k8rV65cuXJl+vF999333//93++8805+W2cIY37PJxwE2rgj9Pl8/cGefuD1egcd0Nra2v9ja2treXl5Pluoe/X19UuXLn3mmWduuumm0Y9cuHBhMplsaWnJT8MMxWazLViwYGi9zMATpLW11WazYTFo2W3ZsuXuu+9OD9CM4qKLLjp9+nR+mmQ0g77nDMMM6qaecBBoIwhXrFixffv2aDRKEMQ777yzZMmS9N3G119/3djYSBDE4sWLu7q60rulnzhx4tixY/1dyZC9xsbGK6+8csOGDYO2hP7iiy/SXztBEPqf3LFjB8MwNTU1eW6kjvV/vOFweOfOnbNnzyYIIh6Pf/zxx+kypRUrVmzbti09Lrh169arr77aYrEo2GD9aWho2LVr1+23397/DM/zH3/8cfq+pP9/kCRJ77///tlnn61MK/VuxYoV7777bvqeb+vWrStWrEg//8UXX6QD8sorrzx48GD6SnHv3r08zy9atCijl5alwifXksnkFVdcMW/evNtvv33SpEm7d+9OP79kyZKf/exn6cfPP/+81+tdvXp1ZWVl/5Mgi+XLl7MsO+9f7r333vTzc+fO/c1vfiNJ0qZNm2bPnv29731v+fLlLpfr97//vaLt1Zvi4uIVK1asWrXK6/Vec801iURCkqT0mX/mzBlJkqLR6MUXX7xw4cJVq1aVlJR89dVXSjdZb9avX3/11VcPfObYsWMEQXR3d0uStHz58qVLl952223nnHNObW1tY2OjQs3Uj4ceemjevHksy9bU1MybNy/9N5/n+QsuuODiiy++5ZZb3G73kSNH0gefe+65r7zySvrxE088UV1dvXr1ao/H8/LLL2f4dprZfUIUxZ07d3Z1dS1evNjj8aSfPH78uNPp7L/5PXLkSHpC/dy5c5VrqQ7V19eHw+H+H51OZ3qK6z//+c/S0tL0qPWBAwcaGhoKCgrmz5+P2cTyOnPmzNdffx2Px6dPn15XV5d+MpFIfPXVV3V1delRkEQi8cknn/T19V122WUD6wVAFvX19S6Xq/8vD0EQ0Wj00KFD8+bNs1gsPT09+/bt6+3t9fl8CxcuxGoe2Tt16lRfX1//j7W1tela0HRHSCgUWrp06aRJk9L/euTIkbKysv6v/YEDB+rr6+fMmTNr1qwM304zQQgAAJAL2hgjBAAAyBEEIQAAGBqCEAAADA1BCAAAhoYgBAAAQ0MQAgCAoSEIAQDA0BCEAABgaAhCAAAwNAQhAAAYGoIQAAAM7f8DIsVJPg9VEpIAAAAASUVORK5CYII=", "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.2" }, "kernelspec": { "name": "julia-1.7", "display_name": "Julia 1.7.2", "language": "julia" } }, "nbformat": 4 }