{ "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.293939920566 NaN 4.29e-01 8.0 \n", " 2 +0.442250635076 1.48e-01 5.35e-01 1.0 \n", " 3 +0.260880117233 -1.81e-01 1.72e-01 2.0 \n", " 4 +0.243061859882 -1.78e-02 2.71e-02 2.0 \n", " 5 +0.242665272662 -3.97e-04 6.15e-03 1.0 \n", " 6 +0.242635962554 -2.93e-05 2.18e-03 1.0 \n", " 7 +0.242634020177 -1.94e-06 1.01e-03 2.0 \n", " 8 +0.242633398802 -6.21e-07 1.80e-04 1.0 \n", " 9 +0.242633376905 -2.19e-08 3.16e-05 2.0 \n", " 10 +0.242633376341 -5.65e-10 5.01e-06 1.0 \n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown:\n Kinetic 0.0304517 \n AtomicLocal 0.0972470 \n PowerNonlinearity 0.1149347 \n\n total 0.242633376341 \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.03869210860187387, 0.0, 0.0]\n [0.038693388861610487, 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+gvaeTAAAgAElEQVR4nOzdd2AUZf4/8M+Undm+m2x6IIUQegcRUUEFRA8PEMVefoJ+sZztzvMs3zsPvfNEOUXPeueBnqJfFftZzoqCSpGmIJ1AIAFSt+/05/fHcDFAgJTdmdnN5/VXstnyJJtn3/N0ihACCCGEUHdFm10AhBBCyEwYhAghhLo1DEKEEELdGgYhQgihbg2DECGEULeGQYgQQqhbwyBECCHUrWEQIoQQ6tYwCBFCCHVrGIQIIYS6ta4G4Y4dO2KxWPvvr2laF18RJRchBLfZsxqsJlajqqrZRUCHSe470tUgnD179vfff9/++8fjcfzYtRRZlmVZNrsU6GeEkHg8bnYp0GHwHbGa5L4j2DWKEEKoW8MgRAgh1K1hECKEEOrWMAgRQgh1axiECCGEujUMQoQQQt0aBiFCCKFuDYMQIYRQt4ZBmGai1Qmi4Y4ECKW3aI2gKViRrYI1uwCoAxL10o9PV9lzuN4zizylTrOLgxDqMFXQ9nx0cP+3TcXjc8rOyze7OAgAgzCdENjxek3plHybi928qLpoXE6Ps3LMLhNC3ciGDRu+++67rjyDHFNqlzY6C3h/H/d7zzYU/Rjgs2zJKl5my8nJufDCC1P05BiEaWP/d01EJUWnBYACT6ljw+O7iscHKIYyu1wIdReffvrpyy+/PGbMmE4/A1EJFEAzQ9XsBFIEzTv30yxW4RNrbGzcsWMHBmF3Jwbl6o/rBt9UDhQAgD3AOfL45i3R7IEes4uGUDcyceLE+fPnm12KbmfdunWzZs1K3fPjZJn0UL8mmDPM58znW27JG+Wv+z5oYpEQQigzYBCmh+at0ez+7ta35AzzBrdFlQQek4YQQl2CQZgGVEmL7k14K1ytb2TtjL+Pu2FD2KxSIYRQZsAgTAOh7TFPqZPhjnyzsHcUIYS6DoMwDTRvjfr7uo++PaufO1EnSiE8Xx4hhDoPgzANBLdGstoKQoqhvGXO8O6E8UVCCKGMgUFodUKTpAqaq9De5k/dJY5oddzgIiGEUCbBILS64Jaov58bjrHo1lPiiFRjixAhhDoPg9DqmrdG2+wX1blLHNF9uA03QqhjotFoy84Aa9eu/frrr9u8W2Nj48svv2xgucyBQWh1kd1xby/XsX7K2hneb4sfEI0sEkIo3YXD4d///vcAoKrq7Nmz8/Ly2rxbIBB48skn161bZ2zpjIZBaGlyVNFUwvuPty2vG3tHEUIAhJDGxkZVPbTJRjQabfm65ZZE4sjPig8//DA3N7dfv37Hetrrr7/+kUceSXppLQWD0NJiNYKrqO1pMi08JU6cL4NQd3b77bffcMMNQ4YMGTFixA8//LB69ephw4addNJJPXv2/POf/wwA27Zt69ev3+jRo/v06TNp0qTGxsaWxy5evHjatGkAcNFFF73++usAkJ2dHQ6HP/roo3POOQcApk6d+vbbb8fjmfwhg5tuW1qstj1B6DjwXZMx5UEIHWF9I6k1NiMmFFE8c9gt8Xj8o48+WrFiRVFRUSQSGThw4KJFiyZMmBAKhcaOHXvqqacOGzbsm2++CQQCmqbddNNNDz/88K233qo/9ptvvrntttsAIBqNiqIIAM3NzYQQWZYjkQgAZGdnFxUVff/99+PGjTP09zQQBqGlxfYLvt7HHCDUuYrsQqOkihrDY/seIaP9Zx/56oBm5CuenMceEYQAcNFFFxUVFQHAN998Y7PZCCGfffYZAAwdOvTzzz8/44wzVq5c+cwzz0Sj0fr6+h07duiPIoTs37//WAOELfLz82tqapL/m1gGBqGlxWqFonGB49+HYihngT26L+GrOEFkIoSS7ndD6d8NNf8aNCfn0DHddXV1iqLoKQgAJSUlI0aMeP311//3f//3zjvv7NGjh8vl+vDDD/WfUhTF87wkSfq3hPw8/5wQQlGHlm0JguBwOAz6TcyAQWhdRCWJeqn10UvH4il1RKsxCBFC0L9/f0EQ7rvvvtbRdcUVV9x6663XXnstAHz//fet719ZWbl79+5+/foVFxfv2bOn5faqqqri4mIAIIRUV1dXVlYa9RuYAIPQuuJ1oj3bRttOfLHp7uFo3hI1oEgIIYs76aSTTjvttBkzZtx66602m23VqlVjxozp37//4sWLBw8evHXr1hdeeCE3N7fl/pMmTVq+fPk555xzzTXXTJs2LRAIAMCiRYvmzZv3wgsvAMDmzZvtdvuAAQPM+o0MgEFoXe2ZKaNz5vG1Xzee+H4IoUw0efLk7Ozslm9ff/31F1544fXXX5dleeDAgQMGDBg7dizDME899VRlZeX//d//rVq1yuPx3HHHHQBw7bXX/vKXv7z//vvHjh37+eefv/TSSwCwc+fODz/8cPjw4QDw6quvXnfddS3dpBkJg9C64u0OQkc+n6gXgcCxdmJDCGWwGTNmtP6WYZjZs2fPnj279Y133XVXy9ejRo0CgAceeAAA+vbte+qpp7755pszZ84cMmTII488Mn/+/D/96U8+nw8AotHoW2+99e233xrxa5gHg9C6orVC0eknmCmjY3ia4WkxJB9/6T1CCB3tueeeUxSlzR/Z7faVK1e63cfc5TEzYBBaV/tbhADgzOcTdSIGIUKooziO4ziu5dvWc0dZls34FATcWcay2rO5WmuOPD5RhzuOIoRQh2EQWlRsv3CsMwjb5Mjj4xiECCHUcRiEFpWokxx5J15B2MKRyyfqpNSVByFkWatWrdq8eXMSn3DZsmWt9yNtbePGjTt37kzia1kBBqFFCQ2SPYc78f3+Sx8jTF15EEKW9c9//vP9999P1rNVVVXdcMMNXq/3WHe44oorWo8jZgCcLGNRiQbR28vZ/vvzfpsSV3HHUYS6oeeeey6JzzZ//vzZs2fbbG1PUBg0aBDP859++unZZ5+dxBc1F35oWlRHW4RAgT2XS9RjoxChbmfu3LkLFy7817/+dffdd1988cXZ2dmjRo3avXv37373u7y8vL59+3733XcAcNddd5WXlxcWFo4ZM2bFihUAsG/fvrPPPrugoOC000576KGH7rvvPkVRXnnllQsuuCCRSIwaNUoUxc8///ymm24CgPvuu2/x4sUAcOGFF+qbzmQMbBFaEgGhWbYHOhKEAM48PlEnuntk8t64CFmNGmzQhJiRr2jL6wH0YcdP1NTU2O12iqKefPLJ//znP6+++ur1119/6qmnPvDAAwcPHnz66advu+22lStXnnfeeXPnzuV5/v3337/44ot37do1a9asYcOGffzxx1VVVePGjZs4ceKGDRvsdntJSUksFluzZo2mac3Nzdu3bweA3bt366vsR48e/dBDDxn5K6caBqEViUGZdTAM17H2uiMP58sgZLTo1+8IP60y8hVzf/Uw7fa3+aNJkyaNHTsWAKZNm/bpp5/OmjULAM4///zf/va3ADB69OhPPvmkqqoqHo83NTVt3br1s88+W7JkCU3TFRUVl156aX19/d69ewsKCo5fgKKiotraWlmWj9V9mnYwCK0o0SA5OtQvCgAAjly+cWM4FeVBCB2Lb+q1vqnXml2KQ/QtswGA5/mW3UftdnsikRAE4ZRTThk4cODo0aOzsrJsNlt1dTXHcS2TYvLy8urr6xmG0bS2j1dsOZhJURSapmk6c0bWMuc3ySRCYwcHCAFA33EUJ44ihNqydu1aSZJefvnlW2655ZJLLolEIrm5uTRN692eAPDDDz8AQHl5+b59+wghTqfT7XYfOHCg5RlajvDdt29faWkpwxx1OnDawiC0IqFB6ugAIQA4czmhQYKMmtWMEEqOgoKCffv2LV26dOfOnXPmzGFZlmGYW2+99aqrrnrvvfceeuihr776CgAGDhzIcdy2bdsoipo5c+Ytt9yyZcuWSCTy5JNPrlu3bsKECQCwYsWKM8880+xfKJmwa9SKhAYxZ5ivo4+iOZp1MmJQ5rMypOMeIdQeo0ePzs/PJ4S09HPm5+dPnDhR/5rjuJkzZ/bq1eupp5667777WJb91a9+lZ2d7ff7H3zwwYULF7733nuDBw++5ppr4vE4RVHXXHPNa6+99oc//OHpp59++umnFy9eXFVV9dNPP3333Xf68OFrr732xBNPmPbbpgLpmnHjxi1durT9949EIpqmdfFFM97aR7ZH9sY78cANf9sV3BHt0ENEURRFsROvhVJE07RIJGJ2KdBhwuEwIeSRRx75zW9+Y3ZZkqmxsVFRFEJIc3PzwIED33rrLUJIXV1d3759Y7GYfp833nhj0qRJLQ9ZtmzZOeecY3A5165dO2zYsNa36O9IsmDXqBV1bowQAOzZNqEJJ44ihNpl2bJlZWVlo0aN6tu376RJk6ZPnw4Aubm5H3zwgaqqbT6krKxMP7w3k2DXqOVIYYW20ay9MwPR9mxObJKTXiSEUEaaNm3aOeecEwqFAoFA68kvFRUVLV+fe+65p556asu3PXr0MLSIhsAgtByhQezE2gkdn20L7TB0bS9CKK3xPK/PBT0Wl8vlcrkMK48psGvUchId3VytFXs2JzZjixAhhDoAg9ByOrzLaCv2ACc04hghQgh1AAah5XRuWxkd52PlqEJUXEuIEELthUFoOUKTZA904Eje1iia4nw27B1FCKH2wyC0HLFJtmd3fkU8rqBAqBtauXLlWWed1a9fv+uuuy4cxj2HOwZnjVqLJmuqoNrcnX9f+GxOwBUUCBlFVmWFtL3kLkXsLE8B1fqWvXv3Tps2bdGiRcOHD7/zzjtnzZq1ZMkSI4uU7jAIrUVslvks2+H/5B1jz7KJ2CJEyCjzVz21tPpbI1/xlanPBRxZrW959dVXzz333HPPPRcAHnvssby8vHA43LLdGjohDEJrEYMyn9XJmTI6PsA1b44mqzwIoeO7+5Tb7j7lNnPLsGfPnrKyMv3rQCDg9XpramowCNsPxwitRWjq6pbZ9mwOW4QIdSuKouzbt0//OhwOh8PhwsJCc4uUXjAIrUUMyry/i0GIk2UQ6nbefPPNVatWiaJ47733nnPOOX5/20fYozZhEFqL2CzxXZgyCgCc16YkNE3BpYQIdSNXX3313LlzKyoqampqFi1aZHZx0gyOEVqL2Czz/i6NEQIFvJ8VmyRHXicXIyKE0k5+fv5jjz1mdinSFbYIrUVs7tIiQh2uoEAIofbDFqGVEJBCMufrahDifBmEupXbbrvN4XCYXYo0hkFoIVJYZp0MzXZhFSEA4HwZhLqZ/v37m12E9EYDwP3339+jR48ePXr88Y9/JIQAQDgcvvzyy3Nzc/v16/fWW2+ZXcjuQmzu6iJCHZ+FhzEhhFB7sW+88caiRYuWLVtG0/SECRP69et3ySWX3HvvveFweOfOnatXrz7//PNPOumknj17ml3UzCc0d3URoY73s2IQgxAhdGINDQ1Op9PpdHbu4TU1NQUFBa1Pt09H9PPPP3/zzTeXl5eXlpbefPPNzz//vCRJ//rXv+677z6v1zthwoSJEye++OKLZpezWxCbZXsygpDz2zAIEepWlixZUltbe5w7vPbaa3V1dUfffvXVV3el269v377Hf91jWbly5apVqzr9uslFb9myZejQofo3Q4cO3bp1a21tbTgcPuJG80rYjYhBKTktQp9NjihEw6WECHUXDzzwwLZt245zh9///ve7du0yrDwn9Oabb77zzjtml+IQtqmpyePx6N94vd6GhobGxka73W6zHfpE9vl89fX1x3p8MBg866yzKOqw+R0333zzAw880Ob9Y7EYIeSI+yNdrC7hKLVFo0nYKZRx0MGDYZvnxP0VkiQBAMclYWwSJQUhJB6Pm10KdJhYLEZRlCiKZhekbS+//PK+ffv++te/vvrqqzNnzpw4ceLHH3/80ksvybI8derUK664YtGiRXV1dfPmzcvLy7vsssvGjx9/9JMsWbLkjTfeoGn6kksumTZtGgBIkvS3v/3t22+/BYCZM2fOmDFjwYIFK1euVFV13LhxN998c0tMtBBF8ZZbbrnmmmsWLFhAUdQNN9wwbtw4ANi9e/cjjzxSXV09cuTI3/72t1u3bv3kk08oimpsbBw4cOAtt9xywt9R07TWn436O9KeP47T6aTpEywUZLOzs1sOrwqHw4FAIBAICIIgSZL+4RgKhXJyco71eL/f/+mnn+q/aguapo/zwi6XC4OwTWr0gK/Q43Lbu/5U9iyOlTi3+8QzqjEIrUa/UnS73WYXBP2MEOJ2u3neoptUjBw50u/3n3HGGUOHDu3Tp8/777//P//zP88++6zb7b7pppsaGxsnTpzodrvPOuus/v37V1RUHP0MixYteuCBB5555hlVVefMmROPxy+99NIZM2YwDHPnnXeqqlpVVaV/Vtxxxx0A8Mc//rGhoeHBBx884nkURfn73/++c+fO+++/f8+ePeeff/7SpUt79Ohx8skn33777VddddXjjz8+derUV199tU+fPgzDzJw5MxAItOd3pGm6daXQ35HO/8kOx1ZWVm7atGnixIkAsHHjxsrKysLCQpfL9dNPPw0bNgwANm3aNHDgwOM8BcMwLIvLMJJA7PKO2y04v00MyR7ApUUIpdaO12safjD0INwRv6vkPId95Pbv39/tdo8cOfKMM84AgFmzZt177716q+7RRx+96aabbr31VqfTedJJJ40ZM6bN51ywYMFf/vKXyZMnA8ADDzywYMGCvn37fvfdd3v37tXn0eitnTvvvLO2tvbAgQNXXXXVgw8+eHQQ6v70pz+NGTNm7Nixa9asefbZZ/v37z9kyJC77roLAJ5//vmioqL9+/eXlZWxLKtHj+nYa665Zu7cuRdffDFFUU8++eQ999zD8/xll1325z//+aWXXtqwYcPHH388b948s8uZ+VRBI4SwjuRMvuL9NgnnyyCUeuXTCsumFBj5iqzzBJ8SVVVVgwYN0r8ePHhwdXW1oigdesiuXbu2b9/et2/f1rNJm5ubp0yZomlaeXk5ABw8ePBYz9ayrnHAgAFvvPEGz/MtrSmn09mrVy9LjVYCAHvZZZf9+OOP+p9g9uzZV1xxBQDMmzdv9uzZeXl52dnZzz77rP5ro5QSmqWkLCLU8T4bLiVEyAAMT4MFekxbjzfl5OS0TOyoq6vz+/0syx5/QOqIh+Tk5OTl5R0RdYsWLSovL1+8eDEALF++/JNPPjnWszU0NPh8PgCor6/PycnJycn56aefWn5aX1+fm5sLAPqydSugKYp66KGHGhoaGhoa5s2bp4/tZWVlvfXWW+FwePfu3ZdddpnZhewWxCQtItTxWTYxhEGIUHeRn5+/adMmPVqmT5/+xBNPxONxRVHmzZt3/vnn63fYuHHjsR4+ffr0v/71r5IkCYLw6KOPnn/++WPGjCGEPP300/odduzYoWlaMBjUNE0QhIcffvg4hfnrX/9KCGlsbFy4cOG0adOmTp36/vvv61n44osvUhQ1cuTIgoKCbdu2ybIlPqZw022rkEJdPYmwNVxKiFC38vvf/37x4sW9e/deuHDh7bffPmjQoN69e5eUlBBC9NCaO3fuP/7xj4qKCr1J1yInJ8flcs2dO9fn85WVlfXq1atnz5733nuvw+F47733Fi9eXFRUVFRU9MQTT8yePTscDhcXFw8cOHDcuHFlZWUAUFxcfPQckR49epSXl/fv33/69OkXXHDBoEGDnnzyyV/84hfFxcWPPfbYm2++6XA4rrjiCkVRBgwYcM011xj1RzomqouN0/Hjx99///1tTsZtUzQaxVmjbdrzcR1FQ8nZeUl5NjEo//DErpP+0PeE98RZo1ZDCInFYjhr1FIikYjH45k/f/6BAwfmz59vdnHahRCiqmqHZjKqqkpR1BFz/iVJstlsLR/aoigeZ/as/q8rCILNZtM07YhXl2X56BUX7bFu3bpZs2atW7eu5Rb9HenEU7UJW4RWIQVlvsvnTrTgvKwcxTX1CHVfFEV1dD4/wzBHr3zjOK5106Wda0homj761TuXggbAILQKMRkHMLWgaIp1sVL4BFPFEEIoiXief/311y0beMeCQWgVUjCZY4SAKygQQoZjWXbmzJkn3MnFatKsuBlMCimcL5n7EvA4XwYhhNoBg9ASVFEjWtJW0+vwMCaEEGoPDEJLkEIyl9R+UcAVFAgh1D4YhJYgBmU+qf2iAMD7ORwjRAihE8IgtAQppCRxyqiO97O4uQxCCJ0QBqEliEndVkbH+3G7UYQQOjE8PskSpKDsLEzCMYSt2TysEleJSigG9/FBKAl69+796KOPfv75523+VJM1oCk62dVNFTTG3t1bLIlEwuFI4aFyGISWIIUVf98kvxcUTdncrBRWkriXN0Ld2fTp03v16qWqaps/3fXO/sBgr6/CldwX3fT33ZWX9jjiAMJuqPWBUEnX3f+4FiEGk3n0RAvOZ5NCKXlmhLqnIUOGHOtH1Jfe3mOK3D2S3HCh+3l7lRZ6SlMYA6i7t7gtIrn7q7XgfayIu6whZIiU1WKbGMRanFoYhOYjKlETKudOfutcbxEm/WkRQkcgKlHiqs2VzD0xdFiLDYBBaD4xJNs8LKRgRgvnY7EKIWQAKaxwHpaik1+NMQgNgEFoPimkJPEAptY4n00MYacKQiknpaZfFHBBsCEwCM2Xiv3VdDxeSyJkCDGkJH1zKB3ns0l4OZtiGITmE0PJ319Nx/lYrEIIGUAKp65FiJsGpxwGofmkYPL3V9Ph6AJCxpBCcnKPUWtxqBaTVDw3OgSD0Hyp2F9Nx3A00KAIba//RQgli5iC7YJ1NEvRPC3HsRanEAah+aSQzHlTtbMBjwMMCKWeFJI5b6p2ruB9NjxJJqUwCM0nhZQUTZYB7B1FyBBSSOb9qbqc5Xw2KYy1OIUwCM1GQArLqdtIkPPh3GuEUk4KKalsEbK4DiqlMAhNpsRVmqdpW6reCOwaRSjVFEGlaIrhU1WLOS+2CFMLg9BkYljmU3YhCXoVwhYhQqkkhZQUTRnVcV5cB5VaGIQmS3kVwqWECKVY6raV0eEYYaphEJospZPNQN+6HluECKWSmLJFhDpsEaYaBqHJpLCSurUTAMD5sWsUodRK3XbBOmwRphoGoclStzOTzuZm5JhKNNyXAqFUSXXXqM3NKHGsxSmEQWiyVLcIKZqyuVk5gv0qCKWKFEp5LWZdWItTCIPQZCk61bo1DhchIZRKYkhOadco4Ky3FMMgNFmqryUBD2NCKMVSd5JaC96Ls95SCIPQVATkaMqDEM+pRyh1iEbkmGpzMyl9Fc7LSmFsEaYKBqGZpKjCOhmKplL6KjyeU49QysgRxeZKeS3GiaMphUFoJin1QwuA+zMhlEpiKjfNb4EtwpTCIDSTAQOEgMPsCKWSQZezONKfShiEZkr1IkIdtggRSp1Ur4DSYYswpTAIzYQtQoTSnRRO7S6JOg73SkwlDEIzGdMiZB2MphJV0lL9Qgh1Q8ZcztqcjCZqmoKby6QEBqGZREOqEABwHlbGfhWEUsCYy1mgcIwjhTAIzWRQFcJ+FYRSxpgWIeAZFKmEQWgmKaTwhlUhbBEilAKicZezLLYIUwSD0DREJaqg2txGBCGPVQihFNAUoomazZnabWV0nNeGLcIUwSA0jRRWbB4WUrsfxSE2rw1bhAgl3aEpo4bUYs6H/TqpgkFoGmOWH+l4XEGBUApIYcVmVC3mvLimPlUwCE2T6sM8W8P5ZgilghSSeZ9RQYgtwpTBIDSNMetwdbimHqFUkEKKcbUYL2dTBoPQNEZ2jWKnCkKpIIUVzsgWIV7OpgYGoWmkkHFViOFpoEEVcHMZhJLJyH4d1s4Qjagi1uLkwyA0jRiWOY9BVQiwXwWhFDDychYAbLggODUwCE1jcBXifSwez4tQchnZIgQAHo/nTQ0MQtNIYZk3sAphixChpBNDspGXs7jLWopgEJpDU4gmaawhG1LoOB+L82UQSiJV0ogGrN3AWoyXs6mBQWgOKWTchhQ6DjeXQSipjFxEqMNNg1MEg9AcRs661mGnCkLJZeQiQh3uu50iGITmMHiMHQA4HGZHKKmksKEDhID9OimDQWgOg6eMArYIEUo2c1qEWItTAIPQHEZuK6M71CIkRr4mQpnMhFqMk2VSA4PQHMZ3jdIsRXO0EleNfFGEMpjxtZjhaQDcIir5MAjNIYWMvpYEAN5nE/FyEqEkEQ0f4AAc7E8N+p133qEOt3bt2gcffLD1LXv37jW7nJlGCht3BlMLHGBAKImkkMwbX4u9uEVU8tHTp08n//Xyyy+Xl5cPHz4cAK6//vqW23v27Gl2OTONKS1CPIMCoSSSIsadytsCW4SpcFjX6MKFC6+99lqKOrTMW1HwuiMlNEnTVMI6jNuQQocHeyKULEpCpRmK4YweXeJx+ncK/PwuVlVVff3111deeaX+7QsvvOD1egsKCh544AFCjjnXUNO06urqLYerr69PecHTmfGr6XWcD1uECCWHFDR6EaHO5sU19cn38xu5aNGiyZMn672gl1xyyY033uj3+9euXXveeecVFRXNnj27zcdHo9F77rnH4XC0vvHKK6+87bbb2rx/LBbTNK2l0dk9RfcLrJuORCIGv65qk2ONiSNeV5IkAOA4zuDCoGMhhMTj8eNcfSLjRaPRI285KNAuyvharHFKvEkw/nWt5uh35FicTifDnKD77VAQapr24osvPvbYY/q3vXr10r8YMWLEnDlzPvzww2MFodfrXbBgwfjx49tZJoqiXC5XNw9CQdEcWbzH4zH4dal8tiEePeJ1MQithhBC07Tb7Ta7IOgwR1ScuKS4sh3G12KSTzevSRj/uhaUxD/Coa7RTz75JBaLTZky5eh7hEIhl8uVrNdDcGj5kRldo37sGkUoOXCAI5MceiMXLlx49dVX8zyvf/vwww+PGDEiJyfnm2++ee655959913zSpiBpLBiM3Ydrs7mZuSYSjRC0d26RY5Q10kh2ZlvN/518crGVNkAACAASURBVACKVGABQFEUnuevu+66lltVVX3wwQej0WjPnj3ffffdiRMnmlfCDCSFFHeRCVWIoimbi5EjivFLGBHKMFJY8fcxoUVI22iapZSEavy08wzGAgDLsi+99FLrW+++++67777bpCJlPlNW0+s4n00KYRAi1FXG76/WQt9AH4MwiXCLNRNIIXMmXgPusoZQkohBc8YIAdfUpwAGoQmksNGnt7TAw5gQSgICctSEzaF0nI8Vcb5MUmEQGk0RVKAO7SJvPJxyhlDXSRHF5mTMmnTGeW14OZtcGIRGk0KK8Rv1tsApZwh1nRQybZgfAHgfbi6TZBiERpNCpo2xA44uIJQMZi0i1OlT3sx69YyEQWg0yYwzzFpwPlYMYhVCqEvMbRHiMTJJh0FoNBPXTgAA78UWIUJdJYZNmykDhybL4OVsMmEQGk004yTCFqyL0WSiyZpZBUAoA5hyJG8LzsvKUYVouC170mAQGs3cKgQ4XwahLpNC5mwXrKNoyuZk5KhqVgEyDwah0aSwaavpdZyXxQEGhLpCCpu8PROug0ouDEKjmb7DGeezYYsQoa4wd7IMAHA+7NdJJgxCY+kbUnjMbBHyuC0FQl2gKUQVNZvTzK0+sUWYXBiEhpKjCmNnKMbMU5BwWwqEukIKyzYPC6YeZcZ78XI2mTAIDSWGZN5v8skPeC2JUFdIQWvUYuwaTR4MQkNJpq6d0PF+VgxiECLUSaLZw/yAl7PJhkFoKHNX0+twfyaEusISLUI8RiapMAgNZe7+arpD243iYlyEOkU07zzRFniwaHJhEBpKNHXHbR3NUoydkWN4OYlQZ5h7gIyOdTKaqOEWUcmCQWgoKazwZl9Lgr6CAocJEeoUMSRzZneNAgWcF+fLJA0GoaFMX4erw2FChDpNCpq8S6IO19QnEQahoawShH4bLkJCqDP0PTHMnvsNOHE0qTAIjWOFDSl0vM8mYdcoQh0nRRTWYfKeGDrcIiqJMAiNI4Vlzmszd0MKHefHKoRQZ0hBCwwQAgDgGGEyYRAaxwqr6XU8jhEi1Cli2MwDmFrjfHiMTNJgEBpHtM61pM+Gs0YR6gQpKPNZ1qnFeDmbHBiExrHChhQ63o/D7Ah1hmiBRYQ63o8j/UmDQWgcK2xIoWN4GihQBDzhGqGOkYKWmPgNuEVUUmEQGkeywNETLXCYEKFOEK2xAgpatoiKYi1OAgxC44hBq3SqgH45if0qCHWQFJJ5vyX6dQCAxwXBSYJBaBzJCjsz/Rfvx/kyCHWYFFJM3y64BYd7JSYJBqFBiEakiFWWTwDOvUao45SESjEUw1vlY5P32yScOJoMVnlHM54UVmxulqItsJweAPRjXHCMEKGOEC0zU0bH+bBrNDkwCA1iqaEFAOBwBQVCHSSFFEvVYlwHlSwYhAax1EwZwGtJhDrOii1CHCNMBgxCg1hqpgwA8D4WZ40i1CGWWgEFALwfR/qTA4PQIJIFzqZvzeZiVTzhGqGOkEJW2WhUx+FIf5JgEBpEtMz+aodQwPlxr0KEOkBolvkszuxS/IzhaJqllBhuEdVVGIQGEYPWmiwDAPYsmxiUzC4FQmlDbLbKjtstcE19UmAQGkQKKZYaIwT9nPpmrEIItZfVxggB58skCQahIch/T+W1Ej6LwyBEqJ2UmErRFlpNr8MVFElhrTc1U0lRhXEwNGuV1fQ63GUNofYTLHMSYWucj8X5Ml2HQWgEKShbahGhDoMQofazzpG8rfG4e34yYBAawXJTRgEAgM/CMUKE2ktslqwYhHg5mwwYhEaQQjJnsSmjgFUIoY6w5uUs7pWYFBiERhBD1tpfTcfwuAgJofYSm60YhHg5mxQYhEaQLLZFYQvejztTINQuosVW0+sYngYKFAEvZ7sEg9AIovWWH+n4LBxpR6hdrFuLfTYJB/u7BoPQCGKTFeebAQCfZZOxRYjQiRCNyBGF81lupB/0WW94Ods1GISpR0AKW/VaEgcYEGoHJaJZ6mDt1nBnjK7DIEw5MSSzToZiLFmF/JyE+24jdCJSSLFmpw4A8Fk2AYOwazAIU04Myny25cbYdXyWDYMQoROSrR2E2CLsIgzClBObLNovCtg1ilD7SEHFsrXYnmUTm/EYmS7BIEw5MSjZsy1ahTgfq8RUohKzC4KQpclh1bJBiGOEXYdBmHLWXIero2iKdTFyBBchIXQ8UlCx4CJCHedj5aiCl7NdgUGYctZch9sC514jdEJWbhFSNGXzsFIYB/s7D4Mw5Sx4qnVrnJ/FNfUIHZ8UtO5kGQDg/TahCYcJOw+DMOWsHoQ+FieOInQcqqgRlbBOxuyCHJM9G4cJuwSDMLWUhEoIYR3WrUKc3yZiECJ0bGKTZMHTY1rjcPp312AQppbFBwgBwJ5tE7FTBaFjE5pkLsvSQWjPxqWEXYJBmFpiULbs2gkdn20Tm7BFiNAxCU2SzdotQt6PSwm7BIMwtay8ml7H6QdQ4NRrhI5BaJT4bGsHIY4Rdg0rCEJtbW3L93l5eW63GwBCodDq1asDgcDw4cPNK17aE4OSlWfKAADNUoyDlsIWPTERIdOJTbKniDe7FMfD+7FrtEvYNWvWnHnmmT179tS/f/TRR6dNm/b9999PmTJl1KhR27ZtGz169OLFi80tZfoSm2VXkcPsUpyAPZsTmjAIEWqb0CQF/C6zS3E8DE9TLCXHVJvLuvPyrIwFgN69e//000+tb73nnntuu+22u+++OxQKDRgw4Ouvvx43bpxJJUxvgrXXTui4bJvQKHnLnWYXBCErEpoki0+Wgf/uOGpzWf2y25poAFBV9ccff9y9e7emaQAQDoc/++yzK664AgB8Pt/UqVPfeustk4uZtiy+iFDHZ7G4GhehNikxlaIoxm716RS442hXsABQV1c3a9asffv2FRYWLlmyRBRFmqaLi4v1e5SUlKxZs+ZYj5dl+csvv6yrq2t9Y58+fQYPHtzm/TVN0zSNoqx4OF/SEZUoMYV10/oVhjVpmsb52Xi1ZOVCdiuEEL2amF0QBAAQbxDsAc767wjnZ4WmblSL2/+O0PSJL2LYESNG1NXV2Ww2WZZnz5590003/eUvf2FZtuXBPM8nEoljPV6SpM8//3zDhg2tb5wyZUplZWWb9xcEgWGYbhKEYpNs87CiJJpdkOORJIlyk0SjKAiC2WVBAACEEEEQWNbqfXHdRORgzOajRVG02SzdtcN4qHhDovvU4va/I3a7/YRZyDoch/qUbTbbnDlzfvGLXxQUFIiiGIlEPB4PANTX1xcWFh7r8S6X6/777x8/fnw7S69pmtPp7CZBKNVE7QHO6bT02BvLskyhrTYYtHg5uw9CCCEE3w6LaIrGnXkOh8Nh8XfEna/U7+tGtVhV1ST+sofl5Pbt2/Py8vLz83v16rV06VL9xqVLl55yyinJer1uRWiU7QFLbyuj43ysHMFjXBBqg9Ak2bPToBbbA5zQiCP9ncQ+8sgjiUSioqJix44dCxYsmD9/PkVRd9xxxy233BIOh1evXl1TU3PJJZeYXc60JDRKaRGEFE1xXpsYTI/YRshIYpMUGOgxuxQn5sAg7AJ2/PjxS5Ys+eSTT/Lz8z/44IOxY8cCwA033BAIBD766KP8/Pxvv/3W5bL0GhrLEhqlwCCv2aVoF3vAli6xjZCR9H4dBayeMYydplhKjio2N44udxg7evTo0aNHH/2Diy666KKLLjK+QJkkjaKFz+aEJpx7jdDhCIhBmffbFNHqQQj/7R3FIOwEqy+OSWtCo2TPSY8gxAEGhI4mhWXWQdNcenxO2rOxFndSerzB6UgRVKKSdNnxCA9jQuhoQpOcFjNldPYAJzRiv05nYBCmitCQNs1BOLTdKAYhQocRGiU+TUY3QA9CrMWdgkGYKmk0QAgAPF5LInQUoUmy+HmireEAR6dhEKaK0JhOnSqcm9UkTRW6y/5MCLWH0JBOl7OOHE5owCDsDAzCVEmjmTIAABTYc7hEg6V3g0PIYIl6yZFr6ZMIW+N8rBxVNAV3xugwDMJUERrTY0OKFo5cPlGHl5MI/SxRLzpy06YWUzTF+XHWW2dgEKZKeo0RAoAjF1uECP1MiamgQXoty8Nhws7BIEwJohEpJKfRMDsAOHJ5oR6rEEKHJOpFe146XcsCrqDoLAzClBCbZZuHpZh0OmTDnssl6rFFiNAh6TVAqMMVFJ2DQZgSQqPkSKOZMgCAY4QIHS5RL6ZdLcau0c7BIEyJtBsgBACbiwEa5KhidkEQsoREQ3q2CDEIOw6DMCWERplPqymjOkcul8BhQoQAIN2mjOrwMKbOwSBMiXRsEYLeO4rDhAgBAAGhQbKnWxAydppmKTmC/Todg0GYEok60ZGXZp0qgDtTIPRfUlimOZq1p8em+a3Z8XK24zAIU4BAIg0nywBWIYT+Kx2njOqceXwcZ711EAZh8gnNks3JMHz6/W1xjBAhXToOEOocebgOqsPS78Pa+hJ1Ujr2i4I+RtggAW5ViLq9RIOUrkGYyycOYhB2DAZh8qXpACEAMDzN2mkxhDtToO4uUSc6ctKyFjvycICjwzAIky9RJzrTbWemFo5cHntHEUrUp22LMIcTm2U8g6JDMAiTL1EvpukwOwA4cjkBLydR90Y0Ijal1TFqrVAMxWfhasKOwSBMvvjBdO0aBQBHPh/HAQbUvQkNEuez0bZ0/Xh05HGJOqzFHZCu77RlqaKmCBrvT6dzJ1pzFdrjBwSzS4GQmWIHRFeh3exSdJ4jj8cg7JB0OmorLSTqREcuB+l07MRhnAX2WC0GYWfIGjQI0CCQRhGaRdIkQlCCoETCEkRkiMoQkkhChbgCcQVEFQQVEioBgJgMknbEk/EAh01Z8nNAUUAD+DgKALJ4YCjwcuBiKRcLbhtk8eDjKD8Hfg6yeCqbhxw75NopT7pekpkpvl9wFqZrpw4AOHL5SFXc7FKkEwzCJEvUic607RcFAM7LAoAcUWwe/N84kkqgNk72RGBfjNTGYV+MHEzAvhipE+BggkRkyOEhx04F7JDFUdk8+Hnwc1QvD3g48NjAx9F2BlwsOFiwM8DT4GQpAHCywLfawIQQEovF3G5365cOSkAIqATCMiEEghIoGkRkiCokJkNUgaAIQYlURyEoQZOgNUtQl4AGgagEcu1UkQvyHVShAwqdVE83FDmpEjeUuikXvsltiR8QAkN8Zpei85x53MFVzWaXIp1gPUiyeH26LiJs4Sywxw4Ifo/7xHfNXBEZtofI9jDZGYadYVIVIbujUBsnOXaqzA09XFSRE3q6qVG5UOykcx2QZ6dyUtmX5v/vvI0c+xG9DSfofEgoUCeQ/XGoS5DaOOyPk+UHoDauVUdhT5Q4GCjzUOUeqpcHKrxUby9V6YOerrTt0EiS2H6x5Ow0rsXYNdpRGIRJljgoBgZ7zS5FlzgL7fH9or+yGwVhUIKNTWRTkGxqJpuDZEsQGkXSx0tV+qgKL4zJoy7rTZe5ocRNpd38CQcLpW6q1A1tRma9AFURsjtCdkZgVT15Zae2NUgiMvTzU/391IAsaqAfBmVT5R6q+2SjphAxKKf15azNzQKAHFX0L9AJ4Z8pyRL1oiNtFxHqXIV8dG/C7FKkVk2MfN9A1jWSdY3wQxNpFMjALGpQFjUgizqvhO7rgxJ3t/joz7VDrp0anXvY7xqWYWuQ/BQkPzWT57aQTc3QJJLB2dSwADU8QI3MoQZlpd8FQfslDoqOAEcx6f3+O/L4RL2EQdhO+GdKKgJCGh7meQRngf3gqqDZpUiymAKr68l3dWRlHVldTxRCRuZQIwLUVb2pYQG6l7dbxF47eW1wUi51Uqt0DEqwoZGsbyTLDpDHN2q7o2RwFjU6jzo5lxqbT5W6M+qPF0vzmTI6Zy4fPyh6y51mFyQ9YBAmk9Assem53XZrzkI+fkAAcsLhJ6trEmHZAW3pfvLNQfJTMxkaoMbkUZf3ph4/hc6wz+5U83MwvpAaX3jojxaVYV0jWVlPllSR36xUaYo6LZ8aV0CdUUT196f9JUX8gOgsSOO1EzpHPi4l7AAMwmSK7xed6bz8SMfaGdbBCM2SPTv9+ngTCnx1gHxeo31eS3aGydh8anwh/ejJ9KhcKg2PlrMotw1OL6BOL6BgMAAwO8Pkm4Pkq/3k0Y1aRCZnFdFnFVJn90jXlmJsv1A4NtvsUnSVs8C+f3mj2aVIGxiEyRSrFVxFaR+E8N/5MmkUhFtD5INq8vE+bUUdGZFDTSiinxxLj86l2PRunKeHCi9V4aWuqgQAqI6Sz2vJF7Xk92tUP0dN7kFNKaHHF1B8+lyFxA8IzoK07xp1FdljNbgguL0wCJMpViuk+5RRnauQj+0Xsgd6zC7I8agEvj1I3tujvbOHCCpM6Un9agD91kTajUvIzVPipq7pQ13TBwgwGxrJx/vI/WvVH5vIxGJ6Wik1pYQOWDtiVEFTYmoaXQIeC++3aRrBiaPthH+jZIrVCiWT88wuRRI4C+zNW6Jml6JtKoEvasmbVdrbe7RiJzWtlH5jAj0skJa9cBmMAhgWoIYFqLuG0o0ifFCtvVdNbv5WHp1HXVBGzyij8xxmF7EtMb05mBH/Ta5Ce6xW8PfpRuugOg2DMGk0WRNDcpoe3XIEZ6G95ssGs0txGALw7UHyyk5tSZVW5qYuLKe/m8r28mTEJ1amC/BwVSV9VSUkFOajfdqbVeTu1fJJudSlFfQF5bTXSi34+H4hA4b5dRiE7YdBmDTx/aIzN+2XH+mc+bzQKGkKoVnzf51dEfLiNu3lHcTBwuW9Mf/SmIOFGWX0jDJIKMwHe7VXdpJfr5An96CvqqQn97BE1YnVCmm93XZrriJ7aGfM7FKkBwzCpInWCq4iS3b3dBzNUo48PlYjeEpN+40EFd6s0p7fqm1qJpf1ppdMZIZj/2emcLBwYTl9YTk0i8xru7T716nXLYOrK6lr+9HmXuVEqhN5I/0mFiCJXEX22mU4cbRdMAiTJlYrOIusPROgI9wljkh13JQg3BEmz2zWXtqujcihfjWA/mUpzeHkzwyVxcP1/enr+9M/Bck/t2pj3lWGB6gbBtC/LKGNbyBqCkkcFF3FGdIidBbwiXqJqMQSbW1rww+YpIllUIsQADw9HQZvtEYAPt5Hzv1YOfV9haNh5TT243PYC8oxBbuFAX7qryczey+1XVVJz/9B6/Wa8tAGrdHYFeGxGsGRx6fvebxHoG00n2WL47L6dsAWYZIQiO8XXOm/M1MLT4ljn1HzZQQVXtyuLdioOVm4ZSD99iQaF793TzwDl/emL+9Nr28kT2zSKl+XLyqnfz2Y7uMzok0TqY67SzLnWhb01YQZNOqZOkYHYb1ALd6rURSln8HmZCkXCx4bZPMQsFOWmj/WIWJQpm1UJi3ZcRbYpbCsJFTWkcJQahbh6c3ak5vU0Xn0c6cx4wqwDwcBAAwLUAvHMX9JMM9sVk//t3JaPv27ofQRm4MnXXRvwlfhSulLGMxVZI/XCjDS7HJ0lqhCo0iaRIjIEJOhWSIAcHYx7Uv23HyjP7hjKqxtBACiH9IdV7SoAhEZmkRoFIigQr6D6uGCQidV4oYyN9XLA5U+qsJLWbx/LFYruIoz6loSKHAXO6L7Eik6j6kuAQs2qn/fok0tpb+Ywvb3YwSiI+U74I8jmDuHMAu3aRd/ofb2wv8OY1q2PE26SHWi+MycFD25KVxF6bHR2v44bAuRnRGyO0L2RKE6Sg4moCZGEirk2CGbpzw2cNsgi6MA4LR8SPsgLHORZ0+lqWNszCuqUC+QvTHYHyd7olAVIZ/WkB1h2BsjPV3UoCxqcDYMC1AjcyirnR2akf0P7hJHtDr5QVgvwCM/qAu3apdU0GvOZ9N0R0pkGCcLvxpAz+lHL96h/c9ytdAJc0ckPw4VQZVCsjM/c0Y3wKobrcUVWNdI1jeS9Y1kYzPZHCQ8A5VeqtJHlbmps4qgp4sudEKhk/IbtSrbWl15PAM9XFQPFxxx8IGswfYw2dRM1jeSf2zRbmggFAWjc+mx+dRp+dSoXPPbi7GaRGCIz+RCJJunp6N+XSiJTxiS4JEf1Gc3a5dW0BtmsMUWu5pBVmaj4f/1oa+spF/ZqV23XC1xwZ9HMSfnJe1fKLpXcBU7KDqj/if1jdaksMJ5Tf6o3x0hyw6Sbw6QFXVkR5gMzKL0Js3VlfSALCrb7MsPawXhsdhoGOCnBvipmeWHbqmOkpX1ZPkBcvN32vYQOSWPOrOIntyDGhYw5xSY8J5E2dQCM145hdwljl3vHkjKUwkq/G2TNv9H9Zcl9LoZrNUa9ChdMBRc2Zu+tBf94nZt5ufqqFzqwVF0v2T0q0dNWiyUap4SZ2RP3JQ9kOsF+LRG+6yGfLmfiCo5LZ8+rYCa3ZceGjC/6XKE9AjCo5W4qRL3oVwMSfDVfu3zWnLpF1pIIr/oSf+ylDq7mHYa9cuJzTJRSQZs1HsEezZH1K5eTmoEXtmp/e/32sgc6qspbFI+s1A3x9Iwuy99eW/6qZ+08R8oM8ro+0YwBV1LsUh1Ind4pnXqAIC31GFwEG5oIu/tIf+u1raFyJlF9MQi6s4hyblYSZ10DcLWfBxMLaWnlgIAVEXIB3vJUz9pVy9VJxTTF5ZT55WkfDPDyJ64tywzT4L2lDgi1fHAoE7Wom8OkttXqAwFi89kTs23dE1IKVXU5Igix1UlrqqCqiRURdA0SdNkoggqEFBFlag/319RFJZtbvmWsdMUTTEcTbEU66AZjmHsNOtgWAfNOlnWxdhcTIb16bWHnYHfDKZn9aH/vF4d/Kb8m8HM7YPoTp/3FN2bKJ+WaZ06AOApc1Z/Up/qVyEAq+rIkirtzd2EpmBaKTVvNHNqPpUuazIzIQhbK/dQvxpA/WoA3SzCe9Xaqzu1G79RJ/egL6+gzu2ZqpWy4d1xT1kGdqoAgKfUGa7qTBDWxMhvV2nLD5CHRtOXVnSDD2kCUkQRGiWxWRKDshRSxGZZiipSUJajCsVQNhfLuhjWybB2hnUyjJ1mbDTrpPlsG0UBw9Mt238QQkRRtNt/nnulJDQgRNWDM66JQVlNaIqgKnFViatyXJWjis3F2tws52M5L8tncbyP5bM4e7aNz+assGFs6mTxMP9k5ob+9G9XaQOWKI+OoaeVdrieC01SRnbqAICnxBmrSaRuf5mtIbJ4h7Z4B+EYmFlOvTOJGZKdfv9vmRaELbJ4uLqSvrqSbhLhjSrtkR+1/1muXt6bnt2XHpDsRnpkT6LsvPzkPqdF+Cpdu96s7dBDZA0WbNQe/kG9vj/9/OmsYR3URlJFLXFQjB0UEnVSok5M1EtCg8g6GD7A2bNsvN9mz+F8vV2ch+V8Npub6dAlGCEkFou53R2ZrEtAjilyRBFDihSRxWY5ujfR8ENYbJLFZsnmZu05nCOXd+Rxzny7s4Dn/Wm7YvcYKrzUWxOZz2vJzd+qf9+iPX4K09vbgWoe2h7zpWaZkOkYO23P5mK1grtnMi/WozL83y5t0TZtV5hc1pt+YwIzIif98q9FJn5KHS6bhzn96Dn96O0h8sJ2bdKHam8vzOlPX5ikvbs0hcT3C56k/pNZh6fEITTKcky1udrV5bT8ALnhG7XEDd9NZTv0SWRxYlCO7kvE9gmx/UKsRpCiijOPd+bzjnw+d4TfkcvZczjGxAkAFNjcrM3NOguP/AnRiBiUhQYpUS/FD4rNm6PxA4ImEWch7yqyu4od7h52V6E9M7ajnFBEbZjBPr5JO+U95eaBzO+GtLenNLg95q/MqKX0rXnKnJE98WQF4Y9N5JnN2mu7tHEF9F1D6XN70Gya9H8eR+YHYYtKH/XnUcwfR8D71dozm7U7Vqpz+jHX96fzu/bvEduXcOTxmbohJkVT3l7O0I5YztAT9I4GJfjtSvXjfeSxMfSF5Wn/19AkLbI3Ea6KR6sTkeo40cDd0+EutueO8JedZ3fkcOlyditFU/Zszp7N+fv8fKMSV2O1QqxWCO+K1X7dKDZJzkK7p8ThKXV6y518Vhq3F2003DGYvrgXdet32tC3lOdOa8eKQwKh7dHSX2TCkdpt8pQ6gttihad16Uk0Au9Va09s1LaGYE5/+ofMWv7UjYJQZ6P1E9HoTc3kb5u0AUvkC8vp33RhM8PwnowdINT5Kl3B7dHjB+GSKu3W77Tzy6hNF7Lpu0+eKmrhXbHQzlhoZyy+X3QW2b2lztyRvl7nF6Z1NhyNdTK+3i5f70NtIFXSYvsSkT2Jhg2hXe/upxnKV+HyVrh8Fa40PWi6p4t6ayLz3h7tyqXquT2peaOZ4yzNjh8UaY7OyAFCnbfMuffTzs+XSSjw4nbt0Y1agIdbB9IXlGfMtuQ/63ZB2GJgFvXsacyfRjFPb9bG/VsZV0DfM4we1vET7yK7E4FBnlSU0CL8le4D31Yf66cHE3DjN+qWIHljAjM2DeeFEo1E9yaat0SD26KxWsFT4vBWuMrOK/CUOjN7jklrDEd7e7m8vVzFAAAgNEj61cDeT+qAAn8fd1Zft7+Pm3Wm2VboU0vpM4rou1apg99UnjmVOa+k7Tc0uD2awf2iAODI5ZWEKkcUm6djH/gRGZ7ZrC3YqJ6cSy8al8kTv7tvEOpy7PCH4fQdg+m/b9F++Yk6PEDdN4Ie2ZFR38ieeKbOlNG5Cu1KQhOb5aNbRa/u1G5foc7qQ79yJtvpaeumUBJq8+Zo00+R5q0R3mfL6ucumZznLXdm4LVux9lzOHsOl39yFgAk6sTgtljdmuCO12udBXz2QE/WAE8a7SbotcHTpzIX9SLXLlNf30U9MbaNpmFwWzRjDuNtGwWeUmd4dwdWE4ZleHyj9rdN6uQe9Ge/YJM+wdBqunsQ6pws3DaIvqE//fxWbfqn6sgc6oGR9OB2TAIWm2VNJfZAxnaqAABQ4OvtCu6I5p+U1XJbgwA3fKNuDpIPJrMdum4wlxxRBpqjoQAAIABJREFUGn4MN/4QjlTHfRWu7IHe8qkFpu8+ZWWOPN6Rxxeelk1UEtoZa9oU2fzPaqAgMMSbM8TrKXGmxVjpGYXUhhns3avVIW8q/zidmdzj50ITjYR3xSsvLjaxeAbwljoi7QvCuAKPb9IWbFTP6UF/m1lT3o4DPwJ+xjNw0wB6dl/6uc3a2R8pk4rpuSPpcs/x/g+C26IpOpzBUvx93KFtsZYg/Hc1mbNcvbw39dIZbFocHKjE1IYNofr1oViNkD3AUzg2e8Cskkyd35QiFEP5+7j9fdy9zi+M1QiNG8PbX6tVRTVniC93uM/6x/i5WHjiFGZ6KZm9TD2nBzX/ZMbFAgBE9wp8li2TDlBrk6+P+4TroGQNnt+q/WmddnoBtew81pgzIC0iw9/+TrAzcOsgenZf+tGN2knvKP+vD33vMCbrGHvCNm+JZg/M5AFCnb/SVf3xQSAQVeA3K9XPashrZzGnWf7sQE0hTRvDdd8Hw7viWQPcxeNz/H3d3WfkL3VcxXZXsb1kcl78oNiwPrR18V4AyB3pzxvlt/iUk7OKqA0z2Fu+VUe8rfxrPHNyHhXcluEDhDpPiUNolo8zTPj2bu2u1Vq5B/49mRne8akS6Q6DsG1uG/xhOD2nH/3HtWq/JfL/DmOu73/k+BHRSHBHtNeMo5ZuZRx7gLO52ZXrY1fu4E/Pp9bPYD3WnkQZrREOrmiqXx9yFzvyTvL3vbInw2P7L/mc+XzJ5LySyXmR6kT9muCGBbucBXz+yVk5Q32WveDw2uCF8cxbu7Xpnyo3DmCmbAj17gZVmKIpf6WreWs0b9SRo6FrG8jtK9SQBH8by5xdbNF3LdUwCI8n3wHPnMr8agD96xXqM5u1Bacc9o8S3ZvgfTaugxOx0pFKYFO+d/PHwb9cUniBhdcIqpJWvzZ04NsmJa7mj/YP/03vzNtCxZo8JQ5PiaN8akHTpsiBlc273tmfN8JfeGq2I8/s83WOYUYZPSaP+s2HiRGNanaO04SjGQzn7+s5IgjrEnDP9+qHe7W5I5jZfelusBHiMWX+h3jXDcyi/nMu++9qctM36uBs6tGT6TIPBQDNW6JZ/TJ/gHBfjFy5VM1xeH4t7BlTWmR2cdqWqJf2L2+sWxv09XKVTcn393GnxSSODEMxVGCINzDEKzbLB1Y0/fh0lSOPLzotkD3IY8E9wYuc1EPe6OoKz8nvK4+PYS6psO4VXlJk9XPv+fAgEAAKFA2e+kn783r16kp6y0xb+q79TRYMwvY6r4SaVMw+ulE76V3l1oHMb4fQzVujpedm8sIJAHh7t3bDN+qtg5jfDXH8UG0L74q3rMK2iOD2WO1XDZHqRP7JWdgEtAg+y1Z6bn7J2XmNP4ZrljZUvXeg8PTsgpOzGbu1wqZxfWjyRcUfudnLvlT/U0P+dgrjztx/H95vs7mYaE1iPWu/8Vu10AFfn4fHoh3CAkBVVdXSpUtDodDw4cPHjx8PAJs3b964cWPLnaZMmeJ0ZuYxQx3CM3D3UPryCurW77Qxb4h/qxV9vTL2z5JQ4Ncr1U9ryLuTWP0c8Jxhvvr1IYsEIdFIw/pwzZf1mkKKxuf0+38llh2U6rYohsoZ5ssZ5ovuTdR81bD6s635o7OKx+dYZL1K/ICoipqn1DmCgjXns7d8q458R3n1zPTePPr47L3dz38cftRlezQj9kFMInbjxo3jx48/77zzcnJyHnvssYkTJ/7zn/98++23Fy5cOGLECP1OZ511FgZhixI39fYk5j+fR9Y7Hc8u0/56MtPF3UotaGMzufQLdUg2tfb8n7dMyxnm3bBgV8WMQnO7uTSF1K1u3vdFA+ezlf4iP6ufB3tBLc7d09H3ip5is1zzVcPah7fnDPH2mJBr+urbhvWhnKE+/Z/HxcI/xzGv7dLO/Y9y11DmtkHW68ntGgLwwjbttYOO60NNP12dn8EN385hy8vLq6urXS4XANx44429e/f+05/+BACTJk165plnzC6edfWqi/U7w7vXDkPekh8cxczqmzk155nN2n1r1EdOZq6uPOya0Z7N2bNtoR0xfx9zRkY1hRz4rqnmywZXkb3PZT285Xhxlk74LFuv6YUlk/JqlzdueHxXVj93z4m5Js6mqV8f6ntZj9a3XNyLPjmXunyp+lmNtmgcm5cpF7hbQ2TOcjWhwLMXeMTH99sVFWzpsP7XQLTL5dJTEABcLhdFURRFAcDevXtfeumlL774QlEUU0toRZqsNW4KFw73zhvNfHYu+/xW7Yx/K1uCxOxydVWTCDM+U5/fqi3/JXtECuryR2ftX95kfMGISvYvb/r+z9tC22P9Z5UMuLYUUzBNsS6mZHLeqHv6OPL4H56s2vbKPqFBMr4YwW1RiqaOPpmozEN9NYUdHqBGvKN8UpP2NVpUYe5a7fT3lQvL6e+mssPzGX8/d8OGsNnlspzDOuvvuuuuCy+8sKCgwG63q6r65ZdfrlixgmXZpUuXZmdnt/l4QRCef/75Tz/9tPWNp5xyysSJE9u8vyiKLMvqWZu+mn6IuHrwhFNFUe3jgi/Ohr9vo05/X7u+L/x2EEmvXTclSQIAQsjXB6lZy6kLy8iLpxKOVkWxjTv7hjh2f3wwtC9iN+pQAqKRxnWRms8anQV87ysLXMV2ABDbLFym0E+ot9kyuveKgrzTvYHR7oPLm9c/vjNrgKtoQoDzGTd2uPeLuoLT/aLU9j/S7wfDuDxq9lfkwjJy/3DC0SCKIsdZeq+Ao31bR924gurjhZVTtCKnJksAAFlDXQe+bsoanvbXke1/RziOO2Hi/HzV/+CDD65evVrvDv31r3/90UcfLVy48McffwwEAvPmzTveU9A0dZT2lC99Na6PBIb/vPSIpuD6vmTVeeSHZurkD+hv69Ls15c1+MM66upl1DOnkIdGkuNsPUbb6LwxvgPLmo0oFoHmjdFNj1c3rA33vrSw8uoiPQVRxmB4umhCYPBvSlkXu/GJPXs/alDiqgGvG98vJg5K2UOP18M/Pp+s+iWpilLjP6a3hgwoVDIFJfjVCurKZdQfh2mvn6EVtUo9X6UzUSeJzbJ5pbOiQ5dgjz322AsvvLB06dJAIND6xwzDnHPOOcuXLz/W4+12+6xZs/S5pu0hyzLP82mdlHJMje4W+l9desSJ5OU8vDsZ3t6tXbVc+0VP6qGTjrkxm6VsC1PXLKeK3fSGC5jcdgRNz/H5a/6yrdd5TEePdOmQcFW86v0Dmqz1ml7UHRZrtkYIURSF59PhvycZeB4qpjp7nplX/Z+6jY9VF58RKDo9kNJjQPZ8U180LsfuPMG/eyEP75wN/9iiTfwE7h7ovG04nxYfW29Uabev0KaVUpsuZHxtNZlyh/tCGxM9J6Z3tZIkKYl1hAaAf/zjH0888cQnn3xSVHRouTQhP3eOL1++vKKiIlmvlwEa1oeyB3iYY7Sbzi+jN13I2mgY+Kb8yk7N4LJ1CAF4YpM24T/UNZXk3UntSkEAsLmY3BG+2mWNKSpVol7avKh668t7C0/NHv7r3t0tBbstzsP2vrBoyM3l0erEmoe2160JQmpG6MSg3LQlUnBK1onvCgAA1/Wjl/+SXVzFnPuxUhOz9KhhVYRM+Y9y/1rtjQnMU2PbTkEAyBvpr18TNLZoVseuW7duzpw5Y8aMufPOO/Wb5s6de+ONN1ZUVGRlZa1cuXLPnj3PPfecuaW0lLo1wZKz845zB68NnhzLXNmbnrNcfWGb9tRY5v+3d+dxUZ33wsCfc2bfd2ZnGxCQHVRENhVjNGhd2iRNP73NTdo0bZJP2txr3ry9adMluU1u27TJbZrlzdpma2K2aoyNggmoICrKIiCyMzA7s+9nzjnvH5MSaqOCDMzAPN+Pf8ycmeH8ZDjnd87z/J7nyU68edxHPeT3juMhHDRvJ3XznDZcXSftenpEWS2J7YCwiA+fOGKxnnepN0pz/k0LxwUmIZaMkfvvqe5R/+gBk6F5OmOXQqCL8bjVib9bFBXieS2bskqAHK0PPzPCKfso8mQF5dtZCTcCL0yA3/cQT/bg/1lI+c/Ca9xO81LZJEF6JwNczUqpi10wqkajeeedd2ZvkkqlTz/9dHt7u8vluv/++xsaGlgs+Pv6gs8YDNnDwlXXPjgrUpCzu6nP9BFVByN356E/KaawE2IYMSBI8Kc+4tHz+P8ppjxQgOLYvGv2mBK6coN4+H1D3h2pMQmJxEnDienJJpu0hF/2UDaNs6zKjaBY42ewi+/PtHa6Bt+e4qiZ6TsUrBgVZzkveV1DvtIHs+b7QSoKflaK7khF7mjG3xkhnq+iqDmJcqHWZCDva8Wz+cjpXdSrLxv3BQSklAtNbY6sm+GJ/QvI7FbQ61BXV/erX/1q7n2EXq83OkhjITuNo4HX9VwNS71JOvePGP1gXzt+wkz+Zh16a2acryV77OQPTuIUBLxUQ4muNxatGp1vRRyJk+efHE7dKpOWCBYUEAls3a6xj81sBSN9h4ItT5aOsasgSdLn83G5sE0YEBHS0DI99blNVipI3ZpCXdgVEh4mzv9mSPeN6+l19ng8PB4PAIAR4PEu4ple/JEyyj15cZ6oetRDPniaOG8j/7Ae/VraPM4tmA/veOJS2b4sumC5FifPfCMxsdQ3KUFzeOhvZoAAKouCUBAKHaWyKTQOhcal0nhUhpAW/ZeYc4UELCHnoC/rlvktZq1kgzc3UY6byB+14X/qI35XQVkni8N/z4OBR8/jfx4kHiunfC93occvQkGyb1X3vzouyOZe9w2ca9g39rGJxEH2reoEmbkNSigoFdFslsorRBOfWjqeuKTeKFXVSK57ReXxQ2aBjrPAXmcaCh4pRW/OQH54Ev/LIPFsFWVNPKZkc4XB/3ThLw4QPy6gvLGRMt/1sWkcinytaPIzW+buBF2CKhLEQ3Ys5MTCLizsiUR8OObHiTCBh4jsb6pjPqXwUidCupimu1mJkEgkiJM4iYeJ6P8waA+7x/xhFxa0YxFfhCGmsxUMdgqDo2Jy1EyWlJEIqVHfZFXVSK5vZbsaBXJ2N/W1QeLrjXiVHHm0HF2yjkOCBK8NEj87S9yoQXr20mI1XwYvjSVfJ+p7aSz/7vT5LlTvnQxM/N3iN4XStqfIyoSJ8OVCCYvGoej2KlU1kvFPzB1PDGq3yOQVIoQyvz8a44lpe5+n5D9iU/eXJ0Q+a6D+ZZDYdQTfqkF+vYaiXKqBeSEcvHCR+HUn3qBFO/dQr7uFVr1Jeu5/BrX1skUt/56jiB/3TgV8hqDfFPKbQwFLiMRJppjOENHofCqNT2NK6FwthcJAKQyUyop918lS/wpQGsLRsK7eNEpEyKAt7DeHAuagtdM1dsiM+SK8VDYvjcXP4PAz2HFZZDVoDzv6PJkPX/8FFIqAO1eh38xEn+4lqg5GGlLRn5agOv4iJgESgAPjxCMdBI8G/rY19peuadvleIjofWEeudA7EdA3WT0TAW29NO+O1PmezqCkxZLRc2/XeicC459aJo/ZNPXSlLWiOVZUGY5PG1qmC+/JiOE5FAHg9mx0bzr660688H3srlz0wSKKeDGb9sMEeHmAeLyTKJEgjTdRC0QLOnboPGpKuXCq2Za+QxGrCOeBBH5T0DXs94z73eN+zBvhqFhcNZOXxpavE7FkdBp3SXPT8ugjjPhxz7jfMx5wjfi8+gBbyRSt4gpzuLw01pJNAH3pzUmGhJ627Wr1onPnCoOnLhB/6sM3q9AHi9DyWKcoggQfjRNPdBFhHDy6Bt2ZesVLh+vrI/wSCUb+ZnQN+7K+oeKlXfGqmMRJe5/H0DIdtIfVdVJFpWhRB4ota7CP8JrcY/7JRqtHH1BuEMsrRFdpKIsEcf2n1ukL7sJ7Mhii629Pu3qPlN5HPnae+GCM+EEeet/q2M/C7wqD5/uJP/YRJWLwi/KYXdGGnNj5J4dK/yNrIb+Z+e7RcdHrvOR1DXqpbKpAx+als3lpbHbKvNv8YttHuDwS4WxEhPSM+R0DXueAN+gIi3J4kgKeKI+3qLeJ0z3usY/Npf+pu+7+ia/kxcDLA8QfLhBKNrg7F70lE114ZaktCN4YIv7UR6SwwINF6K60a1wpLDQRAgBIYDnrHD9s5mewFZVirpY1812EPZHolzXd5WIrmYpKsbSYn4BrtCYUmAjnyG8OGY9P27pcHBVTUijgZ7LZCsYXf10k8FtCjoueyWM2ST4vbbt8gQ2AczntjnjI3/cQbw8Te9LRu3LQ6OJlC3TGSr5wkfhgjNiZiu4rRAvFMT52Jpuszku+gh+kL2L3BAm8k4HpHre9zxN2YcJcniiHK1zFXeDgq2RPhLOFXZi9zzN9weMZ9fN1bGmJQFLAj3lGxLyR878byrsj9Sp3PAuBk+ATPfH/LhInTORNWvTmTGSLCp3vOimuMDisJ94bI48ZiJ2p6A/y0Mq5HYcxSIQAAACIMDHVMu3o8/gMQRqXChAQ8eMAAfx0Nl/HkZUK4JK5cwQT4bwQEdLR77H3eTxj/qADo7EpCAWJ+HEqm8LPYKs3SjmqGMzMN/fTri0IXh4gXrlE0FFwmw7dlYbkz7MNkwSg205+MErsHyXDBPheDnrHKnSR1nojCbLnmVFpqUBVI7n2u+f3o4FnImDrdNq63SgNlRTyxfk8XmrM2vBgIvwKkSDu6PVYO13uYb9wFUdWLhTl8WIzIpsE/a9OsBWMtJsWfTF6axB8MEa8N0qcspClEqROiRSLkVIJkspF/rUdESPAqIfsd5KnLOQJM9k9TdYp0d3pyNfT0SvNKPGVYpUIZ5A4GXJiAAAKk0JjU2AhzHzBRHjdosV3JEFSGGhsO5nme9olAThhIt8bJf42TlIQUK1AquVIsQTJ5iNfOfOiLQh6HWTnNNluJY8ZCB4N2ZWG3JKJrpUt+lCzoC3c9b8jRfdlxGpJLL8pZOlw2s67UBoiLRFIiwVsRez7TmEivJpIAJ/udls6nD5jUFokSFkj5Kezr/9cTIKh9w1+Q7Dw3oylLOvwR8AJM3nSRHTZQbednPKRYgZIYX2RDv0R4AyTjhBQc5AcAVgnQ6vkSKUc4VzXgR/zRAgtEEyECWghp91eB3nSTLaayQsO8pKLRBEgYXyZDh0hYPCTbCrIEyLFEmSNFNmsQtK4S3r9aD7j0H9qLfhBOlN6/eeBsCdiPee0nnViPlxWKpCVC2NyL34ly3sc4WKjsijyCpG8QhRyYtZzrqH9BiJMyMoFKeXC+V7vkAQ59I4haA/n352+xMWNbCrYqka2qr+ocCMBMPmBNUhiBAAAsKhASAdSJhLT/koIghZFvgjJFyHfz/3iqT0EHCHS+Y8JnUQMoGQhrLieieVrRSQOep4dzb87fb6TWuAhYrrHbT3n9IwHJIX8jF1KgY6z7NqBVloinMEQ0jSbpZrNUp8haOlw9jw3RuNSZCUCSRGfJbv2N+03hUY+NAIE5N+VFtsCmeuAAKBkAyV7uf1xQRD0L8QMIE68dSwU60UoDel5djRtu1xRIbpmJsNDhOOi19bpdA74+Dq2fJ0o747U5VsKvmIT4QyOipmhUmTsULhGfLYu14XnxlAaIsrjCTI5vHT2ZZVLJE56xv2Wc67pbrf2BplygxgOdIMgKBmklAvZcsbIRybjyWnNRqkwh3tZPysRIb2TAc+o3zHg9Uz4+elsabEg6xb1YoxwX2IrPxF+AQECHUeg4+j2AK8h6LzoMZ92DO2fAgjCENJoPCqJk5EgHrSEmVK6aDWv/P9mU9nL/tuFIAiaO66GVXRfhq3bbe1wDn9gpPOpNC4VpaN4iAg7sbA3wpYz+OlsZZU4747UuExsskiSJhHOQABXzeT+Y61zzBMJuTDME0EoCIVJYUrocOkDCIKSmbSILy3ikwTpN4UifpzACJSOMoQ0uoC2UhdHS75E+M9oPGoiTLUHQRCUUBAUWdSyz4Sycu5tIQiCIOg6wEQIQRAEJTWYCCEIgqCkBhMhBEEQlNRgIoQgCIKSGiyYhKClFowEMSICAPBhfoIk/rExhBEYAIAkyUAgwA5/sdQJj/7lpKMsKouKUgAAXDoHWXbTWEFQooKJEIKuE4ZjrrDHHfK4Q25P2OsOez1hrzfs82F+P+YPRILesC8QCYRxzI8FApFAhMCjmY9JZdBQGgCATWNTkC9aZWY2AgAIgkDRL7Z7wt6ZPfqxAE7iMxvZNBYFoXBobDqFxqKy2DQWi8pkUZlsGptL53BpHC6dw6Vz+HQel84VMHhCpoBDW5SlxCBoWYOJEIKuyB3y2AJ2s8/qCDqtfpsj6LIF7M6g0xF02YMODMf4DB6fwefTuXwGj0fn8uhcLp0jYirZNDaLyuTRuSwqk06hR1MUFaVGU9fVdzr31SeiedEb9mEEFogEfWF/EA8GsKA/EvCEvb6wz+K3eTGfJ+R1hz2ukNsVcmM4JmDwhUyhlCUWMvkSlljCEklYYilLImNLJCxx9I4TgpIKTIQQBBxBl9FrMnrNRq/Z6DObfBaLz2b2WRgUhoQtlrNlYpZQxpakCTTliiIhUyhg8MVMIZfOiW/YbBoL/HPb6TVhOOYMuaKJ3Bl02wL2KY+py9I7HbCbfTZn0Clg8OWcFAVHJuekKLny6D8FRw4TJLSCwUQIJZdgJDjhntK7p8Zc+kmPYdJjmPIYKShFzVVGT/o54qy61A1ytkzOSWFSY7+gaHzRKDQZWypjS7/yVYIk7AGH0Wcx+ywmn3XAPvT5xEmj12z1T8vYEg1PpeGrUvnqVL4mja+RsmO9pjkExQlMhNBKFiHwMdfEiHN8xDk25poYc+kdQaeGp07lq7V8dZVmnYanUvOU87qpWsFQBJWyJVK2pFCWN3t7hMBNPvOkxzDhnhp1TjRPtI67J0ORUJpAmylMyxCkZQrTdKJ0AYMfr8ghaCFgIoRWlGAkNOgYuWQfHnSMDNpH9O5JJVeeKUzXiTJ2Zm3LEKYqOHIUgfWW80NFKRqeSsNTrVetmdnoCXvHXBOjzokR53iLvnXYOcaisrJFGdnizGyRLleSdaX7TghKNDARQssbQZJjrok+20CfbeDi9OCU15Qu0K4S6/KlObuzt2cI0xgUerxjXJl4dG6hbHWhbPXMFqPXPOQYuWQfOTR85PdnngMA5Iqz86TZ+dLcXEk2LFiFEhZMhNDyE4wE+2yXuq29F6wX+2wDEpYoT5qzWrJq16rtOmEGLOuIl2gna422MvrU4rddnB7ssw38ueevg44RBSelQJZXKMsrSslXcFLiGyoEzQYTIbQ8BCPBHmt/p7mn03JhyDGWLcoolK3em7Pjkap9fAYv3tFBXyGFLU1hS2u1lQAAnMSH7KM91v6Tk6efO/8aDaWVyAtKUgrKFEUwKUJxBxMhlLgIkuifvnTW2HXW1DloH14l1pXKi75X/G+rpTmwwXN5oSCUHElWjiTrG7k7AQAT7qluy4UOU9eLXa8zKPRyRfEaRUm5ohhe00BxARMhlHBsAfupqbNnjOfPmbpTONI1ipJ/K7i5SJa/8gYzJK1UvjqVr96RdSMAYMyl7zB1Hhn97Lftz6TyNetUpRWq8jzJKhSBMyFDSwQmQighECQ5YB9snTzTNnXG7LeuVZZWadb9eO3dIqYw3qFBiytdoE0XaL+eszNC4BesfaeN539/+jmrf7pCVbZBvW6dqgxW2UCLbakToT8SmLQbkX/Ur9NQGpPKYFAYdAoNjuVKQhiOnTf3HJ88dXKynUfnblCv+9Ha7+dLc+HdQBKiopQSeWGJvPD7Jd+x+m1tU2f/PnLst+3P5ElX1WjWV2nWwfEYSSgYCYWJsDfswwncHwkAADIEqfRY94wsdSK0BexPdb4w8xQjsGAkFMRDYTzsDftYVCaHxubRuTwGLzqLlYgpFLOEUpZEyhbLWBIxS7TEAUOLIYyH2w0dzRNtpwxnM4Rp1ZqKZ7b+j4qriHdcUKKQsaVfy972textwUjwjPH8icn2l7vf1PCUtdoNG1OrlFx5vAOEYiAYCVn8Vqt/2uaftged0wG7M+iOTorrCXs9Ya8P89FQGp1C59G5KIJG2wZ+WfNQzP8AEJIkF/L5urq6X/3qV3V1dXN8v9fr5XA4yBVGNAciQW/Y6w373GGvM+hyhlyOoNMecFr907bAtNVv82EBBSdFyZWreQoVVxmdH0TBSYE3ENctHA4DAOj0pag9wXDslKHjs4kT7YaO6ExmNZr18OLmMnOfdDup4CTeab7QPNHaom9L4Ug3pVZvTquWL1XFqcfj4fFgIc/1c4c9evfUhGtyymua8hiNXrPJZ/ZjATlHJmVLZGypmCmUsMQCBk/IEAgYfB6DG11B5Urn9th+I4nVRxhdROYqDSAhPBydGXnKY5zyGtoNHXr3lDPk0vLU6YJUnShdJ0zXiTIk8NyaSAiS6DB1NY41t06e0YnSN6fV/GjN9+F0XNC8UBBKuaK4XFH847V3d1l6j40fv+vwf2j5qvq0uk1p1SKmIN4BQl8K4+FR58SQc3TYMTbmmhh1TYQiIS1fncpXa3iqKs06FVeh4KQkzkVwYiXCa2JQ6NGu9dkbg5HQuEs/6poYdo6d6f9w0DFCQ6k54qwcSVaeZFWuJBuec+NlwD50dPTzpvHjcrZsS3rd90tuh9co0AKhCFoqLyyVF/5ozd0dps7GseaXu9/Il+ZuzdhYrVkPS4vjIkLgw87Ri9ODF6cHB6aHJj0GLV+TJUrPFKZvUK9NF2gTfIr2ZZYIvxKTyogOUZrZYvZZLtmHL9qH3un/6OL0oJglypfmFMryCmWrUwUauLT3YrP5p4+Mff7pyLEwjm3N2PTHGx7X8FTxDgpaaagopUJVXqEqD0ZCJyfbj4x+/tSZF6o1FdsyNxfLC+Bhvtg8YW+Ptf+Ctf+Ctf+SfVjJledJVuVJVu3OvilTmEaj0OId4DyshET4r+ScFDknJTrVE0GS466JC7aLXZa+N3p0P7lrAAAUCklEQVTfC0aCxSkFxSkFZYqiy+4soQUK4+ETk+2HR5ou2gbrUjfsq7ivQJYLz0fQYmNSGfXptfXptY6g8+hY8x87XvKGfdsy67dlboZlNbHlDfvOm3s6LT2d5l6j15QryS5KWX174TeX+1yyiVUsswSsftt584VOc885c3cYD5crStYqS9cqS5J2vFpMimWGHCOHhhubxlqyxZnbM+trtJVw5pfrBotlFm7QMXJ4uKlprEUnSr9Jt6VWW7nAgvtkLpaJEPgFW/8Z4/mzxk69eypfmlumKCpOKciR6ChI3Ob1je03knSJcDaj13zW1HnW2HnO1C3nyCpUZRWq8gJZXlLVoC4kEfowf+NY86Gho66Q+ybdlm2Z9XKOLNYBJh2YCGMFIyInJ9sPDR+9OD24Jb22Qbc1S5RxfT8qCROh2WdtN3ScNp47Z+pO5WvWKkvXKEtWS3NoaEK0I8JEGHsESfTaBtoNHe2GDpPPslZZWqWpqFCWcemceIe26K4vEfZY+z4eOnJisn2tsnSHbmuZohgu8hcrMBHGnMVv+2T46CfDTSKmYEfW1vq0WjaNNa+fkCSJkCDJi9ODJ6fa2ybP2AL2ClX5elX5GmVJAtYbwkS4uKITXbZOne40X8iVZFdr1ldp1q3gG515JUJXyP3pyLGPh48CAHbobrgxc3MCHiHLHUyEi4QgybPG8x8PHzln6q7VVu7I2rpamjPHz67sRBjGwx2m7hOTp05OnhYxBRs066rU63IlqxL56nYljyNMBFKWeEfW1h1ZW0N4ODqlxZ97/irnyGq062u1lWlJWV9DkOR5c/fHQ0dOG89VayoerLh39nKsELQsoAiyTlW2TlXmCDr/PnLs121PUVHqzqytN2Rs5NNXbJK7Ch/mP2XoOK5vO204ly3WVWsqvp1/c3KWF8E7wmsjSKLb0tuiP3Vc38aisWq1lbXaylViXbzjio2r3xHa/NOHR5oODR/l0jg7srbekLFxWdeGLQvwjnBpkIDssvQeGjrSOnVmvWpNg+6GUkXhlYqcV9IdoTvkOTnZ3qJv67L0FqWsrtFWVmsqll3TDmwajRsSkBenB1sm2j6fOAkAqEvdUKutzJOuWtYjBL4yEWI4dnLq9CfDjf22S5vSqhuybsgRZ13hB0AxBhPhEvOEvUdHmw8NH/Fh/u2Z9V9Z87UCEqEj6DyuP9Wsb+23XVqrLK3VVq5Xr1m+17UwESaEIcdIi76teaLVhwVqtOvrUjcUyVYvx3LTyxLhwPTQ30ePNY21ZArTbtJtqUutggMhlhhMhPEyYB86PNzUNN6SLcrclllfq62cmadm+SZCi9/WMtHWom8dcY5XqMrrUjesU5atgPl3YCJMLBPuyZaJtmZ9q8Vnq9ZWVGvWlyuKYr5KyOKJJsLpsKNprOXI2OcRInJjxqYbMzcrlmo6Y+gyMBHGF4ZjJybbPx09dsF6cYNm3Q3pdeWKYp/Xt7wS4bhLf3yy/bi+zeg1V2nW1Wo3rFEUL6/ZXq4OJsIEZfJZjutPndCfGnSMrFGWVKrXVqrXCBkJPRewxW87Nnq8Wd9q9JnqUqu2ZmzMl+bGO6hkBxNhgnAEXU1jLU3jzUavpUq57gZdXVFKfiK3+uAk3mPpa5s6e2KyHSOwKk1FjWZ9ibwgkWO+bjARJjp3yNM6daZt6sxZU2cqX7NetaZCVbZKnJU4tcjjLv3JqdPH9W1THlOlak2dZkOFtjyOk0RAs8FEmGgMXtPhgcY281lbwF6tqajWVpTJE6jVxxF0thvOnTZ0nDF2KrnySvXaam1Ftigz3nEtLpgIlw2MiHRbek8bzrUbOqaDjjJ5UbmiuERekMrXLH0wPszfZblw2nC+3dCBk3iVZl21Zn2JvADHcLBU6xFCcwETYQKKnnYNXtMJ/amTk6cHHSPFKQUVqrJyRbGWr176eLxhX7e177yp+6yp0+qfLlcUV6jK16nKpCzx0gcTFzARLku2gL3D2HnO3N1pvhAmsAJpboEsL0+SnS3WsajMRdqpPeDotV28YL3Ybe0dc+nzJKvWKksrVOWZwrSZ9yzlwrzQXMBEmIAuO+26w54OY1e78VyHsZMAZJm8MF+aWyDLyxCmLlLLCkGSkx7DxenBPttAj7Xf6DXlSVeVphSWKYpzJVkrsvHz6uCA+mVJyhLfmLn5xszNAACzz3rB2t9rG/h84sSIc0LOkWUK0zIEaWkCjZqnVHOV853/Kcob9hm8pgn35JhLP+wYu+QYDkfC+bKc1dKcu0v/PV+Ss5K6yiEojvh03qa06k1p1QCAKY+xy9LbY+374NIhs8+aKUzLFmVmCFNT+RoNTyVjS64jSxEkYfJZDB6T3jM16pwYdU0MOUYEDH6uJDu68uIqcRYVhX0ZMQMTYRzIOTI5R1afXgsAwEl8zKUfc06MuMabxlqmPMYpr5GO0iVssYwl4TN4fAaPSWHwGFwAAAIQDo3txXwAAD8W8GE+b9jnCLqmA3arf5ogCSVXkcpXpwu02zI33y++KzkniYCgpaTmKdU85U26LQAAPxYYcowMOcZGnRPNE62THqMz6JKxJRKWSMgUChl8Lp3DpXFQFKWiVJIkcRIHAAQjoTAedoXc7pDHGXKZfTZn0ClmidU8hZanzhCmbkytyhZn8uiwkWCxwEQYZxSEohOm64Tp9bM2ukMea2Da5p92hTzukCeEhzwhLwCABKTBa+LSOAAAFo0pYqp5dI6IKZSwxFKWmM9YTuXdELTysGmsopT8opT8mS0Yjln8NnvQ4Qi6XCG3N+zzhr0ESWIEhiJotB2VQaULGHwtX82n84RMQQpbKmGJ4Q3fUoKJMBFFbwR1wvR4BwJB0ILQKLToLWO8A4GuJum6WCEIgiBoNpgIIQiCoKQGEyEEQRCU1GAihCAIgpIaTIQQBEFQUqMCAE6fPv3aa68BAG6//faKioroC+++++7hw4flcvl9992n0cRhSjAIgiAIWgJoV1fXli1bcnJy8vLytm7deu7cOQDA888//9BDD9XX1wcCgaqqKr/fH+84IQiCIGhRUILBYHV19SOPPFJRUeF0OhsbG3fv3n3bbbc99dRTN9988/bt29955x0Gg1FaWvqVn3/ttdc2bdqUnp4+x/2Fw2E6nZ6Ec40mLBzHAQAUChy9m0AwDIOzvyaUcDjMYCz7xWxXkth+I2hbW9umTZuiTzZu3Nja2mo2m4eHh2c2btq0qbW1NVb7gyAIgqCEQjWZTFKpNPpEJpMZjUaTycRgMGYm9pbJZN3d3Vf6vM/ne/jhh2d+QtSOHTu+9a1vffX7J0d8Le/DO8LEQRAEAABFYdlUoiBJEsdxPxXO+pRAIpFIAH4jiYG983soXxwIBObYjsVkMq95fqMyGIzoQjwAgFAoxGKxmEwmhmEEQUQ/HAwGWawrLoZAp9NvuOGG/Pz82Ruzs7OZzK9eWggTSalr62EiTByRSAQAQIUHecIgSTIUCl3pCILiIhAIXOU0CC0lhkCEMJgYhs3xGJnLVT5Vo9Ho9froE71er1arVSoVSZIGgyFaLBrdeKXP02i0jRs3zn09QgqHzy6pgYkwccD1CBMNSZKkz8eG6xEmEtzjYcdu9Tto4VAUjWE7Frpnz54333yTJEmSJN966629e/fy+fz6+vo33ngDAOB2uw8ePLh3795Y7Q+CIAiCEgr1nnvueffdd6urq1EUdTqdr776KgDg8ccfb2hoOHny5MDAQG1tbW1tbbzjhCAIgqBFQRWLxR0dHW1tbSRJVlZWRpvI1qxZMzAw0N7enpKScqWBExAEQRC0AqAAABqNVltbW1dXN7ujSCgU3njjjTHPgnq93uv1xvZnQgthtVqtVmu8o4C+5PV6Z7rtoUSA4/jg4GC8o4D+ycjIyEyZ58ItddH8/ffff+LEiSXeKXQVzz333AsvvBDvKKAvtbS0/PjHP453FNCXnE7nli1b4h0F9E/27t07Pj4eq58GR49BEARBSQ0mQgiCICipLXQYdSQS6erqmvv7HQ5HT08Pm81e4H6hWBkfH0dRtLm5Od6BQF+4cOGC3W6H30jicLvdkUgEfiMJJRgMnj592mAwXPOda9as4XA4V38PQpLkQqJ5++23n3nmmbnPS+L1eplMJpzHJHEEg0EAAJzHJHFEIpFgMMiFA+oTBkmSbrdbIBDEOxDoS263m8vlzmVM/csvv5yVlXX19yw0EUIQBEHQsgb7CCEIgqCkBhMhBEEQlNRgIoQgCIKSGkyEEARBUFJb3OpNDMOOHTvmdrs3bdp02eK9AIBTp06NjY2VlZWtWrVqUcOAZrjd7ra2No/HU1hYmJOTM/uljo6OmccymSw1NXXJo0tGvb290cJdAACPx7vsWBgeHj5z5oxWq62qqopHdMlo9jcCABAIBDM1h1NTUyaTaealkpKSOa4NC10fs9k8OTmp0+mEQmF0i8fjaWpqQlG0vr7+skEROI5/9tlndru9pqZGqVTOa0eLWDUaCoU2b95MkmRaWlpjY2NTU1NRUdHMqw888MCBAwdqa2s/+eST3/72t9/5zncWKQxoxsWLFysqKtauXSuTyY4cOXLPPfc8+uijM69SKJSamhoajQYAaGhogLN8LY3Vq1dzOJzocV5cXPy73/1u5qV333333nvvbWhoaG1t3bx58/PPPx+/MJPIHXfcMTk5GX187ty5vXv3vvjii9GnDz/88JtvvpmdnR19+uGHH8JRLounoKBgdHQUw7D9+/fv2rULAGA0GisrKwsLCyORyNDQUGtrq0wmi76ZIIjt27fbbLbVq1cfPnz4wIEDGzZsmMfOyEXz+uuvFxcXYxhGkuRPf/rTb3zjGzMvjYyMsNnsqakpkiQ//fRTlUoVfRu0qBwOh8FgiD4+f/48giBWq3XmVRRFZz+FlkZeXt6JEyf+dTuO4xkZGR9++CFJkhaLhcfj9fX1LXl0Sc3v9wuFwuPHj89s+a//+q99+/bFMaSk0t/fj2FYTk7ORx99FN3y0EMP3XrrrdHHu3bt+vnPfz7z5o8//jgzM9Pv95Mk+eSTT0bvweZuEfsIDxw4sGfPnujY+ZtvvvngwYPkP+4+Dx06VFlZqVKpAABbtmwJBoNnz55dvEigKKFQONNioNVqSZKc3QQEAOjo6GhtbXW5XPGILnn19/e3tLRctgZIT0+P2WxuaGgAAMhkso0bNx48eDBOASap/fv3p6SkXNYobbPZmpqaLl26FK+okkdubu5lU68cOHDglltuiT6O5pSZlw4ePLhz504WixV96bPPPvN4PHPf1yImwqmpKY1GE32s1WpDoZDNZvvXl1AUValUU1NTixcJ9K8ef/zxLVu2zHwLAACxWPyb3/zmgQceSE9Pf//99+MYW1Jhs9mvv/76z372s4yMjD/84Q8z2w0Gg0KhiLZUAwA0Gs1cZpOCYuiVV1757ne/iyDIzBYEQbq6up544omqqqodO3aEQqE4hpeELssps7PG7JfUajWCIPM6XhaxWCYSicz0JEcfYBgWfYph2OypcahU6sxL0BJ45ZVX9u/ff9l6WAaDIXrafeutt+68885t27Zdc4I+aOHa2tqiv/ZTp05t3LixoaEhWi+DYdjsQgwqlRrD1degaxoZGWltbX377bdnb/z5z3/+2GOPAQDcbveGDRv++Mc/7tu3L04BJqPLcsrsI2L2SwiCoCg6r5yyiHeECoViprXHYrFQKJSUlJToU6VSObshyGKxzLfIB7pub7755iOPPHL06FGtVjt7+8zNxze/+c1gMAhXIl0aM7/29evX63S6zs7O6NPoMTLTm2CxWKJdCdDSeOmll7Zv337ZeWnmy+Lz+bt3755daA0tgdmJ47IjYna6sdvtkUhkXjllERNhbW1tY2Nj9HFjY2NVVRWVSg2FQsFgsLa29sSJE9EOqt7eXo/HU1ZWtniRQDPee++9Bx988NNPP50p0w+Hw4FAYPZ7on3UarU6HgEmL7vdrtfr1Wp1JBLx+XwFBQUoikb7zjEMa25urq2tjXeMyQLH8TfeeOPOO++c2eLxeAiCmP2erq6u2T0L0BKora09evRo9HFjY2NdXR0AIBAIYBgWTTfRC8fGxsaCggKJRDL3n7yIwyfsdnthYeGePXt0Ot1jjz32xhtvbN++/Yc//KHP5/vLX/5y4403Igiya9euZ5999mtf+9p///d/L1IY0Ize3t6SkpJt27bl5+dHt9x9992vv/768ePH77333g8++KCsrMzlcr344otf//rXn3766fhGmwy6urp++ctfVlZWkiT55z//OS0t7dChQ3/9619/8pOfjI2NRY+a+++///Dhwy6Xq6WlJd7xJotDhw5997vf1ev10VvAcDjMYDA6Ojp+8YtflJaWSqXSY8eOnTx5sqOj47JmFSiGXnjhhdHR0RdffHHjxo3Z2dn33HOPy+Wqrq7et29fJBJ5+umn29vbc3Jyampqdu7ced999xUVFdXU1JSWlj7++ONPPvnkt771rbnva3FXnzAYDK+++qrT6dyzZ090VEdzczOGYdFK0ZdeemlkZKSiouKWW26Z3SMNLRKj0XhZ5eHOnTsnJyeNRmNFRcX+/ftHR0e5XG5NTc2WLVviFWRS8Xg8+/fv7+/vp9FoZWVle/fuRVF0cHCwra0tOrL2/fffP3nyZFpa2l133QVX8VwybW1tfr+/vr4++pQgiGefffbWW289e/ZsW1ub1+vNysq67bbbRCJRfONc2Q4ePGg0Gmee7t69OyUlpbe396233kJR9Nvf/nZ0SpAPPvggIyOjtLTUarW+8sorFouloaFh8+bN89oXXIYJgiAISmpwrlEIgiAoqcFECEEQBCU1mAghCIKgpPb/ATYFXSnxXD8+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.4" }, "kernelspec": { "name": "julia-1.5", "display_name": "Julia 1.5.4", "language": "julia" } }, "nbformat": 4 }