{ "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 Δtime\n", "--- --------------- --------- --------- ---- ------\n", " 1 -0.143384299049 -0.42 8.0 \n", " 2 -0.156027075014 -1.90 -1.10 1.0 2.37ms\n", " 3 -0.156771942256 -3.13 -1.57 1.0 1.92ms\n", " 4 -0.157045929883 -3.56 -2.32 1.0 1.99ms\n", " 5 -0.157054237732 -5.08 -2.79 1.0 1.92ms\n", " 6 -0.157056383897 -5.67 -3.62 1.0 1.88ms\n", " 7 -0.157056406567 -7.64 -4.49 1.0 1.95ms\n", " 8 -0.157056406872 -9.52 -4.97 1.0 2.09ms\n", " 9 -0.157056406917 -10.35 -5.78 2.0 2.90ms\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380294 \n AtomicLocal -0.3163464\n LocalNonlinearity 0.1212606 \n\n total -0.157056406917" }, "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-5, ρ)\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.05568274704827848, 0.0, 0.0]\n [0.05568162010102383, 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+gvaeTAAAgAElEQVR4nOzdd3hUVd4H8HPL9D7pnfSEEEIJvYeOEVBABAFFZXGxodh1lV3XXRuLimXtvRdAeq8CoSchpPdepveZW94/kpdFIBCSSe6U3+fx2SdhLzOHM3Pv995TMZZlEQAAAOCvcK4LAAAAAHAJghAAAIBfgyAEAADg1yAIAQAA+DUIQgAAAH4NghAAAIBfgyAEAADg1yAIAQAA+DUIQgAAAH4NghAAAIBf87ggfP755202WxcPpigKlojjkMvl4roIfg3qn1tQ/9xyY/17XBB++eWXWq22iwe7XC6GYXq1POA67HY710Xwa1D/3IL655Yb69/jghAAAADoSxCEAAAA/BoEIQAAAL8GQQgAAMCvQRACAADwaxCEAAAA/BoEIQAAAL8GQQgAAMCvkVwXAPQuF4OKDWyJga00IYpBTgYp+ChYiBIV2EA1xoMbIQC4w7CoUM8WG9hGKzI4EYEhDEPREpSkwFKUmBguz30Fato3NdnQD+XM7nrmjyY2TIylKLE4GeITiI+jciM61oyK9EyZkR0UgM3th8/vh/WTYVwXGQB/0WJDG6uZ3yqZEy1sqBgboMLCxEjBRzSLaAadbkUlBqbcxA4OwKZG4HfGYYkKOD17FwShr9lTz75TQP/RzN4Wg9+fjH87EVcJrn2klULHmtlfq5gRv9MD1dhjA4iZURiccAD0nhMt7PoLzO46ZkYU/pcU/Ies652eR5rYnXXMuK10ogJ7MBVfEIcTcH72DszTFq2OjIzMycmJiIjoysE2m43P5xME0dul8go7atm1Z2mzCz05EF8Qh0u6fJPjYtBPFcx/LjAuBr0+nJgR2dWzzWQyyWSybhYX9BjUP7duqv7PtrFPnqSrTOjRAfjyJFzG6+q7uBi0vZZZl8802dBzGfiyRByHOEQIufX7D0HoC4r07GMn6CozeiUTnxvT/fNkczXzzCkmWoLeG0MkyG/8KnAh5hbUP7e6WP9tdrQmh95bz740BL8vuftPdYeb2OdP01YKvT2SGBsKYejO7z8MlvBuLgb94xwzfis1PRLPu528vV+P7hbnxOD5t5Mzo/BRv1Pr8hnas+6RAPA+P5QzA39zBQtR8QLyLyk9atscH4odziafTMcXH6D/cpQ2wh5Q7gNB6MUu6Njhm6mcFubcbeTqAbhbhoCSOFo9AM+ZQ26rYbK2UfUWCEMAusPkQosP0K+cZzZPJd8YQUi73BZ6HRhCd8bjF+aTGEIDf6X2NcDp6R7uCcJvvvkmMTExODj4vvvuu862utu3b8/MzNy6datb3tTPfVLMZG2jHknDt00nIyRubieJk2F7Z5HTIvHMTdSuOjjZALg5uVo2cxMl5aGTc8hhQW4+PeU89OFY4uNxxN2H6L+doaHlpufcEISFhYUPPvjgF198UVxcXF1d/fLLL1/zMIPB8PTTT2u1Wo1G0/M39WdWCi05SG8oYA5nk8uTeuuZHsfQ84PwHyeT9x2h1+XD7scAdNVvVcy0HdTaIfhHYwlRrw3MnxqBnZ1LnmxhJ2+nmjt9+gBd4obL6GeffTZ37twxY8aoVKoXXnjh008/veYAnDVr1qxevTooKKjn7+jPaszsuK0UgaETs8kUZa93mI8PxY7PJr4tY+49TDshDQG4kX+eYx47weycQS6K7/WOp2AR2jGDnBiGjdhMnWmDB8Puc8NHVVRUlJGR0f5zRkZGS0uLTqe74pj9+/eXl5ffe++9PX87f3asmR35O7UkAf9yQi/eaV4hSoIduZXUOVD2LsoE/fMAdIJm0V+O0r/XMDlzyMEBfTSqE8fQ2iHE+pH4rF3UTxVwr9pNbriaarXaS2NY5XI5QqitrU2tVl86wGKxPPLII5s2bcK6MF27tbU1MjLy0q/Z2dnff/99Zwf71fSJTbX4Y6fJj0ZSU8McZnNfv/sXI9Hjp8kJW5hfx7uChB33nhaLpSufKeglUP/curz+bRRafpznYNCWCZSEtvfxGTo1EG2egN1xhF+ida5Oofr0vbnTxe+/UCgkyRsknRuCUK1WG43G9p/bfwgICLj8gLVr1w4bNsxgMJw5c8ZisVRVVVVXV8fExFzz1YKCgro+j5AgCD8JwvUXmPX5zN5biAw1n6syfDoJrT1LzzqI751JtA/PYVlWKpVyVR4A9c+tS/VvdqGFh6hwCfb5eIKHd7JUTC8bKUXHZ7PZu7EGB+/d0YQ/TLp34/ffDU2jiYmJBQUF7T9fuHAhICDg8sdBhBDDMPn5+StXrly5cmVVVdUXX3zx9ddf9/x9/QSL0LOn6E+LmWOziQw1x9/utUOI+5LxCdvoajN0SACAEEIGJ5qxk0pQYF9NILhdxT5Cgh3OJkuN7J37aQfNZUm8jhtWlsnPzx83btzBgweTkpIWLFjQv3//N954AyH08ssvZ2Zmzpw58/KDR4wYsWrVqrvvvruzV4OVZS5Hs2jlUfqijt06nVRzc6N5De9eZNblMwdmEQHIDCubcAhWluGWyWRiBLJpO6gRwdjbozxlHVAng5YepDV2dtNU0i2TFz2WZ60sk56e/sYbb2RnZ4eHh8vl8pdeeqn9z4uLi5ubm684OCIiAk7dLnIyaNF+utbM7pnlQSmIEHqoP/5kOp61na6zesi5DwAHzBQ2cyc10pNSECHEx9H3k4gEOTZlB6VzcF0aLwFrjXooG4Xm76P4OPZDFiHwyH/fhgJmfT51+FZepLun84MugidCDlkoNGWrY0gw+e5oD0rBS1iEnsqh99Szu2aSISKuS9M7POuJELidhULZuymVAPt5soemIELo4TT8/gR66g66BSbzAj9jp9HcPVSSnPXMFEQIYQi9MYK4PRafuI1qsHrW044HgiD0OBYK3bqLipVhX00gSM/+fB5OphbGYVN3UFpogQF+g2bR0oO0nIe9nenyzBS85MXB+H3J+NgtdKUJsvB6PPtC63+MLjR1O5WixD4e5x0DoNcOISaHY7N3U1Z/mbwE/BqL0D2HaDvN/pBFeOjD4J89kY4/moZP3k5XQRZ2DoLQgxicaPoOaqAae2+MV5xiHdaNJFKU2Jw9FKzBBnzeEzl0uZH9IYvkdqbETXl0AP5sBj5xG11uhCy8Nu/5MH2dzoGm7qBGBmMfjPWmFEQIYQj9dwwhIrD7DsM6+MCXvXKe2VfPbp9BSvpqgUN3WZGCPzsIn7wdsvDaIAg9QnsKjgvF1o/0shRsR+Loxyyiysw+mQPzeIFv+rSY+ayY2TmTVHK2uFOPrEzBnx+EZ22nyyALrwJByD2tA03ZQU0Kw9aN8NQRol0gItHvU8mddSzs2QR8z5Ya5m9n6J0ziFBvnoqwIgV/cTCetY0uNUAW/om3PeH7HI0DTd1OTY3AXhvuxSnYTiVAO2cQY7bQYWK0uPf3oAGgbxxvYe8/Qm+dRiYqvLG95k/uS8YxDE3eTu+dRSR5/z/HXSAIuaRzoBk7qCk+kYLtIiXYjhlE1jYqUIhNi4DTDHi9MiM7fy/95QT3bzTPlXuTcBGBJm2jd88k0lQ+8o/qIbht50yLDU3cRk2PxF73lRRs11+J/TSZXHqQytVC8wvwbo1WNG0H/UomPiPSpwJjUTz+6nB8+k76oh5OUoQgCLnSZEOTtlG39cP+melTKdhufCj23mgiexddA5tUAK9lcqFbdlH3JeP3JPngdXJpAv76cHzqdjofblihaZQT9RZ28nZ6WSL+3CAfPMHazY/F6y1o5k766K2kypNWDAegK1wMWrCPGh6EPe+7J+nieJzE0PSd1Lbp5OAAn3rkvVk++xl7rCoTO2EbfX+KL6dgu0cH4DOjsLl7KDtMqQBehUVoxRGaj2PvjfHBBpvL3RGHvzeamLmTymnx6+dCH78We5pSAztxG/3YAPyJdL+o+deHE+ESbOlBmvHrswx4mb+dposMXrOIWg/d1g//bDw5ew91pMl/z1K/uBx7iDwtO2k7/dIQ/MH+/lLtOIa+GE+02dnHYaI98BL/LWR+qmS3TCPFftNxNCsK+34SOX8ftbPOT7PQX67InMtpYaftoNaPxJf7Ysf7dQgItHEqub+BfRMm2gOP93s18/I5ZucMIkjIdVH6VlY4tnkqec8h6pdKfzxP/euizJU99ezsPdTn48kFsf5Y4Uo+2jGdeLeA+a7cH88x4C2ONbMrjtK/TyPiZH7QJHqVkcHY7pnk6hPMx0V+d57643W5j22sYpYdpH6dTM6M8sezq12EBNs+g3j8BL3LX9tegIcrM7IL9tFfTSCHBvrveTpQjR3JJt7IZ5455V99GRCEvevtC8yjx5k9s8ixof57drXrr8R+nUIuO0SdboMsBJ6lzsJO2U6/Nhyf7lsT57shVoYdziZ31bEPH6P9ZzcZCMLewiK09iz9URFz5FZiAKxjhBBCaEwI9uk48tZdVDGs+Qs8hsGJsnfRD6fhSxLgeogQQqEidCSbLDOy8/bSfrLhNnzwvcJGoYX76AMN7NFbyRgppOD/ZEdjrwwjZu6kG6yQhYB7VgrdsouaGoGt8Y8ZTV0k5aHN00g5D03eTjXbuC5N74PP3v2abGjiNkpAoN0zYVGVa7g3CV+Zgk/fQescXBcF+Lf25WMS5Njr3rwDWi/h4+jLicT0SGzU71SBzsdvWyEI3eychh25mbolGv9qIiGAk6sTT2fgMyKxW3ZRFv9oeAEeiGHRPYdoEsc+GecP8+a7A0No7RDi5aF41nZqW60vZyEEoTv9VMHM2Em9OQJ/cTAOp9b1vT6CSFVit++hHP41PA14ioeP0/VW9scsgoSr4HXdlYBvnko+cJR+LddnV4iCr4B7UAx6+iT9zClmz0xyvl9OFrxZGEIfjSPkfGzxAZryu2lLgGPPnqJPtbK/TyOF0GzTBSODsZw5xG9VzJ37abOL69L0Arhku0GzDU3dQeXr2FNzyYFqeBTsKgJD304ibDR772FYjBT0nX+dZ7bWsDtmkHIe10XxHuFi7HA2qeKj4ZupQp/bxRCCsKf21rNDN1ETw/Ct08gAGBpzk/g4+nUyWWthH/jDf+YsAS6tv8B8Wcrsngln600TEOi/Y4mnM/CJ26gvSnyqGQeCsPtcDHr+NL38MP3VBOKlITj0CnaPiERbppEXdeyjx6G3EPSu9y8y7xYw+2YRYWKui+K17k7E988i1+UzSw7SRl9pJoUg7KYiPTv6dypPy569jcwKhwzsESkPbZtO5rSwj5+ALAS95aMi5vU8Zt8sIlICJ2yPpKmwk3NIGQ8N+o067BObN0EQ3jSaResvMOO3UitS8C3TSH9bpb6XKPho10zySBO7BjZsAr3g4yLmlfPMvllEP79cUNvtRCT6YAzx7mhi8QH6iRza5uXzoCAIb06Bjh27hdpSzRyfTf4lBWrPnZR8tGcWebiRXZMD/YXAnT4uYv55ntk/i4iXQwq606wo7PxtZJMNDfyNOtToxWctXMq7ykqh507Rk7ZRy5PwfbeQcEb1BiUf7Z5JHm1iHzkGWQjc492LzL9yIQV7S6AQfTORWD+SWHaIXn6YbrVzXaBugSDskl8qmbRfqRoLypvH+0sKDIvpRSoB2jOLPKdhVx6FORWgp97MZ966wBy8BVKwd2VHYwXzyAABGvCr672LjNfNDIYgvIFTreyErdQr55nPxxPfTCRCRVwXyA/IeWjnDLLMwC49SLu87YwCnuOlM/SnxcyhWwhY+L4PSHnozRHE/lnk5mom4zdqu1ctyQZB2KmLenbeXvr2vfSyRPzMXHJiGJxLfUfKQ9tnkCYXmreXtsPoGXCTWIRWn6C31LCHs8kIGCPah9JU2O6Z5GvDiadO0uO3Uke9ZEwpBOE15GnZO/fTWduo0SFYyQLyvmSYI8gBIYF+nULIeGjGTkrv5Lo0wHu4GHT3QfpMG7v/FhjUzY3saCzvdvL+ZPzuQ/Tk7dQBjx9HA0H4J3vr2Vm7qJk76WFBWPlC3pp0XERyXSY/xsPR1xOJQQHYhK1Uo5Xr0gBvYKHQnD2U3ol2zyCVfK5L48dwDC1LxIsXkHcn4n89So/8nfqpwnP7DiEIEUJI70QbCpgBv1KPnaDnx+IVC8k16bgEItAD4Bh6aySxKB4fs8UHVzgE7tViQ5O3U2Ei7LcpBNzCegISR8sS8YvzyWcz8PcuMnE/Uv84x3jgptx+/WVxMWhPPftNGbOjlpkRhb87mpgQhkEjqAd6JgOPlKBJ26gfs8gJ0FkLrqXYwN6yi16agL84BLoyPAuOoTkx+JwYPFfLfljIpP9KjwzGliTgc2JwsWdEkGeUom9ZKLS3ntlUzW6pZpKV2JIE/N3RPDWswOvZliTgERJs4X7qjRHE0gRoyQB/crCRXbSfenU4cXcifDc8V4Yae38M8eYIYlM183UZ88BRemokflsMNiMK53YNdH8JQgeNTrWyh5rYg41MTgs7Ihi7NRp/eSgJqw56kUlh2IFbyFt304U69p+ZBNz2g3afFjMvnKa/yyInQWuBNxCTaHE8vjge1zrQ79XML5Xsg8dcqUpscjg2PgwfFYzJ+nx7LIxlPau5NjIyMicnJyIioisH22w2Pp9PEFfurUmzqN7ClhlRgY69oGPPtLGFejZNhY0PxSaG4ZPCMej/cwuTySSTyfr4TTUOdPseSi3AvppI9P0J41E4qX+PQjHoyZP09lp26zQiUdHXKQj17y5OBv3RzB5sYA42smc1bIwUGxqIDVBh6WosQY5ipBjvWs/5bqx/7w6EA03YoVYWw2izCxmdyOhCLTa2yYYarGyQEEuQo1QlNigAW56EDwrAYCtq3xAgQHtmkY8co0f9Tm2aSiTAiiH+SuNAC/dRPBzlzIEBot6Nj6NJYdikMAIhRDGoQM+eaWMLdOyefKbciOqtbJAQCxejYBFS8DE5Dz3YH09TufPE9+4gFOBIyUM4jkVKkJyHFHwUJMTDxChcjAkg9nwXH0f/HUt8WMSM3UJ9Mo7MjoYs9Dtn29gF++j5sdi/hhEEfP4+hMRRhhrLUP/vQ6UY1GRjG6yo2cYanMjsQlJ3NwV5dxCODmYnRmIEAd3j/mhlCp6hxhbup3NasbVD4GroRz4tZp47Tb8/mpgXC+e+7yNxFCnBIiUIod46yeFrBLzYyGDs9FzyRAs7dTvlgZOTgNuZXWjZQfqtC8zhbBJSELgLfJOAdwsSol0zyKxwPHMTtbMOstCXndewmZsoPoFOziGT+3xoDPBhEITA6+EYemEw/mMW+cBR+rETsEi3D2IR+k8+M30ntXYI/sk4WDUGuBkEIfAR40Kxc7eRDVY0fDOVq4VHQ99Ra2Gn76A2VjM5s8k74+GSBdwPvlXAd6gE6Mcs4smB+LQd1CvnPXeFX9B1X5YymZuoiWH4wVvIfjJoDgW9ApoYgK9ZmoBPCsPuP0JvqmI+HU8MVMPV0yvVWdgHjtL1VrRnJgkfIuhV8EQIfFCkBNsxg1zVH5+6g/rbGeg19DIMiz4oZIZspEYG4yfnQAqCXgdBCHwThtDyJPz8bbwiPUr/ldpbD72G3iFXy47ZQn1XzhzMJl8YjF9zbS0A3AuaRoEvCxOjnycT22rZvxylhwVh60bgsMy6xzI40Ytn6B8qmFcyifuSYU110Hfgdgv4vluisIJ5ZIoCDd5I/es8Y6O4LhD4M4ZFnxYzqb+47DS6OJ93P6Qg6FsQhMAviEj096HEyTnkeQ2b+gv1XTkDTaUeYl8DO3QT9VUps2Ua+eFYgtt96YB/gqZR4EdiZdhPk4mjTeyaHHpdPvPaMGJKBDx7cOa8hn3mFF1uRK8Ow2G9NMAhCELgd8aGYifmkL9WMg8eoyPE6J+ZxOgQiMM+VaRn155lDjUyLwwm/pICI2IAx+ALCPwRhtD8WLxgHrkkEb/rID1zJ3WsGdpK+0KRnl12kJ6wjRoUgJUt5D3YH1IQcA++g8B/kTi6NwkvXkDe3g9fcpCevJ3aA7Mses3ZNvbO/fTEbVSKEiu9g/dMBi6BBingGeCbCPwdH0crUvDlSfh35czjJ2g+gZ5Ix+fHwpOKe7AI7a5j1+XTRXr0WDr+yTie23dVBaCHIAgBQAghEkfLEvGlifi2Gnb9Bfqpk8yD/fH7kvEgIdcl81oWCn1TxrxzgeET6LEB+KJ4uLcAHgqCEID/wRDKjsayo8lcLbuhgEn+2ZUdha9MxcfAaJqbUaBjPypivi1jJoTh744hJoVB7QGPBkEIwDVkqLFPxhFvDCc+L2FWHKExhO5Nxu9KwENFXJfMgxld6KcK5vMSptqM7k3Czt1ORsE6PsAbQBAC0CmVAD2ejj+ejh9pYj8vYfr/4hodjC1OwOfEwECP/3ExaHc9+10Zs72WyQrHn8nAZ0XhBCQg8B5wNgNwY+NCsXGhxAaK2FjFfFvGPPgHPT0SXxCLzYzCxf56DlEMOtDI/lzJbKpikhTY4nj8ndE8WBcGeCN/PYkBuHkSEi1JwJck4BoH2ljFfFTE3HeEnhSGz47BZkXhIf7Rampwot31zO/V7I5aJlGBLYjFT88lo6XwAAi8GAQhADctQIDuT8bvT8Z1DrS9lvm9hl2T44qXYTOisGkR+KhgjPSt4ZEMi85r2D317K465nQbOzYUmx2NvzacDBdD/gFf4M4gdDgcAsH1WkZcLhePB3OIgO9QCdBdCfhdCYhiiD+a2d31zOMn6BIDOzoEmxCGjwvFMgMxAcF1KbuFZlGelj3SxB5qZA81MsEibFok9sRAYmIY5retwcBXuecbXVdXt2jRory8PD6fv379+iVLllxxwEsvvfT555+3trZKJJKHHnpo7dq1bnlfADwEiaMJYdiEMOKVTKRzoMNNzMFGdvVxpsjAZqix4UHY8CBsSCCWIMc8eYehGjN7VsOebmVPtLCnWtlICTY2FJsXi707mhcm5rpwAPQa9wTh6tWrBwwYcOjQoTNnzmRlZU2aNCkiIuLyA4KCgvbs2ZOcnFxUVDR+/Pi0tLQFCxa45a0B8DQqAZoTg8+JQQghswudbmNzWtifK9lnTzNaO5uuxgaqsXQ1lqLEUpUYh/MxtA5UbGAv6tgCHZunZXO1LImhIYHYsCBsTToxIhhTw8gX4B8wlu3p4op6vT4oKKi0tLRfv34Ioezs7PHjxz/11FOdHT9//vy0tLS///3v1/x/IyMjc3JyrsjRzthsNj6fTxDe2fbk/Uwmk0wm47oU3kTnQLlaNl/LFujYQj1bqGftNEqQY7EyLFaGoqVYtBSFi7FwMQoS3rhNtSv1TzGoxc4221C9BdVa2BozW2VGlSa21MDSLEpWYGkqLFWJDVRjA9UYPPbdFPj+c8uN9e+GJ8Lq6mo+n9+eggih/v37l5eXd3awVqs9evToAw880NkBLMsaDAaxuOOMFAqFIpF/jMYDfkAlQBPDsImXrbSid6IyI1tpYqtMqMzI7m9AjVamwYpabKyAQIFCTC1ASj5S8DERiSQkEpFI+P8B6XSSfD7d/rOLQWYXstPIRiGdkzU4kd6J2uysyYWChChEhEVKUKQEi5Rgt0ajWBmeIMdg9TgA2rkhCPV6vUQiufSrTCarrKy85pEul2vZsmXTp0+fMmVKZ6+m0WhGjRqF4x2j7qZOnfrZZ591djA8EXLLYrFgmAd3eXkDEqEUIUoRIhR05f9ldGEaB9I7kd6JGV3IRiMbjdko5GA66pxknQL0/z8TbD8h4uOsmEQKHlLwkYLHqgVIze+kyYdCZnNv/aP8BHz/udXF+hcKhSR5g6RzQxAGBgYaDIZLv+p0upCQkKsPo2l66dKlCKGPP/74+q/W9aZRgiAgCDnEsqxUKuW6FD5LilD4dQ8wmVwyGTzWcQa+/9xyY/27YbpTv379SJIsKCho//XcuXOpqalXHMMwzPLly7Va7S+//MLn83v+pgAAAIBbuCEIJRLJ4sWLn3/++cbGxu+///7cuXOLFi1CCOXl5WVnZ7cfs2rVqn379j388MNHjx7du3dvYWFhz98XAAAA6Dn3TJ9Yt27dY489Nnr06LCwsM2bN6vVaoQQy7IURbUfYDab+/fv/84777T/OmvWrKufGgEAAIC+54bpE+4F0ye8CAwf5xbUP7eg/rnlxvr3rSURAQAAgJsEiwb6LJ1df6G18EJrUbWxttHcrLHpXAxlp+xyvkwhlAeLA+NVsYmq2KGhgwJEKq4LC4B/MTstZ5vzijWl5fqqBnOzwW40Oo0kzhMQfIVAHiYNiZKH9w9MTg/qHyoJ5rqwvg+C0Ne0WTV7qg4dqT1RY6wbEJSSFpgyO3FGmCQkQKzm4zwhKTQ4jAaHqcnSXK6r+qPu5IYznwSLA8dHjZ4RlxUiuWouGwDAfYwO056qQweqj1ToqwcG908NSM5OmB4pC1cIZAqB3MVQTtqpsxuazM1VxtqjtTnvn/1cIZCPjxo5OWZ8jCKK6+L7LOgj9B1nmnI3lmzLbS4YHz1qYvSYwSEDSfzGNcOwTKGmZE/lof3VR5IDEhal3j4kdGAX3xH6SLgF9c+tm6r/Ml3l9xd/O9FwelT4sGlxEwcHp/OIG2/Fw7BskabkUM2xPVUHo2QRc5JmToweg2PQpYWQW7//EIS+4Fj9ya8v/GR12RakzJ7cb4KI7M4kaxft2lt9+IeLvwlJ4YqMpZlhg274V+BCzC2of251sf5LdRUfn/+6XF+1IGX2rQnTJbzurOhKMfTRuhO/Fm/V2nRL0uZPi5tEYP5+3YMg7ABBeLGt+P2zn9ko+93pd46NHIn3eMEnFrGHa45/kvt1sCTowSH3xSljrnMwXIi5BfXPrRvWf6u17cNzX51tzr07/c5Z8VN5uBu6onJbCr7M/6HVqlk5+O6xkSN6/oLeC4Kwgz8Hodame//s57ktF+7LWDItdlLPI/ByNEtvLdv9ed5302InLR+4uLNHTLgQcwvqn1vXqX+apX8p2vJtwS9zk2Yt6n979xppruNkw9kPzn2uFCoeG/ZAtDzSvS/uLSAIO/hnEBQBVwwAACAASURBVDIsu6lk+5f5P2QnTluadoeQ7K1d4/QOw3/PfnG2Of/pkQ8PDc24+gC4EHML6p9bndV/ub7q38feUgmVq4etjJCF9dK7MyyzsWT7V/k/ZidMvTv9Tj7hd0tXQhB28MMgrDHWv37iHRzDnxjxULS8S7XUQ2eacl87sWFk+NBVQ5YL/3xjCxdibkH9c+vq+mdY5tuCX34t3vLA4OUz4rL6oAxam27DmU/KdBVPjXwkPci/luuCIOzgriCkHYzT4HJZaNrB0A4aIYTzcIKPC1Q8gYqH4R6x0wqL2F+Ktnxz4ed7Bt45J3GWe9tCr8/ism44/XFBW9HfxjyRpI6/9OdwIeYW1D+3rqj/ZkvrP4/9h0/wnh35aKA4oC9LcqT2xFunP5wcM+7+jCWe82jo0LvsGifjYtsvqgQfx/k4T0rypSQpccPTCwRhh24EoctMWZsdtlanvc1ha3XaNU6HzsXSLF9B8qQkzsdJIYEQop0M42TsWqfLRInDhaokqSpVJo/lbAPvFmvbv4+/5aKp50avDpeGclKGfdVHNpz+aEnaHfNSsjGEIbgQcw3qn1uX1/+R2hPrTr63MPW2ham39eVN6iVGh2ndyferjXUvjH4sQRXX9wVoZ220ay+adMVmU42NFODCQD4hwAlBx0WVdjAuM+UyUQzNClU8YQBfGMgXBQlEQXxRkECg5KGbqTkIwg7XCUKGYp0Gl0Pnsuucdo3L3ua0tTnsLU5EIFGQQBwsEAXxhYECYSBfqOKR4k6jlKFYU7VVX2LW5BkxAgsfFxCcqcSIPv2iH6z5461T/12QMmdR/9u5nULUaG5+6ehrYZKQp0Y+LOGJ4ULMLah/brXXP8XQH57/8nDNsbXjnkoNSOK2SLsrD75/9tPFafMXpMzGbipVeohFmgvGhiMae5szMEOuTJLK4ySEoNOLFe1gHFqnTeO0tzltrU5bq8PW4qDsjCiQLwriCwM6rswCFY+v5BH8a78OBGGHpvMaU4kdwzCWYWk7QzsZykZTVtplohgXw5fz2ts2O+47AvmiQEH3H8lZpC811x9sc+hcCXdE9M3ToZ1ybDjz8fnmC38bsyYlILEP3vGGXLRrw5lPzjbnvTz+2UBcBRdiDkEQcstkMrlI6sUjr0n54udGPybne8Rn0WhufvmPdRKe+LnRq1VCZR+8o63FUfZzA+1kIrMCA9Ll3e5Loh2MrdVhb3Pa2jra6hx6l1PvQgjxZCRPQpBighAQhAAPHxcgiRBCEHbQlOsdTRSO4wjDSBGO83BSTPDEBE9GkqLeGkGjyTNWbGoMGCiPnR3aq92HlfrqtUdfTw5IfGzYA24fft1D7TeeKwfcPTN5Ctdl8V8QhNw6WX329bPvzk6cvnTAHX36+HUjNEt/kffD9oq9z41afc3x3m7UcFhTu7c1elpQ2JiAXqoD2sE4jRRlpSkbTTto2sGokqV8BQ+CsANXo0YpO138VS3CsJRlUdd5/O+JLWW7Pjn/zaohy6f3ydizbijTVT5/8JWsfuNWDFoKaz5xAoKQQ+1n6HOjV48IH8p1Wa7tXHP+K8fWT4+bdO/Axb2xDA3LsBUbG42V1v73xQhUN14uzu0gCDtwOH2CZdjy3xpN1db0v8Zep4uxGywu67qc96qNdS+NfapvJkh0W4Om8c3z75M4+eKYJ6R8CdfF8TsQhJygGPqd0x/ltlx4PvOxpNAErotzPXqH4ZVj6+2U48Uxa4LEgW58ZZZmC7+oYWk25e7oXnoYuCHYj5B7GI4lzA9XJkkvflrNUG67mSjRlq/Y8ZiUL/1g+hsenoIIIRlf+kbW2mh5xAO7nqgx1nNdHAB6nd5heHzfCxqb9oPpb4ZJQrguzg0oBYrXJ700MnzoX3auOV5/2m2vy6LSn+oRQv3vj+EqBd3LF/4NHIrNDhUG8Iu+rGEZN2Thb8XbnjqwdkXG0seH/9VzJgNdH4ERDw29/67+8x7Z8+zJhrNcFweAXlSur3pg5xMDg9NeHv+cmCfiujhdgiHsrrT5/xj3zPpTH3xw9nOKoXv+mtU7mq3NjuSlUR4yx7rnIAh7BkMJd0TQDqZ2T2tPXsbkNL9w+N87K/e9P/2NSTFj3VW6PjMzfso/xz/32ol3fircxHVZAOgVR2qPr9n3t78MWnZ/xhJOZgr2RHpQ6icz36ox1j+85+lGc3NPXqot19h63pB2f0xnsxq8ke/8S7iCk1jKsqim41pTta17r5DfWnjf9tVhkuD3pr3O1WT5nhsQlPLBjDf3VB369/G3XLSL6+IA4DYsYr/M/3HDmU9en7Q2K2Yc18XpJrlA9q+Jz0+OmfDXXU8cqD7avRdxGqny3xqSl0TxpD61qTsEoRvwpGT8vPDib2ppB3NTf5FhmS/zf3jxyKuPDXvgwaH3uWWXFg4FiwM3TH3VQTsf3fu8xqbjujgAuIGdsq898vrJxjP/nf7m5esLeiMMYfNTbn190tpP87597cQ7dsp+c3+fRaU/1oePUcuivaNZuOsgCN0jIF0ujxVXbbuJNocmS8sje57LbSn4ZOZboyIye69sfUlICl4a++SoiMwHdq4p0pT23hsxFGttduhLzO3/WZscN3sXArwRQ7G2FoehzNL+uVvq7bSzFz/3JkvLg7ufFvNEb01+RS1S9d4b9aUkdfzHM9ezLLtix2PFmrKu/8XG41rKSkdOCeq9snHFux9BPErcbWFnXi0NHaWShN148vuuiv3vn/18Udrtd6TM9br+huvDELZ0wB1xyn7PHPzHXwe7eR6krcXRlmvU5BmszQ6Bit+xOCGLnEaXXeciRYQsRiTvJ1YmS7vyKQCvYNc69UVmY6XVVGN16FwCFY+v6FgK32WibK0OgYofkC4LHKiQuvVJ5Vxz/j/+eHNJ2vx5ybe68WU9gYgUPjPq0QPVR58++Pf5ybMXp8274VRgykrX7GxJXxXrMwNkLgfzCN2p4YhGe9E0YGW/6xyjdxjW5bxfa2p4YfTjCarYvipar7j+PJ5qQ+3zh/81PGzIg0Pv7fl8XrvWWb2jxVBiDhysaH/+vvqEtGudpiqbscqiKzSzNBuQLg8crJDHiD1p0Q938u15hNYmR+s5vSbP6LLSqhSpIl4iixGLgwVXfposMjfYNXmG1rMGYQC/X3aINNINcfhz0e/fXfz1xTFPDA5J7+wYH6j/Vmvbv4+/baccz41eHSkLv86RFZsaGYpNmH+9Y/oYTKjv4GlByNLs2TfK4uaEqlKv/fEcqjn29ukPp8dm3TtwMY/gYC0G97rhF9HstPzz2Dqry/b3cU/3ZNnD+kNttXtbw8cFREwI7OK8JVuLo+28ofW8gXGywcOUwZlKYYB3zEjpOh+4EF/NZaZazuhbTukpGx04SBGUoZBGibpyK8PSbHOOrmZ3S2CGInZ2aLdXxrdTjjdy3q0x1v1z/LMhkuDrHOkb9c8idmPx9i/yv787feFtSdnXbKCytTrz3ikf8nSiR42RgSDs4GlBiBDSFpiqtjUNfiLhiucVnV3/1qkPKw01T498OC0whaviuVdXvogMy36Z/8O28j1rxz41IOim/+G0kyn7sd6mcabeEy1QdufWwVJvbz6laz1rkIQLQ0epAgbI+3jzkN7jGxfiDizSl5qbjmv1JRb1AFnIcJUiTtKNR3nazpR8V+ey0il3R/FlN33Vrjc1/u3wvxPVcY8PXyW40VxeX6r/OlPDayc2sCzz1MhHrl7Ko/DzGlmMKDLLs3oHIQg7eGAQIoTyNlSEjwsIHKRo/5VF7PbyvR+f/+qW+Kl3p9/pLTPlu6LrX8Tj9adfO/HOXWnz56fc2vXliWkHc+G/leIQYfz8cJzsUXoxFKu9YGw6rrU2O0KGq0JHqTlZHdG9fONCTFno5lO6xmNaUoCHjlIHDVESwp4N4mNR7d7W5hxd+kOxN3XzdKT2+LqT7y8fuHhO4syuHO8b9X8Jw7KbS7d/nvf9vORbF6fNuzSI3dJgL/i4OvP5pB6eg24HQdjBM4NQe9FUs7Nl0OPxCKFKQ81/Tn5AMa41wx/09h7Bq93UF7HJ0vLikVdDJcFPj3xEwrvxJlYMxV78pFqo5iUsiHBjJ5+txdF4TNtyRi/vJw4bE6BKlnpvD6K3X4hNNbbGPzTaCyb1AFnY6ABZjDuHutQfams+oUt/KI7XhZ3XXAz10bkvD9ce//u4p7u+35m31/81NVta3z79YZ2p8fFhDwwKSUcIlXxbJw4TRma5c6lSt4Ag7OCZQYhYdPaN0vBs9a+2TXuqDi4fuPjWhBk+NjS03c1+EV206/1zn5+oP/3i2CdusIUpi4q+qkEIS14a2Ruj1Bgn03rO0PiHlrLRoaPVIcNVXblcehovvRDTTqb1rKHpmJay02Gj1CHDVd3fJfS6qrc364rN6Q/GXn8NlEZz89+PvhEgUj096pGb2lPQS+u/K47Unnj3zCdpgckr4pZVf6Ad9nxyTx/TewEEYQfPDEKGZfZs+6PtrKFxevWKQUsVAjnXJeot3fsiHqk9/p9T/70jZc7C1Ns6uz9oOKxpO29IfzC2t/vzTDW2pmNaTb5RlSoNHaXuXr8UV7zuQmxpsDcd17aeMyjiJKGj1X3wOF7yXR3OwxMWdDrWcV/1kQ2nP16StmBeSvbN7inodfV/Uxy089uCX0y7nPGBMVOWjhJ62JaoCILwEg8MwqN1OZ+c/zpQEDDvjzvSV8RLIzzu2+NG3f4itljbXvnjPziOPzdq9dW7w1ga7Bf+W5WxOk6o7qP+VMpGt5zWNx3XsjQKGaEKHqbsxjiLvuctF2LKTredMzTn6JxGKmSEKmSEqnvjnrqBdjDn1pXF3hoakH7l/ajFZX3r1IfF2rK/jVmTqIrrxot7S/13G2WhT/6r+GjWgZPGM/ekL5wZP6U39jXsNgjCDh4VhGeacj/N/dZBO1YMWjoyPLPuQJutxZG40NO3UuqJnnwRGZb97uKvvxRtXjXkvmmxE//35xR7/j/lkZMDg4d2f7pFt5mqrU0ndJo8ozxWHDxMqU6Te9oAgct5+IWYZVhDqaX5tF530aRMkoSM6ItHwKuZqm2Fn1UPWpPAl//v5uZcc/6rx98eET501ZB7haSgm6/s2fXfc/UH26yNjsRFEUWa0o/Of9VibVs+cPGk6LEe0tEDQdjBQ4LwVOO5L/N/MDhM96TfOSlmXPu3xGWmzvy7dNiLyb6xX9c19fyLWKqreOXY+ihZ+OPD/9o+0bB6Z4utxZGyLMpNZewO2slo8owtp/XmOltAujxosEKRIPHABTU89ELMIlONrfWcvu28QaDkBQ1VBg9R9lIvYBfV7GqxNtpT7olGCNkpx8e5Xx2qOf7UiIeGhw/pyct6aP27z9nXShMWRsj7dQxtO9OU+3nedyanedmAhZNixt5wMZreBkHYgdsgpFn6UM2x7y/+5mKoJWnzs2LGXfHNKPqiRpkiCx3pI0sUXs0tX0QX7foi//tt5XtXDbl3gmrM+XVlg9Yk9FnT2fU5jVT7Bd2ucQWkywMz5Ip4iedMQ/SsCzGLTDXWtjyjJteI87DAQYqgIQpRUDcfttyLodizr5YmLoqokFa8fmLDgKDUhzPvv6lxMdfkWfXvbsZKa9lP9UOevnIM7enG819d+LHVqlmYOndG3ORuP0/3HARhB66C0OgwbSvfs7FkW5gkZGH/20ZFZF6zm11XaKrZ3ZLxqHevWH8dbvwilmjLXzvxztSi6ekJKWmzu9Nh06vsWqcm19iWZ7S1OVQpsoA0mTJZSoo4borwhAsx42QM5RbNBZP2ookU4YEDFQED5ZJwj+sarz3VVLCz8vP0Tx8b/sDIcPesce8J9d97Sr+vF4cLIiZce9ZEQVvR9xc35rdcnBU/ZW7SzOsvwdNLIAg79HEQsojNbSnYWrr7eMOpcVGj5iVnX7+PnWXY0/8sSVsRI/bRBaDdeyHQV1lyPy19J/3t2/rPWpg61zOXoHMaKe1Fk7bAaCi3SMKEqmSpMkkqjRZx0nDK2YWYRZYmu77ErC82Gyut0iiROk0WkCYXBnriYhEsYndW7P/43Ncri1emTY6PHOG2S7YPByFtZ069XDz02RusqdZobv6tZNvOin1pgcm3JswYGTG0L0fTQBB26LMgrDM17K06tKvigIAUZMdPmx43ScaXduUvVu9soe103Nyw3i4hJ9x7Ich/vzJ4mJLt73r3zCeV+pqHM+931517b2Ao1lhu0RWb9SVmh84ljxUr4iXyWLE0StRnbad9eSFmGdba5DBWWA0VFkO5hRDgyiSpKkmqTJJ64AyzS4o1ZW+d/hAhtHrYynBTRNGXNZnPJ7nrA/LhIGw6rtWXmFPuju7KwQ7aub/6yLay3XWmxqn9JkyJnZCsTujtEiIIwkt6OwhrjPVH604cqD6qsWknRo+dHjfpZj9gu9aZ+3bF8JeSPXCoRc+58YtoqrYVf1079LnE9oo62Xh2w+mPQyUhq4Ysj1XGuOUteo/LTBkqrMZyi6HCYm9zSsKFshixLFokjRIJ1fzeGyfZ2xdip8FlrrObaqymGpu52saTk+15r0iQeEgn7nW0Wts+zv3mdOO5FYOWzYjLau+8KPioKnCQImS4e7rtfTgI8zZURE0J6mzzgM7UmRp2Vx7YW3UYQ9jEmDHjI0clBcTf7OzMroMg7NAbQeiknfmthTkNZ47Xn7ZR9rGRIyZEj84ITuv2EKnct8pjbglRJnbpCdK7uPGLWPhFjTJBEjY24NKfUAy9uXT71xd+GhM5Ynn6okBxwHX+uuegHYypxmZuD49aG+1gJOFCSbhQHCYUhwjEoQI39iy690JMOxhbi8Pa7LA22s0NdkuDHbFIGimURotl0SJZjNhbFt+xuKzfX/xtc+mO2Ykz7uo/X8z738pt+lJLxW8NQ55KdMvF2VeD0GlwnXuzbPjalG4/Ohdryw7VHDtSe8JO2UeEDx0RPnRoaMblH4RbQBB2cFcQ2ilHkaYkt+VibsuFQk1JnDKm/cNLUrvhdqb+QJtN4/SofbzcxV1fRFurI29D5bAXkvCrlsIyOy3fXfx1S9muW+Kn3tn/NqVA0fO360suC21psFka7NYmh7XJbmt2YgQSBQmEgXxhAF8YwBeoeAIVjy/ndWPCYvfqn2VYl4mya10Oncuhc9o1Tlub09bqpKy0OEQgChFIQgXtye35j31XsFP2jSXbfyzcOCpi2PKBi4OvWqsBIXR+fXn0tGB1mhu+t74ahA2HNZZGu1vmQNca6483nM5pOHOxrbifInpISHp6cP8BgalSvqTnLw5B2KHbQWin7BX6mnJ9ZYm2vEhTWmOsi1fGpgelZoQMGBQ8wL13Lj7cOuquL2LZzw18ORk9vdNRDG027TcXft5XffjWhOkLUuaohF4Wh5dzmihbi8Oucdo1TofOZde6HDqn00iRIoIvI3nt/4kJUkKQIoIUEYQAJwQ4ISQwHLU/TWIE1j451Ww2S6VShBDjYhkXgxCi7QzLspSVZlwsbacpO0PZaJeZoqy0y0K7jC6niXJZaJ6EEKj5AiVPqOYLA3jCAL4oSCBQ8rxoebkr2Cn75tKdPxRuzAhOW56+KEbR6TzUtvOGhiOagQ+7YWSyrwZh3oaKqKnBqhR3NmI5aWdBW/G55vz81otFmtIgcWBqQGKSOj5eFZugjO1eLrqx/r1gHakeMjiMrVZNi7W10dxcZ2qsMzXUGOr0DkOMPCpe1S9RFT8zbkqiKrb3xigK1XyhimcotyoT3XAT5HtcZqot1zD02eutwR0oUq8etnJx2rxvC35ZtmXV1NiJC1Pnhkg8a3e0LuLLSL6MVMRf+WVwmSinmXKZKKeJoiy0y0rbWp20naYdDO1gKDuNGETZaIQQS7O0g0EIsSyLYRhCCOdhOA9HCBFCHMMwUkzgPIwQEKQQJ8WEQMWTRoh4UoIn5/FlJE9K+NI9mdFp2li8fWPJ1kEh6f/J+scNe5QDBsqrtjWbqm3u3ezCZzgNLluLw+0XKz7BHxySPjgkHSHEsEyloaZIU1qsKdtffaRSXyMihTGKqAhZWIQsLFwaGiwODBYHKoXKPlvCxrufCE/XnS/Sl+I4bnXZaIa203Y75bC6bCan2eAwGRxGvV0v4omCRAEhkqAQSXCELCxSFh6jiAyVhPTlKkF1B9rsvtg66pY7svpDbdZGR+KdXW2H0dp0PxVt3la+Z1jY4IWpc/tmfJpn8tUnki5qMDf9XLR5b+XhsVEjFvWfd/V2sp2pO9Bmb3Uk3NHTpj+frH83tot2XYu1rcZQV2dqqDc1NlqaWyxtLdY2k9OkEMgVArlcIJfzpWKeSEAIJDwxhmGz4qdEysLhibCDk3aanRYMw8Q8EcEjAokAISmQ8MRSnkQhlCsFcqVQeWl7SQ4FZshz36qIvz3Ml+7E3aU5R5ew4CbOOrVI9cDge5YOuGNr2e6/HX41UKSel5w9Pnq0J3zQoA8wLHu68dxvJVsL20qzE6Z+kf1ugOjmRoGGZCrPvFoaOyfMh5c/7La2XEPU1L6eHd/+CJgZNujyP6RZWm836h0Go8NkcBhtlN1BOayUrTcK4N3XjsHB6SMih3K+1ugNCdV8gYpnrLRe3SDm50zVVpZBlxYz7DoJT7wwde6ClNl/1J3cWLJtw5lPZsRlZSdMi5T52mM3uERj0+2o2LetbLeUL7kt6Za/j3tGQHRnCj9PRioSJG3nDSEjfHb5w+5xmihrs/vbRbuHwIgAkepm73K6x7uD0Iuo+8t0hSYIwis05+hCRqi6PUYDx/BxUSPHRY2sMzVsLdv98J5nI6ShM+ImT4we45ZhacATOGnnsfpTOyv2XWgtmhA9eu3Yp5IDetoeHjJCVbevFYLwCrqLJmWS1HNW0+0zEIR9RJUiK/2xvl821+XwJLSTacszDnnqylV9uyFSFv7A4HtWDFqa03BmZ8X+989+NiQ0Y0q/8SPDMzlcFBj0BMXQZ5tz91cdOVqXk6SOnxGXtXbs0+76NFUp0rKf6q3NDnEIfD3+R1dkUvf32Y3ErwOCsI/IokUuM+XQuQQqL5ub1Xvacg2KOMnlu8T1EIERoyOGj44YbnFZD9cc21a+5/UTG4aFDR4fNWpkRKaEd9MNsKDvOWnn6abzR2pP/FF3MkoePjF67IpBy9zePobhWPAwVctJXb9bQ937yt6LpVl9iSX+dn/sXIAg7CsYUqVIdYWm0NFqroviKVrPGnppjyoJTzwzfsrM+ClGh+mPupy9VYfXnXw/JSBxZETmyPDMrg8vBH2m1drWvpzTueb8JHX8uKiRnc2Id5fgIYqCj6v7ZYd67+xJ9zJWWoVBfJ7MH0PBH//NXFGlyFrP6SEI21E22lxtUy3v0qq+3SYXyNoT0U7ZzzTlHa8/9UvR7xjChoUNHhqaMSR0oELgjw1BHsJG2c83XzjblHuq6bzWphsWNnhSzNinRz3S850Cu0IcJiQEuKkGJhR20BWZ1De5uKjPgCDsO6pUafkvDQzFdmMxLd+jyTMqkyXEVWuq9RIhKRwTOXxM5HCEULWh9lTj+d2VB9/IeTdEEjQoZMDAoLT0oFRvWc7UqxkdpgtthXktF3NbCir11amBSUNDMp4Z+UiSOqEvp/a2C8hQtOUZIAjbaS+auz6d18dAEPYdUkSIwwTGcosy2QcX4L5ZbXnGkGFKTt46RhEVo4ian3IrwzIl2vLcloI9VQffOvWhgBQMCExODUxODUhKUMXCKBu3oBi6Ql9VpCm92FZ8UVPSZtWkBiYNDOr/wOB7UgMS+d2a/+AugQPlhZ/VxELrKEIOvctlpmRRfnpPAEHYp1SpMm2hCYKQstGmSmvKsk4XhOwbOIanBCSmBCQuTJ2LEKo11hdqSi62leyrOlxpqAmThiSqYhNVcfGq2DhlP69e4LQvWVzWCn1Vua6qVFdRqq2oNtZFSEOTAxLSglLmp8yOU8Z0eyMXt5OECzECM9fZpP4aAJfoCk2qFKnf3hBAEPYpZaK07Od6rkvBPc0FoyJR4mnrekTJI6LkEdNiJyGEKIauNFSXaitKdRV/1J+q0FXhGB6niomWR8bIo6LlEZHy8GBxUN+35nmaNpu2ztRQa6yvMdZX6WuqDDVml6WfIjpe2S9JHT8rfkq80qOfrQMz5G25RghCfanFbzsIEQRhH5NGCR06l8tM8aR+XfOaPGPQII9+wCJxIlEVl6j63x4FGpuuylBTZaitMtQcrTtRa6zXO4xh0pAIaWi4LDRUEhIqCQ6RBAVLAr1ur6iusLisLZbWRktLs6WlydzSYG5qMDfVmxqFpCBSFhGjiIySRwwJGRirjA6RBPXeXqxuF5ChKPqipl92CNcF4RSLDGWWWD+eSeLXl+O+h+GYPFZsKLcEZvjgtbKLGBdjKLck3RXJdUFuTvtqT0NDMy79iYN2Npga681NDeamRnPzuea8Zktbi7XVTjmCxYEBIlWQOFAtVKpFKrVQqRQqlAKFUqhQCGRCUsjhP+SanLTT6DDpHUadXa93GHR2g9am09h0Gpu2zaZtsbTiGB4sDgyRBIdIgsKkIamBSeHS0AhZmLfPzpRGCBGL/HxmvbXJTghxf57iDEHY1xSJEkOZXwehocwijRSRQk9fIfaGBAQ/Vhlz9b4/DtrZYmnV2HWt1jad3aCxaisNNTq7XmfTGxxGo9PEsKycL5XxpVK+VMITS3giMU8s5UtEpJBP8KV8CR/nCUgBn+ALCD6B4aL/3yBTxv9T77LFajFhlsv/xOqy0iyDELJTDoqhKIayUXYX47JTDpvL5qCdVpfN4rL+//9aTU6z0Wk2OU00wygEMoVArhQq1EKVSqhQi1Rxyhi1UBUoDggWB7p9e3HPoUyR6gpN/hyE+jKLMtGvBy5AEPY1ZYK06Hgt16XgkrbQrErx5d4IAcFv727s7AAH7TQ5TCan2eQ0W1w2q8tqpWwmp9lO2fV2Q72p0UE7nbTTtWISBAAAIABJREFUQTuctItmGZvLhhBiEWt2/in2GIbB8T/1s4p5YgLDEUICUsDDSQInxKSIh/OEpEBECvkkXyaQhkqDRaRIwhNL+GIpTyLjS+UCmcjzHlL7jDpV1nBYEzGxFyfvezhDmSXQs7sqehsEYV+ThAspC+U0uPgKP22I0BWZUnt5Hr2HExB8gTig59MWfXI/vL6nSJQUf1tLOxhPG73VR1hkrLDE+9xuqTfFLz94bmFIHi8xlFtufKQvsrU4WIqVhPrv8wfwNAQfl0WL9aVmrgvCDXO9jScl+X65stolEIQcUCZIDGV+GoS6IrMqVeY9gwqBX1ClSnWFfhqEhlKLwr87CBEEIScUiVJ9qZ8GobZ93i4AnkSdKtMVmbguBTcM5RZFgr9v3glByAFxsIB2MA69i+uC9DXayZiqrcokCELgWUTBAgzHrE0OrgvS51hkrLIq4rx7DkzPQRByAUOyGJGp2sp1OfqasdwijRT56ZAE4NlUyVJdsd89FFqb7DwJ4efreyAIQq7IYyXGSr8LQkOZRen3jTDAM7VP8OW6FH3NWGmVx8IpCUHIEXms2A+DUF8G3fLAQyniJcYKK8uwXBekTxkrrbJ+/t4uiiAIuSKNEtlaHLSD4bogfYey07YWB6xuDDwTT0ryFTxLnZ3rgvQpY5VVHgtBCEHIEZzEJGFCc62N64L0HWO5VRYjhk2JgcdSJkr0/tQ66jJRlJUWB/vv2nKXuC0IrVZrTU0Nw3T6iEPTdHV1td3uXzdc1yHzs9ZRfalZmQi9EcBzKRIkhjI/mk1oqLTKY8UwqRe5KwjfeeediIiIrKyslJSUwsLCqw84ffp0XFzc1KlTw8PDv/rqK7e8qbeT9xMbK/3o9tNQBtOVgEdTxEuMVVaW9pduQlOlBdpF27khCKurq5977rljx46VlZUtXrz40UcfvfqYlStXrlmzpqSkZPfu3atWrdJoND1/X28nj5OYqmx+0jlPWWmH1iWNhA5C4LlIMSEK4PtPhwWMlLnEDUH4ww8/TJo0KTU1FSG0atWq/fv3Nzc3X35AYWHhxYsX77//foRQZmZmRkbGxo0be/6+3o4nIXgywtrsF3N4DWUWWawYI6AVBng0RYLUT7oJGRdjbXLIYPAaQsgtu09UVlYmJSW1/xwcHCyTyaqrq0NCQi4/ICIiQizuuPVITEysqqrq7NVoms7Ly2tqamr/NSwsLDy802XRWYfV1VJD49465EcSROnOVfMoby0/bbU6xV26o9Sep6SBmLO2tLeL5Fe6Xv+giyRKpukcHZKs78rBXl3/pgZWqGap5nKuC3LTMAwnw2Iwwp2LALjhtUwmU2Dg/7bykkgkBoPhigNEItF1DricwWB48sknebyOLYqGDh26fv36zg42/7EDFRzrftE5Zx2grQwkLhzkuhzdxDCMrWt3IXrNHUGyw5rSpt4ukl/pev2DLmJYgblladv372Hoxn0WXl3/estAglZovj/CdUG6QzjtLjJ+oNncpZFNQqHwUqB0xg1BGBwcrNPpLv2q0+kufxxsP0Cv119+QP/+/Tt7NbVavWvXroiITjc1vRyZNY8/YxFBeOte57IaW+nP9WFr7uC6IN3Uxf3waCdT9WJR9DN/g7kT7gX7EfaG5tdKFXetk0TceKcwr65/4ze1oSmy4MxFXBekR9xV/264ncnIyDh58mT7zxcuXMBxPD4+/vID+vfv39LSUl9f3/7ryZMnMzIyev6+PkASIbS3Ommnj0+rN1VbJRFCSEHgFfxk1SdTtU0WDR2EHdwQhAsWLKipqXnzzTdzc3NXr169fPlyiUSCEHrmmWdeffVVhFBISMj8+fMfeuihvLy8l156iWXZmTNn9vx9fQBGYOJQgaXex+dWmiph9QrgNWT9xMYqHw9CykpTVloUBFPpO7ghCCUSyd69e48dO7Zy5crMzMzXX3+9/c8jIyPDwsLaf/7ggw+io6NXrFhRUlKya9cukvT3xc4vkUaLTTU+ftYZq2CUNvAa8n5ik68HoanGJo0UwVT6S9wTSOnp6b/99tsVf/jQQw9d+lkul7/99ttueS8fI4sW+fiOoCwy1diSFkMQAu8gChLQTsahdwmUNxhh4b3MNVYptItexluHPPkMaZTIVOPLE3gtTXaehIQNz4DX8IPtQk21NphBeDkIQo6JgwWUhaYsNNcF6S2mKqsc2kWBV/H58TLmWhs8EV4OgpBrGJJGCk2+u6qTsdIqg5EywKvIfXq8jEPnQizy4YbfboAg5J40SuzDyxsa4YkQeBtplMjW5GB8dF6TCR4HrwJByD1ptMhXB45SFpoy0+IQGKUNvAnOw8WhArOPzmsyQwfhVSAIuSeNFPrqVEJTrU0aBaO0gfeRRot8tZ3G3H5WgstAEHJPqOLTTsZlorguiPuZYJQ28E6yKLGvDue21NskEXBW/gkEoQfAkDRC6JPtMOYaWMYJeCVf7bBw6FwIx/hymM70JxCEHkESIbLU++Dtp7kOGmGAV/LVeU3mOhvsj301CEKPII0Qmut87YkQRmkDL4YhSYTQXOdrt6eWers08sYba/gbCEKPIIkUmX3uidBUY5PFwMQJ4K1k0T7YTWiut0mhg/AqEIQeQRwscJkoyu5T7TDmWhgpA7yYT3YTmuvsEngivAoEoWfAkCTM1yZRmGCkDPBmsmiR2beeCF1minEyQhWf64J4HAhCTyHxsdmELLLU2aERBngvgZKHMOTQu7guiNuY6+zSSCHM670aBKGnkEaIfKln3tri4MkIUkJwXRAAuk8a5VMPhZZ6mwSGjF4LBKGnkESKfGngqLkWRmkDryeNEvnSgvjmers0HDoIrwGC0FNIQgV2rZNx+cg6v+Y6uPcEXk8a6VMTfOGJsDMQhJ4CIzBRIN/a5OC6IO4B05WAD5BGCH1mxVHayTgNlCgIRspcAwShB5GECy0NPtE6yiJLPYyUAV6Pr+BhOOYb42WsDXZRqADDYajMNUAQehCfCUJbm4MUE6QYRsoAryeJFPnGcG5LA3QQdgqC0INIwoWWRp845eph0i7wEVJfWWjN0mAXh8FZeW0QhB7EZ54IzTCDEPgKia/Ma7I02CXwRNgJCEIPwpOSOOELHRLmOhuMlAG+wUf2zWaRpQmCsFMQhJ5FEuELD4Ww8yfwGUK1L+ybbdc6SRFBiqDb/togCD2LJMzrgxB2/gQ+xSf2zbY02CXQQdg5CELP4gPdhLDzJ/AxPrBvNnQQXh8EoWcRe38QwlR64GN8YN9sCMLrgyD0LOJggUPv8uqF1szQQQh8iw/sm21phCC8HghCz4IRmCjIuxdaszTYpRFwygHfIQriu4wUbffW21PawbiMlDAQFlfrFAShx/Hq8TKUjaastFANpxzwHRiOiUIF1iZvPSutjXZRCCyudj0QhB5HEib03lOuY3AanHHAt3j17Sm0i94QBKHHEf9fe3ce2ESZNgD8zZ1M0vtI0hMoLV1A6AHlKHIUwVYsQuWsZXFhBQEvWNzV1Q/wWBFQbhE5CutFFRRXpKi4nMpdBVoKbS2Fpk16pEna3Od8f4wbawtt2hwzSZ7fX5Ppm5knaSZPZt5n3lfM1cq89dIo9MkDn8SP4mq8NhHqGoyYiEN2FJQGiZByMBFH57UjjmqlBj50EAKf49X3NWllcBNhNyARUg4nmGWz4maNV45koa2HM0Lgg/hRXJ3MgHCy4+gVXYMRhtvuGiRCKsJEXG8sHMVtuL7JiIngkAO+hsljMDGGocVEdiA9ZlJbEI6zA2Ckp65AIqQivojjjfMx6ZtN7CAmgwMfKuCDvLSbUCeD2Ze6B99ZVISJvfKMECplgA/jR/N03pgIG4zQQdgtSIRUxBd7Zb2MVmrgR8GYMsA3eWm9jFZmgJLRbkEipCJMzNV6Yc88nBECHybwzkQIl0YdAYmQipg8BoPL8LoZerX1ekiEwFdxw9hmrcWit5IdSE/gSNdoxIRwRtgNSIQUxRd7Wb2MRWu1mXFOMIvsQABwDxrCxFzv6rMwKExMDObj7R4kQoryukNOS1yBgcHVgO/ii7lar6pig1vpHQSJkKL4Im885OAKDPBlfDHHuwpH4VZ6B0EipCjM2wpHdQ1wKz3wcZiI610dFroGKBl1CCRCisIiOfpmE27zmspRrRQqZYCP40dxdTKjF5Vz62RGPvw8dQAkQoqis+nsQKbXDOlEFKfBb0/g05gYg86meUs5N27D9S0mXiRMDto9SITUhYk43jK+jEFhYvKgOA34Pn6U11wdNchNnCAmnQVf8t2D94i6vGjobShOA37Ci8q5tTLotncUJELq4os43jJVvU4GxWnAL/C9Z95sqJRxHCRC6vK2M0I45IDv86IzQijkdhwkQuriCTn6Fu8oHIXxDIGfwIij0uoNRyWcEToMEiF10Zk0ThBL30z1wlGbBTcozVgkHHLA99GZNE4wS99E9Us1uBU3KMy8CDgqHQKJkNK8onBU32jkhrJpDBhdDfgFr+gm1DebOCEsOhOOSodAIqQ0zBvqZbQN0EEI/Agm9oKjUtdg4MN1UYdBIqQ0r6iXgZJR4Ff43jDQmhZGGe0JSISU5hVnhNAnD/wKJvaCDgs4KnsEEiGlYZEcg8JM8RI1uJse+BVuGNustliNNrID6QrcO9EjkAgpjcagcUMpXaJmNdosOis3FMYzBP6CRqfxIjm6RuoelTYLblSaeeFwVDoKEiHVYUJKH3I6mYEn5MB8vMCvYCJKz5KmbzJyw1hQyO04SIRUxxNxqZwItQ0wzwvwOxSvYtM1GjEhHJU9AImQ6jAhpXvmoU8e+CG+mKNrpO4ZIcyJ1lOQCKmO4oWj0CcP/BAmovQ99fDztKcgEVIdxQtHdTDcNvA/nGCWzWSz6ilaOKprgEujPQOJkOpoDBonhKWXU3HEUYvOZrPg7CAW2YEA4Fk0hAk5hiYqHpW4FTcqzbwIKBntAUiEXoCyV0cNjSYYvQL4J0zM1TeZyY7iHnRNRm4olIz2DCRCL0DZehljsxnGMwT+CRNxjM1UTIR6qJTpOUiEXgATUvQOCn2TGSplgH/CRFxDIxUTIdSv9QIkQi9A6Uuj8NsT+CW+mKOnZB+hrtGACeGo7BmXJUK9Xn/lypXa2tr7NZDL5SUlJXV1da7ao//gRXIMLVQsHDU0myERAv/EEjBpNJpJbSE7kI50DUZIhD3lmkRYUlKSkJDwwgsvDB8+/OWXX+7cYM6cOUlJSUuXLk1NTZ0+fbrJRMVfUpT126TYFCscNbWa6QwaS8AkOxAAyMGNZFFtoDViYnouTEzfQ65JhCtXrnz++ed//PHHn3/++f333y8vL+/QYMGCBQ0NDRcvXrx9+3ZZWdm+fftcsl//QcGp6rUyI1cIJdrAf3EiWFQ7KvVNRm4oTEzfYy5IhI2NjadPn164cCFCKDo6Ojs7+/PPP+/QZvLkyWw2GyEUEBAwaNCgxsZG5/frVyjYTahrNHDC4XQQ+C+ekEW1GXp1jXBdtDdc8EUmkUgCAgLCw8OJh3379u2ip/DXX389ffr0mjVr7tfAbDafPXvWvrWoqKiBAwc6H6S3w4TclrI2sqP4A53MyBXBGSHwX9xIdmspxY5KKBntFUcT4dq1a0tLSzusHDZs2IoVK3Q6HXG2R+ByuVqt9p4bUSgUjz/++MqVK1NSUu63I61Wu23bNg7ntx81Q4cOfeONN+7XWK/Xs9lsBoPh4KvwYoFWjVSv0WjIjuN3aqkuNJFLqZD8jVarpdHgIhhpbAKzVqrXqDXUmYasrV4TPIjvJ0elg59/LpfLZHaT6RxNhJmZmQkJCR1WRkdHI4SEQqFKpbLZbHQ6HSHU0tIiEok6b6G1tTU7O/uhhx565ZVXuthRcHDw559/Tmy5WwwGw08SIdYHr1I28Hl8qgwYgSNjszk4PlIgEJAdiv/CcRzefxLhOM7ktbEsHE4IVUYZNMllofFBmMAvTgpd+Pl3NBGOHTv2fn/q27dvcHDwhQsXRo8ejRD66aefli9f3qGNVqudOnVqSkrKO++80+tY/Zm9cJQiHQAGhYmJMRgcuA8V+DW+mKNrMFAkEULJaK+54IuMzWYvXbr0mWee+fbbb19++eXm5uYZM2YghM6dOxcXF0e0mTFjRm1tbXp6+u7du3ft2nXmzBnn9+tvKFU4Cl0RACCEMCGF5mOCktFec03V36pVq8LCwnbs2BEdHX369Gkul4sQEgqFs2fPJhqMGzdu6NChNTU1xEN7LQxwHCYk5gINJDsQhBDSNcDsSwAgTMRprb53SYTnQclor9FwnFrjlcTExFy8eNHBPkI/KpZBqPnn1paytuQ/x5IdCEIIVXxSFzJAwBvACAgIIDsW/6VWq+H9J5FarUZKZvUhacqKjvUTpKj9tgkhFJcdSXYgHuLCzz/08XgNal0alcEU2AAgTMTRNRlxGyVOJ3SNcFT2EiRCr8GL5BhaTFQYcRS34fpmEy8SDjng7xhsOlvANCooMQ0FjDLaa5AIvQadSZWp6g1yEzuIyWDDhwcAhIk5VBhfhigZhZ+nvQPfZd6EIjP0amVGPkxMDwBCCCFMxKXCUamHiemdAInQm/yvcJRkugYDBokQAIQQQnwRJc4IoWTUGZAIvQlFfnvqZAY+9MkDgBBCCBNzqTAgPtza6wxIhN6EInNQaOGQA+B/KDJvNpSMOgMSoTehwiFns+BGpZkXAfNOAIAQQnQmjRvK0jWRfKlGK4NLo70HidCb/FY42kxm4aiuwcCLYEOfPAB2mJhL7lT1v/08hZLR3oJE6GWIQX5JDAC6IgDogPTBLvRNRm4Y/DztPUiEXgYTcbWkHnIwyigAHfDFXHILR2GkJydBIvQymIhD7kUYrQzOCAH4A9LLubUNRijkdgYkQi+DiUiu1dbJDBicEQLQDi+cbdZYrEYbWQHArb1OgkToZXgRbGOrxWYm55CzGKwWvZUbAiWjALRDQ7wIMrsJoefeSZAIvQyNTuOFs3WN5BxyOpkRE3EQdMkD8EcYeVVsVpPNpLZww1ik7N03QCL0PiR2SOgaDDDKKACd8cnrs9A1GLFIDo0Ov097DxKh9+GLSauXgUoZAO4JE3O1UtJ+nkLJqJMgEXofTMQl69KoVmrgR0EiBKAjfhRXK9WTsmvoIHQeJELvg5E32j3crgTAPbEDmYhGM6ktnt813NrrPCbZAbjGypUra2tryY7Cc1RV2qBrmId7BWxmXF2rC6ri29c89thjTzzxhCdjAICy+GKOTmpgDxB4eL/QYeE8H0mEH3/88ZtvvhkUFER2IH6kuLj43LlzkAgBIGBirlZqCPZsIrTorVajlRMMJaNO8ZFEiBDKzc0VCoVkR+FHpFJpZWUl2VEAQBV8Mbe1WuvhnWqlBr6YC3c0OQn6CAEAwAVIGXH0t0QInAOJEAAAXAATc/TNJg9PF6qTGTAo5HYaJEIAAHABOovOCfb0dKFaGZwRugAkQgAAcA2+mKuVevDqKE7cRAj3TjgLEiEAALgGX+zRe3wNLSYWn8nkMTy2R18FiRAAAFwDE3M9OfyhVgpzorkGJEJquXjx4vjx4+0PX3vttS+++KKL9vv27du0aZPbwwIAOIAf5dHCUa0Mhjx0DUiE1GIymeRyObFcV1f373//Ozc3t4v2M2fO3Lhxo1Kp9Eh0AICucEPZFr3Vord6Znc6qJRxEUiEblFaWnrq1KlLly6tXr26tLQUIVRSUrJ+/fotW7bU19cTbVQq1aeffrp69eqtW7dKpdLOG9mzZ8/06dPZbDZC6KuvvpJIJAaDYc+ePQihioqK7777DiEkEAgmTZr04Ycfeu61AQDuh4b4Yq623kMnhRopTEzvGpAI3eL48eN//etfV69eHRwcbLFYdu3aVVBQwOFwlErlyJEja2pqEELHjh0rKSmJioqSSqVpaWkymazDRr755pusrCxiee3atWVlZRqNZsmSJQihc+fO7dixg/hTVlbW0aNHPfjiAAD3JYjhaeo9MQ2F1WQzt1l4EWwP7Mvn+c4Qax08/aO1Wu2hO1uD2bSiLAbjj6Mc0en0b775hsFgaLXa8ePHX716tW/fvgghBoOxefPmLVu2zJ07d+7cuUTj1tbWAwcOrFixov0Wbty4kZiY2O3eBwwYQJx0AgBIx4/y0EBrOpmRB/PxuojPJsLFf6K3eKrTmsdEjE6fxvT0dAaDgRCqrKzU6/WLFi0i1jc2NorFYoTQzZs3n332WYlEYrFYlErl/Pnz2z/dZDIZDAY+n99xu53w+fy2tjaXvBAAgJP40VzpmRYP7AhupXchn02EqWEk/1Dicrn2BQ6Hs3PnThqN1v5PS5YsmTVr1tKlSxFCK1assFr/0MHOZrODgoKUSmV0dDRCiMPhGI2/T8ZrMBg4nN/KppVKZUREhPtfEACge3wxV99isllwOtO9X0Haej0/GhKha0Afodv1798/MjLy4sWL/f6HSIQNDQ3JyckIIbVa/dVXX3V+4rBhw+zXPIcOHXrs2DFi2WazHTt2bOjQocTDa9euZWRkeOKVAAC6Q2PQeOFsD9xNqKk3QCJ0FUiEbsdisYqKilavXj1+/Php06YlJyd/8sknCKFly5bl5+fPmjVrzJgxf/rTnzo/ccaMGd9++y2x/H//93+lpaWZmZlWqzU5Odlqtb7wwgvEn44dOzZjxgyPvRwAQNcE0TyNuwtHcaSTQSJ0GZ+9NEqupUuXtr/UOXz48PLy8qqqKpVKlZSUFB4ejhB69tlnp06dKpFIhgwZwmazbTYbQmjkyJGnTp0inlVQULBhw4aWlpawsLDIyMhz585dvXp12LBhp06dioqKItpIJJKKiopp06Z5+hUCAO6DH83V1usRCnHfLnRNRlYAk8mFwdVcAxKhW9g7CO1YLNbAgQM7rIyPj4+Pj+/QjEiTCCGBQPD222+fP3/+0UcfJdbExMTQaDR7FkQInT9/fvPmzcS9hgAAKuBHc+VXW926C229XhDNc+su/AokQkqbOXNm+4c8Hu/FF19sv2bWrFmejQgA0A1BNE8rMyAcuW/ieE29QRAD10VdBvoIvQmfz3/rrbfIjgIA0BUGl84SMPXNxu6b9pa2Ts+HM0LXgUToFu+99x5xX4Q7bNu2bf369V00KC4ufvrpp920dwBAtwTRXLfWy2ikUCnjSpAI3WLs2LGzZ892x5Y1Gs26deueeuqpLtrk5OScPXu2vLzcHQEAALoliOVp6tw10JpBYaIzaOwA6NhyGUiEbmG1Wi0WC0JILpd/+umnVVVVr7766ltvvaVSqVpaWt5+++3Vq1ffvXuXaFxXV/fee++9+OKL7777blNTk30jN27cWLNmzdq1ayUSCTHWNkKoqKgoMzMzJCQEIfTDDz8Q2W7fvn1arVYqlR46dAghRKPR8vPzd+7c6eFXDQAgCGJ5mlp3JUJtvUEQA9dFXQkSoVucOHHi448/RghJJJJly5a98MIL8fHxFy9efPzxx+fNmxcYGKhUKidMmGA2mxFCR48eValUqampTU1N6enpxHhpZWVl48eP53K5oaGhs2fPJsbaRggVFxdPmDCBWN61a9fp06cRQs8995xKpaqoqHjjjTeIP02YMKG4uNjzLxwAgBASxPE09Xrc5pbhjjX1Bj4kQpfy2ZNrZdFmi6LRM/uisdjhC1cj+r1/VWg0mv3790dEREydOlUkEv33v//NysrCcfzLL78sKytLTU1dvHgx0TI/P7+ysrK4uHjOnDmbNm1atGjRSy+9hBASCoX28tGysjJHeh8TExOrq6v1ej2PBwcMAJ7G5DLYAUx9swkTun4GeW2dXjjCjTcp+iGfTYSCcdNsag9NV0vj8u+XBRFCIpGIGAtUKBTSaLTBgwcjhGg0mlAobGlpQQidOXPmueeeMxgMAoHg7t27Y8aMQQhVVFTk5OQQW0hLS7NvTavVYhjWbUjEaN1qtRoSIQCkEMRimlq9OxKhpk6f8HhU9+2Aw3w2EbLEfZC4D9lRIIQQMQdF54c0Gg3HcYTQk08++cEHH0yaNAkhNHPmTGJImuDg4NbW3+7JValU9qdHREQoFApimcvlGgy/V6a1P/+Ty+VMJjM0NNRNLwoA0LWAOJ5aooscHuzazRpVZoQjTgjLtZv1c9BHSD65XE5MMSGRSL7//nti5SOPPLJr1y6NRoPj+JYtW+yNMzIyrl27RiynpaUVFxcTHY0Ioa+++io1NZVYvn79enp6OpPpsz90AKA4N9XLqO/qA+K7vyYEegQSIflWrlz50EMPTZky5bHHHhsxYgSx8qmnnkpJSUlISOjfv39sbKz9VG/mzJn2KpglS5aEhIT0799fr9dnZmbevn37tddeI/509OhRGHQGABIJornaBqPN4uJ6GU2tThAH/R2uhlNMdHR0XV2dg411Op3FYsFxXCgUNjQ0uDOunjEajTqdDsdxi8WiUqns6xUKhc1mI5ZbW1tNJhOxXFtbe/XqVZPJpNVq9Xq9vb3JZLLZbF9//XVaWhqxxmazpaSklJaW2tu0tLRgGHbt2jX7mtbW1n79+ikUCre9PhzH8c2bNy9dutStuwBda2trIzsEv9bt+//zhiq1ROfanV7ffltZqXbtNr2UCz//cOnMLeyjYDMYjKCgIPt64v4/QmBgoH05NjY2NjYWIcRi/Xbp32QyLVy4cMyYMQqFYvv27Zs2bSLW02i0HTt23Lhxgyi6QQiFhobS6fSwsDD71kpLS99+++32+wIAeB5xddSF9/zhNlxT58oNAgIkQopisVjTp08vLy/HMKy4uNg+DS9CaNSoUaNGjWrfePny5QEBAfaHmZmZngsUAHAfAXE8da1e5LoN6hqMnGAWkwezL7kYJEKKotFoeXl5eXl5jjR+/fXX3R0PAKCnBHGY9GyLCzeortVDB6E7QLEMAAC4BV/MMarMFq21+6aOUd/VQcmoO0AiBAAAt6DRaQFxWNsdnas2qKksf8oLAAAVI0lEQVTVB8AZoRvApVF3MRqNe/bsqaioSEpKWrx4sb0KBgDgPwL7YW01utBBAd037Y7VaDO0mPhimH3J9eCM0C1wHJ8yZcqlS5fGjRt37Nix/Px8siMCAJAgsC+/7bbWJZtS1+r50Vwaw23T3vsxnz0jfP3HDXVqmWf2xWFytjz0Lzrt918VP/74Y21t7ffff0+n0ydPnhwTE1NZWZmUlOSZeAAAFBEYz9PKDDYLTmc6m8DaqrVBCXyXRAU68NlE+NeUeWqTxjP74jA47bMgQujGjRupqal0Oh0hFBAQkJSUdOvWLUiEAPgbOpuOCTmaWl1gP2dzWGu1NvahCJdEBTrw2UQYJXDh3Ts9ptPpNJrf07BGo2l/nx8AwH8E9uW31TibCG0WXFOnD+gDJaNuAX2E7nLy5MmamhqE0E8//UTMuEt2RAAAEgT2xVpvO1s4qr6rw0QcBge+sd3CZ88ISTd69OiZM2fyeLybN2/u2rWr/YBqAAD/EdgPqyqqx204jd77bsJW6CB0J0iE7hIbG/vtt99WV1fHx8dzuVDxDICfYgmYrECmTmbkR/f+e6CtWhc9Pqz7dqBX4ETbjZhM5oABAyALAuDnQgYIlBW9r93Drbi6VhfQFzoI3QUSoVs8+OCDM2fOJDsKAAAlhCQLVBXqXj9dXavnRXKYXBhr213g0qhbDB8+3Jmn22w2q9XqksFojEYjh8NxfjsAgF4LTODf+khiNdp6V+3SWq0NcvruC9AFOCN0lwsXLlRXV3fd5uTJk1KptPP6Q4cOPfLII87HYLVauVxuW1tbr7egVCqLi4udjwQAf8Zg0wNisdbqXg4xoyxXBw8QuDYk0J5rEqHNZjt//vzRo0dVKlUXzYxGY0lJSUuLK+cloaxNmzYdP3686zarVq26cuWKZ+LpnZqammXLlpEdBQBeLzhZoLzVm25Ci9aqazAG9YczQjdywaVRq9Wam5srkUj69eu3cOHC48ePP/DAA/dsuWrVqnfeeaewsHD+/PnO75fKzp49e/369dbW1jt37mRkZOTl5Wm12h07dty8eXPAgAHLli0TCATFxcV37tz58MMPz507N2nSpIkTJ95zU2q1evv27VVVVYMGDVq6dCmP99vY819//fXx48cNBkNWVtbcuXM1Gs2HH374yy+/sFisKVOmTJkypYvwXn755YKCgr1795pMpvnz59sv5F6+fLmoqEiv10+dOjU7OxshtGPHDqVS+dJLLyGEVq1ahWHQXQ9Ab4QMENz6dy1C4p4+UXFTHZTEd36ENtAFF5wRfvPNN7/++uulS5f+85//LFq0aPXq1fdsdunSpbNnz6ampjq/R+oTCoUhISFxcXHp6enx8fE2m23ixInXrl3Ly8u7cePG+PHjrVZrTEyMQCBISEhIT08Xi+99eFgsljFjxvz66695eXmXLl2aPHkyjuMIoVWrVv3zn//MzMzMzc2tqqpCCEml0qampmnTpo0bN27FihVFRUVdhLd+/fqFCxeOHj168ODBOTk5169fRwidOXMmJydnwIABY8eOXbx4cWFhIUIoOTmZzWanp6enp6czmdCjDEAv8cVcYvqInj5RcaPNJZNXgC644Kvt8OHD06dPJ85U8vPzH3jgAYvF0uFL02QyPf3004WFhYsXL3Z+j44o33tX12D0zL6YPMbQF/q1v1s2KSkpNjY2LS2NqB39/vvvGxoafvrpJwaDkZOTk5iY+N133z3yyCPh4eGZmZlTp06935aPHDlCTOdEo9EmT57cp0+f06dPp6SkrFu3rry8PCEhASFEPD0pKWnNmjV6vb6xsXHZsmVFRUVz5szpIuZXXnklNzcXIVRbW7t169Y9e/Zs2LDhxRdfXLRoEUKIy+WuXLlywYIFWVlZ7733HlTAAuAsGgoZEKC8pRFnhjr+JNyKqyq1CXlR7osLIJckwrq6umHDhhHLcXFxFouloaEhJiamfZs333wzNzc3JSWl260Zjcavv/46NPS3z0pMTMzIkSPv19hqtVqt9579OXFWtNVkc/Q1OIfGoHU9ZkRVVdXQoUMZDAZCiMFgpKamVlZWOlIOU1VVlZaWRqPREEJsNnvIkCEVFRV8Pj84OJjIgnYymWzWrFnNzc2xsbEKhYLYVxfs/4vU1NStW7cihCorK59//nli5bBhw2pqakymbn69EtWt3b4K4CZdfP6BB/T0/Q8eyJf9qIgcGeT4U1qrtNwINh2jwT+6MwfffzqdTnyFdsGhRHjt2rUVK1Z0Xl9YWBgfH28ymeyF/sSC0fiHU7Hr168fPnz48uXLjuzLYDAcOXLE3hOWlpbWxdVUo9GI4ziDwSAuGLbHCmBSZybcoKCg9mVESqUyODjYkScGBwd3eGJISEhoaGhbW1uH0+5169ZlZGS8++67CKEDBw5s3Lix6y3bS0lVKhURTPt9qVQqPp/PZrO7/gDZbLYO/2vgSSaTCd5/EvX0/cf6sXQHDeomLTvI0TMQeWlr4AAe/JfvycH3n81md9ut49D/o0+fPqtWreq8PiIiAiEkEonkcjmxprm5GSHUocfrX//6V3R09Jo1axBCEonk4MGDAQEBeXl599xXUFDQ7t27o6OjHQmMRqOx2WwGg9Ftwve88PBwiURCLI8fP37JkiU///xzWlra1atXL126tG/fvg5t7ikrK+tvf/tbeXn5wIEDL1y4cPPmzQcffFAoFCYnJ7/77rv/+Mc/EEJyuTw8PFyr1YaHhyOE9Hr9zp07uw1v586d27ZtMxqN+/btI+YNzs7O3rVr12OPPcZkMrdv304Uy4SFhbW0tOh0unuWyTCZTCifIZHVaoX3n0S9eP/DHgjS3jIFT3Bs5GEctVXU/WlhHIbB6FT34MLPv0PFMkFBQePuhQgiMzPzxIkTRMuTJ0+mpqZ2CG7BggV5eXn9+vXr168fh8OJiIggvrJ926JFi44cOSISiZYvXx4TE7N79+5HH300JSUlJyfngw8+iI+PRwg9//zz27ZtEwqF69evb/9cOp1O/ITp37//9u3bs7KyUlJS8vLy9u3bJxaL6XT6gQMHDh061L9//5SUlBkzZhCb+uijj4YNG5aamjp06FBiOzQajclk3vNXAoPBGDJkSEJCQlxc3FNPPYUQevHFF4krrklJSeXl5Zs2bUIIxcXFFRQUDBgwIDQ0VKlUuvk9A8DHRaQFN//S6mDj1motg0PniyELuh/uNIVCIRaLV6xYsWfPHpFIVFRURKwfPXr0hg0bOjTOyMjYv39/F1uLjo6uq6tzcNc6nc5iseA4LhQKGxoaeh67R1mtVqlUarVae/dEm83WYb1CoWj/qs1m8927d4nLxV2j0+nNzc0qlUqhUHT4k1qt7rzynjZv3rx06VJHWgI3aWtrIzsEv9ab99+GX1xzS9tgcKTtrY8k9WfkPd6F33Dh598Ft0+EhIRcuHCBzWZfuXKlsLBw9uzZxPolS5aMHTu2Q+Nly5b57cx8dDqdOJ/r3RM7n9iFhIQIhUL7QyaTGRcXx2azHdxsUFBQSEhIh5UCgaDzSgCAa9BQREqg3IGTQovWqrypjhzmUDEBcJJr7gyLi4tbu3Zth5UFBQWdW/75z392yR6BM9auXcvnw0AVAJAgIj341n5J7OSIrkvNm0pUoYMCmDwYaNsTYKxRf/T3v//dXpcLAPAkQQyPE8pq/rmbk8KGCwrRSLg24yGQCAEAwKPiHo6UHG/CbR1v+rJr/VWL21BgX7hs4yGQCAEAwKOCEvjsoPueFOI2/PZ/ZPHZkYhyN4X5LEiEAADgafE5wtrv731SKPtRweIzw1N6MAANcBIkQgAA8LTAvhgvnH33aGOH9Sa1RfJDc0JejyepAM7wkfkEOBzOiBEjuh1g02/hVtxmxRksuv1iC27DbWacwab3+vJLa2vrvHnzXBUhAP5mQEHs9W23WYHM6HG/DTBi0VorP5YIR4TwIjnkxuZvfCQRlpSUODMPuz9ovKyS/dgSNzmCHcjSyoyyn1oSZ0fzo3o/aAWNRgsKgqs3APQSE2MMWtzn+rbbFr0tKAFDOKr6rD58aFB8diTZofkdH0mE4eHh/jBsmzP69UOaTH1VUT3OpEVEhwz6W3JgH2eH6VOr1S6JDQD/xAlmDV7cR/aTQvJ9s7HVnPB4VOhAmHqQBD6SCIEjBDG81JX9yY4CAPA7XiSn33ToESQZFMsAAADwa5AIAQAA+DVIhAAAAPyadyfCr7/+ury8nOwo/Ne+ffuamprIjsJ/bdq0yWw2kx2Fn7JYLBs3biQ7Cv8ll8v37t3rqq15dyI8evTo5cuXyY7Cfx04cKCqqorsKPzXtm3bVCoV2VH4KbVavXnzZrKj8F/V1dWffPKJq7bm3YkQAAAAcBIkQgAAAH4NEiEAAAC/RsPx+86JRQqBQBAeHu7gqKHNzc08Hk8gELg7KnBPMpksNDSUw4FxEclRW1sbExNDp8PPWRLYbDaJRBIfH092IH7KZDLJ5fKoqKhuW+bn57/xxhtdt6FcIpRKpQaDwcHGZrOZwWDAFwFZjEYjZEESwftPLnj/yeXg+y8Wi3k8XtdtKJcIAQAAAE+CcykAAAB+DRIhAAAAvwaJEAAAgF+DRAgAAMCvec18hM3NzXv37m1ubp4yZUpWVlbnBlar9cMPPywtLU1OTv7LX/7CYrE8H6Svam5uPnLkyM2bN4ODg2fNmpWYmNihgc1m27Nnj/3h4MGDR48e7dkYfdmpU6cqKyuJZSaTuWDBgs5tpFJpYWGhSqWaNm3amDFjPBugj9u9e3f7osLOH+/KyspTp07ZH06bNi0yEmaZd0pTU1NJSYlEIhk/fnxSUpJ9fV1d3b59+9ra2h5//PGRI0d2fqLJZNq7d29VVVVKSkpBQYGD9xR4xxmhTqcbNWpURUVFfHx8fn7+Z5991rnNsmXLduzYkZiY+PHHHz/55JMej9GXPfPMM999951YLG5qakpJSTl37lyHBhaLZfHixbdu3bp9+/bt27dbWlpIidNX7d+//7PPPiPe25qams4NVCpVRkZGfX19TEzMY489dvToUc8H6cNqampu/8/y5ctLS0s7NDh//vyGDRvsbYxGIylx+pKJEye+/vrrL7/88vnz5+0r5XL58OHDm5ubxWJxTk7O8ePHOz9x7ty5n332WWJi4pYtW1asWOHo/nBvUFhYOHz4cJvNhuP4p59+OmTIkA4NZDIZh8ORSCQ4jisUCi6XW11dTUKgPkqv19uXFy9evGDBgg4NiCNfrVZ7Ni5/MX/+/I0bN3bRYPPmzVlZWcTyzp07MzMzPRKX37ly5QqPx1MqlR3W79+/Pzc3l5SQfJXVasVxfMSIEfv377evXLduXU5ODrG8ZcsW+2fe7ubNmxiGtba24jheU1PD4/Hkcrkju/OOM8IzZ85MmjSJRqMhhCZNmnT9+nWlUtm+wfnz5/v37x8TE4MQCgkJSUtLO3v2LDmx+iIul2tfNhgM9xvK54MPPti6devPP//sqbj8yIULFzZs2HDw4MF7zrt05syZyZMnE8uTJk06f/48TM/kDnv37p05c2ZwcHDnP9XV1W3YsKGwsLC5udnzgfmee17S7PA5P3v2rM1m69AgIyMjMDAQIdSnT5/Y2NiLFy86tDunA/YEmUwWERFBLIeFhTGZTJlMdr8GCCGhUCiVSj0aon84f/784cOHn3322c5/mjhxYktLy82bN7OysjZs2OD52HxYXFxcRESEQqFYu3ZtRkaGTqfr0KD95z8yMtJmszU0NHg8TB+n1+uLioru2UEbFBQ0ePDg1tbWw4cPJycnd752Clyiw+fcbDbL5fL2DRoaGtongsjISAcTgXcUyzCZTIvFQixbrVar1cpmszs0sFqt9odms7lDA+C8ioqKmTNn7t69u3///h3+xGazf/jhB2J57ty5Dz300NKlS/l8vsdj9E2vv/66fSEtLa2wsPCZZ55p36D9AUIswOff5b744ouQkJCxY8d2/tO0adOmTZtGLC9ZsuS11147dOiQZ6PzC91+znudCLzjjDA6Otqe2IkFsVjcoUF9fb39YX19vSODsQLHVVZWTpw48e233541a1bXLUeNGmWxWOrq6jwTmF9hsVgZGRmd62XaHyD19fUsFis8PNzj0fm4wsLChQsXEh00XRg9evTt27c9E5K/6fA5xzCsw2XqXicC70iEubm5R44cIQbjPnToUFZWFnG2cfXq1bt37yKExo0bJ5fLS0pKEEJVVVW3bt2yX0oGzrt79+7DDz+8atWqgoKC9usvX75MfOz0er195TfffINhWJ8+fTwcpA+zv71qtfrUqVODBg1CCJlMphMnThBlSrm5uYcPHyb6BQ8ePDhlyhQH528BDqqpqTl79uy8efPsa3Q63YkTJ4jzEvs/CMfxo0ePDh48mJwofV1ubu6XX35JnPMdPHgwNzeXWH/58mUiQT788MPXrl0jfileuHBBp9NlZmY6tGmXVPi4m8VimTx5cnp6+rx588LCws6dO0esz8rKeuONN4jlzZs3i8XiBQsWxMbG2lcCl8jOzubz+en/8/TTTxPrU1NTt2/fjuP47t27Bw0a9MQTT2RnZwcGBn700UekxutrQkNDc3Nz8/PzxWLxo48+ajabcRwnjvw7d+7gOG4wGMaMGTNq1Kj8/Pzw8PBffvmF7JB9zSuvvDJlypT2a27duoUQamlpwXE8Ozt74sSJBQUFDzzwQGJi4t27d0kK03csX748PT2dz+f36dMnPT2d+M7X6XQjRowYM2bMnDlzIiMjy8rKiMZDhgzZuXMnsfzqq6/Gx8cvWLBAJBK9//77Du7Oa2afsFqtp06dksvl48aNE4lExMqKioqAgAD7yW9ZWRlxQ31qaip5kfqgyspKtVptfxgQEEDc4nrjxo2IiAii17qkpKSmpiYoKGjYsGFwN7Fr3blz5+rVqyaTKSkpKSUlhVhpNpt/+eWXlJQUohfEbDafPHlSpVJNmDChfb0AcInKysrAwED7Nw9CyGAwXL9+PT09ncFgKBSKS5cuKZXK6OjoUaNGwWgezquurlapVPaHiYmJRC0ocSGkra1t4sSJYWFhxF/LysqEQqH9Y19SUlJZWTl06NCBAwc6uDuvSYQAAACAO3hHHyEAAADgJpAIAQAA+DVIhAAAAPwaJEIAAAB+DRIhAAAAvwaJEAAAgF+DRAgAAMCvQSIEAADg1yARAgAA8GuQCAEAAPg1SIQAAAD82v8D0190JFAeWaYAAAAASUVORK5CYII=", "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.5" }, "kernelspec": { "name": "julia-1.8", "display_name": "Julia 1.8.5", "language": "julia" } }, "nbformat": 4 }