{ "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.293818258993 NaN 4.29e-01 7.0 \n", " 2 +0.441704844053 1.48e-01 5.35e-01 2.0 \n", " 3 +0.260910645994 -1.81e-01 1.72e-01 2.0 \n", " 4 +0.243051307531 -1.79e-02 2.68e-02 2.0 \n", " 5 +0.242665843579 -3.85e-04 6.21e-03 1.0 \n", " 6 +0.242635943577 -2.99e-05 2.16e-03 1.0 \n", " 7 +0.242634049067 -1.89e-06 1.04e-03 2.0 \n", " 8 +0.242633425565 -6.24e-07 3.24e-04 1.0 \n", " 9 +0.242633377584 -4.80e-08 3.86e-05 1.0 \n", " 10 +0.242633376344 -1.24e-09 7.55e-06 2.0 \n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown:\n Kinetic 0.0304516 \n AtomicLocal 0.0972473 \n PowerNonlinearity 0.1149345 \n\n total 0.242633376344 \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.038687275293251125, 0.0, 0.0]\n [0.038692961619048116, 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+gvaeTAAAgAElEQVR4nOzdeXxU5dk38Ossc87sM8lkBxJCCPsOIqKCCogWC4riVpdHsA8udWutVfu0Fm2tuKJ1bS1oFX1V3OtSlxZlEVBAFGQnEEiArLPP2e/3j0NjgABZ5iwzub4f/iCTWe5kcs/v3DtFCAGEEEKou6KtLgBCCCFkJQxChBBC3RoGIUIIoW4NgxAhhFC3hkGIEEKoW8MgRAgh1K1hECKEEOrWMAgRQgh1axiECCGEujUMQoQQQt1aV4Nwx44diUSi/ffXNK2Lr4jSixCC2+zZDVYTu1FV1eoioMOk9x3pahDOmTPnm2++af/9k8kkfuzaiizLsixbXQr0I0JIMpm0uhToMPiO2E163xHsGkUIIdStYRAihBDq1jAIEUIIdWsYhAghhLo1DEKEEELdGgYhQgihbg2DECGEULeGQYgQQqhbwyDMMPHqFNFwRwKEMlu8RtAUrMh2wVpdANQBqXrp+6ernHlc31klvjK31cVBCHWYKmh7Pjq4f2VTj4l5vc8rtLo4CACDMJMQ2PF6Tdm0QoeH3byoumRCXs+z8qwuE0LdyIYNG7766quuPIOcUGqXNrqL+GA/73vPNpR8H+JzHOkqXnbLy8u76KKLDHpyDMKMsf+rJqKSktNCQIGvzLXh8V09JoYohrK6XAh1F59++unLL788bty4Tj8DUQkUQTND1ewEUgLNO/fTLFbhE2tsbNyxYwcGYXcnhuXqj+uG3lgOFACAM8S5CvjmLfHcwT6ri4ZQNzJ58uSHH37Y6lJ0O+vXr589e7Zxz4+TZTJD/dpw3oiAu5BvuaVgTLDum7CFRUIIoeyAQZgZmrfGcwd6W9+SN8If3hZXUnhMGkIIdQkGYQZQJS2+N+Wv8LS+kXUywX7ehg1Rq0qFEELZAYMwA0S2J3xlboY78s3C3lGEEOo6DMIM0Lw1HuzvPfr2nAHeVJ0oRfB8eYQQ6jwMwgwQ3hrLaSsIKYby93ZHd6fMLxJCCGUNDEK7E5okVdA8xc42v+stdcWrkyYXCSGEsgkGod2Ft8SDA7xwjEW3vlJXrBpbhAgh1HkYhHbXvDXeZr+ozlvqiu/DbbgRQh0Tj8dbdgZYt27dl19+2ebdGhsbX375ZRPLZQ0MQruL7U76+3iO9V3WyfBBR/KAaGaREEKZLhqN/u53vwMAVVXnzJlTUFDQ5t1CodCTTz65fv16c0tnNgxCW5PjiqYSPni8bXm92DuKEAIghDQ2NqrqoU024vF4y/9bbkmljvys+PDDD/Pz8wcMGHCsp73uuuseeuihtJfWVjAIbS1RI3hK2p4m08JX6sb5Mgh1Z7fddtv1118/bNiwUaNGfffdd19//fWIESNOOumkXr16/elPfwKAbdu2DRgwYOzYsf369ZsyZUpjY2PLYxcvXjxjxgwAuPjii19//XUAyM3NjUajH3300TnnnAMA06dPf/vtt5PJbP6QwU23bS1R254gdB34qsmc8iCEjvBtI6k1NyMmlVA8c9gtyWTyo48+WrVqVUlJSSwWGzx48KJFiyZNmhSJRMaPH3/qqaeOGDFixYoVoVBI07Qbb7zxwQcfvOWWW/THrlix4tZbbwWAeDwuiiIANDc3E0JkWY7FYgCQm5tbUlLyzTffTJgwwdSf00QYhLaW2C8E+h5zgFDnKXEKjZIqagyP7XuEzPavfeSLA5qZr3hyAXtEEALAxRdfXFJSAgArVqxwOByEkM8++wwAhg8f/vnnn59xxhmrV69+5pln4vF4fX39jh079EcRQvbv33+sAcIWhYWFNTU16f9JbAOD0NYStULJhNDx70MxlLvIGd+XClScIDIRQmn3m+H0b4Zbfw2al3fomO66ujpFUfQUBIDS0tJRo0a9/vrr//d//3fHHXf07NnT4/F8+OGH+ncpiuJ5XpIk/UtCfpx/TgihqEPLtgRBcLlcJv0kVsAgtC+iklS91PropWPxlbni1RiECCEYOHCgIAj33HNP6+i64oorbrnllmuvvRYAvvnmm9b3r6ys3L1794ABA3r06LFnz56W26uqqnr06AEAhJDq6urKykqzfgILYBDaV7JOdOY6aMeJLza9PV3NW+ImFAkhZHMnnXTSaaedNnPmzFtuucXhcKxZs2bcuHEDBw5cvHjx0KFDt27d+sILL+Tn57fcf8qUKcuXLz/nnHOuueaaGTNmhEIhAFi0aNH8+fNfeOEFANi8ebPT6Rw0aJBVP5EJMAjtqz0zZXTuAr72y8YT3w8hlI2mTp2am5vb8uXrr7/+wgsvvP7667IsDx48eNCgQePHj2cY5qmnnqqsrPx//+//rVmzxufz3X777QBw7bXX/vSnP7333nvHjx//+eefv/TSSwCwc+fODz/8cOTIkQDw6quv/vznP2/pJs1KGIT2lWx3ELoK+VS9CASOtRMbQiiLzZw5s/WXDMPMmTNnzpw5rW+88847W/4/ZswYALjvvvsAoH///qeeeuqbb745a9asYcOGPfTQQw8//PAf//jHQCAAAPF4/K233lq5cqUZP4Z1MAjtK14rlJx+gpkyOoanGZ4WI/Lxl94jhNDRnnvuOUVR2vyW0+lcvXq113vMXR6zAwahfbW/RQgA7kI+VSdiECKEOorjOI7jWr5sPXeUZdmsT0HAnWVsqz2bq7XmKuBTdbjjKEIIdRgGoU0l9gvHOoOwTa4CPolBiBBCHYdBaFOpOslVcOIVhC1c+XyqTjKuPAgh21qzZs3mzZvT+ITLli1rvR9paxs3bty5c2caX8sOMAhtSmiQnHncie/3X/oYoXHlQQjZ1t///vf3338/Xc9WVVV1/fXX+/3+Y93hiiuuaD2OmAVwsoxNpRpEfx93++/PBx1KUsUdRxHqhp577rk0PtvDDz88Z84ch6PtCQpDhgzhef7TTz89++yz0/ii1sIPTZvqaIsQKHDmc6l6bBQi1O3Mmzdv4cKF//jHP+66665LLrkkNzd3zJgxu3fv/s1vflNQUNC/f/+vvvoKAO68887y8vLi4uJx48atWrUKAPbt23f22WcXFRWddtppDzzwwD333KMoyiuvvHLhhRemUqkxY8aIovj555/feOONAHDPPfcsXrwYAC666CJ905msgS1CWyIgNMvOUEeCEMBdwKfqRG/PbN4bFyG7UcMNmpAw8xUdBT2BPuz4iZqaGqfTSVHUk08++a9//evVV1+97rrrTj311Pvuu+/gwYNPP/30rbfeunr16vPOO2/evHk8z7///vuXXHLJrl27Zs+ePWLEiI8//riqqmrChAmTJ0/esGGD0+ksLS1NJBJr167VNK25uXn79u0AsHv3bn2V/dixYx944AEzf2SjYRDakRiWWRfDcB1rr7sKcL4MQmaLf/mO8MMaM18x/xcP0t5gm9+aMmXK+PHjAWDGjBmffvrp7NmzAeCCCy749a9/DQBjx4795JNPqqqqkslkU1PT1q1bP/vssyVLltA0XVFRcdlll9XX1+/du7eoqOj4BSgpKamtrZVl+VjdpxkHg9COUg2Sq0P9ogAA4MrnGzdGjSgPQuhYAtOvDUy/1upSHKJvmQ0APM+37D7qdDpTqZQgCKeccsrgwYPHjh2bk5PjcDiqq6s5jmuZFFNQUFBfX88wjKa1fbxiy8FMiqLQNE3T2TOylj0/STYRGjs4QAgA+o6jOHEUIdSWdevWSZL08ssv33zzzZdeemksFsvPz6dpWu/2BIDvvvsOAMrLy/ft20cIcbvdXq/3wIEDLc/QcoTvvn37ysrKGOao04EzFgahHQkNUkcHCAHAnc8JDRJk1axmhFB6FBUV7du3b+nSpTt37pw7dy7LsgzD3HLLLVddddV77733wAMPfPHFFwAwePBgjuO2bdtGUdSsWbNuvvnmLVu2xGKxJ598cv369ZMmTQKAVatWnXnmmVb/QOmEXaN2JDSIeSMCHX0UzdGsmxHDMp+TJR33CKH2GDt2bGFhISGkpZ+zsLBw8uTJ+v85jps1a1afPn2eeuqpe+65h2XZX/ziF7m5ucFg8P7771+4cOF77703dOjQa665JplMUhR1zTXXvPbaa7///e+ffvrpp59+evHixVVVVT/88MNXX32lDx++9tprTzzxhGU/rRFI10yYMGHp0qXtv38sFtM0rYsvmvXWPbQ9tjfZiQdu+Muu8I54hx4iiqIoip14LWQQTdNisZjVpUCHiUajhJCHHnroV7/6ldVlSafGxkZFUQghzc3NgwcPfuuttwghdXV1/fv3TyQS+n3eeOONKVOmtDxk2bJl55xzjsnlXLdu3YgRI1rfor8j6YJdo3bUuTFCAHDmOoQmnDiKEGqXZcuW9e7de8yYMf37958yZcr5558PAPn5+R988IGqqm0+pHfv3vrhvdkEu0ZtR4oqtINmnZ0ZiHbmcmKTnPYiIYSy0owZM84555xIJBIKhVpPfqmoqGj5/7nnnnvqqae2fNmzZ09Ti2gKDELbERrETqyd0PG5jsgOU9f2IoQyGs/z+lzQY/F4PB6Px7TyWAK7Rm0n1dHN1Vpx5nJiM7YIEUKoAzAIbafDu4y24gxxQiOOESKEUAdgENpO57aV0XEBVo4rRMW1hAgh1F4YhLYjNEnOUAeO5G2Noiku4MDeUYQQaj8MQtsRm2RnbudXxOMKCoS6odWrV5911lkDBgz4+c9/Ho3insMdg7NG7UWTNVVQHd7Ovy98LifgCgqEzCKrskLaXnJnECfLU0C1vmXv3r0zZsxYtGjRyJEj77jjjtmzZy9ZssTMImU6DEJ7EZtlPsdx+B95xzhzHCK2CBEyy8NrnlpavdLMV3xl+nMhV07rW1599dVzzz333HPPBYDHHnusoKAgGo22bLeGTgiD0F7EsMzndHKmjI4Pcc2b4+kqD0Lo+O465da7TrnV2jLs2bOnd+/e+v9DoZDf76+pqcEgbD8cI7QXoamrW2Y7czlsESLUrSiKsm/fPv3/0Wg0Go0WFxdbW6TMgkFoL2JY5oNdDEKcLINQt/Pmm2+uWbNGFMXf/va355xzTjDY9hH2qE0YhPYiNkt8F6aMAgDndygpTVNwKSFC3cjVV189b968ioqKmpqaRYsWWV2cDINjhPYiNst8sEtjhEABH2TFJslV0MnFiAihjFNYWPjYY49ZXYpMhS1CexGbu7SIUIcrKBBCqP2wRWgnBKSIzAW6GoQ4XwahbuXWW291uVxWlyKDYRDaiBSVWTdDs11YRQgAOF8GoW5m4MCBVhchs9EAcO+99/bs2bNnz57z5s0jhABALBa74oor8vPzBwwY8NZbb1ldyO5CbO7qIkIdn4OHMSGEUHuxb7zxxqJFi5YtW0bT9KRJk/r373/ppZfefffdkUhk586dX3/99QUXXHDSSSf16tXL6qJmP6G5q4sIdXyQFcMYhAihE2toaHC73W63u3MPr6mpKSoqan26fSain3/++Ztuuqm8vLysrOymm256/vnnJUn6xz/+cc899/j9/kmTJk2ePPnFF1+0upzdgtgsO9MRhFzQgUGIULeyZMmS2tra49zhtddeq6urO/r2q6++uivdfv379z/+6x7L6tWr16xZ0+nXTS96y5Ytw4cP178YPnz41q1ba2tro9HoETdaV8JuRAxL6WkRBhxyTCEaLiVEqLu47777tm3bdpw7/O53v9u1a5dp5TmhN99885133rG6FIewTU1NPp9P/8Lv9zc0NDQ2NjqdTofj0CdyIBCor68/1uPD4fBZZ51FUYfN77jpppvuu+++Nu+fSCQIIUfcH+kSdSlXmSMeT8NOoYyLDh+MOnwn7q+QJAkAOC4NY5MoLQghyWTS6lKgwyQSCYqiRFG0uiBte/nll/ft2/fII4+8+uqrs2bNmjx58scff/zSSy/Jsjx9+vQrrrhi0aJFdXV18+fPLygouPzyyydOnHj0kyxZsuSNN96gafrSSy+dMWMGAEiS9Je//GXlypUAMGvWrJkzZy5YsGD16tWqqk6YMOGmm25qiYkWoijefPPN11xzzYIFCyiKuv766ydMmAAAu3fvfuihh6qrq0ePHv3rX/9669atn3zyCUVRjY2NgwcPvvnmm0/4M2qa1vqzUX9H2vPLcbvdNH2ChYJsbm5uy+FV0Wg0FAqFQiFBECRJ0j8cI5FIXl7esR4fDAY//fRT/UdtQdP0cV7Y4/FgELZJjR8IFPs8XmfXn8qZw7ES5/WeeEY1BqHd6FeKXq/X6oKgHxFCvF4vz9t0k4rRo0cHg8Ezzjhj+PDh/fr1e//99//3f//32Wef9Xq9N954Y2Nj4+TJk71e71lnnTVw4MCKioqjn2HRokX33XffM888o6rq3Llzk8nkZZddNnPmTIZh7rjjDlVVq6qq9M+K22+/HQD+8Ic/NDQ03H///Uc8j6Iof/3rX3fu3Hnvvffu2bPnggsuWLp0ac+ePU8++eTbbrvtqquuevzxx6dPn/7qq6/269ePYZhZs2aFQqH2/Iw0TbeuFPo70vlf2eHYysrKTZs2TZ48GQA2btxYWVlZXFzs8Xh++OGHESNGAMCmTZsGDx58nKdgGIZlcRlGGohd3nG7BRd0iBHZB7i0CCFj7Xi9puE7Uw/CHfWbSs532EfuwIEDvV7v6NGjzzjjDACYPXv2b3/7W71V9+ijj95444233HKL2+0+6aSTxo0b1+ZzLliw4M9//vPUqVMB4L777luwYEH//v2/+uqrvXv36vNo9NbOHXfcUVtbe+DAgauuuur+++8/Ogh1f/zjH8eNGzd+/Pi1a9c+++yzAwcOHDZs2J133gkAzz//fElJyf79+3v37s2yrB49lmOvueaaefPmXXLJJRRFPfnkk3fffTfP85dffvmf/vSnl156acOGDR9//PH8+fOtLmf2UwWNEMK60jP5ig86JJwvg5DxymcU955WZOYrsu4TfEpUVVUNGTJE///QoUOrq6sVRenQQ3bt2rV9+/b+/fu3nk3a3Nw8bdo0TdPKy8sB4ODBg8d6tpZ1jYMGDXrjjTd4nm9pTbnd7j59+thqtBIA2Msvv/z777/XfwVz5sy54oorAGD+/Plz5swpKCjIzc199tln9R8bGUpoltKyiFDHBxy4lBAhEzA8DTboMW093pSXl9cysaOuri4YDLIse/wBqSMekpeXV1BQcETULVq0qLy8fPHixQCwfPnyTz755FjP1tDQEAgEAKC+vj4vLy8vL++HH35o+W59fX1+fj4A6MvW7YCmKOqBBx5oaGhoaGiYP3++PraXk5Pz1ltvRaPR3bt3X3755VYXslsQ07SIUMfnOMQIBiFC3UVhYeGmTZv0aDn//POfeOKJZDKpKMr8+fMvuOAC/Q4bN2481sPPP//8Rx55RJIkQRAeffTRCy64YNy4cYSQp59+Wr/Djh07NE0Lh8OapgmC8OCDDx6nMI888gghpLGxceHChTNmzJg+ffr777+vZ+GLL75IUdTo0aOLioq2bdsmy7b4mMJNt+1CinT1JMLWcCkhQt3K7373u8WLF/ft23fhwoW33XbbkCFD+vbtW1paSgjRQ2vevHl/+9vfKioq9CZdi7y8PI/HM2/evEAg0Lt37z59+vTq1eu3v/2ty+V67733Fi9eXFJSUlJS8sQTT8yZMycajfbo0WPw4METJkzo3bs3APTo0ePoOSI9e/YsLy8fOHDg+eeff+GFFw4ZMuTJJ5/8yU9+0qNHj8cee+zNN990uVxXXHGFoiiDBg265pprzPolHRPVxcbpxIkT77333jYn47YpHo/jrNE27fm4jqKh9OyCtDybGJa/e2LXSb/vf8J74qxRuyGEJBIJnDVqK7FYzOfzPfzwwwcOHHj44YetLk67EEJUVe3QTEZVVSmKOmLOvyRJDoej5UNbFMXjzJ7V/3QFQXA4HJqmHfHqsiwfveKiPdavXz979uz169e33KK/I514qjZhi9AupLDMd/nciRacn5XjuKYeoe6LoqiOzudnGObolW8cx7VuurRzDQlN00e/eudS0AQYhHYhpuMAphYUTbEeVoqeYKoYQgilEc/zr7/+um0D71gwCO1CCqdzjBBwBQVCyHQsy86aNeuEO7nYTYYVN4tJEYULpHNfAh7nyyCEUDtgENqCKmpES9tqeh0exoQQQu2BQWgLUkTm0tovCriCAiGE2geD0BbEsMyntV8UAPggh2OECCF0QhiEtiBFlDROGdXxQRY3l0EIoRPCILQFMa3byuj4IG43ihBCJ4bHJ9mCFJbdxWk4hrA1h49VkipRCcXgPj4IpUHfvn0fffTRzz//vM3varIGNEWnu7qpgsY4u3uLJZVKuVwGHiqHQWgLUlQJ9k/ze0HRlMPLSlEljXt5I9SdnX/++X369FFVtc3v7npnf2ioP1DhSe+Lbvrr7srLeh5xAGE31PpAqLTr7r9cmxDD6Tx6ogUXcEgRQ54Zoe5p2LBhx/oW9R9/33El3p5pbrjQA/x9yop9ZQbGAOruLW6bSO/+ai34ACviLmsImcKwWuwQw1iLjYVBaD2iEjWlct70t871FmHanxYhdASiEiWpOjzp3BNDh7XYBBiE1hMjssPHggEzWrgAi1UIIRNIUYXzsRSd/mqMQWgCDELrSREljQcwtcYFHGIEO1UQMpxkTL8o4IJgU2AQWs+I/dV0PF5LImQKMaKkfXMoHRdwSHg5azAMQuuJkfTvr6bjAixWIYRMIEWNaxHipsGGwyC0nhRO//5qOhxdQMgcUkRO7zFqLQ7VYmLEc6NDMAitZ8T+ajqGo4EGRWh7/S9CKF1EA7YL1tEsRfO0nMRabCAMQutJEZnzG7WzAY8DDAgZT4rInN+onSv4gANPkjEUBqH1pIhi0GQZwN5RhEwhRWQ+aNTlLBdwSFGsxQbCILQaASkqG7eRIBfAudcIGU6KKEa2CFlcB2UoDEKLKUmV5mnaYdQbgV2jCBlNEVSKphjeqFrM+bFFaCwMQouJUZk37EIS9CqELUKEjCRFFIOmjOo4P66DMhYGocUMr0K4lBAhgxm3rYwOxwiNhkFoMUMnm4G+dT22CBEykmjYIkIdtgiNhkFoMSmqGLd2AgC4IHaNImQs47YL1mGL0GgYhBYzbmcmncPLyAmVaLgvBUJGMbpr1OFllCTWYgNhEFrM6BYhRVMOLyvHsF8FIaNIEcNrMevBWmwgDEKLGXSqdWscLkJCyEhiRDa0axRw1pvBMAgtZvS1JOBhTAgZzLiT1Frwfpz1ZiAMQksRkOOGByGeU4+QcYhG5ITq8DKGvgrnZ6UotgiNgkFoJSmusG6GoilDX4XHc+oRMowcUxwew2sxThw1FAahlSTjhxYA92dCyEiikZvmt8AWoaEwCK1kwgAh4DA7QkYy6XIWR/qNhEFoJaMXEeqwRYiQcYxeAaXDFqGhMAithC1ChDKdFDV2l0Qdh3slGgmD0ErmtAhZF6OpRJU0o18IoW7InMtZh5vRRE1TcHMZQ2AQWkk0pQoBAOdjZexXQcgA5lzOAoVjHAbCILSSSVUI+1UQMow5LULAMyiMhEFoJSmi8KZVIWwRImQA0bzLWRZbhAbBILQMUYkqqA6vGUHIYxVCyACaQjRRc7iN3VZGx/kd2CI0CAahZaSo4vCxYOx+FIc4/A5sESKUdoemjJpSi7kA9usYBYPQMuYsP9LxuIICIQNIUcVhVi3m/Lim3igYhJYx+jDP1nC+GUJGkCIyHzArCLFFaBgMQsuYsw5Xh2vqETKCFFHMq8V4OWsYDELLmNk1ip0qCBlBiiqcmS1CvJw1BgahZaSIeVWI4WmgQRVwcxmE0snMfh3WyRCNqCLW4vTDILSMGJU5n0lVCLBfBSEDmHk5CwAOXBBsDAxCy5hchfgAi8fzIpReZrYIAYDH43mNgUFoGSkq8yZWIWwRIpR2YkQ283IWd1kzCAahNTSFaJLGmrIhhY4LsDhfBqE0UiWNaMA6TazFeDlrDAxCa0gR8zak0HG4uQxCaWXmIkIdbhpsEAxCa5g561qHnSoIpZeZiwh1uO+2QTAIrWHyGDsAcDjMjlBaSVFTBwgB+3UMg0FoDZOnjAK2CBFKN2tahFiLDYBBaA0zt5XRHWoREjNfE6FsZkEtxskyxsAgtIb5XaM0S9EcrSRVM18UoSxmfi1meBoAt4hKPwxCa0gRs68lAYAPOES8nEQoTUTTBzgAB/uNQb/zzjvU4datW3f//fe3vmXv3r1WlzPbSFHzzmBqgQMMCKWRFJF582uxH7eISj/6/PPPJ//18ssvl5eXjxw5EgCuu+66ltt79epldTmzjSUtQjyDAqE0kmLmncrbAluERjisa3ThwoXXXnstRR1a5q0oeN1hCE3SNJWwLvM2pNDhwZ4IpYuSUmmGYjizR5d4nP5tgB/fxaqqqi+//PLKK6/Uv3zhhRf8fn9RUdF9991HyDHnGmqaVl1dveVw9fX1hhc8k5m/ml7HBbBFiFB6SGGzFxHqHH5cU59+P76RixYtmjp1qt4Leumll95www3BYHDdunXnnXdeSUnJnDlz2nx8PB6/++67XS5X6xuvvPLKW2+9tc37JxIJTdNaGp3dU3y/wHrpWCxm8uuqDjnRmDridSVJAgCO40wuDDoWQkgymTzO1ScyXzweP/KWgwLtocyvxRqnJJsE81/Xbo5+R47F7XYzzAm63w4FoaZpL7744mOPPaZ/2adPH/0/o0aNmjt37ocffnisIPT7/QsWLJg4cWI7y0RRlMfj6eZBKCiaK4f3+Xwmvy5VyDYk40e8Lgah3RBCaJr2er1WFwQd5oiKk5QUT67L/FpMCunmtSnzX9eG0vhLONQ1+sknnyQSiWnTph19j0gk4vF40vV6CA4tP7KiazSIXaMIpQcOcGSTQ2/kwoULr776ap7n9S8ffPDBUaNG5eXlrVix4rnnnnv33XetK2EWkqKKw9x1uDqHl5ETKtEIRXfrFjlCXSdFZHeh0/zXxQMojMACgJxjOu8AACAASURBVKIoPM///Oc/b7lVVdX7778/Ho/36tXr3XffnTx5snUlzEJSRPGWWFCFKJpyeBg5ppi/hBGhLCNFlWA/C1qEtIOmWUpJqeZPO89iLACwLPvSSy+1vvWuu+666667LCpS9rNkNb2OCzikCAYhQl1l/v5qLfQN9DEI0wi3WLOAFLFm4jXgLmsIpYkYtmaMEHBNvQEwCC0gRc0+vaUFHsaEUBoQkOMWbA6l4wKsiPNl0gqD0GyKoAJ1aBd58+GUM4S6ToopDjdj1aQzzu/Ay9n0wiA0mxRRzN+otwVOOUOo66SIZcP8AMAHcHOZNMMgNJsUsWyMHXB0AaF0sGoRoU6f8mbVq2clDEKzSVacYdaCC7BiGKsQQl1ibYsQj5FJOwxCs1m4dgIAeD+2CBHqKjFq2UwZODRZBi9n0wmD0GyiFScRtmA9jCYTTdasKgBCWcCSI3lbcH5WjitEw23Z0waD0GzWViHA+TIIdZkUsWa7YB1FUw43I8dVqwqQfTAIzSZFLVtNr+P8LA4wINQVUtTi7ZlwHVR6YRCazfIdzriAA1uECHWFtZNlAIALYL9OOmEQmkvfkMJnZYuQx20pEOoCTSGqqDncVm71iS3C9MIgNJUcVxgnQzFWnoKE21Ig1BVSVHb4WLD0KDPej5ez6YRBaCoxIvNBi09+wGtJhLpCCtujFmPXaPpgEJpKsnTthI4PsmIYgxChThKtHuYHvJxNNwxCU1m7ml6H+zMh1BW2aBHiMTJphUFoKmv3V9Md2m4UF+Mi1CmideeJtsCDRdMLg9BUoqU7butolmKcjJzAy0mEOsPaA2R0rJvRRA23iEoXDEJTSVGFt/paEvQVFDhMiFCniBGZs7prFCjg/DhfJm0wCE1l+TpcHQ4TItRpUtjiXRJ1uKY+jTAITWWXIAw6cBESQp2h74lh9dxvwImjaYVBaB47bEih4wMOCbtGEeo4KaawLov3xNDhFlFphEFoHikqc36HtRtS6LggViGEOkMK22CAEAAAxwjTCYPQPHZYTa/jcYwQoU4Ro1YewNQaF8BjZNIGg9A8on2uJQMOnDWKUCdIYZnPsU8txsvZ9MAgNI8dNqTQ8UEcZkeoM0QbLCLU8UEc6U8bDELz2GFDCh3D00CBIuAJ1wh1jBS2xcRvwC2i0gqD0DySDY6eaIHDhAh1gmiPFVDQskVUHGtxGmAQmkcM26VTBfTLSexXQaiDpIjMB23RrwMAPC4IThMMQvNIdtiZ6b/4IM6XQajDpIhi+XbBLTjcKzFNMAhNQjQixeyyfAJw7jVCHaekVIqhGN4uH5t80CHhxNF0sMs7mvWkqOLwshRtg+X0AKAf44JjhAh1hGibmTI6LoBdo+mBQWgSWw0tAACHKygQ6iApotiqFuM6qHTBIDSJrWbKAF5LItRxdmwR4hhhOmAQmsRWM2UAgA+wOGsUoQ6x1QooAOCDONKfHhiEJpFscDZ9aw4Pq+IJ1wh1hBSxy0ajOg5H+tMEg9Akom32VzuEAi6IexUi1AFCs8zncFaX4kcMR9MspSRwi6iuwiA0iRi212QZAHDmOMSwZHUpEMoYYrNddtxugWvq0wKD0CRSRLHVGCHo59Q3YxVCqL3sNkYIOF8mTTAITUH+eyqvnfA5HAYhQu2kJFSKttFqeh2uoEgLe72p2UqKK4yLoVm7rKbX4S5rCLWfYJuTCFvjAizOl+k6DEIzSGHZVosIdRiECLWffY7kbY3H3fPTAYPQDLabMgoAAHwOjhEi1F5is2THIMTL2XTAIDSDFJE5m00ZBaxCCHWEPS9nca/EtMAgNIMYsdf+ajqGx0VICLWX2GzHIMTL2bTAIDSDZLMtClvwQdyZAqF2EW22ml7H8DRQoAh4OdslGIRmEO23/EjH5+BIO0LtYt9aHHBIONjfNRiEZhCb7DjfDAD4HIeMLUKEToRoRI4pXMB2I/2gz3rDy9muwSA0HgEpatdrSRxgQKgdlJhmq4O1W8OdMboOg9BwYkRm3QzF2LIKBTkJ991G6ESkiGLPTh0A4HMcAgZh12AQGk4My3yu7cbYdXyOA4MQoROS7R2E2CLsIgxCw4lNNu0XBewaRah9pLBi21rszHGIzXiMTJdgEBpODEvOXJtWIS7AKgmVqMTqgiBka3JUtW0Q4hhh12EQGs6e63B1FE2xHkaO4SIkhI5HCis2XESo4wKsHFfwcrYrMAgNZ891uC1w7jVCJ2TnFiFFUw4fK0VxsL/zMAgNZ8NTrVvjgiyuqUfo+KSwfSfLAAAfdAhNOEzYeRiEhrN7EAZYnDiK0HGookZUwroZqwtyTM5cHCbsEgxCYykplRDCuuxbhbigQ8QgROjYxCbJhqfHtMbh9O+uwSA0ls0HCAHAmesQsVMFoWMTmmQux9ZB6MzFpYRdgkFoLDEs23bthI7PdYhN2CJE6JiEJslh7xYhH8SlhF2CQWgsO6+m13H6ARQ49RqhYxAaJT7X3kGIY4RdwwqCUFtb2/J1QUGB1+sFgEgk8vXXX4dCoZEjR1pXvIwnhiU7z5QBAJqlGBctRW16YiJClhObZF8Jb3UpjocPYtdol7Br164988wze/XqpX/96KOPzpgx45tvvpk2bdqYMWO2bds2duzYxYsXW1vKzCU2y54Sl9WlOAFnLic0YRAi1DahSQoFPVaX4ngYnqZYSk6oDo995+XZGQsAffv2/eGHH1rfevfdd99666133XVXJBIZNGjQl19+OWHCBItKmNkEe6+d0HG5DqFR8pe7rS4IQnYkNEk2nywD/91x1OGx+2W3PdEAoKrq999/v3v3bk3TACAajX722WdXXHEFAAQCgenTp7/11lsWFzNj2XwRoY7PYXE1LkJtUhIqRVGM0+7TKXDH0a5gAaCurm727Nn79u0rLi5esmSJKIo0Tffo0UO/R2lp6dq1a4/1eFmW//Of/9TV1bW+sV+/fkOHDm3z/pqmaZpGUXY8nC/tiEqUhMJ6af0Kw540TeOCbLJasnMhuxVCiF5NrC4IAgBINgjOEGf/d4QLskJTN6rF7X9HaPrEFzHsqFGj6urqHA6HLMtz5sy58cYb//znP7Ms2/JgnudTqdSxHi9J0ueff75hw4bWN06bNq2ysrLN+wuCwDBMNwlCsUl2+FhREq0uyPFIkkR5SapRFATB6rIgAABCiCAILGv3vrhuInYw4QjQoig6HLbu2mF8VLIh1X1qcfvfEafTecIsZF2uQ33KDodj7ty5P/nJT4qKikRRjMViPp8PAOrr64uLi4/1eI/Hc++9906cOLGdpdc0ze12d5MglGrizhDndtt67I1lWabYURsO27yc3QchhBCCb4dNNMWT7gKXy+Wy+TviLVTq93WjWqyqahp/2MNycvv27QUFBYWFhX369Fm6dKl+49KlS0855ZR0vV63IjTKzpCtt5XRcQFWjuExLgi1QWiSnLkZUIudIU5oxJH+TmIfeuihVCpVUVGxY8eOBQsWPPzwwxRF3X777TfffHM0Gv36669ramouvfRSq8uZkYRGKSOCkKIpzu8Qw5kR2wiZSWySQoN9VpfixFwYhF3ATpw4ccmSJZ988klhYeEHH3wwfvx4ALj++utDodBHH31UWFi4cuVKj8fWa2hsS2iUQkP8VpeiXZwhR6bENkJm0vt1FLB7xjBOmmIpOa44vDi63GHs2LFjx44de/Q3Lr744osvvtj8AmWTDIoWPpcTmnDuNUKHIyCGZT7oUES7ByH8t3cUg7AT7L44JqMJjZIzLzOCEAcYEDqaFJVZF01zmfE56czFWtxJmfEGZyJFUIlKMmXHIzyMCaGjCU1yRsyU0TlDnNCI/TqdgUFoFKEhY5qDcGi7UQxChA4jNEp8hoxugB6EWIs7BYPQKBk0QAgAPF5LInQUoUmy+XmireEAR6dhEBpFaMykThXOy2qSpgrdZX8mhNpDaMiky1lXHic0YBB2BgahUTJopgwAAAXOPC7VYOvd4BAyWapecuXb+iTC1rgAK8cVTcGdMToMg9AoQmNmbEjRwpXPp+rwchKhH6XqRVd+xtRiiqa4IM566wwMQqNk1hghALjysUWI0I+UhAoaZNayPBwm7BwMQkMQjUgROYOG2QHAlc8L9ViFEDokVS86CzLpWhZwBUVnYRAaQmyWHT6WYjLpkA1nPpeqxxYhQodk1gChDldQdA4GoSGERsmVQTNlAADHCBE6XKpezLhajF2jnYNBaIiMGyAEAIeHARrkuGJ1QRCyhVRDZrYIMQg7DoPQEEKjzGfUlFGdK59L4TAhQgCgtwgzbYwQD2PqHAxCQ2RiixD03lEcJkQIAEiG7ZKoY5w0zVJyDPt1OgaD0BCpOtFVkGGdKoA7UyD0X1JUZjiadWbGpvmtOfFytuMwCA1AIJWBk2UAqxBC/5Wql5yZNkCocxfwSZz11kEYhOknNEsON8Pwmfe7xTFChHSZtadMa64CXAfVYZn3YW1/qTopE/tFQR8jbJAAtypE3V6qXsrUIMznUwcxCDsGgzD9MnSAEAAYnmadtBjBnSlQd5eqFzNu7YTOVYADHB2GQZh+qTrRnWmzrlu48nnsHUUoVZ+Rw/wA4MrjxGYZz6DoEAzC9Mvca0kAcOVzAl5Oou6NaERsyry1EzqKofgcXE3YMRiE6Zc8mKldowDgKuSTOMCAujehQeICDtqRqR+PrgIuVYe1uAMy9Z22LVXUFEHjg5l07kRrnmJn8oBgdSkQslLigOgpdlpdis5zFfAYhB2SSUdtZYRUnejK5yCTjp04jLvImajFIOwMWYMGARpF0ihAk0iaRAhLEJZIVIK4DDEZojJJKqD/E1UQVEipBAASMkjaEU/GAxw2ZSnIAUUBDRDgKADI4YGhwM+Bh6U8LHgdkMNDgKOCHAQ5yOGpXB7ynJDvpHyZeklmpeR+wV2cqZ06AODK52NVSatLkUkwCNMsVSe6M7ZfFAA4PwsAckxx+PBv40gqgf1JsicO+xKkNgl74+RgCvYlSJ0AB1MkJkMeDyEnFeIhl6dyeQjyEOSoPj7wOsDPgd9Bu1jwsOBiwckAT4ObpQDAzQLfagMTQkgikfB6va1fOiwBIaASiMqEEAhLoGgQkyGukIQMcQXCIoQlUh2HsARNgtYsQV0KGgSiEsh3UiUeKHRRxS4odlO9vFDipkq9UOalPPgmtyV5QAgNC1hdis5zF3AH1zRbXYpMgvUgzZL1mbqIsIW7yJk4IAR93hPfNXvFZdgWIdujZGcUdkZJVYzsjkNtkuQ5qTIv9PRQJW4o9VJj8qGHm853QYGTyjOyLy3433kbec4jehtO0PmQUqBOIPuTUJcitUnYnyTLD0BtUquOw544cTHQ20eV+6g+PqjwU339VGUAenkytkMjTRL7xdKp2DXajWAQplnqoBga6re6FF3iLnYm94vBym4UhGEJNjaRTWGyqZlsCZMtYWgQSaWf6hegKvwwroC6vC/d2wulXirj5k+4WCjzUmVeaDMy6wWoipHdMbIzBmvqySs7ta1hEpNhQJAaGKQG5VCDgzAklyr3Ud0nGzWFiGE5Q1fT6xxeFgDkuKL/B50Q/prSLBOPbjmCp5iP701ZXQpj1STINw1kfSNZ3wjfNZFGgQzOoYbkUINyqGm96AFBKPV2i4/+fCfkO6mx+Yf9rFEZtobJD2HyQzN5bgvZ1AxNIhmaS40IUSND1Og8akhO5l0QtF/qoOgKcRST2e+/q4BP1UsYhO2Ev6a0IiBk4GGeR3AXOQ+uCVtdijRLKPB1PfmqjqyuI1/XE4WQ0XnUqBB1VV9qRIju4+8WsddOfgeclE+d1CodwxJsaCTfNpJlB8jjG7XdcTI0hxpbQJ2cT40vpMq8WfXLS2T4TBmdO59PHhT95W6rC5IZMAjTSWiW2Mzcbrs1dzGfPCAAOeHwk901ibDsgLZ0P1lxkPzQTIaHqHEF1M/6Uo+fQmfZZ7fRghxMLKYmFh/6pcVlWN9IVteTJVXkV6tVmqJOK6QmFFFnlFADgxl/SZE8ILqLMniAUOcqxKWEHYBBmE7J/aI7k5cf6Vgnw7oYoVly5mZeH29KgS8OkM9rtM9ryc4oGV9ITSymHz2ZHpNPZeDRcjbldcDpRdTpRRQMBQBmZ5SsOEi+2E8e3ajFZHJWCX1WMXV2z0xtKSb2C8Xjc60uRVe5i5z7lzdaXYqMgUGYTolawVOS8UEI/50vk0FBuDVCPqgmH+/TVtWRUXnUpBL6yfH02HyKzezGeWao8FMVfuqqSgCA6jj5vJb8u5b8bq0a5KipPalppfTEIorPnKuQ5AHBXZTxXaOeEmeiBhcEtxcGYTolaoVMnzKq8xTzif1C7mCf1QU5HpXAyoPkvT3aO3uIoMK0XtQvBtFvTaa9uITcOqVe6pp+1DX9gACzoZF8vI/cu079volM7kHPKKOmldIhe0eMKmhKUs2gS8Bj4YMOTSM4cbSd8HeUTolaoXRqgdWlSAN3kbN5S9zqUrRNJbB0P3ljl/b2Hq2Hm5pRRr8xiR4RysheuCxGAYwIUSNC1J3D6UYRPqjW3qsmN62UxxZQF/amZ/amC1xWF7EtCb05mBV/TZ5iZ6JWCPbrRuugOg2DMG00WRMjmb38qIW72FnznwarS3EYArDyIHl1p7akSivzUheV019NZ/v4suITK9uFeLiqkr6qElIK89E+7c0qctfX8kn51GUV9IXltN9OLfjkfiELZsroMAjbD4MwbZL7RXd+xi8/0rkLeaFR0hRCs9b/OLti5MVt2ss7iIuFn/WlV2L+ZSwXCzN70zN7Q0phPtirvbKT/HKVPLUnfVUlPbWnLapOolbI6O22W/OUOCM7E1aXIjNgEKZNvFbwlNiyu6fjaJZyFfCJGsFXZtlPJKjwZpX2/FZtUzO5vC+9ZDIzEvs/s4WLhYvK6YvKoVlkXtul3bte/fky+J9+1Jz+tLVXObHqVMHooIUFSCNPibN2GU4cbRcMwrRJ1AruEnvPBOgIb6krVp20JAh3RMkzm7WXtmuj8qhfDKJ/WkZzOPkzS+XwcN1A+rqB9A9h8vet2rh3lZEh6rqB9PQy2vwGoqaQ1EHR0yNLWoTuIj5VLxGV2KKtbW/4AZM2iSxqEQKAr5fL5I3WCMDH+8i5Hyunvq9wNKyewX58DnthOaZgtzAoSD1yMrP3MsfV/ehHvtf6vKb8eYPWaO6K8ESN4CrgM/c83iPQDprPcSRxWX07YIswTQgk9wuezN+ZqYWv1LXPrPkyggr/2K4t2Kg5GbhlCP32FBoXv3dPPAOXV9CXV9DfNpInNmmVr8sXl9O/HEr3C5jRpolVJ72l2XMtC/pqwiwa9TSO2UFYL1CL92oURXkd4KDBzVIeFnwOyOUh5KRsNX+sQ8SwTDuobFqy4y5ySlFZSamsy8BQahbh6c3ak5vUk/Lpp09lzijGPhwEADAiRC2cwPw5xTyzWT39n8pphfRvhtNHbA6edvG9qUCFx9CXMJmnxJmsFWC01eXoLFGFRpE0iRCTISFDs0QA4JyedNqPmzb7gzuhwrpGACBxGWQNkooWVyAmQ5MIjQIRVCh0UT09UOymSr3Q20v18UFlgKrwUzbvH0vUCp4eWXUtCRR4e7ji+1IGncdUl4IFG9W/btGml9GfT2MHBTEC0ZEKXfCHUcwdw5iF27RL/q329cP/jWAmGna1FKtO9Tgzz6Ant4SnJDM2WtufhG0RsjNGdsfInjhUx8nBFNQkSEqFPCfk8pTPAV4HBDmKAhhfQHyONP8NmB2EvT3k2VNp6hgb84oq1Atkb+LQOeBVMfJpDdkRhb0J0stDDcmhhubCiBA1Oo+y29mhWdn/4C11xavTH4T1Ajz0nbpwq3ZpBb32AjZDd6REpnGz8ItB9NwB9OId2v8uV4vdMG9U+uNQEVQpIrsLs2d0A+y60VpSgfWN5NtG8m0j2dhMNocJz0Cln6oMUL291Fkl0MtDF7uh2E0FzVqVba+uPJ6Bnh6qpweOOPhA1mB7lGxqJt82kr9t0a5vIBQFY/Pp8YXUaYXUmHzr24uJmlRoWMDiQqSbr5erfn0kjU8YkeCh79RnN2uXVdAbZrI9bHY1g+zMQcP/9KOvrKRf2an9fLla6oE/jWFOLkjbn1B8r+Dp4aLorPqb1Ddak6IK57f4o353jCw7SFYcIKvqyI4oGZxD6U2aqyvpQTlUrtWXH/YKwmNx0DAoSA0KUrPKD91SHSer68nyA+Smr7TtEXJKAXVmCT21JzUiZM0pMNE9qd7Ti6x4ZQN5S1273j2QlqcSVPjLJu3h79WfltLrZ7J2a9CjTMFQcGVf+rI+9IvbtVmfq2PyqfvH0APS0a8et2ixkNF8pe7YnqQleyDXC/BpjfZZDfnPfiKq5LRC+rQiak5/enjI+qbLETIjCI9W6qVKvYdyMSLBF/u1z2vJZf/WIhL5SS/6p2XU2T1ot1k/nNgsE5VkwUa9R3DmckTt6uWkRuCVndr/faONzqO+mMam5TMLdXMsDXP60z/rSz/1gzbxA2Vmb/qeUUxR11IsVp3KH5ltnToA4C9zmRyEG5rIe3vIP6u1bRFyZgk9uYS6Y1h6LlaMk6lB2FqAg+ll9PQyAICqGPlgL3nqB+3qpeqkHvRF5dR5pYZvZhjbk/T3zs6ToH2lrlh1MjSkk7Vo5UFy6yqVoWDxmcyphbauCYZSRU2OKXJSVZKqKqhKSlUETZM0TSaKoAIBVVSJ+uP9FUVh2eaWLxknTdEUw9EUS7EumuEYxkmzLoZ10aybZT2Mw8NkWZ9eezgZ+NVQenY/+k/fqkPflH81lLltCN3p857ie1PlM7KtUwcAfL3d1Z/UG/0qBGBNHVlSpS3ZTVgKZpRR88cypxZSmbImMxuCsLVyH/WLQdQvBtHNIrxXrb26U7thhTq1J/2zCurcXkatlI3uTvp6Z2GnCgD4ytzRqs4EYW2S/Hq1tuwAeWAsfVlFN/iQJiDFFKFREpslMSxLEUVslqW4IoVlOa5QDOXwsKyHYd0M62RYN8M4acZBs26az3VQFDA83bL9ByFEFEWn88e5V0pKA0JUPTiTmhiW1ZSmCKqSVJWkKidVOa44PKzDy3IBlvOzfA7HB1g+h3PmOvhczg4bxhonh4eHT2auH0j/eo02aIny6Dh6RlmH67nQJGVlpw4A+ErdiZqUcfvLbI2QxTu0xTsIx8CscurdKcyw3Mz7e8u2IGyRw8PVlfTVlXSTCG9UaQ99r/3vcvVnfek5/em0z9SP7Un1Pq8wvc9pE4FKz643azv0EFmDBRu1B79T5w6g/3Y6a1oHtZlUUUsdFBMHhVSdlKoTU/WS0CCyLoYPcc4cBx90OPO4QF8P52O5gMPhZTp0CUYISSQSXm9HJusSkBOKHFPEiCLFZLFZju9NNXwXFZtksVlyeFlnHufK510FnLvQ6S7i+WDGrtg9hgo/9dZk5vNactNK9a9btMdPYfr6O1DNI9sTAWOWCVmOcdLOXC5RK3h7pfNiPS7D/9ulLdqm7YqSy/vSb0xiRuVlXv61yMZPqcPl8jB3AD13AL09Ql7Yrk35UO3rh7kD6YvStHeXppDkfsGX1j8y+/CVuoRGWU6oDk+7upyWHyDXr1BLvfDVdLZDn0Q2J4bl+L5UYp+Q2C8kagQprrgLeHch7yrk80cFXfmcM49jLJwAQIHDyzq8rLv4yO8QjYhhWWiQUvVS8qDYvDmePCBoEnEX854Sp6eHy9vT6Sl2Zsd2lJNKqA0z2cc3aae8p9w0mPnNsPb2lIa3J4KVWbWUvjVfb3dsTzJdQfh9E3lms/baLm1CEX3ncPrcnjSbIf2fx5H9QdiiMkD9aQzzh1HwfrX2zGbt9tXq3AHMdQPpwq79eST2pVwFfLZuiEnRlL+PO7IjkTf8BL2jYQl+vVr9eB95bBx9UXnG/zY0SYvtTUWrkvHqVKw6STTw9nJ5ezjzRwV7n+d05XGZcnYrRVPOXM6ZywX7/XijklQTtUKiVojuStR+2Sg2Se5ip6/U5Stz+8vdfE4GtxcdNNw+lL6kD3XLV9rwt5TnTmvHikMCke3xsp9kw5HabfKVucLbEsWndelJNALvVWtPbNS2RmDuQPq77Fr+1I2CUOeg9RPR6E3N5C+btEFL5IvK6V91YTPD6J6sHSDUBSo94e3x4wfhkirtlq+0C3pTmy5iM3efPFXUorsSkZ2JyM5Ecr/oLnH6y9z5owN9LijO6Gw4GutmAn09gb6H2kCqpCX2pWJ7Ug0bIrve3U8zVKDC46/wBCo8GXrQdC8P9dZk5r092pVL1XN7UfPHMsdZmp08KNIcnZUDhDp/b/feTzs/XyalwIvbtUc3aiEebhlMX1ieNduS/6jbBWGLwTnUs6cxfxzDPL1Zm/BPZUIRffcIekTHT7yL7U6FhviMKKFNBCu9B1ZWH+u7B1Nwwwp1S5i8MYkZn4HzQolG4ntTzVvi4W3xRK3gK3X5Kzy9zyvylbmze45JawxH+/t4/H08PQAAQGiQ9KuBvZ/UAQXBft6c/t5gPy/rzrCt0KeX0WeU0HeuUYe+qTxzKnNeadtvaHh7PIv7RQHAlc8rKVWOKQ5fxz7wYzI8s1lbsFE9OZ9eNCGbJ3533yDU5Tnh9yPp24fSf92i/fQTdWSIumcUPbojo76xPclsnSmj8xQ7lZQmNstHt4pe3andtkqd3Y9+5Uy209PWLaGk1ObN8aYfYs1bY3zAkTPAWzq1wF/uzsJr3Y5z5nHOPK7w5BwASNWJ4W2JurXhHa/Xuov43MG+nEG+DNpN0O+Ap09lLu5Drl2mvr6LemJ8G03D8LZ41hzG2zYKfGXu6O4OrCaMyvD4Ru0vm9SpPenPfpL9WwF39yDUuVm4dQh9/UD6+a3a+Z+qo/Oo+0bTQ9sxCVhsljWVOENZ26kCAEBBoK8nvCNeeFJOy20NAly/Qt0cJh9MZTt03WAtOaY0fB9t/C4aUnkiLQAAIABJREFUq04GKjy5g/3l04ss333KzlwFvKuALz4tl6gksjPRtCm2+e/VQEFoqD9vmN9X5s6IsdIziqkNM9m7vlaHvan87XRmas8fC000Et2VrLykh4XFM4G/zBVrXxAmFXh8k7Zgo3pOT3pldk15Ow78CPgRz8CNg+g5/ennNmtnf6RM6UHPG02X+473dxDeFjfocAZbCfbzRrYlWoLwn9Vk7nL1Z32pl85gM+LgQCWhNmyI1H8bSdQIuYN8xeNzB80uzdb5TQahGCrYzxvs5+1zQXGiVmj8Prr99VpVUPOG+/NHBu1/jJ+HhSdOYc4vI3OWqef0pB4+mfGwAADxvQKf48imA9TaFOjnPeE6KFmD57dqf1yvnV5ELTuPNecMSJvI8re/E/SzYef0px/dqJ30jvI//ejfjmByjrEnbPOWeO7gbB4g1AUrPdUfHwQCcQV+tVr9rIa8dhZzWpHd64mmkKaN0bpvwtFdyZxB3h4T84L9vd1n5M84nhKnp8RZOrUgeVBs+DaydfFeAMgfHSwYE7T5lJOzSqgNM9mbV6qj3lb+MZE5uYAKb8vyAUKdr9QlNMvHGSZ8e7d259dauQ/+OZUZ2fGpEpkOg7BtXgf8fiQ9dwD9h3XqgCXy/41grht45PgR0Uh4R7zPzKOWbmUdZ4hzeNnV3yau3MGfXkh9O5NN+8GY6RWvEQ6uaqr/NuLt4So4Kdj/yl4Mj+2/9HMX8qVTC0qnFsSqU/VrwxsW7HIX8YUn5+QND9j2gsPvgBcmMm/t1s7/VLlhEDNtQ6RvN6jCFE0FKz3NW+MFY44cDV3XQG5bpUYk+Mt45uweNn3XjIZBeDyFLnjmVOYXg+hfrlKf2awtOOWwP5T43hQfcHAdnIiViVQCmwr9mz8O//nS4gttvEZQlbT6dZEDK5uUpFo4NjjyV32zbwsVe/KVunylrvLpRU2bYgdWN+96Z3/BqGDxqbmuAqvP1zmGmb3pcQXUrz5MjWpUc/PcFhzNYLpgf98RQViXgru/UT/cq907mpndj+4GGyEeU/Z/iHfd4BzqX+ey/6wmN65Qh+ZSj55M9/ZRANC8JZ4zIPsHCPclyJVL1TyX75fCnnFlJVYXp22pemn/8sa6deFAH0/vaYXBft6MmMSRZSiGCg3zh4b5xWb5wKqm75+uchXwJaeFcof4bLgneImbesAf/7rCd/L7yuPjmEsr7HuFlxY5A7x7PjwIBIACRYOnftD+9K36P/3oLbMcmbv2N10wCNvrvFJqSg/20Y3aSe8qtwxmfj2Mbt4aLzs3mxdOAMDbu7XrV6i3DGF+M8z1XbUjuivZsgrbJsLbE7VfNMSqU4Un52AT0Cb4HEfZuYWlZxc0fh+tWdpQ9d6B4tNzi07OZZz2CpvGbyNTL+7xkZe9/D/qv2rIX05hvNn758MHHQ4PE69Jfcs6b1ipFrvgy/PwWLRDWACoqqpaunRpJBIZOXLkxIkTAWDz5s0bN25sudO0adPc7uw8ZqhDeAbuGk7/rIK65Stt3BviX2rFQJ+s/bWkFPjlavXTGvLuFFY/BzxvRKD+24hNgpBopGFDtOY/DZqslUzMG/A/pbYdlOq2KIbKGxHIGxGI703VfNHw9WdbC8fm9JiYZ5P1KskDoipqvjL3KArWXsDevFId/Y7y6pmZvXn08Tn7ep//OPqox/FoVuyDmEbsxo0bJ06ceN555+Xl5T322GOTJ0/++9///vbbby9cuHDUqFH6nc466ywMwhalXurtKcy/Po9963Y9u0x75GSmi7uV2tDGZnLZv9VhudS6C37cMi1vhH/Dgl0VM4ut7eYiKjm4pnnfvxu4gKP0nILcgT7sBbU5by9X/yt6ic1yzRcN6x7cnjfM33NSvuWrbxu+jeQND+h/PB4W/j6BeW2Xdu6/lDuHM7cOsV9PbtcQgBe2aa8ddF0Xafrh6sIsbvh2DlteXl5dXe3xeADghhtu6Nu37x//+EcAmDJlyjPPPGN18eyrT11iwBn+vU4Y9pZ8/xhmdv/sqTnPbNbuWas+dDJzdeVh14zOXM6Z64jsSAT7WTMyqink4Kqmff9u8JQ4+13e01+OF2eZhM9x9Dm/uHRKQe3yxg2P78oZ4O01Od/C2TT130b6X96z9S2X9KFPzqd+tlT9rEZbNIEtyJYL3K0RMne5mlLg2Qt94uP7nYoKjkxY/2si2uPx6CkIAB6Ph6IoiqIAYO/evS+99NK///1vRVEsLaEdabLWuClaPNI/fyzz2bns81u1M/6pbAkTq8vVVU0izPxMfX6rtvyn7BEpqCscm7N/eZP5BSMq2b+iae3928LbEgNnlw66tgxTMEOxHqZ0asGYu/u5Cvjvnqza9so+oUEyvxjhbXGKpo4+mai3j/piGjsyRI16R/mkJuNrtKjCvHXa6e8rF5XTX01nRxYywQHehg1Rq8tlO4d11t95550XXXRRUVGR0+lUVfU///nPqlWrWJZdunRpbm5um48XBOH555//9NNPW994yimnTJ48uc37i6LIsqyetZmr6buYpydPOFUU1X4e+PfZ8Ndt1Onva9f1h18PIZm166YkSQBACPnyIDV7OXVRb/LiqYSjVVFs486BYa7dHx+M7Is5zTqUgGikcX2s5rNGdxFfcUWRp4cTAMQ2C5ct9BPqHY6s7r2ioOB0f2is9+Dy5m8f35kzyFMyKcQFzBs73PvvuqLTg6LU9h/S74bChAJqzhfkot7k3pGEo0EURY6z9V4BR1tZR92wiurnh9XTtBK3JksAADnDPQe+bMoZmfHXke1/RziOO2Hi/HjVf//993/99dd6d+gvf/nLjz76aOHChd9//30oFJo/f/7xnoKmqaO0p3yZq/HbWGjkj0uPaAqu60/WnEe+a6ZO/oBeWZdhP76swe/XU1cvo545hTwwmhxn6zHaQReMCxxY1mxGsQg0b4xvery6YV2072XFlVeX6CmIsgbD0yWTQkN/VcZ62I1P7Nn7UYOSVE143eR+MXVQyh1+vB7+iYVkzU9JVZya+DG9JWJCodIpLMEvVlFXLqP+MEJ7/QytpFXqBSrdqTpJbJatK50dHboEe+yxx1544YWlS5eGQqHW32YY5pxzzlm+fPmxHu90OmfPnq3PNW0PWZZ5ns/opJQTany3MPDqsiNOJC/n4d2p8PZu7arl2k96UQ+cdMyN2WxlW5S6ZjnVw0tvuJDJb0fQ9JpYuPbP2/qcx3T0SJcOiVYlq94/oMlan/NLusNizdYIIYqi8Hwm/PWkA89DxXR3rzMLqv9Vt/Gx6h5nhEpODxl6DMieFfUlE/Kc7hP8uRfz8M7Z8PxWbconcOcg922j+Iz42HqjSrttlTajjNp0ERNoq8mUPzIQ2ZjqNTmzq5UkSWmsIzQA/O1vf3viiSc++eSTkpJDy6UJ+bFzfPny5RUVFel6vSzQ8G0kd5CPOUa76YLe9KaLWAcNg9+UX9mpmVy2DiEAT2zSJv2LuqaSvDulXSkIAA4Pkz8qULus0aBSpeqlzYuqt768t/jU3JG/7NvdUrDb4nxs34tKht1UHq9OrX1ge93aMBgzQieG5aYtsaJTck58VwAAuLY/veKn7Cu7mXM/VmoSth41rIqRn/xLuXed9sYk5qnxbacgABSMDtavDZtbNLtj169fP3fu3HHjxt1xxx36TfPmzbvhhhsqKipycnJWr169Z8+e5557ztpS2krd2nDp2QXHuYPfAU+OZ67sS89drr6wTXtqPFNpv33cq2Lk2mWq+P/bu+/wOKqzUeBnZnvv2q4uS7K6ZFuW1WzLGBvZcUmAkCdfuJAQEuAh4Qtcbj4SUuALfElI4IZQLjWhBUyLjXGwJYNkW7Jsy1axJMvqWmm7tteZnZn7xxKhONhI1kq78p7fwx+7s2VetJ55Z855zzkEaNlOZS1w2nBtvbz7qVF1jSy2A8IifmLysNV2zq3dKM/9Dz0cF5iEOApW3v9K9YwFxvabjS0zGbtUoqwYj1ud/IdVVSld0LIpOSLkSAP251Fe+YeRJypp385OuBF4GAn+0Es+0Uv8pIj2k6KvuJ0WpHIpkvJNBfm6a6UudtHoOp3u7bffnrtJLpc/9dRTJ0+e9Hg89957b2NjI4cD/16f85tCYQcmXvXVB2dlCnJmN/3pfrL6QOTOfPSnJTRuQgwjBiQF/txPPnKO+N8ltPsKUQJfcM0eW8ZUb5COvGfMvy01JiFRBGU8PjPVbJeXCssfzGHwVlS5ERRrwgxuyb2Z9m730FvTPC07fYeSo4hNI5jros897C97IHuhH6Sj4GdlaGMqclsL8fYo+Vw1TctLlAu1ZiN1TxuxSoSc3kVPv+KycZ9DQEqF2NzuzL4Rntg/h8xtBb0K9fX1v/71r+ffR+jz+aKDNBaz0zgafM3A13G0m+Tz/4gpAO7vII5bqN+uQ2/OjPO1ZK+D+sEJgoaAF2tp0fXGolWjC62Iowjq3BMjqVsV8lLRogKigL3HPf6Rhatipe9QcZXJ0jF2BRRF+f1+Ph+2CQMyQhlbZ6Y/syvKRKlbU+iLu0IiMPLcb4ezvnE1vc5er1cgEAAAcBI81k0+3Uc8XE67Kz/OE1WPeakHTpHn7NRTVbQdqQsIBfcTnY9fLL8/mylaqcXJs79ITCz3TUrIgg3/3QIQQOfQEBpCY6J0Lo3BozH4dKaQzhQxWBIGS8RIzLlCgtawa8iffdPCFrNWc8Ebm2jHzNSP2ok/95O/r6StU8Thf8+Lg0fOEX8ZIh+toH0vb7HHL0JDcm7WDrwyIcrhX/UNnHvEP/6RmSJAzs3aBJm5DUooKB3RbZYrKyWTn1g7H7+o3SjX1MquekXliYMWURZvkb3ODBQ8XIbemIH88ATx1yHymWramnhMyebGwP90Ey8Mkj8upL2+kbbQ9bEZPJpyrWTqU3vm7gRdgioSIsIOPOzCMTeOeSMRP4EHCBIjiTCZ801tzKcUXu5EyJQysm5UIxQSCREUQREYGf0/DDkwz3gAc+OhGSwSIFhSJlfF4qaweBo2T8vmyFmJkBoNzTZNrezqVrarVSFndtNfHSK/3kRUK5FHKtBl6zgkKfDqEPnzM+T1OqR3LyNW82UI0jjKdZL+F8cL7kxf6EL1vqng5D+sAXM4bXuKolycCD8ulLAYPFrWXrWmVjbxsaXz8SH9FoWyUoLQFvaPxnR8xtHvLf3P2NT95YuRTxvpfx0idx0mtuqQ36yhqZdrYF6YAM9fIH/TRTTq0a499KtuodVukp/9nyF9g2JJy7/nKRIgfNNBvzEUMIeDlnDAFqYiFFvKZEkYTCGdIWSwZUy+nkZjoTQWSufEvutkuf8EKAPh6ThXbholI1TIFg5YsaAlZOtyjx+04P6IIJUrSOMIM3jCDG5cFlkNOTBnvzfzoau/gEIRcPsq9JuZ6FN9ZPWBSGMq+rNSNEu4hEmAAmD/BPlwJylggL9vjf2la9p2JREm+55fQC70TQYNzTbvZFDfIM+/LXWhpzMoaXEUzLxb9b7J4MQn1qmjdl2DPGWtZJ4VVcZjM8bWmaK7MmJ4DkUAuDUH3ZuO/qaLKHoPvyMPfaCYJl3Kpn2MBC8Nko91kaUypOkGeqFkUccOU0BPqRBPt9jTd6hiFeECUCBgDrlHAt6JgGcigPsiPA2Hr2UL0rjKtWJOCovBX9bctDL6CCMBwjsR8E4E3aN+nyHIVbMlq/jiXL4gjbNsE0BffGOKJWOmbbtSvej8uTHw5Hnyz/3EZg36QDFaEesURVLgwwny8W4SI8Aja9CdqZe9dLi6PsIvUGD07yb3iD/7GxpB2mWviimCcvR7ja0zIQemrZerqiRLOlBsRYN9hF/JMx6YarJ5DUH1BqmyUnKFhrJIkDActs6c9xbdlcGSXH172pV7pAx+6tFz5Pvj5A/y0XtWx34WfjcGnhsg/9RPlkrBLytidkUbduHnnhgu+8/sxfxlFrpH5wWf66LPPeSjc+miLK4gnStI43JTFtzmF9s+wpWRCOciI5R3POAc9LkGfSEnJskVyAoFknzBkt4mzvR6xj+ylP0k66r7J76UDwcvDZJ/PE+queDOPPSmTHTxlaX2EHh9mPxzP5nCAQ8Uo7vSvuJKYbGJEABAAesZ18QhizCDq6qS8vWc2d8C80aiP9ZMt5urZquqpPISYQKu0ZpQYCKcp4AlbDo2Y+92c9VsebFQkMHjqVmf/+uiQMAadl7wTh21ywoEaduVi2wAnM9pd9RL/aGXfGuE3JOO3pGLRhcvW6TTNur5C+T74+TOVPT+IrRIGuNjZ6rZ5rroL/xB+hJ2T1DANxWc6fU4+r2YGxfnCSS5fPEq/iIHXyV7IpwLc+OOfu/Mea93LCDM4spLRbJCYcwzIu6LnPv9cP5tqVe441kMggIfG8j/d4E8bqZu0KM3ZiJbNOhC10lxY+CQgXx3nDpqJHemoj/IR6vmdxzGIBECAAAgMXK6dcbZ7/UbQww+HSAgEiAAAoTpXGEWT1EmgkvmzhNMhAsSbWxwDng9Y4GQE6dzaCgdiQQIOpcmzOBqN8p5mhjMzDf/0649BF4aJF++SDJRcEsWuisNKVhgGyYFQI+Den+M3DdGYST4Xi562yp0idZ6o0iq9+kxeZlIUyv76ncv7KuBdzJo73LZezwoA5UVCaUFAkFqzNrwYCL8EpEQ4ezz2rrcnpGAeBVPUSGW5AtiMyKbAgOvTHJVrLQblnwxelsIvD9OvjtGnrRSZTKkXo2USJEyGZLKR/69HREnwZiXGnBRJ63UcQvVM0PVq9Hd6cjX09HLzSjxpWKVCGdRBBV24QAAGpvG4NJgIcxCwUR41QiMjAQIiqBoLDS2nUwLPe1SABw3U++OkX+foGgIqFEhNUqkRIbkCJEvnXnRHgJ9TqprhuqwUUeNpICB7EpDbspE1yqWfKhZyI51/9/R4nsyYrUkVsActna67OfcKAORl4rkJSKuKvZ9pzARXkkkSMz0eKydLr8pJC8WpawRC9O5V38upsDwe8aAMVR0d8ZylnUEIuC4hTphJrsdoMdBTfspKQukcD5Ph4EIcGGUMwy0PCRXBNYp0GolUqVEeFd14Mc8EUKLBBNhAlrMabfPSZ2wUG0W6ryTuuimUATIWIj4n51ijjAwBSguHeSLkRIZskaObNYgafxlvX60nHYaPrEV/iCdLb/68wDmjdjOumxnXLifUJSJFBXimNyLX87KHke41OgcmrJSoqyUhF247ax7eJ+RxEhFhSilQrzQ6x2KpIbfNoYcWMGd6ctc3Milg61aZKv28wo3CgBzANhCFE4CAACHDsRMIGcjMe2vhCBoSRRIkAIJ8v28z586wsAZplz/nNBJwgJqDsKJ65lYuVZCEaD3mbGCO9MXOqkFESZnej22sy7vRFBWJMzYpRZl8VZcO9C1lghnscQM3Wa5brPcbwxZO129z44z+DRFqUhWLJzPdE0Bc3j0AxNAQMEdabEtkLkKCABqLlBzV9o/LgiC/o2UBaSJt46Far0EZSC9z4ylbVeqKiVfmcmIMOm84LN3uVyDfmEWV7lOkn9b6sotBb9mE+EsnoadoVFl7FC5R/32bvf5Z8dRBiLJF4gyeYJ07iWVSxRBeScC1rPumR6P/jqFeoMUDnSDICgZpFSIuUrW6Idm04kZ3Ua5OJd/ST8rGaF8U0HvWMA56PNOBoTpXHmJKPsm7VKMcF9m134i/BwCRFk8URYvaw/wGUOuC17LKefwvmmAICwxgyGgUwQVCREhK8aWMyWrBRX/J4fOXfG/LgRB0PzxdZziezLsPR5bp2vkfRNTSGfw6SgTJcIk5sIxX4SrZAnTuepqaf5tqXGZ2GSJJE0inIUAvpbN/+da57g3EnbjuDeC0BAam8aWMeHSBxAEJTN5sVBeLKRIKmAORwIEiZMoE2WJGUwR41pdHC35EuG/YgjoiTDVHgRBUEJBUGRJyz4TyrVzbwtBEARBVwEmQgiCICipwUQIQRAEJTWYCCEIgqCkBhMhBEEQlNRgwSQELbdQJISTEQCAHw+QFPnPjWGcxAEAFEUFg0Eu9vlSJwLmF5OOcugcOkoDAPCZPGTFTWMFQYkKJkIIuko4gbsxryfs9YQ9XsznwXxezOfD/H48EMADwUjIh/mDkSBG4AE8GIwEIyQRzXxsOouBMgAAXAaXhnzeKjO7EQBAkiSKfr7di/lm9xjAgwRFzG7kMjg0hMZjcJk0BofO4TI4HDqbQ2dzGVw+k8dn8PhMHp/JEzIFAhZfxBKKWEIeY0mWEoOgFQ0mQgi6LE/Yaw86LH6bM+SyBezOkNsedLhCLmfI7Qg5cQIXsgRCllDI5AtZAgGTL2Dy+UyehK3mMrgcOlvA5HPobCaNGU1RdJQeTV1X3un8V5+I5kUf5sdJPBgJ+bFAiAgF8VAgEvRiPj/mtwbsPsznwXwezOsOeTyYFydwEUsoZovlHKmYLZRxpDKORMaRyjkyBVcm40ijd5wQlFRgIoQg4Ay5TT6zyWcx+Swmv8Xit1kDdovfykSZMq5UyVVIOWIFV5Ym0lWoisVssYgllLLFfCYvvmFzGRzwr22nXwkncFfYHU3krpDHHnRMe83d1r6ZoMPit7tCLhFLqOSlqHgKJS9FzVdG/1PxlDBBQtcwmAih5BKKhCY90wbP9ITHYPAYp7zGaa+JhtK0fHX0pJ8rza7Xb0jhKVQ8BZt+rc2swaAxFFy5giv/0ldJinQEnSa/1eK3mv22QcfwZ5MnTD6LLTCj4Mp0Ao1OqEkValOFujShTs6N9ZrmEBQnMBFC17IISUx4DKOu8VHnxLh7csw96Qy5dAJtqlCrF2o36NbqBVqtQL2gm6prGIqgcq5MzpUVKfLnbo+QhNlvmfIaJz3TY67Jlsm2Cc9UOBJOE+kzxWkZorRMcVqWJF3EEsYrcghaDJgIoWtKKBIeco5edIwMOUeHHaOTnik1X5kpTs8Spzdmb80Up6l4ShSB9ZYLQ0dpOoFGJ9Cs16yZ3ejFfOPuyTHX5KhrotXQNuIa59A5OZKMHGlmjiQrT5Z9uftOCEo0MBFCKxtJUePuyX77YL998MLM0LTPnC7Sr5JmFchzd+dszxCnsWjMeMd4bRIw+UWK1UWK1bNbTD7LsHP0omP04MjhP5x+FgCQJ83Jl+cUyPPyZDmwYBVKWDARQitPKBLqt1/stfX32gb67YMyjiRfnrtatmrXqu1Z4gxY1hEv0U7WWn1V9Kk1YL8wM9RvH/xL79+GnKMqXkqhIr9IkV+cUqDipcQ3VAiaCyZCaGUIRUK9toEuS2+X9fywczxHklGoyN+b2/jz6p/ArqnElMKVp3DldfoqAABBEcOOsV7bwImpU8+ee5WBMkqVhaUpheWqYpgUobiDiRBKXCRFDsxcPGPqPmPuGnKMrJJmlSmLv1fyH6vlubDBc2WhIbRcWXauLPsbeTsBAJOe6R7r+U5z9wvdr7FozApVyRpVaYWqRMgSxDtSKBnBRAglHHvQcXL6zGnTubPmnhSefI2q9D8KbyxWFLDprHiHBsVGqlCbKtTuyL4eADDuNnSauw6Pffq7jqdThbp1mrJKTUW+bBWKwJmQoWUCEyGUEEiKGnQMtU2dbp8+bQnY1qrLqnXrfrz2TglbHO/QoKWVLtKni/Rfz90ZIYnztv5TpnN/OPWsLTBTqSnfoF23TlMOq2ygpbbciTAQCU45TMg/69cZKINNZ7FoLCaNAcdyJSGcwM9Zeo9NnTwx1SFg8jdo1/1o7fcL5HnwbiAJ0VFaqbKoVFn0/dLv2AL29ukz/xg9+ruOp/Plq2p166t16+B4jCQUioQxEvNhfoIkApEgACBDlMqMdc/IcidCe9DxZNfzs09xEg9FwiEijBGYD/Nz6Gwegytg8gUsQXQWKwlbLOWI5RyZnCtVcGRSjmSZA4aWAkZgHcbOlsn2k8YzGeK0Gl3l01v/R8NXxTsuKFEouPKv5Wz7Ws62UCR02nTu+FTHSz1v6ATqOv2GjanVar4y3gFCMRCKhK0Bmy0wYw/MOEKumaDDFfK4wx532OPFfF7M58f9DJTBpDEFTD6KoNG2gV/VPhjzfwAIRVGL+Xx9ff2vf/3r+vr6eb7f5/PxeDzkMiOag5GQD/P5ML8H87lCblfY7Qy5HEGXLTBjD87YAnY/HlTxUtR8pVag0vDV0flBVLwUeANx1TAMAwAwmctRe4IT+Elj56eTxzuMnbnS7PrUDbW69fDi5hLzn3Q7qRAU0WU53zLZ1mpoT+HJN6XWbE6rUS5XxanX6xUIYCHP1fNgXoNnetI9Ne0zT3tNJp/F7LcE8KCSp5BzZQquXMoWyzhSEUsgZolELKGAxY+uoHK5c3tsf5HE6iOMLiJzhQaQMIFFZ0ae9pqmfcYOY6fBM+0Ku/UCbbooNUuSniVOz5JkyOC5NZGQFNlp7m4ab2mbOp0lSd+cVvujNd+HYx6gBaEhtApVSYWq5Mdr7+y29h2dOHbHof/UCzUNafWb0mokbFG8A4S+gBHYmGty2DU24hyPTm0YjoT1Qm2qUKsTaKp16zR8lYqXkjgXwYmVCL8Si8aMdq3P3RiKhCfchjH35Ihr/PTAB0POUQZKz5Vm58qy82Wr8mQ58JwbL4OO4SNjnzVPHFNyFVvS679feiu8RoEWCUXQMmVRmbLoR2vu7DR3NY23vNTzeoE8b2vGxhrdelhaHBcRkhhxjV2YGbowMzQ4MzzlNeqFumxJeqY4fYN2bbpIn+BTtK+wRPil2HRWdIjS7BaL33rRMXLBMfz2wIcXZoakHEmBPLdIkV+kWJ0q0sGlvZeaPTBzePyzT0aPYgS+NWPTn657TCfQxDso6FpDR2mVmopKTUUoEj4x1XF47LMnTz9fo6vclrm5RFkID/Ol5sV8vbaB87aB87aBi44RNV+ZL1uVL1u1O+eGTHEag8aId4ALcC0kwn+n5KUoeSnRqZ5IippwT563X+i29r+/SNX9AAAUCUlEQVTe924oEipJKSxJKSxXFV9yZwktEkZgx6c6Do02X7AP1aduuL/ynkJFHjwfQUuNTWc1pNc1pNc5Q64j4y1/6nzRh/m3ZTZsy9wMy2piy4f5z1l6u6y9XZY+k8+cJ8spTll9a9E3V/pcsolVLLMMbAH7Ocv5LkvvWUsPRmAVqtK16rK16tKkHa8Wk2KZYefowZGm5vHWHGnm9syGWn0VnPnlqsFimcUbco4eGmluHm/NkqTfkLWlTl+1yIL7ZC6WiZDEefvAadO5M6Yug2e6QJ5XriouSSnMlWXRkLjN6xvbXyTpEuFcJp/ljLnrjKnrrLlHyVNUasorNRWFivykqkFdTCL044Gm8ZaDw0fcYc8NWVu2ZTYoeYpYB5h0YCKMFZyMnJjqODhy5MLM0Jb0usasrdmSjKv7qiRMhBa/rcPYecp09qy5J1WoW6suW6MuXS3PZaAJ0Y4IE2HskRTZZx/sMHZ2GDvNfutadVm1rrJSXc5n8uId2pK7ukTYa+v/aPjw8amOteqyHVlby1UlcJG/WIGJMOasAfvHI0c+HmmWsEU7src2pNVxGZwFfUOSJEKSoi7MDJ2Y7mifOm0POio1Fes1FWvUpQlYbwgT4dKKTnTZNn2qy3I+T5ZTo1tfrVt3Dd/oLCgRusOeT0aPfjRyBACwI+u66zM3J+ARstLBRLhESIo6Yzr30cjhs+aeOn3Vjuytq+W58/zstZ0IMQLrNPccnzp5YuqUhC3aoFtXrV2XJ1uVyFe31/I4wkQg50h3ZG/dkb01TGDRKS3+0vs3JU9Rq19fp69KS8r6GgpQZ809Hw0fPmU6W6OrfKDy7rnLsULQioAiyDpN+TpNuTPk+sfo0d+0P0lH6Tuzt16XsVHIvGaT3BX48cBJY+cxQ/sp49kcaVaNrvLbBTcmZ3kRvCP8aiRF9lj7Wg0njxnaOQxOnb6qTl+1SpoV77hi48p3hPbAzKHR5o9HmrgMbvSUsaJrw1YEeEe4PChAdVv7Dg4fbps+vV6zpjHrujJV0eWKnK+lO0JP2HtiqqPV0N5t7StOWV2rr6rRVa64ph3YNBo3FKAuzAy1TrZ/NnkCAFCfuqFOX5UvX7WiRwh8aSLECfzE9KmPR5oG7Bc3pdU0Zl+XK82+zBdAMQYT4TLzYr4jYy0HRw778cD2zIYvrfm6BhKhM+Q6ZjjZYmgbsF9cqy6r01et165Zude1MBEmhGHnaKuhvWWyzY8Ha/Xr61M3FCtWr8Ry00sS4eDM8D/GjjaPt2aK027I2lKfWg0HQiwzmAjjZdAxfGikuXmiNUeSuS2zoU5fNTtPzcpNhNaAvXWyvdXQNuqaqNRU1KduWKcuvwbm34GJMLFMeqZaJ9tbDG1Wv71GX1mjW1+hKo75KiFLJ5oIZzBn83jr4fHPImTk+oxN12duVi3XdMbQJWAijC+cwI9PdXwydvS87cIG3brr0usrVCV+n39lJcIJt+HYVMcxQ7vJZ6nWravTb1ijKllZs71cGUyECcrstx4znDxuODnkHF2jLq3Srq3SrhGzEnouYGvAfnTsWIuhzeQ316dWb83YWCDPi3dQyQ4mwgThDLmbx1ubJ1pMPmu1et11WfXFKQWJ3OpDUESvtb99+szxqQ6cxKt1lbW69aXKwkSO+arBRJjoPGFv2/Tp9unTZ8xdqULdes2aSk35Kml24tQiT7gNJ6ZPHTO0T3vNVZo19boNlfqKOE4SAc0FE2GiMfrMhwab2i1n7EFHja6yRl9ZrkygVh9nyNVhPHvK2Hna1KXmK6u0a2v0lTmSzHjHtbRgIlwxcDLSY+07ZTzbYeycCTnLlcUVqpJSZWGqULf8wfjxQLf1/CnjuQ5jJ0ER1bp1Nbr1pcpCAifAcq1HCM0HTIQJKHraNfrMxw0nT0ydGnKOlqQUVmrKK1QleqF2+ePxYf4eW/85c88Zc5ctMFOhKqnUVKzTlMs50uUPJi5gIlyR7EFHp6nrrKWny3IeI/FCeV6hIj9flpMjzeLQ2Uu0U0fQ2We/cN52ocfWN+425MtWrVWXVWoqMsVps+9ZzoV5ofmAiTABXXLa9WDeTlN3h+lsp6mLBFS5sqhAnleoyM8Qpy5RywpJUVNe44WZoX77YK9twOQz58tXlaUUlatK8mTZ12Tj55XBAfUrkpwjvT5z8/WZmwEAFr/tvG2gzz742eTxUdekkqfIFKdliNLSRDqtQK3lqxc6/1OUD/MbfeZJz9S42zDiHL/oHMEiWIEid7U8986y/1Ugy72WusohKI6ETMGmtJpNaTUAgGmvqdva12vrf//iQYvflilOy5FkZohTU4U6nUCj4MquIkuRFGn2W41es8E7PeaaHHNPDjtHRSxhniwnuvLiKmk2HYV9GTEDE2EcKHkKJU/RkF4HACAoYtxtGHdNjronmsdbp72maZ+JiTJlXKmCIxOyBEKWgE1jCVh8AAACEB6D68P9AIAAHvTjfh/md4bcM0GHLTBDUqSar0oVatNF+m2Zm++V3pGck0RA0HLSCtRagfqGrC0AgAAeHHaODjvHx1yTLZNtU16TK+RWcGUyjkTMFotZQj6Tx2fwUBSlo3SKogiKAACEImGMwFxhjyfscYc9Fr/dFXJJOVKtQKUXaDPEqRtTq3OkmQImbCRYKjARxhkNoWWJ07PE6Q1zNnrCXltwxh6YcYe9nrA3TIS9YR8AgAKU0WfmM3gAAA6DLWFrBUyehC2WcaRyjlTIWknl3RB07eEyOMUpBcUpBbNbcAK3BuyOkNMZcrvDHh/m92E+kqJwEkcRNNqOyqIzRSyhXqgVsYRiliiFJ5eyJfCGbznBRJiIojeCWeL0eAcCQdCiMGiM6C1jvAOBriTpulghCIIgaC6YCCEIgqCkBhMhBEEQlNRgIoQgCIKSGkyEEARBUFKjAwBOnTr16quvAgBuvfXWysrK6AvvvPPOoUOHlErlPffco9PFYUowCIIgCFoGaHd395YtW3Jzc/Pz87du3Xr27FkAwHPPPffggw82NDQEg8Hq6upAIBDvOCEIgiBoSdBCoVBNTc3DDz9cWVnpcrmampp27959yy23PPnkkzfeeOP27dvffvttFotVVlb2pZ9/9dVXN23alJ6ePs/9YRjGZDKTcK7RhEUQBACARoOjdxMIjuNw9teEgmEYi7XiF7O9lsT2F0Hb29s3bdoUfbJx48a2tjaLxTIyMjK7cdOmTW1tbbHaHwRBEAQlFLrZbJbL5dEnCoXCZDKZzWYWizU7sbdCoejp6bnc5/1+/0MPPTT7DVE7duz41re+9eXvnxr1t74H7wgTB0mSAAAUhWVTiYKiKIIgAnQ461MCiUQiQfiLJAbuzu+hQmkwGJxnOxabzf7K8xudxWJFF+IBAITDYQ6Hw2azcRwnSTL64VAoxOFcdjEEJpN53XXXFRQUzN2Yk5PDZn/50kK4RE5f2wATYeKIRCIAADo8yBMGRVHhcPhyRxAUF8Fg8AqnQWg5sUQShMXGcXyex8h8rvLpOp3OYDBEnxgMBq1Wq9FoKIoyGo3RYtHoxst9nsFgbNy4cf7rEdJ4Qm5pLUyEiQOuR5hoKIqi/H4uXI8wkRBeLzd2q99Bi4eiaAzbsdA9e/a88cYbFEVRFPXmm2/u3btXKBQ2NDS8/vrrAACPx3PgwIG9e/fGan8QBEEQlFDod9111zvvvFNTU4OiqMvleuWVVwAAjz32WGNj44kTJwYHB+vq6urq6uIdJwRBEAQtCbpUKu3s7Gxvb6coqqqqKtpEtmbNmsHBwY6OjpSUlMsNnIAgCIKgawAKAGAwGHV1dfX19XM7isRi8fXXXx/zLGgwGHw+X2y/E1oMm81ms9niHQX0BZ/PN9ttDyUCgiCGhobiHQX0L0ZHR2fLPBdvuYvm77333uPHjy/zTqErePbZZ59//vl4RwF9obW19cc//nG8o4C+4HK5tmzZEu8ooH+xd+/eiYmJWH0bHD0GQRAEJTWYCCEIgqCktthh1JFIpLu7e/7vdzqdvb29XC53kfuFYmViYgJF0ZaWlngHAn3u/PnzDocD/iKJw+PxRCIR+IsklFAodOrUKaPR+JXvXLNmDY/Hu/J7EIqiFhPNW2+99fTTT89/XhKfz8dms+E8JokjFAoBAOA8JokjEomEQiE+HFCfMCiK8ng8IpEo3oFAX/B4PHw+fz5j6l966aXs7Owrv2exiRCCIAiCVjTYRwhBEAQlNZgIIQiCoKQGEyEEQRCU1GAihCAIgpLa0lZv4jh+9OhRj8ezadOmSxbvBQCcPHlyfHy8vLx81apVSxoGNMvj8bS3t3u93qKiotzc3LkvdXZ2zj5WKBSpqanLHl0y6uvrixbuAgAEAsElx8LIyMjp06f1en11dXU8oktGc38RAIBIJJqtOZyenjabzbMvlZaWznNtWOjqWCyWqamprKwssVgc3eL1epubm1EUbWhouGRQBEEQn376qcPhqK2tVavVC9rRElaNhsPhzZs3UxSVlpbW1NTU3NxcXFw8++p99923f//+urq6jz/++He/+913vvOdJQoDmnXhwoXKysq1a9cqFIrDhw/fddddjzzyyOyrNBqttraWwWAAABobG+EsX8tj9erVPB4vepyXlJT8/ve/n33pnXfeufvuuxsbG9va2jZv3vzcc8/FL8wkctttt01NTUUfnz17du/evS+88EL06UMPPfTGG2/k5OREn37wwQdwlMvSKSwsHBsbw3F83759u3btAgCYTKaqqqqioqJIJDI8PNzW1qZQKKJvJkly+/btdrt99erVhw4d2r9//4YNGxawM2rJvPbaayUlJTiOUxT1s5/97Bvf+MbsS6Ojo1wud3p6mqKoTz75RKPRRN8GLSmn02k0GqOPz507hyCIzWabfRVF0blPoeWRn59//Pjxf99OEERGRsYHH3xAUZTVahUIBP39/cseXVILBAJisfjYsWOzW/7rv/7r/vvvj2NISWVgYADH8dzc3A8//DC65cEHH7z55pujj3ft2vWLX/xi9s0fffRRZmZmIBCgKOqJJ56I3oPN3xL2Ee7fv3/Pnj3RsfM33njjgQMHqH/efR48eLCqqkqj0QAAtmzZEgqFzpw5s3SRQFFisXi2xUCv11MUNbcJCADQ2dnZ1tbmdrvjEV3yGhgYaG1tvWQNkN7eXovF0tjYCABQKBQbN248cOBAnAJMUvv27UtJSbmkUdputzc3N1+8eDFeUSWPvLy8S6Ze2b9//0033RR9HM0psy8dOHBg586dHA4n+tKnn37q9Xrnv68lTITT09M6nS76WK/Xh8Nhu93+7y+hKKrRaKanp5cuEujfPfbYY1u2bJn9FQAAUqn0t7/97X333Zeenv7ee+/FMbakwuVyX3vttZ///OcZGRl//OMfZ7cbjUaVShVtqQYA6HS6+cwmBcXQyy+//N3vfhdBkNktCIJ0d3c//vjj1dXVO3bsCIfDcQwvCV2SU+ZmjbkvabVaBEEWdLwsYbFMJBKZ7UmOPsBxPPoUx/G5U+PQ6fTZl6Bl8PLLL+/bt++S9bCMRmP0tPvmm2/efvvt27Zt+8oJ+qDFa29vj/7ZT548uXHjxsbGxmi9DI7jcwsx6HR6DFdfg77S6OhoW1vbW2+9NXfjL37xi0cffRQA4PF4NmzY8Kc//en++++PU4DJ6JKcMveImPsSgiAoii4opyzhHaFKpZpt7bFarTQaLSUlJfpUrVbPbQiyWq0LLfKBrtobb7zx8MMPHzlyRK/Xz90+e/PxzW9+MxQKwZVIl8fsn339+vVZWVldXV3Rp9FjZLY3wWq1RrsSoOXx4osvbt++/ZLz0uyPJRQKd+/ePbfQGloGcxPHJUfE3HTjcDgikciCcsoSJsK6urqmpqbo46ampurqajqdHg6HQ6FQXV3d8ePHox1UfX19Xq+3vLx86SKBZr377rsPPPDAJ598Mlumj2FYMBic+55oH7VWq41HgMnL4XAYDAatVhuJRPx+f2FhIYqi0b5zHMdbWlrq6uriHWOyIAji9ddfv/3222e3eL1ekiTnvqe7u3tuzwK0DOrq6o4cORJ93NTUVF9fDwAIBoM4jkfTTfTCsampqbCwUCaTzf+bl3D4hMPhKCoq2rNnT1ZW1qOPPvr6669v3779hz/8od/v/+tf/3r99dcjCLJr165nnnnma1/72n//938vURjQrL6+vtLS0m3bthUUFES33Hnnna+99tqxY8fuvvvu999/v7y83O12v/DCC1//+tefeuqp+EabDLq7u3/1q19VVVVRFPWXv/wlLS3t4MGDf/vb337605+Oj49Hj5p777330KFDbre7tbU13vEmi4MHD373u981GAzRW0AMw1gsVmdn5y9/+cuysjK5XH706NETJ050dnZe0qwCxdDzzz8/Njb2wgsvbNy4MScn56677nK73TU1Nffff38kEnnqqac6Ojpyc3Nra2t37tx5zz33FBcX19bWlpWVPfbYY0888cS3vvWt+e9raVefMBqNr7zyisvl2rNnT3RUR0tLC47j0UrRF198cXR0tLKy8qabbprbIw0tEZPJdEnl4c6dO6empkwmU2Vl5b59+8bGxvh8fm1t7ZYtW+IVZFLxer379u0bGBhgMBjl5eV79+5FUXRoaKi9vT06sva99947ceJEWlraHXfcAVfxXDbt7e2BQKChoSH6lCTJZ5555uabbz5z5kx7e7vP58vOzr7lllskEkl847y2HThwwGQyzT7dvXt3SkpKX1/fm2++iaLot7/97eiUIO+//35GRkZZWZnNZnv55ZetVmtjY+PmzZsXtC+4DBMEQRCU1OBcoxAEQVBSg4kQgiAISmowEUIQBEFJ7f8D8v1dGNjaT2sAAAAASUVORK5CYII=", "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 }