{ "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 a custom\n", "potential" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "struct ElementCustomPotential <: DFTK.Element\n", " pot_real::Function # Real potential\n", " pot_fourier::Function # Fourier potential\n", "end" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "We need to extend two methods to access the real and Fourier forms of\n", "the potential during the computations performed by DFTK" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function DFTK.local_potential_fourier(el::ElementCustomPotential, q::Real)\n", " return el.pot_fourier(q)\n", "end\n", "function DFTK.local_potential_real(el::ElementCustomPotential, r::Real)\n", " return el.pot_real(r)\n", "end" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "We set up the lattice. For a 1D case we supply two zero lattice vectors" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "a = 10\n", "lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "In this example, we want to generate two Gaussian potentials generated by\n", "two nuclei localized at positions ``x_1`` and ``x_2``, that are expressed in\n", "``[0,1)`` in fractional coordinates. ``|x_1 - x_2|`` should be different from\n", "``0.5`` to break symmetry and get nonzero forces." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "x1 = 0.2\n", "x2 = 0.8;" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "We define the width of the Gaussian potential generated by one nucleus" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "L = 0.5;" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "We set the potential in its real and Fourier forms" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "pot_real(x) = exp(-(x/L)^2)\n", "pot_fourier(q::T) where {T <: Real} = exp(- (q*L)^2 / 4);" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "And finally we build the elements and set their positions in the `atoms`\n", "array. Note that in this example `pot_real` is not required as all applications\n", "of local potentials are done in the Fourier space." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "nucleus = ElementCustomPotential(pot_real, pot_fourier)\n", "atoms = [nucleus => [x1*[1,0,0], x2*[1,0,0]]];" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "Setup the Gross-Pitaevskii model" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "C = 1.0\n", "α = 2;\n", "n_electrons = 1 # Increase this for fun\n", "terms = [Kinetic(),\n", " AtomicLocal(),\n", " PowerNonlinearity(C, α),\n", "]\n", "model = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,\n", " spin_polarization=:spinless); # use \"spinless electrons\"" ], "metadata": {}, "execution_count": 9 }, { "cell_type": "markdown", "source": [ "We discretize using a moderate Ecut and run a SCF algorithm to compute forces\n", "afterwards. As there is no ionic charge associated to `nucleus` we have to specify\n", "a starting density and we choose to start from a zero density." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n Energy Eₙ-Eₙ₋₁ ρout-ρin Diag\n", "--- --------------- --------- -------- ----\n", " 1 +0.293966091773 NaN 4.29e-01 8.0 \n", " 2 +0.442105922092 1.48e-01 5.35e-01 1.0 \n", " 3 +0.260881007504 -1.81e-01 1.72e-01 2.0 \n", " 4 +0.243058931308 -1.78e-02 2.70e-02 2.0 \n", " 5 +0.242664736888 -3.94e-04 6.09e-03 1.0 \n", " 6 +0.242635951763 -2.88e-05 2.17e-03 1.0 \n", " 7 +0.242634026708 -1.93e-06 1.02e-03 2.0 \n", " 8 +0.242633384460 -6.42e-07 1.01e-04 1.0 \n", " 9 +0.242633377107 -7.35e-09 3.51e-05 1.0 \n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown:\n Kinetic 0.0304517 \n AtomicLocal 0.0972446 \n PowerNonlinearity 0.1149371 \n\n total 0.242633377107 \n" }, "metadata": {}, "execution_count": 10 } ], "cell_type": "code", "source": [ "Ecut = 500\n", "basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\n", "ρ = zeros(complex(eltype(basis)), basis.fft_size)\n", "scfres = self_consistent_field(basis, tol=1e-8, ρ=from_fourier(basis, ρ))\n", "scfres.energies" ], "metadata": {}, "execution_count": 10 }, { "cell_type": "markdown", "source": [ "Computing the forces can then be done as usual:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "2×1 Array{StaticArrays.SArray{Tuple{3},Float64,1,3},2}:\n [-0.038732124389381964, 0.0, 0.0]\n [0.03870220763913139, 0.0, 0.0]" }, "metadata": {}, "execution_count": 11 } ], "cell_type": "code", "source": [ "hcat(compute_forces(scfres)...)" ], "metadata": {}, "execution_count": 11 }, { "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": 12 }, { "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+gvaeTAAAgAElEQVR4nOzdZ4BU5dk38Ov06TPbG9tYlt6bgAoWUGygKLGE6AvoY0sUEzWWJzGoMWJFY00UNIo+KnZjI0YioIDSlN4WFnaBrdNnTr3fD4esCyywZeacM7PX7xPMTrlnZ+/5n7tThBBACCGEuiva7AIghBBCZsIgRAgh1K1hECKEEOrWMAgRQgh1axiECCGEujUMQoQQQt0aBiFCCKFuDYMQIYRQt4ZBiBBCqFvDIEQIIdStdTUId+7cGYlE2n9/TdO6+IoosQghuM2e1WA1sRpVVc0uAjpCYj+Rrgbh7Nmzf/jhh/bfPxqN4teupciyLMuy2aVAPyOERKNRs0uBjoCfiNUk9hPBrlGEEELdGgYhQgihbg2DECGEULeGQYgQQqhbwyBECCHUrWEQIoQQ6tYwCBFCCHVrGIQIIYS6NQzCFBOujhENdyRAKLVFauOaghXZKlizC4A6IFYv/fRclS2b7zW90F3qMLs4CKEOU+Pa3s8PHVjRVHRGdtkFeWYXBwFgEKYSAjvfrim9II9zslsWVheOz+5xVrbZZUKoG9mwYcN3333XlWeQI0rt0kZ7nuDr7fro+YbCH7OEDC5RxUtv2dnZl112WZKeHIMwZRz4romopPC0LKDAXWrf8NTuoglZFEOZXS6EuoslS5a8/vrrY8aM6fQzEJVAPlAMVbsbSCE07zpAs1iFT66xsXHnzp0YhN2d6JerP68bdHM5UAAAtizenis0bw1nDnCbXTSEupGJEyc+9thjZpei21m3bt2sWbOS9/w4WSY11K/xZw/1OvKElltyR/rqfvCbWCSEEEoPGISpoXlbOLOfq/Ut2UM9/u1hJYbHpCGEUJdgEKYAVdLC+2KeCmfrG1kb4+vtatgQNKtUCCGUHjAIU0BgR8Rd6mD4oz8s7B1FCKGuwyBMAc3bwr4+rmNvz+jritWJUgDPl0cIoc7DIEwB/m2hjLaCkGIoT5kjuCdmfJEQQihtYBBaXbxJUuOas8DW5k9dJfZwddTgIiGEUDrBILQ6/9awr68LjrPo1l1iD1VjixAhhDoPg9DqmreF2+wX1blK7OH9uA03QqhjwuFwy84Aa9eu/eabb9q8W2Nj4+uvv25gucyBQWh1oT1RT0/n8X7K2hjBx0UPikYWCSGU6oLB4B/+8AcAUFV19uzZubm5bd4tKyvrmWeeWbdunbGlMxoGoaXJYUVTieA70ba8LuwdRQgBEEIaGxtV9fAmG+FwuOXfLbfEYkd/V3z66ac5OTl9+/Y93tPecMMNjz76aMJLaykYhJYWqYk7C9ueJtPCXeLA+TIIdWe33XbbjTfeOHjw4OHDh//444/ff//90KFDR40aVVxc/Oc//xkAtm/f3rdv39GjR/fu3XvSpEmNjY0tj120aNHUqVMB4Be/+MXbb78NAJmZmcFg8LPPPps8eTIATJky5f33349G0/lLBjfdtrRIbXuC0H7wuyZjyoMQOsr6RlJrbEacXUgJzBG3RKPRzz77bOXKlYWFhaFQaMCAAQsXLjz77LMDgcC4ceNOPfXUoUOHrlixIisrS9O0m2+++ZFHHrn11lv1x65YsWLOnDkAEA6HRVEEgObmZkKILMuhUAgAMjMzCwsLf/jhh/Hjxxv6Pg2EQWhpkQNxb6/jDhDqnIW2eKOkihojYPseIaN9sZ/856Bm5CueksseFYQA8Itf/KKwsBAAVqxYwXEcIeRf//oXAAwZMuSrr74644wzVq1a9fzzz4fD4fr6+p07d+qPIoQcOHDgeAOELfLy8mpqahL/TiwDg9DSIrXxwvFZJ74PxVCOfFt4f8xbcZLIRAgl3O+H0L8fYv41aHb24WO66+rqFEXRUxAASkpKhg8f/vbbb//v//7vnXfe2aNHD6fT+emnn+o/pShKEARJkvT/EvLz/HNCCEUdXrYVj8ftdrtB78QMGITWRVQSq5daH710PO5Se7gagxAhBP369YvH4/fdd1/r6JoxY8att9567bXXAsAPP/zQ+v6VlZV79uzp27dvUVHR3r17W26vqqoqKioCAEJIdXV1ZWWlUe/ABBiE1hWtE22ZHM2d/GLT1cPevDVsQJEQQhY3atSo0047bdq0abfeeivHcatXrx4zZky/fv0WLVo0aNCgbdu2vfLKKzk5OS33nzRp0vLlyydPnjxz5sypU6dmZWUBwMKFC+fNm/fKK68AwJYtW2w2W//+/c16RwbAILSu9syU0TlyhdpvGk9+P4RQOjr33HMzMzNb/vv222+/8sorb7/9tizLAwYM6N+//7hx4xiGefbZZysrK//v//5v9erVbrf79ttvB4Brr732oosuuv/++8eNG/fVV1+99tprALBr165PP/102LBhAPDmm29ed911Ld2kaQmD0Lqi7Q5Ce54QqxeBwPF2YkMIpbFp06a1/i/DMLNnz549e3brG++6666Wf48cORIAHnjgAQDo06fPqaee+u67706fPn3w4MGPPvroY4899uCDD3q9XgAIh8Pvvffet99+a8TbMA8GoXWFa+OFp59kpoyOEWhGoMWAfOKl9wghdKwXX3xRUZQ2f2Sz2VatWuVyHXeXx/SAQWhd7W8RAoAjT4jViRiECKGO4nme5/mW/7aeO8qybNqnIODOMpbVns3VWrPnCrE63HEUIYQ6DIPQoiIH4sc7g7BN9lwhikGIEEIdh0FoUbE6yZ578hWELew5QqxOSl55EEKWtXr16i1btiTwCZctW9Z6P9LWNm7cuGvXrgS+lhVgEFpUvEGyZfMnv99/6WOEySsPQsiyXn755Y8//jhRz1ZVVXXjjTd6PJ7j3WHGjBmtxxHTAE6WsahYg+jp6Wj//QUfp0RV3HEUoW7oxRdfTOCzPfbYY7Nnz+a4ticoDBw4UBCEJUuWnHPOOQl8UXPhl6ZFdbRFCBTYcvhYPTYKEep25s6du2DBgn/84x9333335ZdfnpmZOXLkyD179vz+97/Pzc3t06fPd999BwB33XVXeXl5QUHBmDFjVq5cCQD79+8/55xz8vPzTzvttIcffvi+++5TFOWNN9649NJLY7HYyJEjRVH86quvbr75ZgC47777Fi1aBACXXXaZvulM2sAWoSURiDfLtqyOBCGAI1eI1YmuHum8Ny5CVqP6G7R4xMhX5HJ7AH3E8RM1NTU2m42iqGeeeeaLL7548803b7jhhlNPPfWBBx44dOjQc889N2fOnFWrVl144YVz584VBOHjjz++/PLLd+/ePWvWrKFDh37++edVVVXjx4+fOHHihg0bbDZbSUlJJBJZs2aNpmnNzc07duwAgD179uir7EePHv3www8b+ZaTDYPQikS/zNoZhu9Ye92ei/NlEDJa+JsP4ptXG/mKOb9+hHb52vzRpEmTxo0bBwBTp05dsmTJrFmzAOCSSy654447AGD06NFffvllVVVVNBptamratm3bv/71r8WLF9M0XVFRceWVV9bX1+/bty8/P//EBSgsLKytrZVl+XjdpykHg9CKYg2SvUP9ogAAYM8RGjcGk1EehNDxeKdc651yrdmlOEzfMhsABEFo2X3UZrPFYrF4PD527NgBAwaMHj06IyOD47jq6mqe51smxeTm5tbX1zMMo2ltH6/YcjCToig0TdN0+oyspc87SSfxxg4OEAKAvuMoThxFCLVl7dq1kiS9/vrrt9xyyxVXXBEKhXJycmia1rs9AeDHH38EgPLy8v379xNCHA6Hy+U6ePBgyzO0HOG7f//+0tJShjnmdOCUhUFoRfEGqaMDhADgyOHjDRKk1axmhFBi5Ofn79+/f+nSpbt27br++utZlmUY5tZbb7366qs/+uijhx9++D//+Q8ADBgwgOf57du3UxQ1ffr0W265ZevWraFQ6Jlnnlm3bt3ZZ58NACtXrjzzzDPNfkOJhF2jVhRvELOHejv6KJqnWQcj+mUhI0067hFC7TF69Oi8vDxCSEs/Z15e3sSJE/V/8zw/ffr0nj17Pvvss/fddx/Lsr/+9a8zMzN9Pt9DDz20YMGCjz76aNCgQTNnzoxGoxRFzZw586233vrjH//43HPPPffcc4sWLaqqqtq8efN3332nDx++9dZbTz/9tGnvNhlI14wfP37p0qXtv38oFNI0rYsvmvbWProjtC/aiQdu+Otu/85whx4iiqIoip14LZQkmqaFQiGzS4GOEAwGCSGPPvro7373O7PLkkiNjY2KohBCmpubBwwY8N577xFC6urq+vTpE4lE9Pu88847kyZNannIsmXLJk+ebHA5165dO3To0Na36J9IomDXqBV1bowQAGyZXLwJJ44ihNpl2bJlZWVlI0eO7NOnz6RJky6++GIAyMnJ+ec//6mqapsPKSsr0w/vTSfYNWo5UlChOZq1dWYg2pbJi01ywouEEEpLU6dOnTx5ciAQyMrKaj35paKiouXf55133qmnntry3x49ehhaRENgEFpOvEHsxNoJnZDJBXYaurYXIZTSBEHQ54Iej9PpdDqdhpXHFNg1ajmxjm6u1ootkxebsUWIEEIdgEFoOR3eZbQVWxYfb8QxQoQQ6gAMQsvp3LYyOt7LymGFqLiWECGE2guD0HLiTZItqwNH8rZG0RTv5bB3FCGE2g+D0HLEJtmW2fkV8biCAqFuaNWqVWeddVbfvn2vu+66YBD3HO4YnDVqLZqsqXGVc3X+cxEy+TiuoEDIKLIqK6TtJXdJYmMFCqjWt+zbt2/q1KkLFy4cNmzYnXfeOWvWrMWLFxtZpFSHQWgtYrMsZHBH/pF3jC2DE7FFiJBRHlv97NLqb418xTemvJhlz2h9y5tvvnneeeedd955APDkk0/m5uYGg8GW7dbQSWEQWovol4WMTs6U0QlZfPOWcKLKgxA6sbvHzrl77Bxzy7B3796ysjL931lZWR6Pp6amBoOw/XCM0FriTV3dMtuWyWOLEKFuRVGU/fv36/8OBoPBYLCgoMDcIqUWDEJrEf2y4OtiEOJkGYS6nXfffXf16tWiKN57772TJ0/2+do+wh61CYPQWsRmSejClFEA4D2cEtM0BZcSItSNXHPNNXPnzq2oqKipqVm4cKHZxUkxOEZoLWKzLPi6NEYIFAg+VmyS7LmdXIyIEEo5eXl5Tz75pNmlSFXYIrQWsblLiwh1uIICIYTaD1uEVkJACsi8t6tBiPNlEOpW5syZY7fbzS5FCsMgtBApKLMOhma7sIoQAHC+DELdTL9+/cwuQmqjAeD+++/v0aNHjx495s6dSwgBgFAoNGPGjJycnL59+7733ntmF7K7EJu7uohQJ2TgYUwIIdRe7DvvvLNw4cJly5bRNH3WWWf16dPniiuuuOeeewKBwK5du77//vuLL7541KhRxcXFZhc1/cWbu7qIUCf4WNGPQYgQOrmGhgaHw+FwODr38Jqamvz8/Nan26ci+qWXXvr1r39dXl5eWlp6yy23vPTSS5Ik/eMf/7jvvvs8Hs/ZZ589adKkV1991exydgtis2xLRBDyPg6DEKFuZfHixbW1tSe4w1tvvVVXV3fs7ddcc01Xuv369Olz4tc9nlWrVq1evbrTr5tY9NatW4cOHar/Z8iQIdu2bautrQ0Gg0OGDGl9o3kl7EZEv5SYFqGXk0MK0XApIULdxQMPPLB9+/YT3OEPf/jD7t27DSvPSb377rsffPCB2aU4jG1qanK73fp/PB5PQ0NDY2OjzWbjuMPfyF6vt76+/niP9/v9Z511FkUdMb/jN7/5zQMPPNDm/SORCCHkqPsjXaQuZi/lwuEE7BTK2Gn/oSDnPnl/hSRJAMDzCRibRAlBCIlGo2aXAh0hEolQFCWKotkFadvrr7++f//+xx9//M0335w+ffrEiRM///zz1157TZblKVOmzJgxY+HChXV1dfPmzcvNzb3qqqsmTJhw7JMsXrz4nXfeoWn6iiuumDp1KgBIkvTXv/7122+/BYDp06dPmzZt/vz5q1atUlV1/Pjxv/nNb1piooUoirfccsvMmTPnz59PUdSNN944fvx4ANizZ8+jjz5aXV09YsSIO+64Y9u2bV9++SVFUY2NjQMGDLjllltO+h41TWv93ah/Iu355TgcDpo+yUJBNjMzs+XwqmAwmJWVlZWVFY/HJUnSvxwDgUB2dvbxHu/z+ZYsWaK/1RY0TZ/ghZ1OJwZhm9TwQW+B2+mydf2pbBk8K/Eu18lnVGMQWo1+pehyucwuCPoZIcTlcgmCRTepGDFihM/nO+OMM4YMGdK7d++PP/74f/7nf1544QWXy3XzzTc3NjZOnDjR5XKdddZZ/fr1q6ioOPYZFi5c+MADDzz//POqql5//fXRaPTKK6+cNm0awzB33nmnqqpVVVX6d8Xtt98OAH/6058aGhoeeuiho55HUZS//e1vu3btuv/++/fu3XvJJZcsXbq0R48ep5xyym233Xb11Vc/9dRTU6ZMefPNN3v37s0wzPTp07OystrzHmmabl0p9E+k87+yI7GVlZWbNm2aOHEiAGzcuLGysrKgoMDpdG7evFnvMt20adOAAQNO8BQMw7AsLsNIALHLO2634H2cGJDdgEuLEEqunW/XNPxo6EG4w39fybuP+Mrt16+fy+UaMWLEGWecAQCzZs2699579VbdE088cfPNN996660Oh2PUqFFjxoxp8znnz5//l7/85dxzzwWABx54YP78+X369Pnuu+/27dunz6PRWzt33nlnbW3twYMHr7766oceeujYINQ9+OCDY8aMGTdu3Jo1a1544YV+/foNHjz4rrvuAoCXXnqpsLDwwIEDZWVlLMvq0WM6dubMmXPnzr388sspinrmmWfuueceQRCuuuqqP//5z6+99tqGDRs+//zzefPmmV3O9KfGNUIIa0/M5CvBx0k4Xwah5CufWlB2Qb6Rr8g6TvItUVVVNXDgQP3fgwYNqq6uVhSlQw/ZvXv3jh07+vTp03o2aXNz8wUXXKBpWnl5OQAcOnToeM/Wsq6xf//+77zzjiAILa0ph8PRs2dPS41WAgB71VVX/fTTT/qvYPbs2TNmzACAefPmzZ49Ozc3NzMz84UXXtDfNkqqeLOUkEWEOsHL4VJChAzACDRYoMe09XhTdnZ2y8SOuro6n8/HsuyJB6SOekh2dnZubu5RUbdw4cLy8vJFixYBwPLly7/88svjPVtDQ4PX6wWA+vr67Ozs7OzszZs3t/y0vr4+JycHAPRl61ZAUxT18MMPNzQ0NDQ0zJs3Tx/by8jIeO+994LB4J49e6666iqzC9ktiAlaRKgTMjgxgEGIUHeRl5e3adMmPVouvvjip59+OhqNKooyb968Sy65RL/Dxo0bj/fwiy+++PHHH5ckKR6PP/HEE5dccsmYMWMIIc8995x+h507d2qa5vf7NU2Lx+OPPPLICQrz+OOPE0IaGxsXLFgwderUKVOmfPzxx3oWvvrqqxRFjRgxIj8/f/v27bJsia8p3HTbKqRAV08ibA2XEiLUrfzhD39YtGhRr169FixYcNtttw0cOLBXr14lJSWEED205s6d+/e//72iokJv0rXIzs52Op1z5871er1lZWU9e/YsLi6+99577Xb7Rx99tGjRosLCwsLCwqeffnr27NnBYLCoqGjAgAHjx48vKysDgKKiomPniPTo0aO8vLxfv34XX3zxpZdeOnDgwGeeeeb8888vKip68skn3333XbvdPmPGDEVR+vfvP3PmTKN+ScdFdbFxOmHChPvvv7/NybhtCofDOGu0TXs/r6NoKDknNyHPJvrlH5/ePeqPfU56T5w1ajWEkEgkgrNGLSUUCrnd7scee+zgwYOPPfaY2cVpF0KIqqodmsmoqipFUUfN+ZckieO4li9tURRPMHtW/9ONx+Mcx2madtSry7J87IqL9li3bt2sWbPWrVvXcov+iXTiqdqELUKrkPyy0OVzJ1rwHlYO45p6hLoviqI6Op+fYZhjV77xPN+66dLONSQ0TR/76p1LQQNgEFqFmIgDmFpQNMU6WSl4kqliCCGUQIIgvP3225YNvOPBILQKyZ/IMULAFRQIIcOxLDt9+vST7uRiNSlW3DQmBRTem8h9CQScL4MQQu2AQWgJqqgRLWGr6XV4GBNCCLUHBqElSAGZT2i/KOAKCoQQah8MQksQ/bKQ0H5RABB8PI4RIoTQSWEQWoIUUBI4ZVQn+FjcXAYhhE4Kg9ASxIRuK6MTfLjdKEIInRwen2QJkl92FCTgGMLWODerRFWiEorBfXwQSoBevXo98cQTX331VZs/1WQNaIpOdHVT4xpj6+4tllgsZrcn8VA5DEJLkIKKr0+CPwuKpjgXKwWVBO7ljVB3dvHFF/fs2VNV1TZ/uvuDA1mDPN4KZ2JfdNPf9lRe2eOoAwi7odYHQiVcd//lWoToT+TREy14LycFkvLMCHVPgwcPPt6PqK89vcYUunokuOFC9/X0LC1wlyYxBlB3b3FbRGL3V2sheFkRd1lDyBBJq8Wc6MdanFwYhOYjKlFjKu9KfOtcbxEm/GkRQkchKlGiKudM5J4YOqzFBsAgNJ8YkDk3C0mY0cJ7WaxCCBlACiq8m6XoxFdjDEIDYBCaTwooCTyAqTXey4kB7FRBKOmk5PSLAi4INgQGofmSsb+aTsBrSYQMIQaUhG8OpeO9nISXs0mGQWg+MZD4/dV0vJfFKoSQAaRg8lqEuGlw0mEQmk/yJ35/NR2OLiBkDCkgJ/YYtRaHazFJxnOjwzAIzZeM/dV0DE8DDUq87fW/CKFEEZOwXbCOZilaoOUo1uIkwiA0nxSQeU+ydjYQcIABoeSTAjLvSdbOFYKXw5NkkgqD0HxSQEnSZBnA3lGEDCEFZMGXrMtZ3stJQazFSYRBaDYCUlBO3kaCvBfnXiOUdFJASV6LkPewuA4qqTAITaZEVVqgaS5ZHwR2jSKUbEpcpWiKEZJZi7FFmEwYhCYTg7KQtAtJAOA92DWKUHJJASVJU0Z1vAfXQSUXBqHJkl6FcCkhQkmWvG1ldDhGmGwYhCZL6mQz0LeuxxYhQskkJm0RoQ5bhMmGQWgyKagkb+0EAPA+7BpFKLmSt12wDluEyYZBaLLk7cyk41yMHFGJhvtSIJQsye4axVqcbBiEJkt2i5CiKc7FyiHsV0EoWaRA8muxk8FanDwYhCZL0qnWrfFeXISEUBKJATmpXaNweIwDa3GyYBCaLNnXkoCHMSGUZMk7Sa2F4OFEHCZMGgxCUxGQw0kPQjynHqHkIRqRIyrnYpL6KjhxNKkwCM0khRXWwVA0ldRXEfCceoSSRg4pnDPptRgnjiYVBqGZpOQPLYC+uQxWIYSSQ0zmpvkteA8rBfFyNlkwCM1kwAAh4OYyCCWTUZezOMCRRBiEZkr2IkIdtggRSp5kr4DS8V4OW4TJg0FoJmwRIpTqpGByd0nU4cGiSYVBaCZjWoSsndFUokpasl8IoW7ImMtZzsGooqYpuLlMUmAQmkk0pAoBAO9mZexXQSgJjLmcBQrHOJIIg9BMBlUhAB7PoEAoOYxpEQIuJUwmDEIzSQFFMKwKYYsQoSQQjbucZbFFmCQYhKYhKlHjKucyIggFrEIIJYGmEE3UOEdyt5XR8R7cbjRZMAhNIwUVzs1CcvejOIzz4NxrhBLv8JRRQ2ox78V+nWTBIDSNMcuPdAKuoEAoCaSgwhlVi3kPrqBIFgxC0yT7MM/WcL4ZQskgBWTBa1QQYoswaTAITWPMOlwdrqlHKBmkgGJcLcbL2aTBIDSNkV2j2KmCUDJIQYU3skWIl7PJgUFoGilgXBViBBpoUOO4uQxCiWRkvw5rY4hGVBFrceJhEJpGDMq826AqBNivglASGHk5CwAcLghODgxC0xhchQQvi8fzIpRYRrYIAUDA43mTA4PQNFJQFgysQjhMiFDCGXw5i7usJQkGoTk0hWiSxhqyIYUO92dCKLFUSdM0wtoMrMU4wJEcGITmkALGbUih43FzGYQSyshFhDrcNDhJMAjNYeSsax12qiCUWEYuItThmvokwSA0h8Fj7KCfcI2dKggljhSUDb+cxVqcFBiE5jB4jB2wRYhQopnTIsRanAQYhOYwclsZ3eEWITHyNRFKZybUYmwRJgcGoTmM7xqlWYrmaSWqGvmiCKUxKWR0LWYEGgC3iEo8DEJzSAGjryUBQPByIl5OIpQgot/oAQ7Awf7koD/44APqSGvXrn3ooYda37Jv3z6zy5lupKBxZzC1wAEGhBJICsiC8bXYg1tEJR598cUXk/96/fXXy8vLhw0bBgA33HBDy+3FxcVmlzPdmNIixM1lEEogKWTcqbwtsEWYDEd0jS5YsODaa6+lqMPLvBUFrzuSQpM0TSWs3bgNKXS4CAmhRFFiKs1QDG/06JKAa+qT4OdPsaqq6ptvvvnVr36l//eVV15xu935+fkPPPAAIceda6hpWnV19dYj1dfXJ73gqcz41fQ63ostQoQSQ/IbvYhQx3lYGVuEifbzB7lw4cJzzz1X7wW94oorbrrpJp/Pt2bNmosuuqiwsHD27NltPj4cDt9zzz12u731jb/61a/mzJnT5v0jkYimaS2Nzu4pfCDOuuhQKGTw66qcHGmMHfW6kiQBAM/zBhcGHQ8hJBqNnuDqExkvHA4ffcuhOO2gjK/FGq9EGuPGv67VHPuJHI/D4WCYk3S/HQ5CTdNeffXVJ598Uv9vz5499X+MGDHi+uuv//TTT48XhB6PZ/78+RMmTGhnmSiKcjqd3TwI44pmzxDcbrfBr0vlsQ3R8FGvi0FoNYQQmqZdLpfZBUFHOKriRGXFmWU3vhaTPLp5Tcz417WgBP4SDneNfvnll5FI5IILLjj2HoFAwOl0Jur1EBxeRGhG16gPu0YRSgzjN4fS4ZS3ZDj8QS5YsOCaa64RBEH/7yOPPDJ8+PDs7OwVK1a8+OKLH374oXklTENSUOGMXYer41yMHFGJRii6W7fIEeo6KSA78kCiUtcAACAASURBVGzGvy5OeUsGFgAURREE4brrrmu5VVXVhx56KBwOFxcXf/jhhxMnTjSvhGlICiiuQhOqEEVTnJORQ4rxSxgRSjNSUPH1NqFFSHM0zVJKTDV+2nkaYwGAZdnXXnut9a1333333XffbVKR0p8pq+l1vJeTAhiECHWV8bskttA30McgTCDcYs0EUsCcideAu6whlCCm7K+mwzX1CYdBaAIpaPTpLS3wMCaEEoCAHDZhcygd72VFnC+TUBiERlPiKlCHd5E3Hq6pR6jrpJDCORizJp3xHg4vZxMLg9BoUkAxfqPeFjzuz4RQl5k4zA8AgpfFrtHEwiA0mhQwbYwdcHQBoUQwaxGhTp/yZtarpyUMQqOZXYVY0Y9VCKEuMflyFtfUJxoGodFM7lTxYIsQoa4STdo3X4dr6hMOg9BoohknEbZgnYwmE03WzCoAQmnAlCN5W/AeVgopRMNt2RMGg9BoUtDMKgQ4XwahLpMC5mwXrKNoinMwclg1qwDpB4PQaCauptfxHhYHGBDqCilo8vZMOOstsTAIjWb6Dme8l8MWIUJdIQXMHOkHfZgQJ44mDgahsfQNKdxmtggF3JYCoS7QFKKKGucwc6tP3BkjsTAIDSWHFcbGUIyZpyDhthQIdYUUlDk3C6YeZSZ48HI2kTAIDSUGZMFn8skPeC2JUFdIfmvUYhzgSBwMQkNJpq6d0Ak+VvRjECLUSaLZw/yAl7OJhkFoKHNX0+twfyaEusISLUI8RiahMAgNZe7+arrDE69xMS5CnSKavQIK8GDRRMMgNJRo6haFOpqlGBsjR/ByEqHOMPcAGR3rYDRRwy2iEgWD0FBSUBHMvpYEfQUFDhMi1CliQObN7hoFCngPzpdJGAxCQ5m+DleHw4QIdZq5G422wK23EwiD0FDmblHYgvdxuAgJoc4gIIfMn/sNOHE0oTAIjXN4Qwqn+VVI8HISdo0i1HFSSGHtJu+JoeNxTX3iYBAaRwrKvIczd0MKHe/DKoRQZ0h+CwwQAgCuqU8oDELjWGE1vU7AMUKEOkUMWmJ0AwAELx4jkzAYhMYRrXQtibNGEeoEyS8LGdapxXg5mxgYhMaxwoYUOsGHw+wIdYZogUWEOsGHI/0Jg0FoHCtsSKFjBBooUOJ4wjVCHWORFVCAW0QlFAahcSQLHD3RAocJEeoE0W+VIDy8RVQYa3ECYBAaR/RbpVMF9MtJ7FdBqIOkgCz4LNGvAwACLghOEAxC40hW2JnpvwQfzpdBqMOkgGL6dsEteNwrMUEwCA1CNCJZY0MKHY9zrxHqICWmUgzFCFb52hR8nIQTRxPBKp9o2pOCCudiKdoCy+kBQD/GBccIEeoI6wwQ6ngvdo0mBgahQSw1tAAAPK6gQKiDpIBiqVqM66ASBYPQIJaaKQN4LYlQx1mxRYhjhImAQWgQS82UAX1/JqxCCHWERQ5gaiH4cKQ/MTAIDSJZ4Gz61jgnq+IJ1wh1hGSZPTF0PI70JwgGoUFEy+yvdhgFvA/3KkSoA+LNspDBm12KnzE8TbOUEsEtoroKg9Agot9ak2UAwJbBiX7J7FIglDKss+N2C1xTnxAYhAaRAoqlxghBP6e+GasQQu1luX4dnC+TIBiEhiD/PZXXSoQMHoMQoXZSIipFW2g1vQ5XUCSEtT7UdCWFFcbO0KxVVtPrcJc1hNovbr1+UdB3WcP5Ml2GQWgEyW+tWdc6DEKE2s+CA4SgHyODtbjLMAiNYMGhBQAQMnCMEKH2EpslKwYhXs4mAgahEaSAzFtsyihgFUKoI6x5OYt7JSYEBqERxIC19lfTMQIuQkKovcRmKwYhXs4mBAahESSLbVHYQvDhzhQItYtosdX0OkaggQIljpezXYJBaAQxYMVrSQAQMnCkHaF2sW4t9nISDvZ3DQahEcQmK843AwAhg5OxRYjQyRCNyCHFUhuNthAysHe0qzAIk4+AFLTqtSQOMCDUDkpIs9TB2q3hzhhdh0GYdGJAZh0MxViyCvl4CffdRuhkpIBizU4dABAyuDgGYddgECad6JeFTMuNseuEDA6DEKGTkq0dhNgi7CIMwqQTmyzaLwrYNYpQ+0h+xbK12JbBic14jEyXYBAmneiXbJkWrUK8l1UiKlGJ2QVByNLkoGrZIMQxwq7DIEw6a67D1VE0xToZOYSLkBA6EcmvWHARoY73snJYwcvZrsAgTDprrsNtgXOvETopK7cIKZri3KwUxMH+zsMgTDqx2aKLCHW8j8U19QidmOS37mQZABB8XLwJhwk7D4Mw6awehF4WJ44idAKqqBGVsA7G7IIcl5CJw4RdgkGYXEpMJYSwdutWId7HiRiECB2f2CRZ8PSY1nD6dxdhECaXxQcIAcCWyYnYqYLQ8cWbZD7D0kFow6WEXYNBmFyiX7bs2gmdkMmJTdgiROi44k0SZ/EWIS4l7BoMwuSy8mp6Ha8fQIFTrxE6jnijJGRaOwhxjLBr2Hg8Xltb2/L/3Nxcl8sFAIFA4Pvvv8/Kyho2bJh5xUt5ol+y8kwZAKBZirHTUtCiJyYiZDqxSXYXCmaX4kQEH3aNdgm7Zs2aM888s7i4WP//E088MXXq1DVr1px//vkjR47cvn37Kaec8vrrr5tbytQlNsvOQrvZpTgJWyYfb8IgRKht8SYpy+c0uxQnwgg0xVJyROWc1p2XZ2UsAPTq1Wvz5s2tb73nnnvmzJlz9913BwKB/v37f/PNN+PHjzephKktbu21Ezo+k4s3Sp5yh9kFQciK4k2SxSfLwH93HOWcVr/stiYaAFRV/emnn/bs2aNpGgAEg8ElS5bMmDEDALxe75QpU9577z2Ti5myLL6IUCdksLgaF6E2KRGVoijGZvXpFLjjaFewAFBXVzdr1qz9+/cXFBQsXrxYFEWapouKivR7lJSUrFmz5niPl2X566+/rqura31j7969Bw0a1Ob9NU3TNI2irHg4X8IRlSgRhXXR+hWGNWmaxvvYaLVk5UJ2K4QQvZqYXRAEABBtiNuyeOt/IryPjTd1o1rc/k+Epk9+EcMOHz68rq6O4zhZlmfPnn3zzTf/5S9/YVm25cGCIMRiseM9XpKkr776asOGDa1vvOCCCyorK9u8fzweZximmwSh2CRzblaURLMLciKSJFEuEmsU4/G42WVBAACEkHg8zrJW74vrJkKHIpyXFkWR4yzdtcO4qWhDrPvU4vZ/Ijab7aRZyNrth/uUOY67/vrrzz///Pz8fFEUQ6GQ2+0GgPr6+oKCguM93ul03n///RMmTGhn6TVNczgc3SQIpZqwLYt3OCw99sayLFPA1fr9Fi9n90EIIYTgx2ERTeGoI9dut9st/om48pT6/d2oFquqmsA3e0RO7tixIzc3Ny8vr2fPnkuXLtVvXLp06dixYxP1et1KvFG2ZVl6Wxkd72XlEB7jglAb4k2SLTMFarEti4834kh/J7GPPvpoLBarqKjYuXPn/PnzH330UYqibr/99ltvvTUYDP7www81NTVXXHGF2eVMSfFGKSWCkKIp3sOJ/tSIbYSMJDZJWQPcZpfi5OwYhF3ATpgwYfHixV988UVeXt4nn3xy6qmnAsCNN96YmZn52Wef5efnr1ixwum09Boay4o3SlkDPWaXol1sWVyqxDZCRtL7dRSwesYwNppiKTmscC4cXe4wdvTo0aNHjz72B5dffvnll19ufIHSSQpFi5DJx5tw7jVCRyIg+mXBxymi1YMQ/ts7ikHYCVZfHJPS4o2SLTs1ghAHGBA6lhSUWTtN86nxPWnLxFrcSanxAaciJa4SlaTKjkd4GBNCx4o3ySkxU0Zny+Ljjdiv0xkYhMkSb0iZ5iAc3m4UgxChI8QbJSFFRjdAD0KsxZ2CQZgsKTRACAACXksidIx4k2Tx80RbwwGOTsMgTJZ4Yyp1qvAuVpM0Nd5d9mdCqD3iDal0OWvP5uMNGISdgUGYLCk0UwYAgAJbNh9rsPRucAgZLFYv2XMsfRJha7yXlcOKpuDOGB2GQZgs8cbU2JCihT1HiNXh5SRCP4vVi/aclKnFFE3xPpz11hkYhMmSWmOEAGDPwRYhQj9TIipokFrL8nCYsHMwCJOCaEQKyCk0zA4A9hwhXo9VCKHDYvWiLTeVrmUBV1B0FgZhUojNMudmKSaVDtmw5fCxemwRInRYag0Q6nAFRedgECZFvFGyp9BMGQDAMUKEjhSrF1OuFmPXaOdgECZFyg0QAgDnZIAGOayYXRCELCHWkJotQgzCjsMgTIp4oyyk1JRRnT2Hj+EwIUIAkGpTRnV4GFPnYBAmRSq2CEHvHcVhQoQAgEC8QbKlWhAyNppmKTmE/Todg0GYFLE60Z6bYp0qgDtTIPRfUlCmeZq1pcam+a3Z8HK24zAIk4BALAUnywCAPRerEEIAqTllVOfIFaI4662DMAgTL94scQ6GEVLvd2vLxjFChABSc4BQZ8/FdVAdlnpf1tYXq5NSsV8U9DHCBglwq0LU7cUapFQNwhwhdgiDsGMwCBMvRQcIAYARaNZGiwHcmQJ1d7E60Z6dkrUYBzg6AYMw8WJ1oiPVdmZqgRNHEYLDY4QpWYvt2bzYLOMZFB2CQZh4sXoxRYfZAcCew+OOo6ibIxoRm1LqGLVWKIYSMnA1YcdgECZe9FCqdo0CgD1PiOIAA+re4g0S7+VoLlW/Hu05fKwOa3EHpOonbVmqqClxTfCl0rkTrTkLbNGDcbNLgZCZIgdFZ4HN7FJ0nj1PwCDskFQ6aislxOpEew4PqXTsxBEc+bZILQZhZ0gaNMRJQxwaRWgWSZMIfgn8IglIEFYgIkOzROIqxBQIyyBrEFWIqAEABCVQjx7QEQCOmLLk5oClAQAyeAoAvDywNHh5cLKUkwUXBxkCeDjKx4NPgEyByhQg2wY5NsqdqpdkZooeiDsKUrVTBwDsOUKoKmp2KVIJBmGCxepER8r2iwIA72GBouSQwrnxb+Noiga1UVIdhuoIORCF/RFyKAY1EVIXh4NRElEg2wbZNipLgAyByhTAx4OPp/K84ObAyYGPp20MOFhwssAzYGcofd8SDw+tD+wihEQiEZfL1fqlQzIoGgBAs0QAwC+CQiAoQVghERnCCvhF8EtkfwT8EjSJWpMI9XFoiBNFgxwbVeiEPDtVYIcCB1XsgkIHVeKCUhflxA+5LdGD8azBXrNL0XmOXP7Q6mazS5FKsB4kWLQ+VRcRtnDkC5GDcZ/bdfK7pq+gDDsCZEeA7ArBriCpCpGqEByMkRwbVeqCYhdV6IBSFzU6B4qcdI4Ncu1UVjI/9paGXYZAAQC4W35yks6HuAr1cVIbgUMxUhuFA1Gy/CDURrXqMOwNEzsDZW6q3E31dEOFh+rloSq9UOxM2Q6NBIkcEEvOSeFabM/FrtGOwSBMsNghMWuQx+xSdIkj3xY9IPoqu1EQNouwsZls9pONTWSLn2zxg18ivb1UpYeq8MC4PGpGL7rMDcVOKuXmT9gYKHZSxU5oMzLr41AVIntCZFcIVteTN3Zp2/wkJENfH9XPR/XPoAb4YGAmVe6muk82agoR/XJKX85yLhYA5LCi/wOdFP6aEixWL9pTdhGhzlkghPfFzC5Fcu2PkB8ayLoGsq4RfmwiTSIZmEENzKQG+KgppXQfL5S4usU3f44NcmzU6Jwj3mxQhm1+stlPNjeTF7eSTc3QJJJBmdTQLGpYFjUimxqYkXoXBO0XOyTas3iKSe0/AHuOEKuXMAjbCX9NCUUgnoKHeR7FkW9LvwGGsAzfN5DvDpFV9eT7ek0jMCKbGp5N/b/e1NAsulu1eE7Kw8GoHGpUq3T0S7ChkaxvJMsOkqc2anvCZFAGNTqXOiWHGpdHlabXRUMkxWfK6By5QvSQ6Cl3mF2Q1IBBmEjxZolNze22W3MUCNGDIpCTDj9ZXZMI3xzUlh4gKw6SLX4yNIsak0td3Yv661i2mzT4EsXHw4QCakLB4V9aWIZ1jWRVPVlcRX63SqUp6rQ8anw+dUYh1c+X8pcU0YOiIz+F107o7Hm4lLADMAgTKXpAdKTy8iMda2NYOxNvlmyZqdfHG1XgPwfIV7XaV7WkKkTG5VET8un5Y+mR2ZSQekfLWZSLg9PzqdPzKRgEAMyuIFlxiPznAHlioxaSyVmF9FkF1Dk9UrWlGD0Yzx+baXYpusqRbzuwvNHsUqQMDMJEitTGnYUpH4QA4CiwRQ+IKRSEW/3kn/vI5/u1VXVkRDZ1diH93Kn0qGyKTe3GeWqo8FAVHurqSgCA6jD5qpb8u5b8YY2aIVDn9qDOL6Yn5KfSVUjkQNyRn/Jdo85CW6QGFwS3FwZhIkVq46k+ZVTnLBAiB+KZA9wnv6t5VAIrDpGP9mof7iWiCheUULcMoM+cSLtwCbl5SlzUzN7UzN5AgFnfSD7fT+5fq/7URCYW0VNLqQtK6KQuMuk6Na4pETWFLgGPR/BxmkZw4mg74e8okSK18ZJzc80uRQI48m3NW8Nml6Jtigb/PkDerdLe36MVu6ippfQ7Z9NDs1KyFy6NUQDDsqhhWdTdQ+hGEf5ZrX1UTW75Th6VQ11aRl9aTudYsuskcjDuyBdSfXRc5yywRWrjvt7daB1Up2EQJowma2JATtGjW47iKLDVfN1gdimOQABWHCRv7NLe3aOVu6nLyulVU9lyd1p8Y6W7LAGurqSvroSowny+X3u3itz9vTw6h7qygp5WTnus1IKPHoinwTC/DoOw/TAIEyZ6QHTkpPzyI50jT4g3SppCaNb8t7MrSF7dob2+kzhZ+GUveuUUzL9U5WBhWhk9rQxiCvPJPu2NneS2lfLkYvrqXvQ5PSxRdSIH4im93XZrzkJbYFfE7FKkBgzChAnXxp2FdrNLkRg0S9lzhUhN3F1q2juKq7C4Sntpm7bFT66qoN+byGD/Z9qwszC9nJ5eDk0i89Zu7U9r1euWwzWV1LV9aHOvckJ7Y7nDfSYWIIGchbZanDjaPhiECROpjTsKrT0ToCNcJfZQddSUINwRIM9v0V7bqY3Mpm4ZQF9YQvM4+TNNZQpwYz/6xn70pmby8jbtlA+V4dnUDf3oi0po4xuImkJih0RnUZq0CB35QqxOIiqxRFvb2vALJmEiadQiBAB3sd3gjdY0Ap/uI5M/V07/RLEx8P1U9rPJ7LQyTMFuYUAG9cQYZt+V3K960Y/+qPV8S3l4g9Zk7IrwSE3cniuk7nm8R6E5Wsjgorisvh2wRZggBKIH4s7U35mphbvEvt+o+TIxBV7doc3fqLk4uHUg/eEkOoWWnaEEEhj4ZS/6l73odY3k6U1ar7fly3vSvx1IV3qNaNOE90VdJelzLQsAzkJbtDZ9Rj2Tx+ggrI9Ti/ZpFEW5OOBocLGUgwU3B5kCZKXyIaKiX6Y5Kp2W7DjybVJQVmIqa09iKDWL8Oxm7dnN6im59N9PZ07Pxz4cBAAwLItaOJ45FGOe26ye9olyej5952D6qM3BEy5UHfNWOJP6EgZzFtoitfGcEWaXo7MkDRrj0CSSkAxhGQIS0QAm96ATnhRGf3FHVFjbCABEP6Q7rGgRGUIyNInQKBJJhTw71cMJBQ6qxAVlLqqnByo9VE8PZfH+sUht3FmUVteSQIGryB7eH0vSeUx1MXhio/rSVm1qKf31BWxfH0YgOlqeHeaOYH4/hFm4Xbv832qlB+4dyrRseZpwoepY0ZnZSXpyUzgLU2OjtYMx2B4gu4Jkb5hUhWB/hByMQk2URBXIEiBToNw8uDnwcBRDwdhc4uYS/DdgdBCWOckLp9LUcTbmjatQFyM1UaiNkL1h2BUiX9SQnUHYHyGlLmpABjUogxqaBSOyqR4WOzs0ko79D64Se7g68UFYH4dHNqgLt2tXVtBrL8H9r9FJOFi4uT/9P33p13dq1y1XixwwdwQzPtGdB0pclQKyIy99RjfAqhutxRRY30TWN5L1jWRjM9ncTDgaKr1UpYcqc8OZBVSxiy5wQIGdyjDq07BWV56NgRIXVeKCow4+kDXYHiCbmsn6RvLiVrK2gTAUNSqHOjWPOjWPGpljfnsxUhPLGuw1uRCJ5i62168LJPAJ/RI8+qP64hbtygr6x0vZQgdGIGovjoaZvemrK+lFO7XZ36ilLnhoFJPAztLwvrizyE7RafU3qW+0JoUU3m3yV/3eMFl2kKw4RFbVkW0B0t9HDcumhmZSv+xFD8igTN94z1pBeDwcDQMyqAEZ1C96Hr5lb5isqiMrDpHffKftCJCxudRZhfS5PaghWeacAhPcGyubkm/GKyeRq8S++8ODCXmquApPb9Ie/0mdUkKvm8YWW6xBj1IFQ8HVlfRVFfQrO7TL/qWOzKEeGkknpF89bNJioWRzlzhCe6Km7IHcEIclNdq/asjXB0hMIafl06flUTN700OzzG+6HCU1gvBYpS6q1HU4FwMS/OeA9lUtufzfWlAi5xfTF5VS5xTRDqPenNgsE5WkwUa9R7Fl8kQlUlDhPZ3/VWoEFu3S/vcHbVQO9c2FbB9Dpv+h9MbScG0fekYv+tnN2oR/KtPK6D8NZ/K6lmKh6ljOsHTr1AEAT6k9tNfQIPyxiXy0l3yyT9vqJ2cW0BOLqNsH0/2sPQkgVYOwNS8PU0rpKaUAAFUh8s995NnN2jVL1bOL6MvKqYtKEj/F6CihvVFPWXqeBO0usYeqo1kDO1mLVhwic75TORrePJMZl2fpmpBUSkyVw6oSVZWoosQ1Na4qMVWVCJE1Ja4BgBJVW+5MNCKLCm9vbrmF5imapWmOolmaddA0R7N2hhFo1sFwDoZ1spyLSbM+vfawMfC7QfSs3vSf16sD35VvH8TMGdj5hTfhfbHyqenWqQMA7jLHviX1BrzQ6nqyuEpbXEUogKml1F9GMaflUamyJjMdgrC1cjf16/7Ur/vTzSJ8VK29sVO7aYU6uQf9y17U5B7JWikb3BN1l6VhpwoAuEsdwarOBGFNhNyxWltxiDw8ir6ioht8SROQgnK8SY43SZJfFv2y6JflkCIFFSmk0CzFuVnOwbAOhrExrJ1h7QzNU4yLtWXTAMA6Wn1/UyApIs/8PGyiSURTNE0mmqwpUU2VFD1KlaiqRFU5qioRlXUwnIvhPRzvYYUMTvByQiZvy+SETN4KG8YmT4YAj53C3NCPvmOVNuBd5YlT6CmlHa7n6dqpAwDuEkd4fyx5+8tsD5BFO7VFuwhLwfSe1PuTmCGZqff3lm5B2CJDgGsq6Wsq6UYR3tmtzdugXbdMndGLnt0n8Y300N5Y2YV5iX1Oi/BWOne/W9uhh8gaPLlRe/RH9YZ+9Euns4Z1UBtJFbXoITF6IB6rF2N1UqxejDdKrIOxZfJCJi9kcI48IaOvm3OxvJflXGyHoogQEolEXK6OTNYlIEeUw7kbVMRmKVwTb/gpKDbJYrPEuVhbNm/PEey5vCPP5sgThIyUXbF7HL081PuTmH/VkFu+U1/cqj01lunl6cDv3L897E3OMiHTMTbalslHauOu4kRerEcUeGu3tmCbtitIrqyg3zqLGZGdevnXIh2/pY6UJcAN/egb+tHbA+SV7drZnyq9vdT1felLyxOzd5emkOiBuDuhf2TW4S6xxxtlOaJyznZ1OS0/SG5YoZa5YOUUtqIj30QWJzbL4f2xcE08UhOLHhDlsGLPExx5giNPyBnps+fw9izezL3gKOBcLOdiHQXH/IiA6JdjDWKsXoodEpu3hKOHRFVUnQU2Z6HNWWR39bA5C2zpsR3lxCJqwzR2/kZt7EfKLQOYOwe3t6fUvyPiq0yrpfStuUvtob3RRAXhxmby/Bbt/3Zp+i4H5xfTbIr0f55A+gdhi95e6qFRzNwRzEd7tee3aL9bpd7Qj7mhL53btT+PyP6YPVdI1w0xKZry9HQEdoazh5xkHkGzCHesVr/YT+aPoS8tT/nfhippob2x0J5oqDoaqo5RAM4edlexPXdUhqvIZsvkU+bsVgqEDE7I4HyVP9+mxNRITTxSGw9WRWq/aYw3Sc4CwV3icJfaPeXOlG4vcjTcMZi+ooK65VttyHvK305vx4pDAoEd4dLz0+FI7Ta5yxz+7ZGC07r0JBqBj6u1pzdpW/3wP33pH6exRWk097sbBaGOo+HScvrScnpTM3l6k9Z3sTy9nP7tILrTsxmDe9N2gFDnrXT6d0ROHITvVGlzvtOmlVGbL2NTd588Na4FdkcCuyKBXZHYQdFZZHOXOfJGZVRcWij4UvZdtYW1M95eTm+vw20gTdLCNfHQ3mjDj8GqDw9SDOWtcHoqnN4KZ4oeNF3spN6fxHy4V/vl1+oFxdS80Yz3+O8jekikeTotBwh1nq7Nl4mr8OoO7YmftAwBbh1AX1aeNtuS/6zbBWGLARnUi6cxD45knt2sTvhEGZ9P3zOU7sSJd6E9sayB7mSU0CJ8la6D31Yf76cHY3DTCnV7gCyeyIzNTb0rRKKRUHWseWvYvz0cPRB3lzq8Fc6eUwtcxfb0nmPSGs3TnnKHp9xRBAAA8QZJvxrY92UdUODr7cro4/L1dh0xoycVTC2lzyig7/peHfSu8vxpzAXFbX+g/h3hNO4XBQB7jqDEVDmkcB1cVh+S4YUt2pMb1dE59MunM6el71bA3TcIdTk2+NNw5s7BzN+2ahd9qQ7Lov40nB7ekVHf0N5ous6U0TkLbEpME5vlY3vM3til/XalOrsP/eaZbGqdF6FE1eYtoabNoeZtYVsG5+vrLj0vz1Pu6D7hdwK2bN6WzeedkgEAsTrRvz1St8a/8+1aR76QOcCd0d+dQrsJenl4/lTm657k2m/Ut/Opp8YyvmMafv7t4dwRaXIYb9socJc6gns7MP07JMNTG7W/blYnFdFLzmMHZKR5vejuQahzsDBnIH1jP/qlbdrUJeqIbOrBkfTAdnz2YrOsqcSWlbadKgAAFHh7Of07w3mjMlpuq4/DjSvUbX7yz3PZFJotQq52JQAAIABJREFUJgWVxh+DjT8FQ/uivkpXRj93+dSCrmwXkPbsuYI9Vyg4LZOoJLAr0rQptOXlaqAga7Ane7DHXeJIibHSMwuoHy9l71qtDn5X+fvpzLk9fi400Uhwd7Ty8iITi2cAT6k9tKddQRhV4OlN2pMb1XOL6OUXssYcgGU6/Ar4mcDAzf3p2X3oF7ZoEz9Vzimi546gy90n+jvwbw8n6XAGS/H1dgW2R1qC8ONq7Ybl2i97UYvOSI2GoBxWGjYEG9YHIgfimf3dBadl9u9bkoYDHclEMZSvt8vX29XzkoJITbxxY3DHW7WqqGYP9uYM81r/GD8nC38dx1xcRmZ/o55fTD12CqMv7AnviwsZXDodoNYmb2/X7ndr4cIT3UfW4OVt2oPrtVPzut0mUGn+8XeCjYE5A+nZfegnN2qjPlD+X2/63qHM8TZBb94azhyQzgOEOl+ls/rzQ0AgrMBvV6lf1ZC3zkqBAQNN1ho3hup+8If2RDP6u4vOzM7o40qPdQLmchbZnEW2knNzo4fEhvWBbYv2AUDOCF/uSJ/Fp5ycXUhtmMbe8p067H3ltTOY0TmUf3uaDxDq3CX2eLN8gmHCD/Zqd63WSl3w0SSmQ2ND6QGDsG1uDv44jL6+L/2ntWrfxfL/DmVu6Hd0E4JoxL8z3HPasUu30o0ti+dc7Kp14Rk7bRMKqPXTrD41NLw/dnBlc8P6gLvEnjvS1++a4nRd32IuR55Qcm5uybm5oepY/Rr/hvm7HflC3ikZ2UO8lh1t9fLw6gTm3Spt6pfKzf2Z8zcEenWDKkzRlK/S2bwtnDvy6NHQdY3ktpVqswhPjT2i07hbwSA8kTw7PH8q8+v+9G9Xqs9v0eaPZc4p+vkPJbwvJng50883MYBKYGOeZ+sXgXlXOqaVWTdRVFGrX+s/+G2TEtfyTskYfkcv3mvtxE4X7hK7u8RePiW/aVPo4Krm3R8cyB3uKzg1055r9vk6x3FpOT02j/rdp7FhjWpmtsOEoxkM5+vjPioI62Jw7w/qJ9Xa/SOYWX3o7txXkv5f4l03IIP64jz2k2py8wp1UCb1xCl0mZsCgOat4Yy+6T9AuC9CfrVUzbW7b4vvHVNSaHZx2harF2uXNdWv9Xt7OcsuyvdVulJiEkeaoRgqa7Ana7BHbJYPrmz66bkqe65QeFpW5kC3BfcEL3RQ8zzh1RXuUz5WnhrDXFFh3Su8hMjo49r76SEgABQoGjy3RXtwnXpNJb11OneCRZbdBAZhe11YQk0qYh//SRv1oTJnIHP7ILp5W7j0vHReOAEA7+/RblyhzhnI3DnY/mM1F9wdbVmFbRH+7eGa/zSG98fysQloGUIGV3peXsk5uY0/BWuWNlR9dLDg9Mz8UzIZm7XCpmF94JxfFH3mYq/6Wv2ihvx1LONK3z8fIYPjnEy4JraBs920Qs2zwzcXsgk5yjENsABQVVW1dOnSQCAwbNiwCRMmAMCWLVs2btzYcqcLLrjA4UjPY4Y6RGDgnqH0L3tRc77Txrwj/rVW9PZM219LTIHfrVKX1JCPzmH1c8Czh3rr1wcsEoREIw3rA/u/biAqKZqQ3W9miWUHpbotiqGyh3qzh3rD+2I1Sxu+/9e2/NEZhROyLbJeJXpQVEXNU+oYTsGaS9hbvlVHfKC8eWY6TxURKlwvfx58zMk9MYaenvr7ICYQu3HjxgkTJlx44YXZ2dlPPvnkxIkTX3755ffff3/BggXDhw/X73TWWWdhELYodVHvT2K++Cq03mF/YZn2+CldPRHUgjY2kyv/rQ7JotZcwnr+e42cPdSzYf7uimkF5nZzaQo5tLq55usGwceVnZ+X0deNvaAW5yq29/lVsdgs1/ynYe0jO7IHe3qcnWP66tuG9YHsIV79j8fJwsvjmbd2a+d9odw1hJkz0Ho9uV1DAF7drr1Vb7/e37T5mjyLT3YzHlteXl5dXe10OgHgpptu6tWr14MPPggAkyZNev75580unnX1rIv0PcOzzwaD35MfGsnM6pM+Nef5Ldp9a9RHT2GuqTzimtGWyduyuMDOiK+3OSOjmkIOftdU83WDs9DW55c93Gl6GHK6EjK4nhcXFE/KPbCsccNTuzP6uoon5dhzTJtNU78+0PuqHq1vubwnfUoOdeXX6pIa7ZXxbBe347eO7QFy/XI1osCL09ziUwfsigpcKqz/NRDtdDr1FAQAp9NJURRFUQCwb9++11577d///reiKKaW0Io0WWvcFCwY5pk3mvnXeexL27QzPlG2+onZ5eqqRhEuWaK+vE1bfhF7VArq8kZnHFjeZHzBNIUcWN74w5+3B3ZG+s0u7X9tKaZgiuKcTMnk3BH3VNpzhB//WrX9jf3xRsn4Yvi3hymaOvb0tDI3texCdlgWNex95cualK/Rkgb3r9NO+1iZVkavnMIOy2N8fV0NG4Jml8tyjuisv+uuuy677LL8/Hybzaaq6tdff71y5UqWZZcuXZqZmdnm4+Px+EsvvbRkyZLWN44dO3bixIlt3l8URZZl9axNXU0/hpw9BMKroqj2dsK/z4G/badO/1i7oQ/cMZCkxGYrLSRJAgBCyH8OUbOXU9PLyD9OIzytimIbd/YOsu/97FBgf8hm1KEERCONa0M1XzU6CoReVxc4CwUAENssXLoghIiiyHFp3XtFQe54T9Zo18EVzevn78oY4Co6O5MzcOxw37/r8k73ilLbf0h/HAQTcqnZ/yHTy8jcYYSnQRRFnk+xuZXf1lE3r6R6eWDlBVqRQ5MlAICMwc6Dy5oyhqX8dWT7PxGe50+aOD9f9T/00EPff/+93h3629/+9rPPPluwYMFPP/2UlZU1b968Ez0FTVPHaE/5Ulfj+lDWsJ+XHtEU3NCHrL6Q/NhMnfJP+tu6FHv7kgZ/WEfNXE69MI78ZQQ5wdJzmqNzxngPLms2olgEmjaGN86vblwf7HVVQeXVhXoKorTB2Oiis7MG/a6UdTA/PbV332cNSlQ14HWjB8TYISlryIn2hJqQR1ZfRHaFqAmf09sCBhQqkfwS/GYVNWMZ9cch2jtnaEWtUs/b2xGrkyQ/9vMd4fAl2JNPPvnKK68sXbo0Kyur9Y8Zhpk8efLy5cuP93ibzTZr1ix9rml7yLIsCEJKJ6UcUcN74v2uKWWOTIxyAT48F97fo129XDu/mHp41HE3ZrOU7UHq/y2nil30+mlMTjsOFSiZkPfDX7b3vJDp6JEuHRLcHan6+BBRSa9LC80akjQLIURRFEFIhb+eRBAEqJjiKD4jt/qLuo1PVhedkV04Piupc4D3rqgvHJ9tc5zkz71AgA/Phb9v1SZ+CXcPcMwZJqTE19Y7VdptK7WppdTmy9o+hTFnmNf/U7R4Yo7hRUskSZISWEdoAPj73//+9NNPf/nll4WFh5dLE/Jz5/jy5csrKioS9XppoGF9ILO/mzlOu+mSMnrTZSxHw4B35Td2aQaXrUMIwNObtLO/oK6tJB9MalcKAgDrZHKG+2qXNSapVLF6ccvC6u1v1BSenjn0toruloLdFu9he00vHPzr8lB1dM1fttev8UNyRuhEv9y0NZQ/NuPkdwUAgOv60ssvYhdVMed/rtRGLT1qWBUiF3yh3L9We+ds5tlxxz2LOHeEr36N39iiWR27bt2666+/fsyYMXfeead+09y5c2+66aaKioqMjIxVq1bt3bv3xRdfNLeUllK3xl9yTu4J7uDh4JlxzK960dcvV1/Zrj07jrHgUSZVIXLtMvX/t3fn8W1Vd6LAj672XdZmrbbl3fHuJF7iLRtZSAJJ2ApvXhno8KDAY8pMGF6HDrRAh2605RUoFMpStkJY2oSQkjiBOI4dJ3Ecx3HixLEsW7J2a1/v1b13/lBq3ECCHcuWbJ3vJ598pCtZOh/L9/50zvmd30EJcPh6UjfDWKNulfQ9p1c2SRK7IAwL4mOf252nvOpV0qL/rYXrAtMQW84s+ecs30hoZJdlvH1Cd4NCmJfgdatjn9sVdWIaawYz+YVCyv416O/03OpPYr+uo/6v/JRbgYcR4Nf9xK/68X8vp/57+bdsrMLP4hA4GTCFeZrFkhc7azSNRvP+++9PPSSVSp977rmjR4/6fL6HHnpo06ZNbDb8fV0StESiLlRU+O0nZ52ccmIr7fmzROPu2L0lyA8rL237knQECV44SzzVi/9HJfXhMgTHZpyzx5IwlI3i4Q/NJXdnJaZJMdJyeML0hVNaJax5tIDOXVDpRlCiCXScyofyHKe8Q++Nc9WsnM0KdoKSszwXAt6hYPUj+TP9QRoCHq9GNmspd7Xjf9YTLzVS1dxU+aJ20Ew+2InnCyjHbqRdfdu4SyhAvkxkPerOvxle2C+hTB0FvQatra1PPvnk9OcIA4FAfJHGbN40ic6/ZeRp2OpV0un/iCUEdnTjHTbyF7XIbblJ/i7Z7yLvO4JTKeDVZmqhkAL+njU604w4Eid7nx3OWieTVgln1SASOPu8hj02rpKVszkzZWs0zyeSJIPBII8Hx4QBESPN7c7xLydkNcKs6+S02X1DwlGi9xcX825WXUOJYL/fz+fzAQAYAZ7pI54fwJ+ooX6/BEnu8mGDn3zkGNHjJJ9rQLZkzeDaggXxnp9dqNmxgEsSTn4iCTHfnZSIDb34VxtCpyB0BKFREAZC41DpHCqdR2MIaAwhnZlBZwrpqVkrJOyIeoaC+bfObDNrJQe8s4p62Er+axf+wlniV3XUeMWyeebHwFO9+JtDxNNLqf9SPNvzl0KlFHxHfe61UWEB75o7cN7hoOFTK4mDgtvUKVK5DUopCI2iWS3LrBOPfW7r+fmQulWiapFc847Ko5/ZBHmcWRbKpyPg8WrkFh3lvg78T0PEC43UZckoyebDwM/78D8MEv9aSn1rJXUmA70AAEDnUjOXZ5i+cOZuTdEtqPAIEXGhUQ+G+WOoL4YFY7EQTqBELEIUfEfNFCU4fs93IGSI6Xm3KMkYIDCCwEgCI7AgHgvhERfqGwmhPiwygcZCOEvC4GQy2XImV8XiqtlsKSMVQqOpzalqllCZ13IeNisoJ7bS3hgibmrDGzMpTy1F5m3ikCDB6xeIx3uI9RpK/3Z6oupl8LPYmXUZZ181lN6bM6MZFwBAwBQe3WsP26LZ18tl1aJU+HChlEXnUvO2q1TNktHPbD3PDGmvk2XWZsx0g2VLh8t1xl/1b4nJ+ysRUb7cTHvzAnHjPnydhvLfy6jK+VqYF8XBHwaJn57Cr9cip7bRrnmEVr1SevIXQ9q1Mjov+XM2sTAeMEWC5nDYFg3ZomF7lMBIppjBzKAzBDQGn8bMoPM0bISO0NgIjZ34qZP5/hUgdApXw7760CgRI8OOaNgeDVmjjl6v4VNbLITzstj8HI5Qx+HrOFfK2JxTERfqOutb+ljhNb8CQgF3FyLfyUWeGyAad8c2ZSE/qkLyBHMYBEgA/jpKPN5DCBlg1zrq0kR/dc3ekIlHiYGXDKX35kzzr9M/FjYdcASMYc1qmeLuLLhfPDRNbBmz+M4s/1h47G820wGnZq1Mvkw0zYwqS8fE+JcT5Q/oaJyEXUMpAPxzIXKTDvnvU3j5R9g9xcgjFVTxXA7towR47TzxTB9RIQb7N9LKxbM6dxgCmqxGNP7lRM7mZGyhQ4KQLeodDvoNIf9YCPXFuCoWT8PmadmypSKOnDmnq7O+bmHMEcaCuG805B8L+4aDAVOYq2SJiniiIh4/iz1vBaAvvGNiiumJ2nfJi4LfniFeOIuvViGPVCAJD1EECT4xED/rI2IkeGopdXPWFV//2uYIv0IC/V8t3ovB/JtVVyl7RuKka8A/3j6BejBVq0TRIIZJoVcC5wi/lc8QMrU5AsawolGsqMu4ykRXLIIbP3dMnPGV369jZlz7eNrVZ6SMQfLpXuJjA3FfCfLgksRX4fei4OVB4v8PEJVi8OMa6vIEza1EPVjvsxer/y1/Nr+Zmb6j53zAfT7gvRikshBhHleQw+HncDhy5kyHhRI7R7gwAuFURIz0jYQ8FwKeQX/Eg2UU8yWl/IwS/rWNWE7TRL/P8Km16t/zE9sZDWDgj+eJ35whlBxwbzFyay4y+8xSZwS8fZF44SyRyQY7KpAbs7/lm8JsAyEAgAT2Hs/oZzZ+Dke5QszTsic/C9QX8xtC7vOBidNejpKlXCGWVAhScI/WlAID4TSFbFHz4Qlnn5enYknKhYJcDkfBvPTXRYKwI+oeDBgPOMRL+DnXZ86yhzGdy67eT/66n3hvmNieg9xTjCQkFeC4g/zDIPGRgdisRXZUIBWz6wV+nemAw3MhWHZfzpxOTwRM4Yl+n/usP+LBMop4oiJeRiFvlnk66R4Ip4p6MPc5/8QZv38kJMznSisF4jJBwiMi5o/1Pnux5K4sfvaczAPgJPjMSPxhkOiwktdrkVtyKWtVyEw3CPWiYK+R+NBAHjQTW7KQ+0qQBvm0fskJCIQAAAAIlDAfnnCd9QfHI3QeDVBALIQDBAiyOcI8rrRamPD57cUKBsIZIWKk+5zfddbvN4QibozOoVKolFgIp7GpglyOeqWUq5peqYirmv5l1xkBr54nXrtAsKjg9jzkxmzKkhlufksC0O8iPzYQH+hJlAD/UoTcVYjM0V5vJEH2Pz8irRaqmiXf/uwZ8o+Fnb3eiX4fhQok5QJxqYCfnbAxPBgIv0EsjLvO+p29Xp8+JCriympEGSX8xAy+keDsa6NcJSv7+jkfSXdEwMcG4sMR4qidrJZQWpWUSjGlWkLJ4lG+niiHEUDvJwc95FE72WEjT0+QrUpkaw7lphzkShUlvlGiAuEkEiejHgwAQGNTaWwqTISZKRgIrxmOErEgThIklYEkdpJpppddEoDDVvLDEeKvoySNApoVlMZMSqWEUiCgfGPlxYkoOOMi+1xkt508aCa4dMrWbMotOqRWPudLzcIO9PTv9BUP6hK1eClkizp6PI5eL4VKkVUJpZUCjjIB30UuAwPh1cTCuLPP5+jxBK0RaaVQvlQkyOFc+7WYBBc/MofMkfIHdPOZ1hGKgQ4becRK9LnAaRc5HiTFTCBnXwqHwRjwREkPCtRcSpEQ1MqQxkxKQyaFe00nfsIDITRLMBCmoNlcdgfcZIeN7LKRZ9zkBS+JUICE+VU49KJgPESyqKBERKkUU5bLKKtVlGzevH5/tB1zG/c5yr6fM5vdklF/zNnrtZ/woP6YrFooWyriqRMf/yYt7HWEc43GpirqMxT1GVEP5ujxXNw5TmCkbKlIXiOc6fcdkiAvvm+OuNDSe3PmObmRQwPr1JR16ksZbiQA1hBwREiMuPSoiEGRsEAykmchCJqZ0gxKaQbl3uJLd11R4I6Snr8XdBLQgYZLYSf1SpxZm0ESoP+FkbL7cmZ6ncRRwtXvs5/0+g0hSZkgZ4tClM9dcONAiy0QTmKK6Jo1Ms0aWXA8Yu/x9P/eQOdRZVVCSYVgOptih2xR/ccWgIDSe7KRZAccCgBKDlByFtofFwRBXyNmAnHq7WOhqM9AaJTTL4xkb8xU1GV8ayTDo4R7MODs83oGA4JcjnyZqOSftddc6yDpFm0gnMRVs3RqhW6LwqsPOvu8/S8aqAxKRglfmMfl53AY/ziRQBKkfzTs6PE4T/u018mUjWKY4ghBUDqQLxNxFEz9X6zWIxPqVVJRIe+ytfYkTgaMYd9oyHM+4DOE+NkcWZUw/2ZVAldnJsviD4SXUIAwjyvM4+ZtAwFzxDPot3W7hz4Yp1AoTBGdIaARMRKP4mEbypIyMkp4Sx8tmGVhQwiCoIWFp2FXPKBznvY5ejzDH1kYAhqdT6MyEDxKRD0Y6o9xZAy+jpNZLy6+M2tOV6zNs7QJhJMogKdmTc7iov4Y6sFQfwyhUahMKktCT4WCQxAEQclBAdJKgbRSQBJkyBqNhXAcJagMhCmiM0T0xVoHI90v+gw+jTG/tXwgCIJSHwWhJGQJ5oKwePq2EARBEHQNYCCEIAiC0hoMhBAEQVBag4EQgiAISmswEEIQBEFpDSZMQtB8C2FhnMRJQAbQ4OTBSCyKERgAgCTJcDjMQS9tdcJnfFV0lEllMqh0AACXzkXmvBozBKULGAgh6BqhOOqL+r2o3x/1+9CAL+r3o4EAFgyioVAsHMbCQSwUwsJRPBqJRePBL4AGSUBy6GwqhQr+McixaEw6cmmnKoIgEOTSaI0fDUw+J4JHMRwDAASxIEGSbBqLhtA4dA6DSufQ2Gw6m0VlcuhsLp3DY3B5dC6PweUxuAIGn8fgCZl8EUvIpc/JVmIQtKDBQAhBV+SN+pxhlz3onAi7JsIuV8QzEXa5Ix5X2OOOeGIELmTyBUy+gMHnM/kCBo/P4PEYXIlAzKGz4wGJTWMzaQw2jRUPftPpyU1/94l4cA1hoSiOhWPhMBYOx6KhWDiEhfxoIIgG7SFnAAv6owEf6vdGfd6oD8MxIVMgYomkbLGIJZCwxRJ2hoQtlrIlMo5EwhbTEFhQCUo7MBBCEHBHvJaA1RywWgI2a9BuDdrtQactaGfSmFK2OJMrE7MypByJTpi1TFElYglFTIGYnZH03hWHzgb/2K38VhiOeaJed8Trirg9EZ8z7Br3W/vsAxNhly3o9EQ8QqYgkytXcGWZXLmSlxn/p+BmwgAJLWIwEELpJRyLjPlMYz7TmHfc5Deb/OZxv4WKUNU8pYqXqeRlFksKWrNWZHJkmVw5i5aYrUpTB51Kl3GkMo70Gx8lSMIVdluCdlvQbg06zrsufjl2xBKwOUITMo5Ew1dpBKosgTpLoMkWaKScxO9pDkFJAQMhtJjFCNzgHdN7RvUew4hnzOAd80S9Wr46S6DWCjSNmloNX6XmK2fUqVrEEAoi5UikHEm5rGTq8RiBW4M2k9885hsf8YwdGusc9ZmisWiOMEsnysoVZeuE2fkZOgEzYRulQtB8goEQWlQiseiQW3/BdXHIpR9yjxh9JhVPoRNl52XobijYoBNlKbiZMN9ypmgIVcNXafiqetWyyYM+1G/wjBm8xmGP4cuxzmH3CIfOKcjQFYhzCzLyiiX5V+p3QlCqgYEQWtgIkhjxjp1znj/rvHBu4oI5YNMJswrFeWWy4m2Fm3SiLAaVkew2Lk4CBr9CXlohL508YgnYLrr1F1z6PcP7fn389wCAYnFBibSgVFpcLClI+pQqBF0JDITQwhOJRQac5/sdZ/sd5845L0jY4iXSwiXSoq2F1+eKcmBaR7LEM2uatQ3xu/aQc3Bi6Kzz/Jv9fx5y6xVceblsSZmspFK+JJMrT25TIWgqGAihhSESi5x2nO219vfZB4Y9hoIMXYW89KaiLWWNxXBqKjXJOVI5R9qibQAA4CQ+5NKfcQweMXX/vvd1OkKvyiyrkpfVKCoUMChCyQYDIZS6CJI467zQYz113HLqoltfJM6vVpT/n+rvLpEUwgHPhYVKoRZLCoolBTcXbwEAjPnGT9vP9Fj7Xul7i0llLFVULlNULVVUwu80UFLAQAilHGdo4qi557il96T1tJwrW66surP8tnLZksW3mCFtZQnUWQL15vz1AACD19hjPbVv5Itfdj+fLdQsV9bUq5YWSwoQCqyEDM0TGAihlECQ5ODE0JHx7qPjPfaQY7myulFT+4Pl92awRMluGjS3coTaHKH2pqItGBEbcJzrNp/8VfcLzrCrTlWzQl1bq6qBWTbQXJvvQBiORUwuC+Xv+evxIsIsKpNOpcO1XGkIw7FeW/9h09EOU7eQwV+hqf3B8ntLpUWwN5CG6AitKrO8KrP83uo77SHn0fETf9Mf/GX380ukRc3a+kZ1LVzCn4aiOIriaLzibvx/nTDxqeDzHQjtYedzp/4weTeKR1Eci+BRFEcDaJBDZ3PoHD6DJ2DwBEyBmCXKYInEbFG8EKKMI4H9g8UhiqPd5p5DY53d5h6dKLtJU/fCup+reIpktwtKFXKO9IaCDTcUbIjEIscsvR2m7lf73tbwVa1ZK1ZmNcL8msUhiqO2oMMZmnCGJybC7omwyxPxeVGfN+rzRwMBNBjAgjSExqQy4t2k+P9Ptvy/hP8BUEiSnM3Pt7a2Pvnkk62trdN8fiAQ4HK5lCusaA5h4QAWDKBBH+r3RnyeqDde3dgRmnCEnI7QRCgWVvIyFVy5mq9Q8ZRZArVWoFZw5bADcc1QFAUAMBjzkXuC4mi3ueeL0Y5uy8liSUGrdkWTpk7MzpiHt15Apl90O63ECPyUrf+QsfOwsSuTK1+V3bQqqymTK5ufd/f7/Xw+TOS5dn40YPSNj/nGTX6z2W8xB2zWoC2IhTO5MilbLONIxSyRlCMRMvlCpkDIFPAZPD6Dx2Vw4vu0fMMLJvQTSa05wnjNfvmVC1JEYlFr0G4J2Mb9lvGA+aj5hNFn9ka9WQJNjlCbK8rJy8jJF+ngtTWlECRxwnqqzdDeaTpWIM5dnd38r8vvFTIFyW4XtJDQEOoyZdUyZdXDy+87ZTtzcPTwPXsf1gpUa7JbV2U3ZbCEyW4g9BUUR0c8Y8Mew7BnJF7aMBKLagSqLL5GK1A1aJareUoFV546F+rUCoTfikVjxqfWpx6MxCKjXtOId2zYYzh+9pMht56B0AvFeUWS/BJJYbGkAF5zk+W86+K+kS8OGA4refK1Oa33Vd2ZOn/60AKFUJAaRUWNouIHy+87Ye09YGj/4+m3S6XF63QrmzT1MLU4KXAS13tGzzrPD04MXXANG33jWoE6L0OXK8quVy3LEWpTvN7eAguE34hFYxVJ8osk+ZNHbEH7BdfwoOvi++f+MjgxJGZnlEqLymUl5bIlWUINBcBSk3PLGZrYZ/jyc/1BjIit0618Yd3P1XxlshsFLTY0hFqvWlavWhaJRTtMR/ePHPrt8ZebtPUbcldXykvhaT7X/Gig33HujOPcGce5C65hBS+zRFJQIincWnB9riibTqUnu4EzsBgC4ddlcuWZXHnvBZdgAAAUCklEQVS81BNBkqPesTPOwT772bcHPozEIpXyskp5WY2i4rKeJTRLKI52mLr36g8MOodas1bsqHuwTFYMr0fQXGPRmGtzWtfmtLrC7rbR9t+deCWABjfkrtmQu1rJy0x26xaVIBbqtfWfsvX32s5YAtZiSUGFfMl3y24rkRYu6FUuqZUsMw8cIWev7cwpW/9J22kUR5cqqpYrq5crq9I2HzUhyTJDbv1nw/sPGA4XivM25K5u1jYwYeWXawWTZWZvyK3fO3zggKE9LyPn+ry1LdqGWSbcp3OyTHx953HrqR5L36jPWCotrlFUVMnLiiT5V8pkmQeJ/UTSLhBOZQnYTlhPnbCcOmk9reDJa1U19aqlpdLitMpBnU0gDKDBNsOhPcP7/WhgY+6ajXlrr5LoBE0TDISJghGxI6buPcP7ByeG1ua0bMpbl5+hu7aXSsNAaA85u8093eaek9bTWoF6mbJqubK6VFpMR1JiHBEGwsQjSGLAOdhtPnl0/IQt5KhV1qzQ1NYpa3gMbrKbNueuLRCetg98Orz/iKm7VlmzOX9ddWYF3OQvUWAgTDh7yPnZ8P7Phg9ksISb89etyW7h0NkzeoU0CYQESZ53DR0xHesaP+4ITdSpaupUS5crq1Mw3xAGwrnlDE10mU90mo712QeKJQVNmvpGTe28LVeafzMKhN6o73P9wd0X9yEUyqb8dRt0q2GV5ISDgXCOECR5wtL76fC+k9bTLdqGzfnrlkiLpvmzizsQYjjWY+vrMHYfGT8mZAoa1ctXaGpLJIWpPDa2mNcRpgIpR7Ilf/2W/PWRWPSE9VSHqfvN/j9ncmXN2voWbUN2WubXECTZazv96cV9xywnmzR1/1H/f8tlJcluFATNDEKh1KpqalU17ojnb/qD/931WzpC25y/bp1uVXrWdwxh4aPmE+3GruOW3jxRTpO2/o7Sm9KzwBPsEX47giRO2wfajUcPG7vYdHaLtqFF21Aozkt2uxLj6j1CZ2hir/7AnuH9PDp3c/6663QrF3Ru2IIAe4TzgwRkn+3Mpxf3d5mP16uWbc5fV5VZdqUk58XUI/RF/UdM3e3Grj77QIV8SZOmvklbJ2IusIoEcGg0aUhADk4MtY91fTl2BADQmrWiRdtQIi1c0CsEvjEQYjh2ZPzYZ8Nt55wXVmU3bcq/rkicf4UXgBIMBsJ55kcD+0a+/Gx4fwgLb8xbsyF3zddzvhZBIHRHPIeNRw8ZO885LyxXVrdoG+rVyxbu91oYCFPCRbe+3dh1aKwziIWbtfWtWSsqZEtSeUj9Si4LhOcnLv5t5OABQ3uuKPv6vLWtWY1wIcQ8g4EwWS64hj8bbjsw2l6Qkbshd02LtmGyTs3CDYT2kPOwsat9rOuiZ6ROtbRVu6JOtXQR1N+BgTC1jPlM7WNdh4yd9qCzSVvXrK2vyaxYQPunxwPhBOo+YGjfZ/gyRsTW61atz10NC/wnCwyEyYXi6BHTsc9HDp5xDK7Q1F6X07pUURkMBBdWIBz1Gg+bug8buywB2wpNbYu2YbmiamFVe7k6GAhTlDVoP2w82mE8OuTWL1dWr1Avr1MvTfGRd3vIeXDk8CFjpyVobc1qXKdbWSotTnaj0h0MhCnCHfEeMLQfGD1kCdgblbXX5a9M8VEfgiT6Hee6xo93mLpRHG3U1DVr6qsyy1K5zdcMBsJU5436usZPdI4f67H2ZQk09apldaqaQnF+6qy0G/Uaj4wfO2zsGvdbG1TLWjUr6rRLk1gkApoKBsJUYw5Y955v67KdcIZdzdr6Rk1tSo36uCOeY+aT3ZaTxy29Cq68Qb28SVO3aLL5rgQGwgUDI2Kn7QPHzCe7zT0TEXdNZsVSRWVVZlmWQDP/jQlioT77mW7zyWPmkziJN2pqmzT1VZllOIaD+dqPEJoOGAhTUPyyaw5YO4xHj5iODbn1lfKyOlXNUkWlVqCe//YEsdBp+9mTttM9llP2kHOporJWVVOnrJFyJPPfmKSAgXBBcoZdPZZTJ22nT9nOoARWJi0uk5WUSAoKxHlsGmuO3tQVdg84B884BvvsA6M+Y4mksFZVU6usyRVlTz5nPjfmhaYDBsIUdNll14f6eyx93ZaTPZZTBCBrMsvLZCVlshKdMGuOxiEJkjT5zYMTQ2ed5/sd5+IFr6szy5cqqool+Yty8PPq4IL6BUnKFq/PXb0+dzUAwBZ0nHGcG3Ce/3KsQ+8Zy+TKckXZOmF2tlCj5ivVPOVM6z/FBdCgOWAd85lGvGN6t+GCaxglsFJp0RJp0fdr7loiKVxMU+UQlEQCBn9VdtOq7CYAwLjf0mcfOO04+9H5T21BR64ouyAjVyfKyhZo1XyljCO9hjkRgiRsQYc5YDX6xvWe0RHv2EW3XsgUFEsK4jsvForzaQicy0gY2CNMMpzEDV6jwTOm946OeU3jfst4wEJH6FKORMoWC5kCAZPPojL5zEv9AyaVGcWjAIAQFg5iwQAadEe8E2GXIzRBkISSp8gSqHOEWTpRVqE4bzqZn7BHmGpgjzAFTbP/EcLCQ279sNtg8I6N+kwmv9kb8ck4EglbLGaLBAw+j8Hl0bkIgjCpDAAo8XM5EouiOOqN+nxRvyfqtQed7ohHzM5Q8RRagVonzM4RagvEuelZ/uZKYI9wUaFSqHminDxRzpopB31RvzPscoYmvKjPG/FH8ag/Gog/5MAnWFQmAIBNZ2Ww1HwGN4MlkrDFUrYYlv2EoOTi0NmV8tJKeenkEQzH7CGnK+KeCLv9aMCPBoJokCBJJ46SgIyfy0waQ8gUaAVqAYMvYgnlHKmELYYdvvkEA2EqEjD5AiZ/6kweBEELEZ1KV/OVar4y2Q2BribtplghCIIgaCoYCCEIgqC0BgMhBEEQlNZgIIQgCILSGgyEEARBUFqjAQCOHTv2xhtvAADuvPPOurq6+AMffPDB3r17MzMzH3zwQY0mCSXBIAiCIGgeIH19fWvXri0qKiopKVm3bt3JkycBAC+99NKjjz66Zs2acDjc2NgYCoWS3U4IgiAImhPUSCTS1NT0+OOP19XVeTyetra2rVu33n777b/97W9vueWWjRs3vv/++0wms7q6+ht//o033li1alVOTs403w9FUQaDASvLpA4cxwEAVCpcvZtCMAyDtX5SCoqiTOaC38x2MUnsJ4J0dXWtWrUqfmflypWdnZ02m214eHjy4KpVqzo7OxP1fhAEQRCUUmhWq1UqlcbvyGQyi8VitVqZTOZkGTeZTHb69Okr/XwwGHzssccmXyFu8+bNd9xxxzc/36QPtn8Ee4SpgyAIAACCwLSpVEGSJI7jIRqs+pRCYrFYGH4iqYGz5V8QgTgcDk9zHIvFYn3r9Y3GZDLjZZcBANFolM1ms1gsDMMIgoj/cCQSYbOvuBkCg8G47rrrSktLpx4sKChgsb55ayEsQ0pbvgYGwtQRi8UAADR4kqcMkiSj0eiVziAoKcLh8FUug9B8YgozKEwWhmHTPEem8y2fptFojEZj/I7RaFSr1SqViiRJs9kcTxaNH7zSz9Pp9JUrV05/9wkqV8CpaoaBMHXA3SdSDUmSZDDIgbtPpBLc7+ckbq8DaPYQBEngOBaybdu2d955hyRJkiTffffd7du3CwSCNWvWvP322wAAn8+3e/fu7du3J+r9IAiCICil0O6///4PPvigqakJQRCPx/P6668DAJ555plNmzYdOXLk/PnzLS0tLS0tyW4nBEEQBM0Jmlgs7unp6erqIkmyoaEhPkS2bNmy8+fPd3d3y+XyKy2cgCAIgqBFAAEA0On0lpaW1tbWqRNFIpFo/fr1CY+CRqMxEAgk9jWh2XA4HA6HI9mtgL4SCAQmp+2hVIDj+NDQULJbAf0DvV4/meY5e/OdNP/QQw91dHTM85tCV/H73//+5ZdfTnYroK+0t7f/4Ac/SHYroK94PJ61a9cmuxXQP9i+ffvo6GiiXg2uHoMgCILSGgyEEARBUFqb7TLqWCzW19c3/ee73e7+/n4OhzPL94USZXR0FEGQQ4cOJbsh0CVnzpxxuVzwE0kdPp8vFovBTySlRCKRY8eOmc3mb33msmXLuFzu1Z9DIUlyNq157733nn/++enXJQkEAiwWC9YxSR2RSAQAAOuYpI5YLBaJRHhwQX3KIEnS5/MJhcJkNwT6is/n4/F401lT/8c//jE/P//qz5ltIIQgCIKgBQ3OEUIQBEFpDQZCCIIgKK3BQAhBEASlNRgIIQiCoLQ2t9mbGIYdPHjQ5/OtWrXqss17AQBHjx41GAw1NTWFhYVz2gxoks/n6+rq8vv95eXlRUVFUx/q6emZvC2TybKysua9deloYGAgnrgLAODz+ZedC8PDw8ePH9dqtY2NjcloXTqa+okAAIRC4WTO4fj4uNVqnXyoqqpqmnvDQtfGZrOZTKa8vDyRSBQ/4vf7Dxw4gCDImjVrLlsUgeP4F1984XK5mpublUrljN5oDrNGo9Ho6tWrSZLMzs5ua2s7cOBARUXF5KMPP/zwrl27WlpaPvvss1/+8pff/e5356gZ0KTBwcG6urrly5fLZLJ9+/bdf//9Tz311OSjVCq1ubmZTqcDADZt2gSrfM2PJUuWcLnc+HleWVn5q1/9avKhDz744IEHHti0aVNnZ+fq1atfeuml5DUzjdx1110mkyl+++TJk9u3b3/llVfidx977LF33nmnoKAgfveTTz6Bq1zmTllZ2cjICIZhO3fuvPHGGwEAFouloaGhvLw8FotdvHixs7NTJpPFn0wQxMaNG51O55IlS/bu3btr164VK1bM4M3IOfPWW29VVlZiGEaS5I9+9KObb7558iG9Xs/hcMbHx0mS/Pzzz1UqVfxp0Jxyu91mszl+u7e3l0KhOByOyUcRBJl6F5ofJSUlHR0dXz+O47hOp/vkk09IkrTb7Xw+/+zZs/PeurQWCoVEItHhw4cnj/znf/7njh07ktiktHLu3DkMw4qKiv7yl7/Ejzz66KO33XZb/PaNN974xBNPTD75008/zc3NDYVCJEk+++yz8T7Y9M3hHOGuXbu2bdsWXzt/yy237N69m/x773PPnj0NDQ0qlQoAsHbt2kgkcuLEiblrCRQnEokmRwy0Wi1JklOHgAAAPT09nZ2dXq83Ga1LX+fOnWtvb79sD5D+/n6bzbZp0yYAgEwmW7ly5e7du5PUwDS1c+dOuVx+2aC00+k8cODAhQsXktWq9FFcXHxZ6ZVdu3bdeuut8dvxmDL50O7du7ds2cJms+MPffHFF36/f/rvNYeBcHx8XKPRxG9rtdpoNOp0Or/+EIIgKpVqfHx87loCfd0zzzyzdu3ayU8BACAWi3/xi188/PDDOTk5H330URLbllY4HM5bb731X//1Xzqd7je/+c3kcbPZrFAo4iPVAACNRjOdalJQAr322mvf+973KBTK5BEKhdLX1/ezn/2ssbFx8+bN0Wg0ic1LQ5fFlKlRY+pDarWaQqHM6HyZw2SZWCw2OZMcv4FhWPwuhmFTS+PQaLTJh6B58Nprr+3cufOy/bDMZnP8svvuu+/efffdGzZs+NYCfdDsdXV1xX/tR48eXbly5aZNm+L5MhiGTU3EoNFoCdx9DfpWer2+s7Pzvffem3rwiSeeePrppwEAPp9vxYoVv/vd73bs2JGkBqajy2LK1DNi6kMUCgVBkBnFlDnsESoUisnRHrvdTqVS5XJ5/K5SqZw6EGS322ea5ANds3feeefxxx/fv3+/Vqudenyy8/Gd73wnEonAnUjnx+Svvb6+Pi8v79SpU/G78XNkcjbBbrfHpxKg+fHqq69u3LjxsuvS5IclEAi2bt06NdEamgdTA8dlZ8TUcONyuWKx2IxiyhwGwpaWlra2tvjttra2xsZGGo0WjUYjkUhLS0tHR0d8gmpgYMDv99fU1MxdS6BJH3744SOPPPL5559PpumjKBoOh6c+Jz5HrVark9HA9OVyuYxGo1qtjsViwWCwrKwMQZD43DmGYYcOHWppaUl2G9MFjuNvv/323XffPXnE7/cTBDH1OX19fVNnFqB50NLSsn///vjttra21tZWAEA4HMYwLB5u4l8c29raysrKJBLJ9F95DpdPuFyu8vLybdu25eXlPf3002+//fbGjRu///3vB4PBP/3pT+vXr6dQKDfeeOOLL754ww03/PSnP52jZkCTBgYGqqqqNmzYUFpaGj9y7733vvXWW4cPH37ggQc+/vjjmpoar9f7yiuv3HTTTc8991xyW5sO+vr6fvKTnzQ0NJAk+eabb2ZnZ+/Zs+fPf/7zD3/4Q4PBED9rHnroob1793q93vb29mS3N13s2bPne9/7ntFojHcBURRlMpk9PT0//vGPq6urpVLpwYMHjxw50tPTc9mwCpRAL7/88sjIyCuvvLJy5cqCgoL777/f6/U2NTXt2LEjFos999xz3d3dRUVFzc3NW7ZsefDBBysqKpqbm6urq5955plnn332jjvumP57ze3uE2az+fXXX/d4PNu2bYuv6jh06BCGYfFM0VdffVWv19fV1d16661TZ6ShOWKxWC7LPNyyZYvJZLJYLHV1dTt37hwZGeHxeM3NzWvXrk1WI9OK3+/fuXPnuXPn6HR6TU3N9u3bEQQZGhrq6uqKr6z96KOPjhw5kp2dfc8998BdPOdNV1dXKBRas2ZN/C5BEC+++OJtt9124sSJrq6uQCCQn59/++23Z2RkJLedi9vu3bstFsvk3a1bt8rl8oGBgXfffRdBkH/6p3+KlwT5+OOPdTpddXW1w+F47bXX7Hb7pk2bVq9ePaP3gtswQRAEQWkN1hqFIAiC0hoMhBAEQVBag4EQgiAISmv/A2EIY1z4Jv+/AAAAAElFTkSuQmCC", "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": 13 } ], "cell_type": "code", "source": [ "ρ = real(scfres.ρ.real)[:, 1, 1] # converged density\n", "ψ_fourier = scfres.ψ[1][:, 1]; # first kpoint, 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": 13 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.3" }, "kernelspec": { "name": "julia-1.5", "display_name": "Julia 1.5.3", "language": "julia" } }, "nbformat": 4 }