{ "cells": [ { "cell_type": "markdown", "source": [ "# Introduction to periodic problems and plane-wave discretisations" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "In this example we want to show how DFTK can be used to solve simple one-dimensional\n", "periodic problems. Along the lines this notebook serves as a concise introduction into\n", "the underlying theory and jargon for solving periodic problems using plane-wave\n", "discretizations." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Periodicity and lattices\n", "A periodic problem is characterized by being invariant to certain translations.\n", "For example the ``\\sin`` function is periodic with periodicity ``2π``, i.e.\n", "$$\n", " \\sin(x) = \\sin(x + 2πm) \\quad ∀ m ∈ \\mathbb{Z},\n", "$$\n", "This is nothing else than saying that any translation by an integer multiple of ``2π``\n", "keeps the ``\\sin`` function invariant. In a more formal one can use the\n", "translation operator ``T_{-2πm}`` to write this as:\n", "$$\n", " T_{-2πm} \\, \\sin(x) = \\sin(x + 2πm) = \\sin(x).\n", "$$\n", "\n", "Whenever such periodicity exists one can exploit it to save computational work.\n", "Consider a problem in which we want to find a function ``f : \\mathbb{R} → \\mathbb{R}``,\n", "but *a priori* the solution is known to be periodic with periodicity ``a``. As a consequence\n", "of said periodicity it is sufficient to determine the values of ``f`` for all ``x`` from the\n", "interval ``[-a/2, a/2)`` to uniquely define ``f`` on the full real axis. Naturally exploiting\n", "periodicity in a computational procedure thus greatly reduces the required amount of work.\n", "\n", "Let us introduce some jargon: The periodicity of our problem implies that we may define\n", "a **lattice**\n", "```\n", " -3a/2 -a/2 +a/2 +3a/2\n", " ... |---------|---------|---------| ...\n", " a a a\n", "```\n", "with lattice constant ``a``. Each cell of the lattice is an identical periodic image of\n", "any of its neighbors. For finding ``f`` it is thus sufficient to consider only the\n", "problem inside a **unit cell** ``[-a/2, a/2)``. In passing we note that the definition\n", "of the unit cell is itself only unique up to translations. A choice ``[0, a)``,\n", "for example, would have done just as well.\n", "\n", "## Periodic operators and the Bloch transform\n", "Not only functions, but also operators can feature periodicity.\n", "Consider for example the **free-electron Hamiltonian**\n", "$$\n", " H = -\\frac12 Δ.\n", "$$\n", "In free-electron model (which gives rise to this Hamiltonian) electron motion is only\n", "by their own kinetic energy. As this model features no potential which could make one point\n", "in space more preferred than another, we would expect this model to be periodic.\n", "If an operator is periodic with respect to a lattice such as the one defined above,\n", "than it commutes with all lattice translations. For the free-electron case ``H``\n", "one can easily show exactly that, i.e.\n", "$$\n", " T_{ma} H = H T_{ma} \\quad ∀ m ∈ \\mathbb{Z}.\n", "$$\n", "We note in passing that the free-electron model is actually very special in the sense that\n", "the choice of ``a`` is completely arbitrary here. In other words ``H`` is periodic\n", "with respect to any translation. In general, however, periodicity is only\n", "attained with respect to a finite number of translations ``a`` and we will take this\n", "viewpoint here.\n", "\n", "**Bloch's theorem** now tells us that for periodic operators,\n", "the solutions to the eigenproblem\n", "$$\n", " H ψ_{kn} = ε_{kn} ψ_{kn}\n", "$$\n", "satisfy a factorization\n", "$$\n", " ψ_{kn}(x) = e^{i k⋅x} u_{kn}(x)\n", "$$\n", "into a plane wave ``e^{i k⋅x}`` and a lattice-periodic function\n", "$$\n", " T_{ma} u_{kn}(x) = u_{kn}(x - ma) = u_{kn}(x) \\quad ∀ m ∈ \\mathbb{Z}.\n", "$$\n", "In this ``n`` is a labeling integer index and ``k`` is a real number,\n", "whose details will be clarified in the next section.\n", "The index ``n`` is sometimes also called the **band index** and\n", "functions ``ψ_{kn}`` satisfying this factorization are also known as\n", "**Bloch functions** or **Bloch states**.\n", "\n", "Consider the application of ``2H = -Δ = - \\frac{d^2}{d x^2}``\n", "to such a Bloch wave. First we notice for any function ``f``\n", "$$\n", " -i∇ \\left( e^{i k⋅x} f \\right)\n", " = -i\\frac{d}{dx} \\left( e^{i k⋅x} f \\right)\n", " = k e^{i k⋅x} f -i∇ e^{i k⋅x} f = (-i∇ + k) e^{i k⋅x} f.\n", "$$\n", "Using this result twice one shows that applying ``-Δ`` yields\n", "$$\n", "\\begin{aligned}\n", " -\\Delta \\left(e^{i k⋅x} u_{kn}(x)\\right)\n", " &= -i∇ ⋅ \\left[-i∇ \\left(u_{kn}(x) e^{i k⋅x} \\right) \\right] \\\\\n", " &= -i∇ ⋅ \\left[(-i∇ + k) u_{kn}(x) e^{i k⋅x} \\right] \\\\\n", " &= (-i∇ + k)^2 u_{kn}(x) e^{i k⋅x} \\\\\n", " &= e^{i k⋅x} 2H_k u_{kn}(x),\n", "\\end{aligned}\n", "$$\n", "where we defined\n", "$$\n", " H_k = \\frac12 (-i∇ + k)^2.\n", "$$\n", "The action of this operator on a function ``u_{kn}`` is given by\n", "$$\n", " H_k u_{kn} = e^{-i k⋅x} H e^{i k⋅x} u_{kn},\n", "$$\n", "which in particular implies that\n", "$$\n", " H_k u_{kn} = ε_{kn} u_{kn} \\quad ⇔ \\quad H (e^{i k⋅x} u_{kn}) = ε_{kn} (e^{i k⋅x} u_{kn}).\n", "$$\n", "To seek the eigenpairs of ``H`` we may thus equivalently\n", "find the eigenpairs of *all* ``H_k``.\n", "The point of this is that the eigenfunctions ``u_{kn}`` of ``H_k``\n", "are periodic (unlike the eigenfunctions ``ψ_{kn}`` of ``H``).\n", "In contrast to ``ψ_{kn}`` the functions ``u_{kn}`` can thus be fully\n", "represented considering the eigenproblem only on the unit cell.\n", "\n", "A detailed mathematical analysis shows that the transformation from ``H``\n", "to the set of all ``H_k`` for a suitable set of values for ``k`` (details below)\n", "is actually a unitary transformation, the so-called **Bloch transform**.\n", "This transform brings the Hamiltonian into the symmetry-adapted basis for\n", "translational symmetry, which are exactly the Bloch functions.\n", "Similar to the case of choosing a symmetry-adapted basis for other kinds of symmetries\n", "(like the point group symmetry in molecules), the Bloch transform also makes\n", "the Hamiltonian ``H`` block-diagonal[^1]:\n", "$$\n", " T_B H T_B^{-1} ⟶ \\left( \\begin{array}{cccc} H_1&&&0 \\\\ &H_2\\\\&&H_3\\\\0&&&\\ddots \\end{array} \\right)\n", "$$\n", "with each block ``H_k`` taking care of one value ``k``.\n", "This block-diagonal structure under the basis of Bloch functions lets us\n", "completely describe the spectrum of ``H`` by looking only at the spectrum\n", "of all ``H_k`` blocks.\n", "\n", "[^1]: Notice that block-diagonal is a bit an abuse of terms here, since the Hamiltonian\n", " is not a matrix but an operator and the number of blocks is essentially infinite.\n", " The mathematically precise term is that the Bloch transform reveals the fibers\n", " of the Hamiltonian.\n", "\n", "## The Brillouin zone\n", "\n", "We now consider the parameter ``k`` of the Hamiltonian blocks in detail.\n", "\n", "- As discussed ``k`` is a real number. It turns out, however, that some of\n", " these ``k∈\\mathbb{R}`` give rise to operators related by unitary transformations\n", " (again due to translational symmetry).\n", "- Since such operators have the same eigenspectrum, only one version needs to be considered.\n", "- The smallest subset from which ``k`` is chosen is the **Brillouin zone** (BZ).\n", "\n", "- The BZ is the unit cell of the **reciprocal lattice**, which may be constructed from\n", " the **real-space lattice** by a Fourier transform.\n", "- In our simple 1D case the reciprocal lattice is just\n", " ```\n", " ... |--------|--------|--------| ...\n", " 2π/a 2π/a 2π/a\n", " ```\n", " i.e. like the real-space lattice, but just with a different lattice constant\n", " ``b = 2π / a``.\n", "- The BZ in our example is thus ``B = [-π/a, π/a)``. The members of ``B``\n", " are typically called ``k``-points.\n", "\n", "## Discretization and plane-wave basis sets\n", "\n", "With what we discussed so far the strategy to find all eigenpairs of a periodic\n", "Hamiltonian ``H`` thus reduces to finding the eigenpairs of all ``H_k`` with ``k ∈ B``.\n", "This requires *two* discretisations:\n", "\n", " - ``B`` is infinite (and not countable). To discretize we first only pick a finite number\n", " of ``k``-points. Usually this **``k``-point sampling** is done by picking ``k``-points\n", " along a regular grid inside the BZ, the **``k``-grid**.\n", " - Each ``H_k`` is still an infinite-dimensional operator.\n", " Following a standard Ritz-Galerkin ansatz we project the operator into a finite basis\n", " and diagonalize the resulting matrix.\n", "\n", "For the second step multiple types of bases are used in practice (finite differences,\n", "finite elements, Gaussians, ...). In DFTK we currently support only plane-wave\n", "discretizations.\n", "\n", "For our 1D example normalized plane waves are defined as the functions\n", "$$\n", "e_{G}(x) = \\frac{e^{i G x}}{\\sqrt{a}} \\qquad G \\in b\\mathbb{Z}\n", "$$\n", "and typically one forms basis sets from these by specifying a\n", "**kinetic energy cutoff** ``E_\\text{cut}``:\n", "$$\n", "\\left\\{ e_{G} \\, \\big| \\, (G + k)^2 \\leq 2E_\\text{cut} \\right\\}\n", "$$\n", "\n", "## Correspondence of theory to DFTK code\n", "\n", "Before solving a few example problems numerically in DFTK, a short overview\n", "of the correspondence of the introduced quantities to data structures inside DFTK.\n", "\n", "- ``H`` is represented by a `Hamiltonian` object and variables for hamiltonians are usually called `ham`.\n", "- ``H_k`` by a `HamiltonianBlock` and variables are `hamk`.\n", "- ``ψ_{kn}`` is usually just called `ψ`.\n", " ``u_{kn}`` is not stored (in favor of ``ψ_{kn}``).\n", "- ``ε_{kn}`` is called `eigenvalues`.\n", "- ``k``-points are represented by `Kpoint` and respective variables called `kpt`.\n", "- The basis of plane waves is managed by `PlaneWaveBasis` and variables usually just called `basis`.\n", "\n", "## Solving the free-electron Hamiltonian\n", "\n", "One typical approach to get physical insight into a Hamiltonian ``H`` is to plot\n", "a so-called **band structure**, that is the eigenvalues of ``H_k`` versus ``k``.\n", "In DFTK we achieve this using the following steps:\n", "\n", "Step 1: Build the 1D lattice. DFTK is mostly tailored for 3D problems.\n", "Therefore quantities related to the problem space are have a fixed\n", "dimension 3. The convention is that for 1D / 2D problems the\n", "trailing entries are always zero and ignored in the computation.\n", "For the lattice we therefore construct a 3x3 matrix with only one entry." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using DFTK\n", "\n", "lattice = zeros(3, 3)\n", "lattice[1, 1] = 20.;" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "Step 2: Select a model. In this case we choose a free-electron model,\n", "which is the same as saying that there is only a Kinetic term\n", "(and no potential) in the model. The `n_electrons` is dummy here." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Model(custom, 1D):\n lattice (in Bohr) : [20 , 0 , 0 ]\n [0 , 0 , 0 ]\n [0 , 0 , 0 ]\n unit cell volume : 20 Bohr³\n\n num. electrons : 0\n spin polarization : none\n temperature : 0 Ha\n\n terms : Kinetic(1)" }, "metadata": {}, "execution_count": 2 } ], "cell_type": "code", "source": [ "model = Model(lattice; n_electrons=0, terms=[Kinetic()])" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Step 3: Define a plane-wave basis using this model and a cutoff ``E_\\text{cut}``\n", "of 300 Hartree. The ``k``-point grid is given as a regular grid in the BZ\n", "(a so-called **Monkhorst-Pack** grid). Here we select only one ``k``-point (1x1x1),\n", "see the note below for some details on this choice." ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "PlaneWaveBasis discretization:\n Ecut : 300.0 Ha\n fft_size : (320, 1, 1)\n kgrid type : Monkhorst-Pack\n kgrid : [1, 1, 1]\n num. irred. kpoints : 1\n\n Discretized Model(custom, 1D):\n lattice (in Bohr) : [20 , 0 , 0 ]\n [0 , 0 , 0 ]\n [0 , 0 , 0 ]\n unit cell volume : 20 Bohr³\n \n num. electrons : 0\n spin polarization : none\n temperature : 0 Ha\n \n terms : Kinetic(1)" }, "metadata": {}, "execution_count": 3 } ], "cell_type": "code", "source": [ "basis = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1))" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "Step 4: Plot the bands! Select a density of ``k``-points for the ``k``-grid to use\n", "for the bandstructure calculation, discretize the problem and diagonalize it.\n", "Afterwards plot the bands." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computing bands along kpath:\n", " -½ -> ½\n", "\rDiagonalising Hamiltonian kblocks: 55%|████████▊ | ETA: 0:00:00\u001b[K\rDiagonalising Hamiltonian kblocks: 100%|████████████████| Time: 0:00:00\u001b[K\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=6}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeViTV/YH8C+b7JuAbIqKCgi4gYoKLq1BXIh1C9aNVu1A22lFWyt0mZaZ/mYKdlPbmRa0VbFqBbUq2LqAWgXcAAUNBUQUkEVF2UTWJL8/3teAGBSU5E3I+Tw+fSCGcCzhPe+999xzNSQSCQghhBB1pcl1AIQQQgiXKBESQghRa5QICSGEqDVKhIQQQtQaJUJCCCFqjRIhIYQQtUaJkBBCiFqjREgIIUStUSIkhBCi1igREkIIUWsqkAhLSkrKysq4joIQQojq6UwbURVIhN99993WrVu5joIQQoiKEYvFLS0tz3yaCiRCQgghRH4oERJCCFFrlAgJIYSoNUqEhBBC1BolQkIIIWqNEiEhhBC1RomQEEKIWqNESIgKi4s7NHTolC+//IHrQAhRYZQICVFV2dkIClqfk7Pv009/PHuW62gIUVmUCAlRMQ0N2LEDEyfC1xcTJgTb2MwcNy5g+XIMG4bvvkNVFdfxEaJqtLkOgBDSWbm52LoVW7di+HCsWoU5c6CjswxYxvxtejqio+HoCB4PQUGYOhUaGtzGS4hqoBEhIcqusRFxcfD1xUsvAcC5czh+HAIBdHQee5qnJ6KiUFAAHg/vvYehQxEZiYoKTkImRJVQIiREeeXkICwM/fohOhpBQSgsREQEBg582peYmSEoCFlZ2LkTBQVwckJAABIT0YkW/ISoKUqEhCidhgZ2CDh1KgBcuCB7CPh0bQeI779PA0RCOkSJkBAlwgwBHRweGwIOGPD8L8gMEDMzaYBISIcoERLCvY6GgNrdV83GDBBv3GAHiC4uiIzE3bvd9vqEqC5KhIRw6a+/unkI+HSmpuwAcdcuFBTA2ZkGiIRQIiSEC9IhII8HABcvdv8Q8OnaDhDXrqUBIlFrlAgJUShmCNiuELR/f26CYQaIly/TAJGoNUqEhCiCdAjo6wsAaWmKHgI+XdsB4gcfwNkZkZG4c4frsAhRCEqEhMhXdvZjQ8CbN7kcAj4dM0C8dAm7d6OgAC4uNEAkaoESISFyIR0CTpsGAOnpyjUEfDpmgHjzJg0QiVqgREhINxMKZawCOjhwHVbXmZg8NkAcOpQGiKRnokRISPeorUV0NHx84OcHtBkCamlxHdkLa7uCuG4dnJxogEh6FEqEhLyo9HQEB6NfP8TFISREhYeAT8cMEDMy8OuvNEAkPQolQkKeEzME9PTEnDkwN8eVKz1nCPh0bQeIoaFwcEBYGG7d4josQp4XJUJCuqztEDAsjC0E7deP67AUixkgpqfjwAFUVsLdHb6+iIuDSMR1ZIR0ESVCQjqrpkZNh4BPxwwQi4shELD94cLCUFzMdViEdJoqlHITwp2jR4+npGRMmfLGnj0We/di6lRERtLh7zIYGyMoiB0jRkdj2DCMGYPXX6+/dWuLra1lYOAirgMkpEOUCAnp0PXr9+bP//jhw6UbNvzzk082/fUX+vThOialxwwQv/oKu3cjLOznkpJiPb3E3r0d/f29uA6NENloapSQ9sRiJCZi8WKMHm2koSExMYn/9NOh69ZRFuwCZoC4ZYuTmVmKhkZRYKDdrFnYtw9NTVxHRsgTaERISKuSEvzyC6KjoaeHwEBs2qRrapp6584de3t7rkNTSX5+vrm5Iw0MDLS1DQ8dwo8/IigICxbgrbcwciTXwRHyCI0ICUFjI+LjERCAESNQUIC4OAiFCA2FpSV0dHQoC74IKysrQ0NDXV0IBDh+HJmZcHTE3LkYPRobN+L+fa7jI4QSIVFzTEfs/v0RGQkeD0VFiIqChwfXYfVcffsiNBTXryMiAunpGDyYduUT7tHUKFFH1dXYswcxMbh5E0uXIjUVjo5cx6RONDXB44HHQ1UVYmPxwQeoq8PixVi+XEnP5SA9G40IiXph9sI7OiIxEaGhbDs0yoJcMTNr7etdWYkxY9hd+c3NXEdG1Im8EmFDQ0NkZOTChQvDw8Nramo6eppEIvnPf/6zZcsWOYVBCKO0FJGRGDwYgYFwdERODmJjweer+1545eHpiY0bUVSEoCBER8PWFsHByMriOiyiHuSVCIOCgo4dO7Zo0aIrV64sWLCgo6dFRUVt2LBhz549cgqDqLmmJrYKZvhwFBRg2za2CsbKiuvIiCx6emxNzaVLcHTE7NkYPRrR0ait5Toy0qPJJRGWlJTs2bNn165dc+bM2blz57lz5y5fvvzk04qKir7//vvVq1fLIwai5nJyEBYGBwe2CqawEFFR8PHhOizSOf36ITQUBQWIiEBiIvr2pZoaIkdySYTp6emDBw+2trYGoKenN3bs2HPnzj35tODg4C+++MLIyEgeMRD1VFODmBj4+oLHA4CUFCQnIygIhoZcR0a6jqmpiY3FzZvg8bB2LVxcEB6OoiKuIyM9i1yqRsvLyy0sLKSfWlpalpWVtXvOzz//bGFhwefzN23a9PRXE4lEP/30U2JiIvOpl5fXxx9/3L0Bkx7g0iXNrVt1DhzQ9vISBQa28Pkt2toA8OAB15GRF6ajg8WLsXgxLl3S3L1bx8NDe/hw8fLlzf7+LTo6XAdHlJhYLNbV1X3m0+SSCPX19RsbG6WfNjQ0GBgYtH1CWVnZF198kZqa2plX09TUnDx58qJFbNNeS0vLdq9G1FlZGXbs0Pj5Zw1tbQQGSrKzJX36aAK9gF5ch0a6n7c3vL3x1VeShATN6Gjd0FC9BQskK1dKhg3jOjKilMRisaQT8+lySYR9+/YtLi6WSCQaGhoAioqKAgIC2j4hJSWlvLx83LhxAKqrqx8+fOjm5iYUCmW+moaGxuDBg319feURKlFRIhFOnkR0NE6cwMyZ+N//mLlQDYBOhej5DAwQEICAAOTlYdcujdmzNaysEBSExYtBKy2kHVEnTsiUyxqht7e3RCI5duwYgKysrLy8vJkzZwK4du3a77//DmD27NlFRUVpaWlpaWkffPDB2LFjU1JS5BEJ6XlycxEWBnt7hIWxVTAxMeyKIFE3Tk4ID8eNG2xNjb091dSQ5yGXEWGvXr02bNiwZMmSsWPHXrx4MTIy0szMDEBiYmJUVNTMmTN79erVqxc7c6Wvr6+jo8M8gZCO1NcjIQHR0cjOxrJlSE7G4MFcx0SUg7RPTXk59uzB+++jsRHLl+P112FtzXVwRBVodGb+9PncuXNHKBQ6OTlJexbX1dU9ePDA+vH3JjM1amtr29HrhIWFmZiYfPTRR3KKkyg55qDXPXswZgyCgjB3LrSpMyB5qnbvmTlzQDU16kksFotEIp1n/fjlmAi7CyVC9cTc3f/0E5qb8frrWL6cjgMkXdNuFmHlSgwZwnVMRLE6mQip1yhRLiIREhMREAB3d6Sn45tvkJ2N0FDKgqTL9PXZPjUnTgDAxIlsnxraUUPaoURIlEVeHsLD4ejIVsHcvMlWwWhQHSh5Mc7OiIhASQlbUzNgAAID8WhnMiGUCAmnhEKhg8M4W9tJXl53XnoJjY04fhxpaQgKojp40s20tNg+NVlZcHXF3/8Od/fmAQP4trZjDh8+ynV0hEuUCAk3xGKcOIEVK04XFy+oqBg+c+alwkJ88QWcnLiOjPR0dnYIC0NODj77rLS8vKm8/JM33/z9wAE0NXEdGeEIld8RRSsqwu7diI6Gnh4WLFhoaRlmbW384YcvUS0oUSQNDQgE/bOypqamxs6Z88///hcrV2LBAixbRs3Z1Q5de4iCVFfj4EHs2IFLlzB/PuLi4OEBoDcQzXVoRH19/vk65oN338WtW9i5E8uXQ0cHAQFYvhz9+3MbHVEQmhol8iUWIzERgYEYOBBxcQgKQlkZoqKYLEiIEunbF6GhuHYNO3agshJjxsDXFzExqKvjOjIiZ5QIibz89VdrFainJ3JzER8PgYC2NhNl5+mJjRtRVISgIMTFoX9/tspU6Tddk+dEU6Okm1VVITYWMTG4eRPz5+PgQYwYwXVMhHSdnh4EAggEKC1FXBzWrkVlJRYtQlAQHB25Do50KxoRku7BbIQPDISjIxITERqKmzexcSNlQaLy7OwQEoLLl7F/P+rrMX48fHxoY36PQomQvCihsPU4CE9P5OcjNhZ8PnUEJT2NdMo0NJQOu+hR6FpFntP9+9i7F1FRuHMHS5bQcRBEXejqgs8Hn8/+CoSGsr8C1MtUddGIkHRNYyPi4xEQgIEDkZiIyEgUFSEigrIgUTu9eyMoCOnpOHIEeNTLdONG3LvHdWSkiygRks5KT0dICBwcEBkJHg8lJYiNpV6ghMDNrbWXaXo6hgxBQADi49HSwnVkpHNoapQ8A1Myt3Ur6uuxaBHOncPAgVzHRIjyYXqZ8nhs4XRkJN58EwsWYMUKKhlTdjQiJLI1NCAuDnw+hg9Hejq++go5OQgPpyxIyDOYmSEoCMnJSEqCuTnmzIGbGyIjcfcu15GRDlAiJO1Jp0CjoyEQoKiIjkMi5Hm4uCA8HNevY+NGCIVwcQGfj7g4NDdzHRl5HE2NElZxMXbtwpYtbKPFixep0SIh3UBTk50yZdrtRkfjrbcwfz6Cg6nRoLKgEaG6q69HXBx8fTFyJAoKsHUrsrMRHk5ZkJBuZmqKwEAcP46MDDg6QiBgp0xv3+Y6MrVHiVBNicVITkZwMGxtER3d2gubDqAhRN4cHNju3lFRKCiAszPb3bu+nuvI1BUlQrWTm4vwcAwejOBgODoiLw/Hj0MgQK9eXEdGiDrR1ISPD6KiUFbGdvfu1w/BwUhO5joy9UOJsOcrLS2VSCTV1YiJga8vXnoJlZXYtw9CIUJD0acP1/ERot709SEQID4ely/D0RErVmDoUISH4+ZNACgpKeE4PjVAxTI9XFBQWGxslqamWCI5Mm0a1qyBnx+0tLgOixDyBOZAxHXrkJKC7dvh6YlevYIfPiwbOVL3zz/juI6uJ6MRYY914QJCQrBt25Xq6jdaWkoLCrBnD2bOpCxIiFLT0ICPDzZvxq1b0NT8q6YmODn5xvLlSEyESMR1cD0UJcKepqgIkZFwccFrr8HcHEePfr9mTebRo9Hm5lxHRgjpCn19JCVFh4ScP3Fiy/jxCA+HgwNCQmgRsftpSJT+BJGwsDATE5OPPvqI60CUWmUl4uOxYwcuXcL8+Vi2DN7etAWekB6lsBC//oqffoK2NgICsGwZBg3iOiblJhaLRSKRjo7O059GI0LV1tDAngUxYAASErBqVesuCMqChPQw/fsjNBR5edixA5WVmDCBPe+Cmre9IEqEKqntLsC2Z0Hw+XjWrQ8hROUxRwSXlrLnXTDN22Ji8PAh15GpJqoaVTFCIeLisH07DAwQGIi//oKNDdcxEUK4ID3vor4eCQmIicF772HmTAgEVBbXNTQiVA23bmHjRnh4YOZMNDTgyBF2FyBlQUKIdCdiVhY8PREZiQEDEBKCjAyuI1MRlAiVWlUVuwt+5Eikp2P9ety8iYgIODtzHRkhRPnY2bFlpceOwdwcAQFwc0N4OG7c4Doy5UaJUBk1NiI+HoGBcHREXByCglBaSmchEUI6i+lNk5+PmBhUVmL8ePj4YONG3LvHdWRKiRKhEmFKYJizACMj4eODmzcRH0+NQAkhz4kpqyksRGgo0tPh5MSeidjUxHVkyoSKZZRCdjZiY7FjB/T0IBDg/HkMGMB1TISQnkJXF3w++Hz2TMSYGLz9NmbMQGAgpk6leSZKhJwqKcHevYiLQ2Eh5s3D3r0YNYrrmAghPRdzJmJgIHvxCQ3F3buYOxfLl2PkSK6D4w5NjXKAOQiCz2dLYMLDUVSEjRspCxJCFMTeHiEhSE/HH3/A3Bzz5rGnBJeVcR0ZF2hEqDgiEU6eREwMDh/GhAkQCLBnDwwMuA6LEKLG3Nzg5oZPP0VqKuLiMHIkhgxBYCAWLYKxMdfBKQqNCBUhPR0hIbC3R3g4PD2Rm8sWhVIWJIQoA+aU4I0bUVSE0FAkJqJvX7asprmZ6+Dkj0aEcpSTg19/xc6d6NULAgFSU+HoyHVMhBDSMWlZTVUVDh1CdDTeeqvn9/GnRNj97t3Dvn2IicHNm5g/H1u3wseH65gIIaQrzMzYspriYuzfj3ffRWUlXn0VK1bAyYnr4LobTY12g+bmZgD19YiLA58PZ2ckJyM0FDdvYuNGyoKEEBXWrx9CQnDpEg4fBoApU9iymtu3gUdXP1VHifCFiMXiMWNmWllNGDt2h50ddu5EYCBu3WKLQrVpvE0I6Snc3BARgeJifPMNhEIMHYrx48/06TNh8OAJlZWVXEf3QigRPiexGCdPYuXKB+npFdXVn9XXJ12/jgMHIBBAT4/r4AghRD60tODnh5gY3LoFG5vk6uplBQV9X301/+BBNDZyHdzzojFLlwmF2LEDO3bAzAwCgcnnn7+bnh6/fv0nvXtzHRkhhCiKgQE2b/6bkdHntrZj3NxGb9mClSvZQ6CmT1exg1E1JBIJ1zE8Q1hYmImJyUcffcRtGMxBgLt2QVsbAQFYtIiOgCCEkFYVFdi/HzExyMlRlmMRxWKxSCTSeVZapqnRZ7h5E5GRcHXFzJmorMTPPyM7G+HhlAUJIeQxlpYICkJyMjIyHjsWMTkZSj7gokQoW3ExW/A5ZgwKChAdTSWghBDSKQ4ObP47ehTm5njjDQwcyHZ0U06UCB9z795jB+GGhqKsDFFR8PHpsTtJCSFETlxdER6OnBwcPgxzcyxcyB4UnJfHdWSPo2IZAGwPhbg4nD4NPz+sWgU/PzoCkBBCugfT0TQ8HOnpiInBlCkwN4dAwB4/zjm1HhHW1yM+HgEB6N8fcXEQCFBaithY8PmUBQkhpPsxBwXfuoWoKFRWYsIEtsdpeTmXUaljImxsZHte29pi0yb4+6O4mH3E0JDr4AghpKeT9vguKWGHiUOHso9UVHARDwffkyMiEZKTERwMW1tERsLTE9eu4fhxBAbCxITr4AghRP1oaYHHQ0wMysoQGor0dDg7g89HTAwePFBcGD1/jVAsZs/Z+vVX9OuHZcsgFMLWluuwCCGEPKKnx5568fAhDh9GTAzefReTJkEgwIIFcj+xricnwsdbwCA1FYMGcR0TIYSQjhkYQCCAQIDKSsTHIy4O770n94Y1PXBqVChEeDicnMDnA0BSEvsIZUFCCFEV5uYIDER8PDIz2e35trYIDERiIsTibv5eT0uEdXV1hYWFJSUljarQS1XaAmbGDLYFTEEBIiLg4sJ1ZIQQQp6XvT27PZ9pWBMejv79u7lhjYxEmJqaGhwcPGTIECMjowEDBvTt29fQ0NDT0/PDDz+8du1a93zb7nPrVvsWMIWF1AKGEEJ6mo4a1mRkvOgrP9Z0Oykpad26dRkZGQMGDJgwYYKTk1Pv3r1bWlru3buXlZV19uzZe/fu+fv7r1+/3kWB4yyZTbfv3cPhw9ixAxkZmDULAgFmzKDz/wghRI0wZyH88gt0dSEQYPFiODk99oRONt1uTR379+9ftmzZG2+8sWXLllGjRj35VJFIlJSUtGXLlhEjRuTl5fXv3787/iFdQy1gCCGEMNo1rJk8Gb17QyDAa69h4MAuvE7riPDatWvGxsY2NjbP/JqrV6/a29ubm5s/X+hdtXr12ri4E5s2pe7Zo3f0KFtQO38+bX4nhBDSSrpZbs8eDB4MgQAJCcuDg5cvWDDp6V+oAucRBgV9uXnzTVtbr3//O3DBAhgbcx0QIYQQJdbYiD/+QHh4embmx6tX7/32W6OnP//Zq2oSieTq1au3b99mPuXxeN0QZldoatZoaCTY2S157z0cO8YuB+rrKzgKQgghyq6lBSdPIi4OBw6gb9/+Wlo55ubpwOSnf9WzE+GcOXMqKirs7e2ZTzufCJOSkg4cOGBiYhIcHOzg4NDub69fv75///7CwkJDQ8MZM2ZMmTKlo9cxM2v+17/+9sknE+7fR0ICtm7F8uXsBOncuTRAJIQQdScS4exZdlLUwgICAVJSMGSIpVhcUF9f/8wvf8aG+urq6paWlpSUlNhHOhnWgQMHFi5c6OrqWl9f7+Xldf/+/XZPSE9Pv3Pnjru7u4mJyfz587dt2/a0KDU1AfTuze6vvHEDAgHi4tC3L9uVrra2k3ERQgjpIZgO0iEhsLNDcDDMzZGczHZQGTKEfU6vTpRTPmNEaGJiYmZm9hzxRUREREREvPHGGwCys7O3b9++Zs2atk8ICAgICAhgPtbS0oqNjX399dc7+eJMRgwMbO3BI+1KR2NEQgjp2aTjv19/haVlN3TQfFoi3LhxY0NDQ01NzaJFi0aOHMk8GBoa+swXbWpqunDhwu7du5lPX3755TNnzrRLhFJ1dXWpqamenp5djBx41INHZkacM4fOlCCEkJ6j7fxn377w98fZs91zru+z1wh9ut6j5fbt2xKJxNLSkvm0T58+ZWVlTz7t/PnzM2bMqK2tnTJlyj/+8Y+OXk0kEu3atSs9PZ351MPD48mcqquLBQuwYAGqqzUOH9b69Vetd97R9PERz5sn4vNFxsbKXhlLCCFEpsZGJCVp/fGHVkKC1qBBkrlzW5KTRXZ27FX94cOnfa1YLH7mbno8PRGGhIR0JdpWzJxsS0sL82lTU5Ouru6TT/P09Lx+/Xp5efmqVatWrVr1448/ynw1TU3NYcOGzZkzh/nUzs5O5qsx+vTB8uVYvhxVVeL4eI29e3Xef7/XxImSBQsks2dLTE2f7x9ECCFEoRobcfy4xt69GocPawwdKlmwQBIeLrazA6AFaHXyRcRicWe2CD62jzAmJsbNzU06S3nlyhUHBwfTR9kjIyNj+/btGzdu7Mz31tfXT0tLGzZsGIBPP/30+vXrO3fu7Oj5p06dCggIuHPnjsy/ldlirfOqq3HsGOLjcfgwJkyAQIBXXgFlREIIUUINDTh+HHFxSEiAqyt7JJOd3XO+WidbrD1WNfqPf/zj2LFj0q8fPnz4H3/8If3bnJycTZs2deZ7a2pq8vl8Zo2wqalp//79r7zyCoD6+vrff/+dOcuipqZG+vzU1NSBXeqH0xWmphAIEBOD4mIEBSExEQ4O8PHBxo24e1dO35MQQkgXNDQgPh6BgbC1RWQkPD0hFLZWhMqbvNpUh4eHT506NScnp7Cw0MbGZu7cuQBu3749a9as0tJSW1vbJUuWVFRU9OvXr7CwsLi4eP/+/XKKRMrAoPUE5KQkxMXhs8/g7g6BAIsWoU8feX9/Qgghj5GO/+Lj4eYGgYA9d1DB5JUI3d3dc3Nzz5w5Y2pq6u3traWlBcDe3j4zM9PKygrAvn370tPTS0tLra2tR48eraenJ6dIniTNiPX1SEykjEgIIQolvfZK89/69ehEo2t5kePBRWZmZnzmkPhHdHR0hg8fznzcq1ev8ePHy++7d4a+focZ8dVXYW3NbXSEENKjSK+0hw6xV9ovv1SKKy2d4AfIyojh4ex9ysKFXN6nEEKIqlPa/CfVPhEeO3bswYMH0k/37Nlz5coV5mOhUKi4uDgizYjSmWvKiIQQ8hyezH9ffaWka0+PbZ/o379/UVHR079A8cc2veD2iRf05FpuQAAHa7mEEKISpNWI0vzHYe1Fl0+oB5Ceni4SieQZlerR0+twjEgZkRBCGNL8d/gwPDzg74+vv4aVFddhdc5jiVDaFI08qW1GPHIEe/ciPBwjR0pKSxdWVRV8803YkiULuI6REEIU59y58/Pnv6upaTxs2G8pKSZTpmDBAnz3nep1LHnGMUzkSXp6mDMHv/yCsjIEB9cWFhbdubP+3XcPff01Cgq4Do4QQuTv9m1ER2PJkqTS0r/dvm07blxOcTEOHsSyZaqXBdE2EcbFxYWEhMjsjt1WVlbW/Pnzr127JufAVICeHl591eRf/wp86aXNGzaEFhVh4kS4uSEsDMnJUPhaKiGEyFdhITZuhK8vnJwQH4+QkOXTp58NDh74ySejVfq0n9ZimYKCgpUrV6akpPj5+c2fP9/Ly8vZ2Zk5EbepqSkzMzMlJWXPnj0XLlxYunTpf//7XyMjI8WEyG2xTJeIxbh0CfHx+PVXNDZi9mzw+ZgyBdq0S4UQorKEQiQkID4eubmYMQMCAaZNQ8dnHyiRThbLaLSrAj18+PB33313/PhxsVgMwMzMrKWlhdlQoaenJxAI1qxZM2rUKPnF/SQVSoRtPfnW8fNDJ45KJoQQpSAUsof/1dfjlVdU8rb+ORMho6ys7NSpU1euXLl9+7aurm6fPn3GjBkzceJEEy5GvyqaCKVu3sTBg4iLw19/YepU+Ptj7lwYG3MdFiGEPEF6+O2+fTA1BZ8Pf394e0NDg+vInssLJUKlouqJUKqiAr//jrg4nD6NSZPA52POHCXdXkoIUSvM5veEBBw4AEtLts2kiwvXYb2w59lHSOTK0hKBgQgMRGUlEhMRH4/QUHZL4oIFsLfnOj5CiJqRXoukDUM+/VQdr0WUCDlgbs6eNim9C/Pw6FF3YYQQZfbk7JTSNj9TDEqEXJK2Nv3f/9h5eR6vJ8zLE0KUULt6BYEAu3ZRvQJAiVBJaGnBxwc+Pti4ka3U+tvfVLhSixCiPKQV7Hl5mD4doaFUwd4eXV+Vjpsb3NwQHs6+fcPDVW/vDiGEW0/uaQ4Pp1vqDtH/FeXFZMTQUBQW4sABbNqEpUsxaRIEAsyZA5Xu40AIkQfp5oe9e2FmBj4fW7bQIsuzUSJUAf37IyQEISGtS9zvvouxY+Hvj1dfVa7zLQkhiic9+e/wYQwcCH9/nDgBZ2euw1IdlAhViXQDRl0dTpx47Eyo+fPRty/X8RFCFEi6+UF68l9EBOzsuA5LBVEiVEmGho+dkpiQgC++gIUFBAIsXIihQ7mOjxAiN3fv4o8/Wjc/CATYtAlmZlyHpcooEao26SmJ0g0Y06ahVy/4+0MgYNcGmpubn9lYgRCinKS/vzdu4NAhxMUhJwcvvwyBALt3Q1FnH/RwlAh7COkGjA0bcOECfvsNK1agqQk2Ns/VBjwAACAASURBVF/n5cWOG+f8++8xXMdICOmaVas+27nzqK3taOD7ykrMmUPFn3Ih+39nSkqKhobGhAkTANTX13/22Wfnz58fM2bM559/rq+vr9gISddoaMDLC15eiIjA1avw9T1aWfnd0aMrX38dfD78/OgWkhBl19iIU6dw6BCio4+1tEQ3Nb12/Di8vKj4U15kn1C/bNmytLQ05uPw8PCvvvpKLBZv3rw5KChIgbGRF+XujoSESH//H6Oj//3yy9i2DTY27Lb94mKugyOEPO7+fcTFITAQtrb4+GNYWeGHH76eNWtDbOz6ceMoC8qRjNMnHjx4YGxsnJ6e7uHhIRKJbG1tly5d+s033xw7dmzWrFkVFRWmpqaKDLHHnD6hDB4+RFISEhJw6BB696ZeboRwr6AA8fFISMDFi5g4EXw+Zs+GjQ3XYfUInTx9QsaIsKamBoCFhQWAS5cu3b17VyAQAJg0aVJLS8vNmze7P1iiKAYG4PMRFYWSEsTEQE8Pq1fDxgaBgYiPR2Mj1/ERoh7EYqSnIzwco0dj/Hikp2PVKty+jfh4BAVRFlQ0GYnQyspKU1MzPz8fwN69e01MTEaPHg2AOadeS0tLwSESedDUhKcnwsORloazZ+HpiU2bYG0NPh/R0Sgv5zo+Qnqi+nrExyM4GH37IjAQDQ3YsAHl5YiJAZ9PDRQ5I6NYRkdHZ/r06W+++aZAIPjxxx/nzZvHjCuvXLmiqanZr18/hQdJ5MvRke1cc/8+kpIQH4+wMDg6wt8fAQFwdeU6PkJU3J07OHIECQk4dozd+f7JJ6BLqfKQXTUaHR0dHBy8bds2Hx+fiIgI5sFt27aNGDFCwQuERJF692YPSmxpwblziIvD9OnQ0YG/P/h8TJ4M2o5ISOdJj31gdv75+2PzZtAVVAnJKJZRNlQswy3pL3NuLl56Cf7+1PKbkA5JbyJ/+41uIrnXyWKZp23LbGlpKSwsvH37NrOhkKgn6SEYzPROXBzefhsjR0IgwLx5NL1DCIDWZYWEBHZZ4bff4OnJdVikc2QnQpFIFB4e/u2339bV1dnb29+6dQtASEhIQ0NDVFSUYiMkyqJPH7blt3QPRkQE7cEgau3JnQ/r11PNp+qRvaH+s88++/LLL999913pAiEAX1/fXbt2NVKJvdqjPRhEndHOh55HRiJsbm7etGlTRETEF198MW7cOOnjI0eOfPDgATM6JAS0B4Ook4cPaedDjyUjEd69e7e2tnbatGntHmfqRe/fv6+IuIiqYfZgHD+OggIEBiI5Ga6uGD0a4eFIT+c6OEKe1507iIlBQADs7BAZCVdXnD8PoRAREfDxoeWAHkLGGqGxsbGmpmZZWZnr4zvIrly5AsDW1lZBoRHV9OQejLlzqXyOqBja+aBWZIwIjY2NJ06cGB4eXltbq/Hohqeqqio0NHTkyJF96Rx00jna2myD76IiHDoEOzuEh8PODgEBiIlBdTXX8RHyuJYWJCcjJAQODpg9G6WlCA9HWRliYxEYSFmwJ5NdNbpx48bJkye7uLi4urrW1NS89tprR48era6uTkxMVHB8pGdotwcjIQGrVsHdHdOmlR48+LaVldG+fVGGhoZch0nUi1gsXrly7YULV197bf3VqyNp54Pakl01OmLEiLS0tKlTp169erW2tvbgwYPjxo07e/ast7e3guMjPQyzByM2FrduYc0aHD58ICNjWmKizcqVZ1JSIBJxHR9RG5cvIzS0aOfOvOzssO+//5nHQ14e0tIQHk5ZUO10uKF+8ODBMTF0pjmRFyMjzJ8PD48ZfP4bGhoGLi6fhIay6zE8Hvh80GI06XYPHyI1FfHxbNuXqVP7jRxpWVHxaUzMeh8froMj3HlaZxlC5G3gwIFXryYxH4eHt06cMl2/eTzaqk+6QUEBEhMRH4+TJzFyJPh86eSnFrCN6+gI92Qnwg8++IA5lfBJ1FmGyI+0eY1IhLNnkZCA1atRVIQpU+Dvj9mzYWbGdYhERTQ0IDmZzX8VFfDzQ2Agdu6kNrlEBtmJ8PTp0xUVFdJPa2pqKioqDAwMbKhrAlEILS34+ICZrbpxA8ePt9bX8Png8WgVh8h2+zaOHkVCAo4fx6BB8PdHTAw8PGhSgTyN7ER4/vz5do9kZ2cvXrx43bp18g+JkMcMHIigIAQFob4eKSmIj8f8+dDQwLRp4PEwfTqMjbkOkXBKJMLly2zPz8JC9oyUqCiYm3MdGVERXTiG6fTp0/7+/uXl5QYGBnKNqR06hok8Sdrs+MIFjB0LHg+vvAIXF67DIgpUUYGTJ9m3ga0tO1UwZQq0qfKBPNINxzC14+zsXFtbm5ubO2rUqBeLjZAXxXR0CwlBXR1OnEBCAng86Oqy/WsmTUKvXlyHSORD2vMlKwuTJ4PPxxdfwN6e67CIKutCIoyPjwe1WCNKxtAQfD74fODRJTI8vPUSOWsWXSJ7AuntzuHD7O1OeDjd7pBu06mq0ZaWlvz8/DNnzvj6+lK9DFFa0v410kmzsDCaNFNhT06AJyXB2ZnrsEiP06mqUW1t7b59+0ZERLzzzjuKCoyQ52dpyTb+lpZRhIWxZRTMaqK1Ndchkg5IS6IOHmRLooKC8NtvMDLiOjLSc3W2apQQVaSlBU9P9tBEaWE9s1ufWU2kwnolwWySSUzEsWPsJpl9+2iTDFEQmioi6sLamt2tz5wPlZCAwEB2qzWfDz8/2mqtaNIfRGJia9uE6Ghqm0AUrTUR3r17Nz8//5lfMH78eHnGQ4jcMedD+fggIoJtvhUTg+XL2eZbtFtf3qSN9BIT2UZ6GzZQIz3CpdZEmJCQsGLFimd+Qef3HRKi/Bwd2d360nbMzDHCTJtTX19UVZVfu3bN29tbU1P2US3k6dLS0szNzQcMGCTd837zJttafeNGaq1OlELrhvrbt2/n5uY+8wsmTZok55Daow31RMEuX8Yff+DwYWRlPWhq8tbU9Fi6dGB09Kdcx6V6fvhhzwcf7GxpKdfX3zFggPPMmZg5E+PGQUuL68iIeujyhnpra2trqqUjBBg5EiNH4sMPcfOmaMQIzbo6y50760+fxrRp8PPDlCmgI4SfoqkJqak4ehTHjiEnp76pyUxXt3rv3qapU7mOjJAOULEMIR0aMMD0/PndV65cnT2bX1KC+Hhs2ICFC+HiQkWn7TGrrcwfZuUvMhI+PoGJiVZWVpZeXsO4DpCQDnXYa/T8+fOxsbEFBQUPHjxo+/jx48cVElgrmholSoVZTWSu+MXFmDyZXVC0s+M6MoWrq8PZs+xRR/fuYdIk+PvD3x+9e3MdGSEAXrDX6O7du5cuXdqvX7+mpiY9PT1TU9Ps7OxevXqNGzdODqESokoMDMDjgccDgLIy9oio0FDY2bFFpz2+9RfTyi4xERcvYswY8Hh01BFRbbIL4T755BOBQHD9+nU/P7/FixdfunQpJyfHyclJ8ZUyhCgzW1sEBiI2FhUViImBuTkiI9GnD3x9ERmJ7Gyu4+s+d+8iLg7BwbC3x+zZKChAUBCKi3H8OEJD4elJWZCoMBkjwrq6uhs3buzatUtLSwtAU1MTgIEDB0ZHR3t7e4eEhJjQxmNCHidtYRMainv3cOIEEhMxfTq7E4PHw7RpMDXlOsouamlBZia75yEvj+1j/tFH6N+f68gI6VYyEmF9fb1EIjE1NQVgYWEhbTrq6ura2NiYn5/v4eGh0BgJUSkWFmynU+CxDfvSEptRo6DMmxKfLHuJiMDkyXjWOgshqkpGIrSwsDAyMiosLHRxcXFycvq///u/Bw8eGBkZJSUlAbC0tFR4kISoKumGfaaXdGIigoPZdmI8nhKdEsWUvcTH49AhNDRg4kQ65J2oEdlVo3PmzLGxsfnxxx9ra2sHDhxoaGjo5OR05syZcePGnTx5UqNzqwHFxcW//PJLU1OTQCBwdXVt97dNTU1JSUlpaWlaWlo8Hm/s2LEdvQ5VjZIeprwcx46x9Sbm5uwwceJE6OoqNAyxGJcusSM/adkLj0dlL6Tn6GTVqOwJmp9//vmTTz4BYGxsfOrUqWnTpkkkknfffffAgQOdzIIlJSUeHh7l5eVisXj8+PEZGRntnrBhw4bPP/+8oaGhpqbGz88vKiqqMy9LSA9gY8OW2Ny9i9hY2NkhMhLW1myJTXq6fL/7nTts2UvfvggIQEEBVq1CeTmVvRD11eE+whf0ySef5Ofn//rrrwD+8Y9/5Ofn7969u+0TmOlW5uOffvrpm2++EQqFMl+KRoREHTx4gJMnkZCAI0egpQVfX/B48PXtnqMY2p7zID2XccYM9OvXDS9OiNJ6oRHhRx99dPbs2Rf59idOnPDz82M+9vPzO3HiRLsnGLU5Z1Mikejp6b3ItyNE1RkZgc9HVBQKC5GYCE9PxMVh4ECMHo2wMCQnQyzu8msWFCA6GgEBsLLC6tUAEBGBsjLExiIoiLIgISzZI8IhQ4bk5+c7Ozu//vrry5Yts+/6gr6jo+N33303a9YsAMxLNTQ0yEzL9+/fHzlyZGRk5KJFi2S+1Lp1606ePDlixAjm02HDhgUHB3c1HkJUUX09zp7VPHFC88QJzeJijUmTxC+/LJ4+XWxv3+FEzoMHuHBB8/ffteLjNRsbMXWqeNYs8csvi83M6NwYonbEYrGWlpbus5bfZXeWEQqFR48e3bFjx6effvrxxx+PHz8+MDBwyZIlhp1uNqylpSV+dAfb0tKioaEh8xSbBw8e8Pn8uXPndpQFGfb29p6PzohzcHDQot71RD0YGcHXF76+EkB0+zaSkzUTErQ+/VTH1FQydapk1iyJgUHy559/+dZbSwYNEiQlaZw4oZmWpjF6tOTll8WxsSIPD2nyU+LtGoTIjYaGRmfqWp6xRnj79u1du3Zt3749MzPTzMxs4cKFP/74Y2e+/eTJkxcvXswM3U6cOLFs2bKSkpJ2z3n48OGsWbOGDBkSFRX1lFhpjZCQtkQiXLyIY8dw8CAyMqYCwcC/nZ0zX3kF06bBx0fR1aeEKK0XWiOUsra2XrNmzeXLl0+fPm1qatr52s5Zs2b99ttvzMf79+/39/dnPs7IyLh//z6A+vr62bNn9+/f/8cff+xkJSoh5M4dxMZiyxZs3YrychgYNGho/Kyt3SASYds2/PQTfvkFN29yHSUhKuUZxzCJxeKkpKTt27f/9ttvDx8+9Pb27uTr/u1vf/vpp5/8/f3NzMwSExOTk5OZx1955ZVvvvlGIBB8/fXXJ0+enDx5MlNTo6urm5CQ8CL/EkJ6qrbnXeTlwcsLPB7eeovZ8Jdy+vTpCRMmaGujrAzJyUhMRHg4evUCjwdvb/B46ngsBiFd0uHUaG5u7u7du2NiYm7cuGFvb7906dIVK1Y4OTl1/qUfPHjwxx9/NDU1TZ8+3cLCgnnw/PnzgwYNsrS0vH79+o0bN6RP1tTUfPnll2W+Dk2NEjUkEuHyZTb5nT0LFxd2t3snW51J26QdP47evdmv9fMD9QkmaqWTU6OyE6Gfn9+xY8f09fXnzZv32muvTZ06VWapi2JQIiTqo22fT1tb+PiAx8P06TA2fs4X7Cih9vizogjBC55HaGhouHnz5oCAADpoghB5k05pHj4MXV3weBAI8MMPeDSN8kLaHoshnWINC2udYqWeaoTIToT79+9XcByEqJWKCpw8ieRkpKS0HnP/4YcYMECO37TtkcJMAImJCAjAgwdsANOmyTcAQpTT04pl7t27V1JS0tzc3PZB6X4+QkiXyKx5iYriZkBmadl6VlRpKXsyxr/+xQ5JeTxMnYrevRUdFSGckJ0Is7Oz33zzzTNnzjz5V3LqTUpIjyRziS4iQrmW6Ozs2h+gyHTlZg4j9PbGlCnPv0hJiPKTnQhfffXVO3fufPvtt87Ozs9cZiSEtPNk0WZQEPbuVYFD6qUHKDLH0ycmYtMmvPoqhg7tWtkqISpERiKsqam5evXqvn375s6dq/iACFFR5eU4cwaJifjjD+josHv4NmxQ1W182tqtVTbMsb1UZUN6qg7XCPtRa3pCnqW2FufPs4M/6bnzYWEYOJDryLqVoWFrlc3duzh1ComJEAjw8CEmTWI3eDg4cB0lIc9LRiI0MTGZMWNGfHz86NGjFR8QIUquvh4pKWzBp/Rg96gojBoF7nbbKo6VVfsFxcREfPQRTE3ZZMnjwdyc6ygJ6QrZI8LVq1evXLmyqqpqxowZVlZWbf+KqkaJWrl48WJaWubixa/m5xu1q3kJDcXEiWrd4Vq6oIg2VTZBQRg0iM2IXl7N+/fHWltbTp/ux3WwhHRIdmcZGxub27dvy/wCxVeNUmcZwgmxGKmplTNm8B4+nKOtXe3i8tXUqZg6FZMno82p0qS9xkZ2QTEpCZcvRzc1XdbVzf/hh4glSzy0n9HbmJBu9kKdZWJjY5uamuQQFSFKTSzG1as4dQqnTuH0aVhY6EkkIgODy2++Oe7LL7kOTkXo6mLKFEyZgv/7P+zcafPWW7kiUeX69WarVmHCBEyejClTMHo0KCkS5SH7zThp0iQFx0EIh6RrXSdPwsQE3t6YPh0bNsDBQb+29kxhYaG7uzvXMaqkJUtme3m5mpiY9OnTR1pY9M47bOmptzd8fJRrSyVRT0+7K7tx44ZQKKyurl6yZAmAmpoaHR0dfX19RcVGiBw9mfx4PHz1VfvqR2NjY8qCL2Lw4MHMB8bGraWn0qQo3Y9BSZFwSHYirKure+211/bt2wfA3t6eSYRr1qwpKSk5cuSIQgMkpPt0MvkReaOkSJSK7ET4zjvvnDp1KiYmRkdHZ+3atcyDS5cunT59+sOHDw0MDBQYISEvRGby+/JL9O/PdWQEACVFogRkJMKGhobdu3dv2bJl6dKlf/75p/RxV1fXpqam4uJiZ2dnBUZISJdR8lNRlBQJJ2Qkwnv37jU2Nj65m15XVxdAbW2tIuIipIso+fUwlBSJwshIhL1799bR0cnJyXFxcWn7+Llz5zQ0NPrTdYUoDUp+aoKSIpErGYlQX1/f39//ww8/HDFihMajlrp5eXlr1qx56aWX2jWaIUTBKPmpOUqKpNvJLpb57rvvJk2a5OLiMmDAgPv373t5eV26dMnCwuLQoUMKjo8QPEp+yck4cQK6upT8CIuSIukWshOhvb19RkbG//73v2PHjgHQ0tJ6//33V69ebW1trdjwiPp6Mvn5+ODzzyn5EdkoKZLnJrvXqFKhXqM9W11dnba2NlOKJTP5+flR8iPPr+1RWU8mRbFYXFNTY2ZmxnWYRC5eqNcoIYqRknJ29uxVIpFowoT9aWkDzMwweTL8/PDFF7C35zo40iO0HSneu4fTp3HqFN5/H4WFGDeuPi1tikikExn5VlDQEq4jJZyRnQgDAgIqKyuffNzY2HjgwIFz58718fGRc2Ckx2psRFoakpOZkZ+wvv4lHZ2SUaMKNm8eQMmPyJWFBebOxdy5AHDvHg4cuP/nn70aGl59551Le/YsmTgR3t4YP55OF1E7HY4IL168WFdX5+rqamVlVVpampOTY21tPWTIkOTk5G+//TYiImLdunWKDJSotAcPcO4ce5jtuXNwcICPDwQCfPPN0l9+2WRqOnDNmpceVSgToggWFli50t7QcPX581khIaH5+UhOxvr1re9Pb2+89BL69eM6UCJ/stcIv/322+3bt+/fv9/R0ZF5JDMz85VXXvn666/9/f1DQkK2bt1aVFSkmNoZWiNUUeXluHiRPcz96lW4u7MXFx8fOsGcKK+WFmRmsjdtJ07A1JR903p7w9UVdLumWjq5RigjETY3N1tZWR08eHDy5MltH9+2bdvXX3995cqV+vp6c3Pz3bt3z2WmGOSMEqEKKShgryDJybh1C2PHshcRNT/JnaguaQHXmTOor8eYMWxS9PLCs66uhHvPXyxTUVFRXV1taWnZ7nErK6tr164B0NfX79evH/VaIwBEIuTksJnv1Cm0tLCXiaAgjBoFTU2u4yPkxTg6IigIQUEAUFrKvtVXr36sANXbG3Q8nUqT3WLN0NBw27ZtX7Y5k1sikWzbtk3aX62iouLJTEnUxMOHyMhgrwjJyejdm70c/POfGDiQ6+AIkRs7OwgEEAgAoKYGFy4gMRHh4bh8GS4u7G/Byy/DwoLrQEkXyUiEurq6a9eu/de//iUUCv39/ZlimdjY2NTU1K1btwL4888/q6urPT09FR4t4Yz01z45ufXXPjAQ27bRrz1RRyYmrbsypLeGMTEICmq9NfT1pVtD1SC7WEYikWzcuHH9+vVlZWXMI4MHD/7nP/+5ePFiAFVVVVVVVQMGDFBMiLRGyBXpRFBKCk0EEdIpHS0W+PjQYgEHnr9Ypq2Kiory8vK+ffty2HmBEqEitS0NaGjA6NFUGkDI86PyMW51T2cZS0tLWgvs2doWiyclwcyM/UUNDaVicUJelKMjHB0RGAi02VAUHk4bipRLayK8fv16SkqKh4eHu7t7XFxcfX29zC8IZH6kRJXJ3N7u74+vv6btw4TIi40N+Hzw+UCb38FNm7B4cesW/ilT4ODAdaDqpzURnj59esWKFf/5z3/c3d3feeedO3fuyPwCSoSqpbS0dOvW3f7+PGvrEczdKNN6ePhw+Phg1SrExYEaDhOiYEZGrbU20lmZhAS89x569YKPD3g8eHvDwODmzp17Fyzwb3dMOulerWuE9fX1VVVVJiYmhoaGd+7cEYlEMr/A1tZWgeEBtEb4vEQiXL2KuXNn37w5X1Pzm969MydMwKRJ8PaGpye0qd06IcpHLIZQiDNnkJKC06dRXv5SS0uQufm3hw9f8PSkM6S6rMtrhPr6+vqPagH79Okjx9CI3Ny7h3PncO4cUlORlgZ7e2hpDTQy+t3GpnduLi34EaLsNDUxbBiGDcPbbwPA1KkOaWmHdXWt33kHubkYMQLjxmHCBIwbR8ezdKcOxwUSiSQ1NVUoFLa0tLz99tsASkpKdHV1qXZG2bQrSxs2DD4+WL0a48fD0hLAxtzc3MGDB1MWJETlJCZuy83NdXJy0tREXR0uXUJ6OvbswdtvQ0sLnp7syuLo0dDT4zpWVSZ7+0RFRcUrr7ySmpoKwN7e/tatWwDefffdrKysP//8U8Eh0tRoO7W1yMxkM19qKnR1aaMSIWpHus03PR2XLmHoUHbVY9IkKGqPtwp4oe0TwcHBxcXFx48fl0gky5cvZx5ctGjRDz/8UFtba2xs3M3Bkmdhhn3p6eze9uHD4emJwED8/DOsrLgOjhCicG37vUkHiwkJWLsW2tqtg8UxY2jD4rPJSIQPHz48dOhQbGwsj8drO/5zcnISiUTFxcWurq4KjFBNtR32nT3LFpJ5e0MgwNixtGZOCGllaAgfH/j4ICQEaDNYDAt7bLBIezM6IiMRVlVVtbS0dFSt29H+QvLiOhr2/fQTqHqJENJJbQeLDx7g8mWkpCAuDu+/Dx0d9pba05NuqVvJSIQWFhb6+vqXL18eOnRo28dPnTqlqakpPaqXvDjpezQ5GefOtb5HadhHCOkWRkbsYDE0FGhzt71jR+vdto8PJk9W67tt2adPLFiwIDQ0dNCgQRqPag2Tk5PXrFnD5/PNqRfQi2nbyVoohJsbe4zDli2wtuY6OEJIj9a25Zt0/SUmBm+91br+ooaDRdlVo/fv3/fz80tLS+vdu3dtbW2fPn1KSkqcnJxOnjxpZ2en4BBVvWpUuo6dkoJTp2gdmxCijGQuzfj4YMoUFa7Ie9HTJ5qamnbv3n3s2LE7d+6YmJhMnjx55cqVhoaGcgj1GVQxEXZU2Tx5Mh6dbUwIIUqqoz1anp4qdhBN9xzDpAyUPBHu3Xvw2LHU9977e0WFA3Mz9eeftNeVENJztO3aUVzMdir29oabW/X33387YIDdqlVBXMcoW/ccw0Q6IhLhr79w5kzte+993tDw/s8/h3t5/TxuHBYuxIYNUPj8MSGEyEvblcX793H2LM6dw4YNSEmJamrS69Xr99JSzwULPEeMUKXBYluUCLugsBAXLuDCBVy8iEuXYGuL0aMNzc2NGhr+99FHy9au5To+QgiRs969MWsWZs0CgCNHPJYs+VhDQ1xS0v+NN5CfD3d3jBmDsWMxZgycnVWmvzElwqeprsaVK+zq8enT0NSEpyc8PfH++9JOnprAqerqalNTU66DJYQQhZo+nVdYOE5fX19LSwtAczOyspCcjOPHERHBtj5mrpmennBz4zrcjlEifIz0B5mejvT01h+kvz8++6zDHyRlQUKIejIyMpJ+rKPD5jxGTQ2ystjGb6Gh0NBoTYqPBhLKQt0ToUiEnBw27aWn4/Jl9O/P1rmEhmLoUOphTQghz8PEpH3jN+YyGx2N11+Hrm5rRaGnJx6dAcgNdUyE0p9HSgpSU2Fjw/4kBAKq8CSEELmws4OdHfh89lPptsWwMGRmwsGhdbyo+O38apEIq6qQlsb+Tz9/nt3S7umJVasQGwtqlUMIIQrWthK1uRl5eezgJDqa3aEhzYuurnIvuumZiZDp4Smd8Cwpgbs728CaOpkRQohS0dGBmxvc3B7r/ZaejsREREayF3BmEtXLSy49UXtIImxpQW5u64SntIcnj4fQUEXcUBBCCOkWxsbs4iKjrAxpaezi4ooVrSU5np6YOBFmZt3wHVU4EbZtY9Z2ipmObiCEkB7D1hZ8fuviIlPkkZKCyEgsXvzY4uJzd29WgRZra9euu3Il6+jRI9Iil/R09qxa5h/v44MJE2BgwHWghBBCFKjtXCDT2HnAgNZK1KFDsX371lGjRo8cOezpr6MCifDNN/8dFXXGwGCtkRHPy6u1bUHv3lxHRgghRGk8eID0dFy8iAsXcO4c7t0rfPhw3vvv7//qq2ecdaASU6MtQGnfvpV37uDSJWhoQCxGUxM8PamlJyGEEAC4dw8ZGeyfS5dQWQlHx0ahsEpHp/aZX6sCI8J169ZVnyRF/gAAIABJREFUVlZu3rwZbbYAMmm/sRFubgqtsiWEEKIMKishFLZmhHYd3VxcoKWFrKwsExOTAQMGPP2lVCARPuUYprZ5MTsb9+6x2ySYP9QXhhBCegzpBT87G0Jh6764p1zwuT+GqaSk5LffftPU1Jw/f761rL17NTU1GRkZd+/enTNnzjMDlaldqwLpDQKz+0TmDQIhhBCV0Haok5aGhgZ2CpBp/tyNU4DyGhFev37dy8tLIBA0NzcnJCSkpaX17du37ROEQqGHh4ejo2NOTk5VVdVT+lY/98G80rMjmD8FBXB07IZCW0IIIfLQNvO17QLm5gZX1+c5v4LjE+r//ve/t7S0REVFAVi2bJm9vX1ERETbJzQ1NYlEosrKSnt7ezklwnak3Qpk5kXOu74SQohaYc48YOY509ORmsp24paOVWxsXvRbcDw1evTo0Q0bNjAf8/n8L774ol0i7NWrF4DKyko5BfCkdt0KHjxAbi77A4iLaz13gvkzahQMDRUWGiGE9Hzttv1dvgwTE/aSGxSEn3+GlRU3gckrEZaVldna2jIf29ralpaWPvdLicXiI0eOSFPm0KFDlyxZ8uIR6ujA3R3u7li4EACam3HtmsalS5oZGRp79mhcuaLZr5/E1VUydKhk1Cjx+PGS3r2VvaqIEEKUStvrakaGxpUrmlZWkqFDJR4ekvfek3FdbWzs5gDEYrFWJ2pD5JUINTRaJ10lEonGi61p6unpmT3qKGcon5Gajg5cXSWuriImyba0IC+P/fl99ZVWZqamiYnEw0Pi4SEZNUrs5SWxtKS8SAghj2lqQn5+a+bLytJ0cJCMGiX28JDMny8eMaJZOWfa5JUIbWxsysvLmY/Ly8ulo8PnoKmpOWXKlBdfI+wSXV2MGoVRo7BiBQC0tCAnRyMjQyM9Hd9+q5WZCSsreHhg1CjJ3r2LS0vzv/32w0WL5ikyQkII4daFCxfnzXtHX9905cp9QqFxRgZu3oSbGzw84OGB11/HsGHQ1dUAOKvXZ9YIn/k0eSVCPz+/+Ph4f39/APHx8X5+fszjubm5Dg4O+qpWl6Ktzc6jMqeEiMXIy0NGBs6erbl69WZz8zfLl2/evn3eyJEYPhwjRsDZGdoq0bSHEEI6TSxGQQEyM5GVhcxM/Pnn8aqqN7W1T2ZkZM+Y4bV2LVxd8Vxb4Tgmr6v1e++95+XlpaGh0dzcnJiYePHiReZxd3f3pKSkSZMmicXiV199tb6+HsDrr79uaGj4yy+/yCmYbqepCRcXuLhg8WJTa+tXjx7d+Omn/9DQgFCII0fw73+3lqQyVb9jx9IhiIQQ1VNbi7y81gYumZkwNmYva/PmISTktS++CHV07Pvf/45W6V3acuwsU1ZWxmyonzdvXp9HZykeOnTI29vbwsJCIpHs3btX+mQtLa1582RPLXbX9glFamrCtWutTRAyM9HSAlfX1j0x7u60i5EQonSYnXzSLQ1tt5m5uWHECM4KO58Px/sIu5EqJsInlZa2vrfavr2Yeys5HbtMCCFP0e6W/fJl9rB46V27qjfk4r7FGmmL6QbH47Gftn3/JSbi0iX06tWj3n+EECX0lDtyHk9978gpEXKDSXtubmz1DdrMSDCNUlV9RoIQwrmnrNHweAgJoTUaFiVCZdGugXjbNeq4uMfWqOlsDUKITE8f8FHVXkcoESopY2M24XVyyDhqFCwsOI2YEKJYzc3Iy2utbblwAc3NNOB7HpQIVUa7IWNNDa5de6xXKtO1TzpqZIaMEonkypUrLi4uTHNXQogKEYvFV65ccXNz09bWxuNH0WZnIycH/fuzv/JBQYiOxgt0LlFrVDXaQ4hEyM9nd7ky212rqjBsGG7ffru8/IGFxY3Ll8886lJHCFEBdXXg8RYLhfoGBqXu7n9cvgwtLYwYgREj2MYdQ4eq5O51RaKqUfWipQVnZzg7QyBgH6msRGYmVqworaub2dSU5eAAU1O4usLdvfW/xsacBk0IeaShAX/9hexsXL0KoRBCIcrLIZHcbmhYpqPzw7p1GD68G44lIjJRIuyxzM0xZQrOno3auTNuxozdQ4e2LqRfuICtWx8rwGH+O3IkjIy4jpsQNdDcjOJiCIWPbV23s2PXNZYtg6srXFxQVrZ9z57f5szZNWgQ1xH3aDQ1qtakqZH5L7PQ2DY1enjAwIDrKAlRcS0tKCpq/UXLzkZuLvr0af1FYzZT6elxHWiPQ1Oj5NnabfNHm9rU5GRERyM7G2ZmrTU4zH9VrWU6IYrWtlFZdvZjv0dMPSf9HikVSoTkMe1qU9veySYmYtMmupMlpL2nzKz4+CAoiGZWlB0lQvI02tpwdISjY2tqbLu2Id3RyKxtSHduUH840oO1S3vtml0sW0Zr7aqHEiHpGh0dGakxL4+9KCQktG72bzubSn1wiIpitu5JS1qysqCt3frGXrYMI0ZQ9bXKo0RIXhTTrt7NrXXnBtPhkLl2xMVBKHxaaiwuLtbX17e0tOTwn0DUVl1dXWlp6ZAhQ5hP26a97GxkZbG9Wpg3LZ+P4cPVtC11z0aJkHQ/aUtxaWqsrWU3SGVnY8sWCIWorcXQoejdO/H06X9qa9ft37/7pZecNTQ4jZuomWvXaidMmFhfbz9ihL+OzltXr0JDo3Wj7bx5cHenzoVqgRIhUQRjY3h5wcur9ZGqKgiF2Ly5tLHRrbGxbPHiu7W1zk5OcHaGiwtcXNj+AFRiQLoFM0uRk4O8POTkICcHubno1auuulpXJBqnpVX86adwd6fRnpqiREi4YWYGb294eS329paYm5ssWOAjnVAtKMDvv+Orr5CTg1692IkpR0f2g4EDQQNH8nSVlSgoYKc3mQ+Yki7mXeTtjWXLMGwYrK1t/vzzy4yMq3/7WyCVt6gz2lBPlBpTodf2olZayi43SlPj0KE0cFRfzA6fdmmvvh6DBrW+Q5g3DG3yUUO0oZ70BMy+xraqq5Gfz17vEhKwaRO7W7ltamTqWknP0/anz6S9v/6CqWnrT18goJ8+6TJKhETFmJqypzBKK3HajQkSEiAUoqGBvSDSmEB1PX0+wN8fbm5wcYGhIdeBEhVHiZCoPOmu/7a94tquEkm3cEhXiZgrKZXCK4+2K8TSbpw6Oq0/LB4Pbm4YMIA2pJLuR4mQ9Ezm5uzAUaptT5z0dMTE4MoVaGk9lhqZlv/UFkfepNv1pGmvtLS1PxHTlow2qhOFoURI1MWTPXEAlJSwlfQ5OUhKQm4uKiog3cXh7Cz+/vvXiovzo6L+NWOGL3exqyqhUPjKK2/q6RmHhf1aVGQi3begr4+hQ9kdMlOmwNkZ/fvTUI9whhIhUWv29rC3x9SprY/U1SE3l02Nv/5658KF0paWyHnzfhk71nfwYAwejCFDwHxABfftNDTg+nXk57N/rl1DevqRqqrXNTXToqMveXtPfvll/P3vcHaGmRnXsRLSBiVCQh5jaAgPD3h4MJ/ZrFkz8cyZyK+//reRETuJFxuLggLk5bWuTUonVwcPhqkpt+ErSFMTbt1CQUH7ShZmFZb58/LLsLR89bPPQmxtrWJixunqch00IR2gfYSEPCemHqdtJmiXHZkEOWyYamdHac6T/jOfzHnMfUD//rS8SpQL7SMkRL6k9TjSjRxokx0LCtjDja9dY0ty2v5xd4eNDXehd6CxESUlz8h5zKYFynmkJ6FESEh3erJaFbKyY34+NDXbZ0c3N9jaKihOJue1TXjtcp50czrlPNLjUSIkRO7+v707j2rqzhcA/k2AkI2EJYBhBwXROlXHunEUjgtTW+20Vh2sbUXPadFWbR077bT2nefR6jzX1seMPGvVDm1dpjh9bV0oZam4oKL0tVYsopUlEpJAViCEQHLfH7/rJSRhEY1J5Ps5Off8+N2bm5uQ/L73t9zf7Tc6kss5rl1zEh3vf56Ujg64datHwOst5uFVemhowkCIkHsMJDrm5cG1a2Ay0TNn2j7I5OMqlSo7O3vp0qVJSUlwd9xmbzGP3EIdYx5CdnCwDEKerrm5+4IE5uKEri6IiIBbt2Z3dU1lsb4dPfrnxkYwm8HuGo8RI+wna0Vo6MDBMgg9CpqboaYGZDKQy0GlAo0G9Hpoa4PAQPDzA4BWAAqgw9cXeDxoawODATQaUKlAIAB/f/D1BRbr4XU9IuSNMBAi5BFMJpDLezRp3r5Nj6mRSrv782bOhIQESEqipx9Tqb7Nzs7OzDyemAhwdxo5uRwaG+H2bThzBv75TycjPxMSuveJEMJAiNBDZXspOglRJGjZxSrSkzd8eD+TsISFhW3evJn5k5lGzg5zXQR55OXRifZ2+wCZkACJiSASueCdI+SpMBAi5Cq2I18ch66QOtm0aT0Gv7iIv7/zANneTodh2wBpNy0AU3fEGyCjRxUGQoTul2PAk8uhpga43O5wwgzX9KjL8ng85wHS9h2RWzw6viPm4VHvCKFBwECIUF9KS89cvPh/K1ZkBgYGarXdLZnMw3HS0dmzvb7+5PTSDnB2dQdTx7XtyCSPuDjo6jJ//PE/pdLQhQvnu+N9IDQgGAgR6sFqBYUC6urgzh24cUOzZcu6jo6M7ds3WiwfsVgQFwfx8RAfD3FxMG0anfbegHevnAbIjg6oq4OaGqithZoaKC6G/fuhthba2kAkOqBSVfv5fVNaGjVlysTYWIiOhogI6G80O0IPFQZCNERpNHDnDtTXQ309yGRw5w7U1dFDLkNCICYGoqNh2DA+j2f18zu9ZMmTmzZBUJC7D9oj+ftDUhIkJdnnt7ZCbm7s+vVHKKqlszPs5En601YqITQUYmMhKgqioyEmBmJi6HR4uDveABry8IJ69CgjoyVt2zNJ+tYt6OjoMWiFGRgSFwcCQfceTCZTQ0PD8OHD3fcmvFtDQ4NQKBT3vAEH08RqO2i2sdG+lZVJ9Dt6FiGn8IJ6NIQMsGBlOvAGXrByuVyMgvcjMjLSMbO3PkjHE5fz5+nrKQd44oLQIGAgRF7DdqyKbVlZXw8BAT2KyEWL6DTOqOldervMA3pe6SGX00N15HKoraWf5ViPxOGsaIAwECJ3oijq9OnTQqFw4sSJJIdMsOLYmHnzJpjNPaoCzAUJ0dE4+OLR19uVHnC3PYD5zlRUdLcHBAX1+M7YRkqiurq6urr6ySef7Lf1DD3CMBAiNyB1O7kcvv762IED/2Kx1JMm7dTrJ8hk0N7eYxjF5Ml0OjZ2CA3ORPekj4ZWMgyKDNK5fh2+/54eEkVREBMDUqn6/PnFXV1TZ8/+afXq/4iKgvBwHLAzFGEgRC7R0gJ37oBCAQ0N9LKxkV7K5eDvDxEREBEBXV1+LJbZx6fz+ed9UlMhOhokEncfOnpU+PvT999wpNeDTAaVlexLl6xWq0ku992zB2QyUCpBp4PwcIiKgmHD6GVkJEil9DI4+KG/DeR6GAjRIJlMoNHQgY1ZkhYquRxMpu5RKmQ5fjydjooCmyGEz128OCwgIOCxxx5z43tBQ41YDGIxjBkTNHXqt7du3UpLS2N6E81maG7u8cUuK+tOa7V0c6vtd9t2ibwRBkLUK8cSwXbpWCJMmADz5tHdMAO/5G7KlCmufBMI9SUmJiYmJsY2h8Ohmysc21oBoKMD1OoePwTbLknmR8H0RDLLnud/yLNgIBzqHIdiDiTUSaUglbpwkmiEPBPTqu80TNoO9SLLyko63dBAXwHitCoZE0PfVwu5BQbCR9zOnZ/s3Llv/Pipr7+eTX6ZTHedQgFqNYSFdXeHhIfD1KndHSShoe4+eoS8CjMpuVM6HR0mmUh55Ur3T1IopHsimb7JqCg4fnzXiRNHnntuzscfb3a+U/QgYCD0bm1toFCAUglNTaBUdicUCmhqIjc0/xdFHSoqWujrS5/JTprU3fMfHo61OoQeksBACAyE0aOdr21uBoWCHmJ25w7cuAFFRfDdd192dHz9ySfzjh7dLJVCaCiEhoJtIiwMQkNh2DBsd70vGAg9GrmImLRSOibImJSgoO4GzKAgiIyEJ57ozqmq2rB+/duZmW+99pq73wxCqHcSCUgkMGZMj8zjxzdt2bJq9er/WLDAvhCoroZLl7pzSEeGbVFgm46IwMtt+4JzjboN6XW3/XLbfdFlMvDx6fFVtkuQNEIIkVHcTs+YSaKpCUizUG/lSWTkIzihK841+pDodLrq6uqJEyeyejYytrf39aXUakGpBLHY/rtIRqOQnKgo4HDc9bYQQt6Ey6X7Pvq4EMmxhYnM5kpyyHCePuqUpDPFbta6q1evSiSSCC8/JcdAOBgaDajVoFaDUmlevjyto2Pc6NExY8Z8QLrlSP+cry/dlB8WRk9XkZhIJ0JDITwcr8xFCD1UfUxTR7S0QGMjNDVBUxM0NoJKBTIZXLnSPfigtbVHP6VKdejs2aM+Pg27d3+dnBwTHAwhIRAS4n0T/GIg7EbaKpkgp1ZDc3N32jY/IABCQyEkBMTirvZ26OxM6OpST59Od1yTrwiP5+73gxBC9yIgAAICnNxakmE29xiOd+yYzmyO8vHRf/650bb8DAqiI2JICDDRkTwkku58z5k0cUgEQtJKyTxs++RsHwoFfZEQ0yxAHtHRPZoIJBLbFkt+ZeXh8vKKRYveEgrd+BYRQsjlOByIjATmzlovvrgiL++r2NiXU1KSbTdjilzbwraqqkd5azvWz+5hVwiHhYGviyOVFwyWeeedd5RKZW5uruMq2yZvxwdZpVYDm93PB00ew4Z5X40eIYS8lF0VpbeKikoFfn59Fd0kPzIS/P3tX+LSpUuhoaEJfTQHA4BX1AhbWgSffXb9xx8LJJIn29uhtRUMBtDpwGh0Xvv+3e/s83HICUIIeRoeD3i8/oe+WyzOu6iuXbPP9PWFoCAQiUAoBB4PurrqL1zI/Mtf/nf79n5ewgsCIYvFBggwGHR+fmA2g8kEFgsAgMUCnZ1AQiMZyGSxQEcHGI2g14NaTV++GhgIYjG9xIiIEEJewWIBnQ70etBqQa+n00yC+bO1FYxG6OgAiwX8/OgowGKB2QxmcxtF+VmtDvVEB14QCPl8/TPPCL/9NsNxlWPN2mSiM+vr7VdpNMBiAY8HXK7zhmnysN3AcaywUyaTicvlPvh3jhBCnq2jo4PD4bAGMENVe3t3+dzHg9lGqQQOx3kRLRZDXJx9fkiIY9PoqLy8/xw9ur3fY/OCQOjr69vbDQoGWLNmGI32ZxM6Hf24ccPJGYfRCGIxBAXRFUqmZskkRCLre+/NUasNa9asePXV5QIBBAbipGUIoUecXg9tbVBU9MO6dX/lcNjbt5+yWILtilBSmWPSfn5OitCgIBCLITm5RyZZPpBZyBcsWGAhTYh98oJA+ADx+cDn30PgtFjofyqpm9v+g+vrQa+H5ubWujptZ+f6Xbu+OXBgudEIWi3weMDng1gMQiHw+SAUglhMv3RgIAgEwOdDQACIRMDnA4mdtpkDqYMihNCgURTodNDWBkYjtLSAwQBGI5DiiyQMBmhpAaMR2troKgHpciLtkK2tdJlmMl3S6Rb7+JTt3387NjaYBLCICBg9mo5nTOeUh/dMuTAQ/vrrr2VlZVFRUenp6WxnwzE1Gk1+fj6bzX766afFHjllrI8PBAf3feW7KDf3rYKCki1bNsTHd+cy1X+mmm/bJtDURLfcOq7VaKCjg26btW2kZdJ9Z/J49jcCrK2tzcn5bMGCOZMnT3LNh4QQcpXq6upPPjn60kvPjh071m6VXdFhW8I4Ztqt1euBw+mnMJFKnZcwXG739QwaTdY77/xXQsL0996b4NUtYa66fOLLL79ctWrVokWLysrKkpOTjx49areBTCabPHlyampqZ2dnRUXFpUuXwsPDne7qUZ1rtDcmE92EO7jzNYAeddDq6rkGwzJ//w9efPGqry9d4yRLsRjYbHpJWnTJkoTSoCA6ByF0/wwGsFjopV4PViu91Ono+hlFgVYLAKDVduccPZra0vIXLnfT449fsS0EOjv7anMSCPppiPKcK9ldzZ1zjVIU9f777+/du3fBggUGgyEhIeHKlStPPPGE7Ta7d++ePXv2Z599BgALFy7MycnZuHGjKw7G63C5wOUOfgK2zk5obQW9no6OGzY8dv78ZyEh0pQU+99hbW3/v0Odjg6HTMgUi7tDqUgEvr4QEEAv/fxAKKSXHA4IBPTS3x/4fHrJ5dInlWTmnebm5uXL3w4MFO3fv93f8SIghFzsr3/dcunSz/v2bU5KSgIAsxna2uglGYJOlqRGRZZMgpyzMpuRJ7a20r/Bzk5oaYGuLnpJfiYDPxMdPhxYLLh6Nbmq6tPY2Njs7O6eFPI7Qg+QS2qEVVVVY8eONRgMpGhbvHjxqFGjNmzYYLtNcnLytm3bnn32WQA4cuTIzp07KyoqnO5tqNUIHziZTBYVFTWQYV1OMaeuen13KDUYun/hLS3dv/zW1n6KEttChMcDFut/jEYTm307LGyuRDKHhEISPgF6VElJOQLQHURJgUKQ0gSAjr6D24lTOTmflpZe2bbt7bi4uMF9gC5lNpsvXLgwYcIEoUfObKTT6d56a/OwYZIPPnjHaf8IQb4wcDcOwd1OLIJ85eDusEMA+lt39yWAlGGD2IlWW9vQsNJiWc3hFHI4/93aeg+ncUyCy+3ewN+ffiJzRujn1/2tG5z6+vqYmJjBP39oc2eNUC6XSyQS5gQ/MjJSLpfbbdPQ0BB5d6KeiIgIxw0YVqu1tLSU+TMxMfG555570If8KBs2bFhXV9f97IGUsX0HjEFob4dffpn+wgtZ/v6cL754n8PpcizI9Ho6fpOgCz0LMrWaXvvbb05KQ6ORBT1LQ2YnJhPLsUhloimfT5Fo6uOjqq090NX1ZmHhtvj4HObT6ONnJRRSfZR6zEs4Zde/a4vNBpHIyTnrp58ubWiQCoXvr11b6ri2rY3V2dnrPpnY4BQTYBzZ/l8ckfMkQqn8VC4PY7N/OXToHJs9HQC6uqC1lWX3EranL3w+Bb2fvnC5FPQ8fWHGafP5FMl0eg7E/F9sz4H4/PClS7l37mzavXvTk092uuhcgqKgj/9Cv6RSaef9PH9os1qtA6kDuCQQWiwW29dms9mOBbHt8fn4+PRdUhuNRo1GQ9IGg2Egw2GR5+NwYMKEEdXVJXcz7ita3z+dDgBYYBNN29uFixdzm5u3r17956eeogsjUvftTUsLq4+vJ6lb934Avf5irVa4+wvowWhsM5ujOzrOazROopZQSPXRGxQVRfUxRLmPC4FYLBCLe21Jso1SlZUT3313LZfre+DA8IgIM9jX4yl3j7BgnTvHDF/AcuURZLVafQYwEN8lgVAqlarVaovFQo5AqVTGxsY6bqNSqUhaoVD0cTsrNpv91FNPYdMocrVhwxzz/Kuri8xmM8dTh36vX5974sTJmTO/ioryxEuhpk6dlpl5wdfXd9At8wjdD9I02u9mLplkOjk5OTg4mLRnms3mkpKSGTNmkHRTUxPZZsaMGd999x1JFxQUzJw50xVHgtD989goCADBwcFLl74cFRXl7gPplZ+fH0ZB5OFcchbp6+v77rvvLl++fPXq1SUlJQkJCSQQFhcXL1myRKvVAsCf//znlJQUPp/f2dn51VdfXb582RVHghBCCPXNVbcdWrNmzf79+/V6/fz58wsKCsgp4eOPP75nzx6ywejRo69cuSIUCkNCQioqKvq9TQZCCCHkCl5wP8Jp06YFBgaeOHHC3QeCEELIm5w8efLQoUOHDx/uezMvuBGtxWK5z9H/CCGEhiCz2dxGLsbqkxcEQoQQQsh1MBAihBAa0jzx2iM7Uqn00qVL6enp7j4QhBBC3kSr1TreuMORFwyWsVqtJSUl/W+HEEII9TR58uSA/m7y6wWBECGEEHId7CNECCE0pGEgRAghNKRhIEQIITSkYSBECCE0pHl0ILxz505+fn5dXZ1tZkFBgVKpdNchIYQQ8mQymezUqVMymcw289SpU8y9jxx5dCDMz8/ft2/fmTNnmJzq6uoPPvggNDTUjUeFEELIY508eXLv3r1lZWVMzvXr17dt2yaRSHp7ikdfUP/qq6+q1WrbnHXr1u3atYvN9uj4jRBCyF1WrlzJ3PWdWLt27e7du/u4L6Y3RZRjx45FRUVNnjzZ3QeCEELIOxw5cmTkyJHjx4/vYxuPrhGqVCqNRsPhcNRqNZfL/dvf/nbs2LG9e/c2NjbOnTt30qRJ7j5AhBBCnkWpVGq1WhI+/Pz8tm/f/tVXX+Xk5KhUqmeeeWbChAmOT/GsQNjc3FxaWgoAAoFgzpw5JSUlbDZbpVJdvHjx3Llzr732GofDGTNmzMSJE5csWXLx4sWgoCB3HzJCCCEPUlRU5O/v39DQUF5eXlRU9Oabb/r4+IwdO5bD4SxevPjHH390nHHNs6ZYu3HjxtatWwEgNDR0+/bttvmvvPJKaWkp6R2UyWR//OMfy8vL/fz83HasCCGEPFhlZeWqVat++OEH0jtYU1OzcOHC8vJyHx8fuy09KxD2Zu7cuRs2bCBtoQqFIjMzc8+ePSNGjHD3cSGEEPJQ6enpO3bsGDduHADIZLLMzMz9+/cnJCQ4bukFg2Xy8vJiYmJIFJTL5enp6c8//3xtbW1ra6u7Dw0hhJAnOnz48KhRo5goOGfOnCVLlty+fdvpDeu9IBAqlcotW7aQtNVqXbNmDYvFun37dkdHh3sPDCGEkGdqbm7etGkTSVut1jfffNNqtd6+fbuzs9NxY+9oGkUIIYRcxAtqhAghhJDrYCBECCE0pGEgRAghNKRhIEQIITSkYSBECCE0pGEgRAghNKRhIEQIITSkYSBECNm7cePGvn372tvb3X0gCD0MGAgRQvbOnTu3YsWKlpYWdx8IQg8DBkKEXMhkMqnV6kE/V6FQPJCpBI1Go0KhcDq5FAC0t7f3sfZej6qjo0OhUFgslkEeK0IPHQajqkD3AAAHRUlEQVRChO7BY4899u6775K0yWSSSqVxcXFMob9q1aopU6aQ9EcffZSUlMTj8SQSSVBQUFZWltFoJKsWLFgwbdo02912dXUlJSWtWbOG/KlQKDIyMgIDA6VSqVgsXr58udMp5imKGjdu3LJly2wzVSpVaGjo7t27yZ+//vrrH/7wB5FIJJVKQ0JC3nvvPdsQVVVV9fTTT5O1XC43LS3NYDBkZ2e/8cYbADBy5Mjg4ODg4GCFQgEAGo3mhRdeIEcVGBi4bNkyg8FA9mMwGIKDg/fs2fPqq6+KxWKpVFpdXT3Yzxihh45CCA3YSy+9NGbMGJIuLi729fVlsVjl5eUkJzY29vXXXyfp9evXHzhw4PLly9euXcvOzhYIBCtXriSrcnNzAeCnn35idnvixAkAKCkpoSiqpaVl5MiRI0aMyMvLq6ysPHToUGho6LPPPuv0eN5//30ul6vRaJicDz/8kM1m19XVURRVV1cXEhIyadKk7777rrKyMjs7m8fjvf3222TLuro6iUQSGxt7+PDhysrKkpKSN954o7m5uba29q233gKAvLy8wsLCwsJCk8lksVhSUlKEQuG+ffuuXr3697//ncvlpqenk13pdDoACAsLmz9/fkFBQWFhoVqtfkAfOUIuh4EQoXtw8OBBFovV2NhIUdT69eunTZv2+OOPb926laKomzdvAsC///1vp0/cuHGjSCQi6dbW1oCAgHXr1jFrFy1aFBsba7FYKIratWuXj49PVVUVs/bIkSMAcOPGDcfdVldXs1isvXv3Mjljx45l4lNWVpZEIrENkxs3buRyuUajkaKolStXcjicW7duOe52//79AKBUKpmcgoICALB9oR07dgDA2bNnqbuBcMyYMeQtIORdsGkUoXswa9YsiqJKSkoAoLi4eNasWbNnzy4uLiZ/stns1NRUZuPTp09v2rTp9ddfX7FixdmzZw0Gg0qlAgCBQLBo0aIvvviCdMvp9foTJ04sW7aMzWYDQEFBQVRUlEwmK7qL3F/72rVrjseTmJiYkpJCqphkm59//jkzM5P8+f333yclJVVUVDC7EggEJpPp1q1bAFBUVJSWljZ8+PCBvPGffvoJADIyMpgckj5z5gyTM3/+fPIWEPIuvu4+AIS8SUxMzIgRI4qLi5955pmKioodO3a0tLTk5OS0t7cXFxePHz9eIpEAAEVR8+fPz8/PT09Pj4+PF4lEQqEQAPR6fVhYGABkZmYePHiwoKBg3rx5hw8fNplML7/8MnkJpVIpl8v/9Kc/2b5uUFBQb4NuMjMzs7KyqqqqkpOTDx48KBKJ5s+fz+xKoVA47qq5uRkAmpqapk+fPsA3Xl9fz+PxAgMDmRypVMpisWyPKjw8fIB7Q8ijYCBE6N7MmjUrPz//hx9+8Pf3nzx5stlstlgsZWVlp0+fXr58Odnm8uXL33zzTV5e3sKFC0lOTk7O119/zexk+vTpw4cPz83NnTdvXm5ubmpqKlMzE4vFo0ePJjWwgcjIyFi7du3nn3++cePGI0eOZGRk8Pl8skokEqWmpn755ZdOnxgYGKhUKgf4KgKBoL29va2tTSAQkJzm5maKosRiMbMNqbki5HWwHQOhezNr1qz6+vqPP/44LS2Nw+EIhcJJkyZ9+OGHTU1Ns2bNItvU1tYCwIQJE5hnnTp1ynYnLBZr6dKl33777YULF8rLy5nGTABITU2trKy8fv36AI+HVAFzc3NPnDihUChsd5WWllZSUtJbVTItLe3MmTNkRKidgIAAALC9oJ6MhiWDeojjx48DwNSpUwd4nAh5Ljf3USLkbZqamkhP2M6dO0nOhg0bAIDD4bS2tpKc69evs9nsV155RavVKpXK9evXk6bR6upqZj81NTVsNjsuLk4gEBgMBiZfpVJFREQkJiaeOnVKp9M1NTWdOXMmKyvLdsyLncLCQgCIi4tLTEy0Wq1M/rVr1/h8/pQpU86dO9fW1iaXy/Pz87OyspiD5PP5TzzxRFlZWVtbW21t7UcffaTVaimK+uWXX1gs1tq1a8+fP3/lyhWz2Ww2m0eNGiWVSgsLCw0Gw/Hjx0NCQn7/+9+T0TFksMyePXse2KeM0EOEgRChezZu3Diwuf6htLQUANLS0my32bFjB4fDIaebKSkp//jHP+wCIUVRM2bMAIClS5fa7f/mzZszZ85kzlZ9fX1nzpzZ0tLS2/FYLJaYmBgA2Lx5s92qCxcukKMleDxeRkYGs/bs2bPJycnM2oSEBBIIKYraunVrdHS0j48PADQ0NFAUVVNTk5KSwmw8Y8YMkk9hIERejkVRlGuqmggNdTqd7rfffgsODo6Pjx/E05VKZV1dnUAgiI2NJRXKQauvr1coFCKRKD4+3t/f327tzZs3tVptaGjoQI6ztrZWpVJJpdLo6Oj7OSSEPAcGQoQQQkMaDpZBCCE0pGEgRAghNKRhIEQIITSkYSBECCE0pGEgRAghNKRhIEQIITSk/T+ZcdYwH8IHDAAAAABJRU5ErkJggg==", "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", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\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", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\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": 4 } ], "cell_type": "code", "source": [ "using Unitful\n", "using UnitfulAtomic\n", "using Plots\n", "\n", "plot_bandstructure(basis; n_bands=6, kline_density=30)" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "!!! note \"Selection of k-point grids in `PlaneWaveBasis` construction\"\n", " You might wonder why we only selected a single ``k``-point (clearly a very crude\n", " and inaccurate approximation). In this example the `kgrid` parameter specified\n", " in the construction of the `PlaneWaveBasis`\n", " is not actually used for plotting the bands. It is only used when solving more\n", " involved models like density-functional theory (DFT) where the Hamiltonian is\n", " non-linear. In these cases before plotting the bands the self-consistent field\n", " equations (SCF) need to be solved first. This is typically done on\n", " a different ``k``-point grid than the grid used for the bands later on.\n", " In our case we don't need this extra step and therefore the `kgrid` value passed\n", " to `PlaneWaveBasis` is actually arbitrary." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Adding potentials\n", "So far so good. But free electrons are actually a little boring,\n", "so let's add a potential interacting with the electrons.\n", "\n", "- The modified problem we will look at consists of diagonalizing\n", " $$\n", " H_k = \\frac12 (-i \\nabla + k)^2 + V\n", " $$\n", " for all ``k \\in B`` with a periodic potential ``V`` interacting with the electrons.\n", "\n", "- A number of \"standard\" potentials are readily implemented in DFTK and\n", " can be assembled using the `terms` kwarg of the model.\n", " This allows to seamlessly construct\n", "\n", " * density-functial theory (DFT) models for treating electronic structures\n", " (see the Tutorial).\n", " * Gross-Pitaevskii models for bosonic systems\n", " (see Gross-Pitaevskii equation in one dimension)\n", " * even some more unusual cases like anyonic models.\n", "\n", "In this tutorial we will go a little more low-level and directly provide\n", "an analytic potential describing the interaction with the electrons to DFTK.\n", "\n", "First we define a custom Gaussian potential as a new \"element\" inside DFTK:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "struct ElementGaussian <: DFTK.Element\n", " α # Prefactor\n", " L # Extend\n", "end\n", "\n", "# Some default values\n", "ElementGaussian() = ElementGaussian(0.3, 10.0)\n", "\n", "# Real-space representation of a Gaussian\n", "function DFTK.local_potential_real(el::ElementGaussian, r::Real)\n", " -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)\n", "end\n", "\n", "# Fourier-space representation of the Gaussian\n", "function DFTK.local_potential_fourier(el::ElementGaussian, q::Real)\n", " # = ∫ -α exp(-(r/L)^2 exp(-ir⋅q) dr\n", " -el.α * exp(- (q * el.L)^2 / 2)\n", "end" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "A single potential looks like:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=1}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXwTdf4/8PdM0ruFttArvWhLKRQEuQqVo/TgLggiKKKiICq4Lt547P5WWTxQF7/oekBBEBVUwAO5KS0gVzmF5T5KzyQthZY2aXpk5vP7I2u2lhZammSSzOv58I9k+CR5x5nOK5/PZw6OMUYAAAByxUtdAAAAgJQQhAAAIGsIQgAAkDUEIQAAyBqCEAAAZA1BCAAAsoYgBAAAWUMQAgCArCEIAQBA1hCEAAAga04YhGVlZXl5eVJXIXeCIEhdgtxhFUhOFEVcw1Jyoijeto0TBuFPP/20YMECqauQu+rqaqlLkDusAsnV1dXh54i0RFE0GAy3beaEQQgAANByCEIAAJA1BCEAAMgaghAAAGRNaZF30el0P/zwQ0VFxejRo7t169Zkmy1btpw8ebJXr16jRo0yLywtLV23bp0gCBMnTgwLCzMvz87OPnz4cJcuXcaPH8/z/03r8vLydevWVVdXp6enx8TEWKRyAACQOQv0CGtqahITE3/66SetVnvPPfdkZ2ff3Oall1568cUXdTrd888//8orr5gWqtXqXr16HTly5Ny5c7169bp48aJp+XvvvTdjxoyqqqr58+fPnDnTtLCioqJfv36ZmZn5+fn9+vU7evRo2ysHAAAg1mYrV67s06ePIAiMsY8//jg5OblRA61W6+7unpeXxxi7cuWKh4dHaWkpY+z1119/4IEHTG2eeeaZp59+mjGm0+nat29/9OhRxlhZWZmnp+fFixcZY4sWLUpLSzM1fvPNNydNmtRcPUuXLp05c2bbvxe0RWVlpdQlyB1WgeQMBkN9fb3UVciaIAg6ne62zSzQI9y2bVt6erppADM9PX3Xrl21tbUNG+zatSsuLi4yMpKIOnXqFBsbu2vXLtMLx40bZ2ozbty4bdu2EVFOTo63t3efPn2IqEOHDgMHDtyxY4f5Uxo1BgAAaCMLzBGq1eohQ4aYHqtUKsaYRqPp1KlTwwYhISHmpyEhIWq1moiKi4vNy00LGWPNNVar1cHBweaFOp2usrKyXbt2TZa0hes14P1M02NPT89+/fopOPKxzHzon3i7kEsrf0vwHLVzYUTkqeTc+IZPyU1BPFE7FyIiLxdydeQjmWpra11dXaWuQtawCiRXW1srCIKDnlN/o45ERkSkM5KRERHp66lOJCKqEbgagYioVqRqIyOiepH0Ru7Wb1hRR9a7zI6Rkc7YxHLG2FOdWazXbV5ugXDgOI798f1MDziOa66BqY2pQaMX3rzwFo1v/pSG3ASDr8t/H3u5cr4uJDCqqG/b92xKUTXV3/7yPX8iMqqs54lIb2R1YsOnVCeQSFRZT/THNsdz1N6ViMhTQa58w6fMVUEKjtq5cEQU5MEivCjSi0V4U6Q319ENV3UCgD+5UUdF1VSg5wr1VKSnQj0V6KmkhhPE/+12iKj8j+G89q7Ec0RE3kpSckQNfp27K5i7gojIjSdPJU9ELjx5KW+z2/F1peb32W2l5Mi8z2+IMVK24FMtEIQhISFardb0WKPRcBxn7rrd3ICItFqtSqUiIpVK1fCFISEhHMfd3Hjw4ME3f4q3t7ePj09zJaXwF5Y9/0rbv5rkBEaVdUR/pKYg/hGTRqoTSGBUWc8YI62B8qpYzjXK17H8KlYjcJ18uE7eFOnDRXpznbwp0puL9OGCPWxXeV1dnZubm+0+D26CVSA5xphSqVQqrTAY1Yx6kYr0rFBPF6/qv3jtydq6+hqBDAKZOnDuij/956Hk4vj/hpPyj7Bo7RCXiUBUebs2t23QEnPmzBk2bFjL24uiaDA01VX8MwusoZEjR3788cdvvvkmz/MbN24cNmyY6c+voKDA29vb398/OTn5sccey8/Pj4yMzMvLu3jxYlJSkumFv/7667Rp04ho48aNI0eOJKIBAwbodLpjx4716dPn2rVrBw8eXLJkianxxo0bn3vuuYaNnZ6CIz83IiI/t+Z+1DSxvEYgdTXLraTcKpZbxU5dJ41BzK0kdTVTeXLR7SjEg1N5UbQPF+3DRbejSG9OYbVfagBgceW1pK5mmur//o2r9aQxsNxKytexdq4U7cOF1Fy7uH97xpLPpa7Ukr7++usjR460KghbyAJBOGXKlH/961/jx4+Pj49ftmzZ+vXrTcunTp06ceLEl156KTAw8Jlnnhk1atSkSZPWr1//7LPPBgQEENGcOXP69ev3+OOPe3h4/PDDDwcOHCAiLy+vV199ddKkSdOmTdu0adMDDzzQuXNnInr88cc//fTTKVOmhIaGrly5MjMzs+2VOyt3hSnkqFFMVhspT8fyqihfx/J1LEvNVujEfB27XksR3lyUN6Wo+Hsjua6+SEUAO3Klim3IZyevs0I9K9RTvo55Kijcm4vw5iK9KdyLu9ufwr35Tt4U7PnfH7VFRYpEL4/JkydLXbslHTp0yErvbIEg9PDwOHDgwLp168rLyw8cOBAXF2da/t5775kPe/nwww+HDx/+n//8Z/HixSNGjDAtDA0NPXny5I8//mg0Gk+ePGkaLyWiV199NTEx8fDhw/PnzzcfKerr63vkyJH169dXV1cfO3YsKiqq7ZXLjaeS4n25eF9qFJA1AuXr2KVK2lIojtgieijp3kju3kg+MZDjkYkAUmBEx8rYz/nihnxWamDjIvl7grhwbz7ciyK9OU/bjbbKwp+OTHEOGRkZOTk5y5Ytk7oQR3W6nG0sYL8WiOdvsNFh/LhIblQY79PURPQtVFVV3WISF2wAq0ByNTU1rZ0jFBgdKGFrr4g/5TEXntIjuHER/LAQTtnKqbuioqLExMTCwsLWvcy+vfzyy0FBQS+99FLLX2K6DZOX120OG8XvCmisux/X3Y+b14sv0LGtRWzVRXHWb0L/AC49nJ8SzYd4Sl0fgNPRGylLLa7NZZsKxSgfLj2C2zpaEY9JCltBEEKzIry5J7tyT3bly2tpS5H4Sz5781h9T3/u//VRpKrwJwpgAcV69s/j4ne54tBgfnwk9+EAl0AbHt0NJghCuD0/N3oohn8ohupExc954tN7hbj29P4A/GIFuHPXa+m9E8KX58UnuvK5D7j441QX6TjyxUvA5lx5mhLNn7lfOT6ST91kfHSXoDVIXROAo9EbaeEJseva+ht1dHKS8r3+CqSgtNAjhFZz4enJrvzkKH7hSaHHuvonuvJ/u1vh3cqjaQBkqF6kFRfEt46Jg4K4feOUse1lPaZy6dKlvXv3XrhwYdCgQWPHjpWwEgQh3CE/N3qvv+KZeP6Nw2K3dca/9+ZnxvE4MR+gSSKj9Xnia4fFKB/aOFLRuwP+VOijjz5Sq9WXL1+uq6tDEIIDC/fiVg1THL7KXswRPj4tvp+gGBOOv3CAP8ksZi/lCJ5KWjZEMSxEdn8g+/fv5zguMTHR9PTgwYNGo3Hw4MGffvopET355JOSVkeEOUKwiP4B3O505Zt9+L8eEMZtN56tcLaTUwHuTKGeDdlofPmQ8G5/xf7xShmmIBHpdLoZM2aYz1mfPXt2eXm5tCU1gh4hWAZHdH8UPz6S//SMOGyTcUKocvEQMl2iHkCedmq4WQfoxbv4F+6yi2s0Td8tqKut/iN1dBj/wl1/6mINHz5cFMW9e/cOGTJk//79ZWVlY8aMsXYZrYIgBEty5en5Hvz0WP6p3ULSRuOPaYpQL3vYAwDYFCN6+7j4+Vn+2yQuOdReBt6e7Mq34E4MbRV601VcOI6bNWvW0qVLhwwZkpGR8cQTTygU9vUbGUEIlufvRl8m1n+Z79rvZ+P3qcqhwchCkJHKenpst1BWw/aOFsKtcUPwOzUoSLK/xJkzZ8bExFy6dGn9+vWnTp2Sqozm2MtPFXA+c3vwK5KUU3Yav7zQypsXAzis49dY7x+NAe60c4wyBNeI+YOfn9+4ceMmTZqUnJwcEREhdTmNIQjBikaFcbvGKt8/IT53UDAiDcHZfXVRHLXV+H4Cv2Sw4s7ucOvEZs+effLkyaeeesq8ZNGiRTExMd99992yZctiYmKWLl0qVW121G0Hp9TVlzt4r/KhbOPYbcafhitx+xhwVkvPie/8Lu5JV8bJ+zT55tTU1ERGRja8p/rs2bMff/xx81MPD8l60PjRAlbn60q/jlCGeHITdhhrBKmrAbCCjHPiO7+LWWMVSMGb1dfXf/3113/961/nzZvX8DAZDw8Pvwbc3d2lqhBBCLag4Gj5UEWQB3fvdmQhOJsVF8QFv4uZYxTRPkjBJhiNxtOnT8+dO/fpp5+WupamYaAKbETB0cokxSO7hIk7jD8PV7rZ1+HTAHdo5QXxb0fEnWMUndshBZvm4eHx3nvvSV3FraBHCLaj4OjrYYr2rtyDWUI9jp0Bx7fygvj6ESFzjKIrbknmyBCEYFOmLBSJpmbjOFJwbN9dFt84ImaOUXZDCjo4BCHYmgtPP6Qoqo1sxh4B1yQFB7XmsvjSITFrLG5P7QwQhCABNwX9mKY8d4P96z/oFYLjOVbG5h4Qto3CMaJOAgfLgDTcFbQ+TZHws/EuP25kGPYm4DDKa2lKlvDpIEV3P/vdbt3d3ZVKZb9+/aQuxJKKiopeffVVa7wzghAkE+7FrUtT3p9p3DdeiePOwSGIjB7eZZwQyU2OsuvhtI4dO+7atausrEzqQiwsLi7OGm+LIAQpDQri5vVS3LdD2D8eF50BBzD/uFBRR+/2d4CzfyIjIyMjI6WuwjHY9Y8akIPnevC9O3JP7sVp9mDvMotZxjm2NhXXEXU2WJ8gvc8HKc5VsE9O48AZsF8FOvbILuN3KQqVJ4bxnQ2CEKRnOnDm3RPCbg3OpwB7VCPQpEzh1V6KIbi5pjNCEIJdiPTmVgxVPrJLuFEndSkAN3l2vxDTjpvbAztM54T1CvZiZBg3NoJ75RAmC8G+fHVR3F/Clg1xgANk4M4gCMGOLExQbC1iWWoMkIK9KDHQK4eENSkKbxepSwGrQRCCHWnnQl8OVTzxm6Crl7oUACIiejFHeCyW7+mPqUFnhiAE+5Kq4oYGc387igFSkN4uDdurZX/vjUFRJ4cgBLuzOFHx4xW2V4sBUpBSrUBP7xU+uYfHoKjTQxCC3WnvSp8O4p/4TcC97EFCC34Xevpz4yKwk3R+WMdgj8ZF8D39ufnHkIQgjfM32JKz4kcDsYeUBaxmsFOfDVKsvCgeKcMAKdgaI5q9V/hHH0WoF46RkQUEIdipju60MEHxzD7cvBdsbcUFUW+k2d2we5QLrGmwXw935hUcrc3FNUjBdq7V0huHhS8GK3j0BmUDQQj2iyN6L0Hx6mGxFnOFYCsvHBSmdeZ7d0AMygiCEOza0GAu3peWnkOnEGxhl4bt1rA3++DEQXlBEIK9ey9B8c7vQiWuNQPW99phYWECThyUHQQh2LseftyYcP79ExgeBevaVMh09TQ5CntF2cEqBwfwz378F2fFQj0OIAVrYURvHRP+2Y/HMTIyhCAEB6Dy5J7sys8/hplCsJaf8kSB0b2R2CXKEdY6OIbX7lZsKhRPl6NTCJYnMpp/TPxnXwV6g/KEIATH4ONCr/RUvHYYnUKwvB9yRVcFjQ5HDsoUghAcxpx4/kwFbtsLFiYwmn9cfLc/uoPyhSAEh+HK0z/78q8fwUXXwJK+vSQGelCqCjkoXwhCcCQPxvAGI20vQhSCZdSLNP+4OL8vzqCXNQQhOBKO6OWe/If/wTmFYBkrL4jRPjQ0GN1BWUMQgoN5IJq/eIOO4fZM0GZ1Ir17At1BQBCCo3HhaW4P/sP/4PBRaKuMc2J3PxoYiO6g3CEIwfE8EcdnFot5VegUwp2rFejdE+JbuL42IAjBEfm40BNx/Een0CmEO7fmsniXH/XpiO4gIAjBMc3tofjmklhWI3Ud4LA+PSPO7YHuIBAhCMFBBXnQfZ34L86iUwh3Yo+W3ainEaHoDgIRghAc10s9+U/PCAaj1HWAA/r4tDi3O240Af+FIARHFdeeGxDIr7qETiG0ToGOZanFR2Kx94P/wqYADuzlnvz7J0Rccg1a5dMz4mOxfDvchh7+gCAEBzYoiAv2pF/y0SmEljIYacUF8Zl47Prgf5QWeZf6+vpVq1ZdvHixb9++999/P8c1MfSenZ2dmZkZEBDw+OOPt2/f3rSwqqpqxYoVJSUlycnJaWlppoXnz5/fvn27RqPp3Lnz1KlTPTw8iEin061evdr8bgkJCXfffbdFigeH9tJd/LsnxPs6Yb8GLfL1JTExiI9ph+lB+B/L7D4eeuihlStXBgUFzZ8//9VXX725wTfffDN16lQ/P7/9+/cPHjy4vr6eiARBSEpKys7O9vf3nz59+vLly4nIaDQOGjTo9OnTvr6+q1atGjhwYHV1NRFdu3btL3/5y9E/aLVai1QOju7eSP5GHe3VYngUbo8RfXJafBbdQWiEtdmZM2c8PT0rKioYYxcuXPD09Lx+/XrDBqIodu3ade3atYwxQRC6d+/+/fffM8Z++eWXzp07G41GxtjGjRujo6MFQRBFUa/Xm15YW1sbFBS0efNmxlheXp6Pj09L6lm6dOnMmTPb/r2gLSorK232WZ+cFqZmGW32cY7ClqvAUewsFruvqxdt9XEGg6G+vt5WnwZNEARBp9PdtpkFfhnt3r17wIABptHO2NjYkJCQnJychg00Gs25c+dGjhxJRDzPDx8+fNeuXUSUnZ2dlpamUCiIKC0tLS8vLz8/n+M4T09P0wt5nq+vrzc/NRqNn3766bJlyy5dutT2ssFpPNKZ31KIk+vh9hafFv/aHSdNQGMWmCPUaDSBgYHmp0FBQRqNpmEDrVbr4eHh4+NjbnDw4EHTC+Pi4kwL3dzc/Pz81Gp1VFSU+YVvvfVWly5dhgwZQkQ8zw8ZMuTKlSulpaXPPffckiVLpk2b1lxJOTk5TzzxhOlxSEjIG2+80favCa1SU1Pj4mKjw/LciEaH8ivP1v6lGwZI/8eWq8Ah5Oton5b/MlGssdVvppqaGqVSqVRa5lAMuAOiKIri7Q+ma+kacnd3v3nhxx9//OSTTyqVyoafZDQaG614pVIpCELDBqa/z5tf2PDvdvny5atWrdqzZw/P80QUHh6+bds20z+NGjVq7ty5twhCf3//vn37mh+bOp1gSwqFwpb/25+Io2cOcHN72OwDHYCNV4H9++wCzehC7dxs9/9E8QebfSI0wnGcJYNQp9PdvNC0glUq1datW80L1Wq1SqVq2CwkJKSurq6srKxjx46mBiEhIaYXFhcXm9+/srLS/MKvv/76H//4R1ZWVmRk5M2fO3To0LKysoqKCl9f3yarjY2NnT17dgu/GliDi4uLLbsjKWHEccYj5YpE3FLnDzZeBXZOV0/fXK4/MkHp4mK7LUQQBPQIpSWKotF4+6tPtXSOUNkU02kSI0eOPHbsWGFhIRHl5ORUV1ffc889RFRQUHD27FkiCggIGDBgwPr164nIYDBs3rw5PT2diMaOHbtt2za9Xk9EP/74Y69evcLCwoho3bp18+bN27ZtW5cuXcwF1NXVmR9v2bJFpVI1l4IgTzPi+IxzOKEQmrbyopii4qN88DsJmmCBnyphYWF/+ctfkpOTR4wY8csvv8yfP9905t+SJUtOnjz566+/EtGCBQseeOCBI0eOnDhxolu3bsOHDyeipKSkfv36DR06tG/fvj/++OOqVauIqLKycurUqVFRUc8995zp/f/617+OGzfugw8+2LRpU9euXbVa7b59+77++uu2Vw7OZHosH7e2/qM6RXtXqUsB+5NxTlyciCFKaBrHmGWOL9i/f//58+f79u3bs2dP05L8/Hy9Xh8fH29+umvXruDgYPORokQkiuLOnTuLi4uHDh0aHR1NRPX19bt37274zl27dg0LC9PpdDk5OQUFBf7+/oMGDTKNsjYpIyMjJydn2bJlFvlecGeqqqrMh0fZzJSdQoqKe7obzhIjkmgV2Kfj19h9mULuA0ob9wdxsIzkRFE0GAxeXl63bmaxILQfCEJ7IMleeHsxeyVH+P0+7HeIEIQNPHdQ8HWlN21+M3oEoeRaGIT47QzOY3gopzPSsTJn+20HbWEU6fvL4iOdsa+DZmHjAOfBET0Wyy87j0Nm4H82FYqx7TlcXBRuAUEITmVmHP9DrqjH3XrhD19dZNNx60G4JWwf4FRCPOmeIH5tLjqFQER0vZay1OKkKOzo4FawfYCzmdWVy8DoKBAR0erLYnoE74szauCWEITgbMaE80V6Ol2OQ2aAvrogYlwUbgubCDgbBUePxnI4ZAbOVDCNgVJUOEwGbgNBCE7o8S78msuiEVEob19dEB/pzCmQg3A7CEJwQtE+XCcfLkuD0VH5EhmtucwexbgotAC2EnBOD0bz311Gl1C+thezEE/q5ov+INweghCc0wPR/M/5Yo1w+5bglL66iMNkoKWwoYBzCvGkXv7ctiJ0CuWosp62FopTorF/gxbBhgJO68EY/rtcTBPK0feXxbRQvqO71HWAg0AQgtOaHMVvLRR19VLXATaHcVFoFWwr4LT83SgxiNtYgNFReblUyS5VspFhOEwGWgpBCM7swWiMjsrO97lsSjTvgn0btBg2FnBmEzvxuzTi9Vqp6wAbWpsr3o+rbENrYHMBZ+bjQmmh/M/5GB2Vi9wqVmJgg4IwLgqtgCAEJ/dgNIcz6+Xjh1x2XxSPy6pBqyAIwcmlR/BHy1iJQeo6wCbWXxEndcJuDVoHWww4OXcFjQnn111Bp9D55VWxPB0bGoz+ILQOghCc34Mx/He4Z70MrM9jEzvxSuzVoJWwyYDzGxHKnatg+TqcR+HkMC4KdwYbDTg/F54mduLXXkEQOrNiPTt/gyXjNrzQeghCkAXclcnprc9j4yN5V+zSoPWw1YAsJIVw6mp2qRKdQqe1/op4Xyd0B+FOIAhBFhQcjY/gN+QjCJ1TiYFOXmfDQ7FDgzuB7QbkYnwkvwEX4HZS66+IYyN4d4XUdYBjQhCCXKSquBPXWFmN1HWAFazPEydhXBTuFIIQ5MJNQamh/OZCdAqdTVkNHS1jI8OwN4M7hE0HZGR8BLehANOEzuanPHFUGO+plLoOcFgIQpCR9Ah+Z7FYI0hdB1jU+jxxUhTGReHOIQhBRvzdqKc/l6VGp9B5VNTRwVI2GuOi0AbYekBe7o3kN+D2hE7k5zwxVcV7u0hdBzgyBCHIy72R3K8FDF1Cp4FxUWg7BCHIS0w7zteVDl9FFDqDaiPt0bD0COzHoE2wAYHs3BvJYXTUOWRrWN+OXDuMi0LbIAhBdsZH8jiJwjlsLRRHhWMnBm2FbQhkJyGAK6thuVXIQoe3tYiNCsMEIbQVghBkh+dobDguwO3wrlQxvZHd5Y8ghLZCEIIc4SQKJ7CpkI0O4xGD0HYIQpCjVBV3tIxdq5W6DmiDbUXiqHDkIFgAghDkyENJKSp+Cy7A7bBqBfpNy1JV2IOBBWAzApkaH8lhmtBx7dGyHn6cv5vUdYBTQBCCTI2P5HcUi7W4ALdj2lYk4r5LYCnYkkCmOrhRD39ulwadQoe0tYiNxgQhWAiCEORrfAS/oQDThI6nSM9KDaxPBwQhWAaCEORrYiful3xcgNvxbC5kI8N4nDkBloIgBPnq3I7zcaFjZYhCB7OtiI3EBWXAchCEIGvjI7hfMTrqUOpFytaII0Kx7wKLwcYEsjY6nN9WhB6hIzlQymJ8uEAPqesAJ4IgBFm7J4g7W4FLzDiSrYUijhcFy0IQgqy58jQ0mN9ZjNFRh7G1iOEMQrAsbE8gdyPDOIyOOgqtgfJ1bEAAeoRgSQhCkLsx4dy2YpxE4Ri2FYmpKl6J/RZYFDYokLsoH85dQafLEYUOYFsRwx0nwOIQhAAYHXUMIqPMYnF4KIIQLAxBCEAjQ7ltRThext4duspCPLlwLwQhWBiCEIBSQ/mcUlZtlLoOuKWtReIoXFAGrABBCEBeSurTkdujxeioXdtWxEaFY5cFlqe0yLvk5+d//vnnlZWV6enpY8aMublBbW3t559/fvbs2fj4+NmzZ7u6upqWnzx5csWKFYIgTJs2bcCAAaaF//73v/V6velxly5dJk6caHp84cKFjIwMg8Fw//33Dxs2zCKVA5iMDOO3FYmjwhRSFwJNK6+lsxXsnkD0CMHyLPDzqqKiIjExsa6urk+fPjNmzPj+++9vbjN9+vQNGzYMGTLk559/fvzxx00Lz507N3jw4MDAwJiYmBEjRhw+fNi0fMGCBZcuXSovLy8vLzcnolqtTkxMdHd379Gjx3333bd9+/a2Vw5gNjKM24rjZezYtiJxWAjvhh8qYAUW6BF+9dVX3bp1W7RoERG5u7svXLjwgQceaNggNzf3559/VqvV/v7+Y8aMUalUb7/9dqdOnT755JNp06a99tprRFRWVrZo0aI1a9aYXvLCCy9069at4ZssWbIkJSXln//8JxHV1dW9//77I0aMaHvxACZ3d+Bu1LErVSzKB30Oe7QVd5wAq7FAj3Dfvn0pKSmmxykpKcePH6+urm7Y4MCBAz179vT39ycif3//u+6668CBAze/cN++feaXLF269K233tq8eXNzn9KwMUDbcUSpKj6zGJ1Ce8SItheLCEKwkhb1CEVRNBgMNy/39PTkOE6r1aamppqWBAQEEJFGo4mJiTE302q1HTt2ND8NDAzUaDSmZublAQEBWq2WMcZxXFpamp+fX319/dNPP52Wlvbll182epOAgICampry8nI/P78mC/7tt98mTZpkehwWFvbuu++25GuCBRkMBoXCwYaxkgP5Xwv4aRFOcvCoI66C5pwo57wVimCF4c+/se1dTU2NUqlUKi1zKAbcAVEUWQsuG9WiNbR///709PSblx88eLBr165ubm51dXWmJaYH7u7uDZu5ubnV19ebn9bW1poaNFxeV1fn5ubGcRwRffPNN6aFs2bNiomJmTdvXlxc3M2f4ubm1lzBkVlMORgAACAASURBVJGR5uHZ9u3bN6oHbKC+vt7h/rend6IXjwgKV3cXpzgy0RFXQXOyr7IxEcwRvw6CUFqiKNbW3v7mMi1aQ4MHD66oqGjuX0NDQ4uKikyPCwsLlUplUFBQcw2IqKioKDQ0lIjCwsIavtC0sKGIiIjAwMCCgoK4uLhGn+Lv7+/p6dlcSREREVOmTGnJVwMr4Xme5x0sTwI9KbqdeOQaNyjIGYbgHHEVNGdbkfG1uxU872Drhf+D1IXAbVhgDU2YMOGnn34yjZ2uWbMmPT3d9Ato7969ly9fJqLU1NTi4uJjx44R0dGjRzUajWkodcKECWvWrDH1W9esWTNhwgQi0uv1RuN/x6b27NlTVlbWvXt3U+O1a9ea/sncGMCycIkZO1RZT79fY0ODHSwFwYFYoM8+bty4pUuXDhgwIDo6+uDBgzt27DAtnzdv3sSJE1966aV27dq98847Y8aMSUpK2r1797vvvuvt7U1ETz311LfffjtkyBB3d/eCgoLFixcT0f79+x999NG+ffvW1NQcOnRo0aJFKpWKiB566KEVK1YMHDgwODj4999/37NnT9srB2hkZBj/yiFhfl+p64AGMovFe4I4T4wvgtVwLZlIvC1RFA8dOnTt2rXBgwe3b9/etDAvL8/Hx6dDhw6mp5cuXTKdUN/wOJq6urp9+/YZjUZTHJoWXrx48eLFi25ubj179jQdfWMiCML+/fv1ev3gwYNNUdqkjIyMnJycZcuWtf17wR2rqqry8fGRuopWqxcp6Nv685NdAhxvNqoxB10FN3tqrxDvy83t4XgDjDhYRnKmIz29vLxu3cwya4jn+YEDBzZa2KlTp4ZPO3fu3Llz50ZtXF1dk5OTGy2MjY2NjY29+VMUCsWQIUPaWitA81z+uGH9gzGOt9t1VtuK2At3YXWAFWHzAviTkWHcNpxNaDdOlzOOo7j2mCAEK0IQAvzJqDBua6GIJLQTW4vYaJxHD1aGIAT4kygfztuFO3UdUWgXthXhgjJgdQhCgMZGYXTUPlQbKaeUpaiwmwLrwhYG0JjplkxSVwGUrWH9AjgfF6nrAGeHIARobFgIl1PKDE5yzVEHtrVQHBmGfRRYHTYygMa8XainP7evBKOjEttaxEZhghCsD0EI0IRUFZelxuiolK5UMb2R3eWPIASrQxACNCFFxe9Uo0copU2FbEy4o11mGxwTghCgCYlB3LkKVlEndR0ytrUQJ06AjSAIAZrgytPAQG6PBqOj0qgVaG8JS8WJE2AT2M4Ampai4rM0GB2Vxh4t6+HH+Td7720AS0IQAjQtRcXtxGn1EtlWJI7CiRNgK9jUAJrWpyOnrmZag9R1yNLWIoYJQrAZBCFA0xQcDQ3ms3EShc2VGkhdzfp0RBCCjSAIAZqVouKycBKFzWWpxaRgXoEcBFtBEAI0KzWU24FpQpvL1rBkFWIQbAdBCNCseF/OyOhKFbLQprI1LDkEQQi2gyAEuJVhIRgdtSl1Nbtew7r7IQjBdhCEALeCaUIb26lmySoel1YDW0IQAtxKmorbqRaRhDaTrca4KNgaghDgViK8OS8ld6YcUWgjOFIGbA9BCHAbqaEc7kRhGwU6ZjCyrr4IQrApBCHAbaTgeBlb2almKSrMD4KtIQgBbiNFxe/WiAKi0PowQQiSQBAC3EagB4V7c0fLkIRWt1uLCUKQAIIQ4PZScRKF9V2qZEaROrdDEIKtIQgBbi9FxWXh6ttWlq1mqegOghQQhAC3NyyEP1jKagSp63BqOHECpIIgBLg9HxeK9+MOlmJ01Ip24xKjIBEEIUCLpGJ01JrOVTAlT518EIQgAQQhQIukqHicVm892RpMEIJkEIQALTIoiPvPdVZZL3UdTipbjQlCkAyCEKBF3BXUP4D7TYtOoeUxot1acRgmCEEiCEKAlhoWwu/WYJrQ8k5dZ+1cuHAvBCFIA0EI0FLJIVw2pgmtACdOgLQQhAAtlRDInb/BKuqkrsPp4BKjIC0EIUBLufI0IJDbi2lCixIZ/aYVkxCEIB0EIUArDAvhd2Ga0KJOXGeBHpzKE0EIkkEQArRCcgiXrUGP0JJw4gRIDkEI0AoJAdzlSlZeK3UdTiRbI2KCEKSFIARoBSVPAwK437QYHbUMgdG+EjY0GDsikBK2P4DWGabid2F01EKOlbEwLy7QQ+o6QN4QhACtg2lCC8rGHSfADiAIAVqnX0fuShW7hmlCS8hWizhSBiSHIARoHSVPiYGYJrQAo0gHStkQTBCC1LAJArTasBBME1rA4TIW7cN1cJO6DpA9BCFAqyWrcNFRC8jCGYRgHxCEAK3WpwOXr2NlNVLX4eCy1WJyCHZBID1shQCtpuRpUBC3B9OEbVAn0qGrbHAweoQgPQQhwJ3ANGEbHSxlXX05X1ep6wBAEALcmWQVhyBsi2w1S8EEIdgHBCHAnejdgSvSs1KD1HU4rGwNJgjBXmBDBLgTCg7ThHeuRqBjZeyeIPQIwS4gCAHuEKYJ79j+EnaXP+fjInUdAESEIAS4Y8khmCa8Q7j1EtgVBCHAHbq7A6epZiWYJmy9bDVLVmHnA/YC2yLAHeI5GhzM79ZgmrB1qo104jpLDESPEOwFghDgzg3D6Gjr7S1hfTpwnkqp6wD4A4IQ4M5hmvAO4NZLYG8s86vs+vXrCxcuvHz58t133/3iiy96eDS+4TRjLCMjY8eOHQEBAS+99FJ0dLRpeV5e3gcffFBSUpKamvrUU0/xPE9Ef/vb34xGo/m1AwcOnDBhQkVFxXvvvWdeOGrUqGHDhlmkeIA71tOfKzUwdTVTeWLP3lLZGrawv0LqKgD+xzI9wvHjxxcWFs6YMWPPnj1PPfXUzQ0+/PDDjz766JFHHvH19R08eLBOpyOimpqapKQkDw+P6dOnf/HFF++8846pcVRUVPQfPvvss4qKCiK6cePG4sWLzct9fX0tUjlAW/AcDQnm96BT2GJV9XS6nA3ABCHYFdZmBw8e9PX1ra2tZYxptVpXV9fi4uKGDYxGY2ho6Pbt201PExMTly5dyhhbtWpV7969TQt/++23wMDAurq6hi/cu3evt7d3VVUVYywvL8/Hx6cl9SxdunTmzJlt/lrQJpWVlVKXYCOLTwlP/maUuoom2Ocq2FQgpmyql7oKGzEYDPX1cvmy9kkQBJ1Od9tmFugRHj58OCEhwdXVlYiCgoKio6OPHTvWsIFarS4uLh4yZIjp6ZAhQw4fPkxEhw4dGjx4sGlhYmLi9evX8/LyGr7wyy+/nDp1qre3t+lpXV3d888///LLL2/fvr3tZQNYBI6XaZUs3HoJ7E9L5wjz8/NvXtihQwdvb++SkpIOHTqYF3bs2FGr1TZsptVqvby83N3dza86c+YMEZWUlHTv3t20UKFQ+Pn5abXa2NhY0xKdTrd27Vpz5rm5uT3xxBPdunUrKSmZNm3aCy+88NprrzVXbVZWVkpKiulxRETEv//97xZ+TbAUvV7PcbIY/opypes1rhdKdSpPqUv5M/tcBTuLXT/oY9TpZHHOSU1NjVKpVCpxgKxkRFFsyV9BS9eQOVcamj9//rRp0zw9PWtra80LDQaDl5dXw2ZeXl6mgVNTQQaDwdTJu/ULv//++9DQ0IEDB5qeBgcHm/MsISFh0qRJ8+bNMx1cc7P4+PjnnnvO9Njb29vcpwSbYYzJ5397kko4XOk5LdC+Ojp2uAoq6uhyVf3QCE9X+/pfZS3KP0hdiHyJomgw3P6aFy1dQ5cvX27un8LCwtatW2f+1IKCgvDw8IYNVCqVKIrFxcVhYWFElJ+fb3oQFhaWm5tralNWVqbX603LTZYvXz5r1qwmPzE+Pl6v1+t0unbt2jXZIDg4OC0trYVfDaCNkkO4bA2b1lnqOuzebo14TxAnkxQEB2KBTTI9Pf3cuXMnTpwgok2bNrm5uZm6cceOHcvMzCQiX1/ftLS0FStWEFFZWdnGjRsnT55MRJMnT966datpHHXlypVDhw4NDAw0vef58+ePHj368MMPmz9Fq9UKgkBEjLGlS5d27dq1uRQEsLFkFZetxjTh7WVrGCYIwQ5ZoM/u5+f3/vvvp6am9urV68SJE8uXLzcNBaxfv/7kyZOmntn7778/ZsyYzMzMS5cuTZkyJSEhgYh69er1+OOP9+7dOy4u7sKFCxs2bDC/57Jly9LT0825SERff/31hx9+GBcXp1ariWj16tVtrxzAIrr5cnojK9CxCG+7m5OzK7vUbMlgBCHYHY4xy/yS1Wq1V65c6dq1q5+fn2mJXq83Go3t27c3PTUYDL///ntwcHBUVFTDF+bn56vV6rvvvrvhafhVVVUuLi7m42tMCgoKioqKOnToEBMTc4th94yMjJycnGXLllnke8Gdqaqq8vHxkboK23kwSxgTzj0aa0d7eXtbBddqKeb7+rKHXZR29D/JunCwjORMc4SNDlu5mcXWUHBwcHBwcMMljT7bw8MjMTHx5hdGRkZGRkY2WtjkH3BERERERESbKwWwvGEhXLaGPRordR12LFstDg7i5JOC4ECwVQJYQLKKy8I04S1la3DrJbBT2C4BLCCuPScwulKFLGxWtprhZrxgnxCEAJaRFMxl4xIzzSg1kNbAenVAEII9QhACWAZOoriFLLWYFMwrkINglxCEAJaRHIJpwmZlaxjuQQh2C0EIYBkx7TgXni7eQBY2IQsThGDHEIQAFpMUgmnCJqir2Y061sMfQQh2CkEIYDHJuCVTU3aq2bAQHjEIdgtBCGAxKSouSy0iCRvJVmOCEOwaghDAYiK8OQ8ld74CUfgn2RpMEIJdQxACWFIypgn/7EoVqxVYnC+CEOwXghDAknA2YSOmWy8hBsGeIQgBLClVxe3SYJrwfzBBCPYPQQhgSSpPrr0rd7ocUfhfuzBBCHYPQQhgYRgdNbtwg3FEMe0QhGDXEIQAFoazCc2yNSwF46Jg9xCEABY2LITfrRExT0iYIAQHgSAEsLAQTwrw4E5cl3sSMqJdGnEYJgjB7iEIASwvRcXtlP004Zly5qXkIr0RhGDvEIQAlpeq4nYWi1JXITFMEIKjQBACWF6Kit9fwmoFqeuQFCYIwVEgCAEsz9eV4ny5g6XyHR1lRHu0YlIwghAcAIIQwCrSVNxOtXxHR09cYx3cuFAvBCE4AAQhgFWkhvJyPl5mRzEbHooUBMeAIASwikFB3H+us8p6qeuQSGaxmIogBAeBIASwCncFDQjkdmvkODpaJ9LBUpYUjN0LOAZsqQDWkqqS6ejoXi3r7sf5uUldB0DLIAgBrCVVxWUWyzEIM4vFNIyLguNAEAJYS5+OnLaaaaqlrsPmMtUsLRT7FnAY2FgBrEXBUVIIL7eTKMpr6XwFGxCAHiE4DAQhgBWlyu+iozvV4uBgzk0hdR0ALYYgBLCitFBuh8ymCTOLMS4KDgbbK4AVdWnPKTk6f0NGWZipZmm4xCg4FAQhgHWlqLidsukU5uuYrp718EcQgiNBEAJYl6ymCbcXseGhPGIQHAuCEMC60kL5bI0oyCMKM9UsFeOi4GgQhADWFeRBoZ7csTLnT0KRUbYalxgFx4MgBLC61FAuUwajo79fYx3duXDcegkcDYIQwOpSVdzOYuc/rX5HMcOV1cARIQgBrG5YCH/oKjMYpa7DynaqRZw4AY4IQQhgdT4udJc/t6/EmUdHawQ6WMqSQrBLAceDrRbAFtJUnHNfdHSvlt3lz7V3lboOgNZDEALYQloo79zXWsO4KDguBCGALSQGcrlVrMQgdR1WswOXGAWHhQ0XwBaUPKWp+K1Fzjk6eq2WLlWygYHoEYJDQhAC2MjocG5LoXOOjmapxcFBnAt2J+CYsOUC2MioMH5HsWh0xj7htiI2Igw7E3BU2HYBbCTEkyK9uYOlztYpZERbCtmYcIyLgqNCEALYzphwbovTTRMeK2PtXKlzOwQhOCoEIYDtjA7nnW+acFMhG4vuIDgyBCGA7QwM5Ap0TF3tVFm4uVAcE449CTgwbL4AtqPgaHgYv7XIeYLwag2dr2CDg9EjBAeGIASwqdFhTnUSxZZCMTWUd8WOBBwZtl8AmxodzmcWi/XOcsTMZhwvCo4PQQhgUwHuFNue2+8Ud6IQGGUWi6PCEITg2BCEALY2OsxJTqLYV8I6+XAqTwQhODYEIYCtjQ7nNzvFNOHmQhHjouAEEIQAtpYQwJUYWIHO4bNwUwEbixMnwPFhIwawNZ6j4aH8Ngc/iaJQz0oMrH8AeoTg8BCEABIYHcZtcfAg3FjARoXxPHIQHJ9lgvDgwYOTJk0aNmzY//3f/zHWxJ/31atXZ8+enZSUNGfOnGvXrpkWFhcXZ2RkPPvsswsWLGjYuLKy8vnnn09KSpo5c2ZxcbF5eVZW1vjx41NTUzMyMixSNoBURofzWWqxVpC6jjbYXCiOjUAMgjOwQBCq1eqRI0eOGDFiwYIFS5cu/eSTT25uM2XKFIPB8MEHH+j1+ilTppgW7tmzZ/v27RqN5scff2zYeObMmXl5eR988IGXl9eYMWNMyXr+/PkJEyZMmTLl73//+9tvv/3NN9+0vXIAqfi7UTdfbp/DnkRRI9AeDRuOW9KDc2BtNn/+/AkTJpge//LLLzExMY0anDhxwsvLq7q6mjFWXV3t6el56tQp87+uXr26d+/e5qeFhYWurq4lJSWMMUEQgoKCsrOzGWPPP//8jBkzTG1WrFjRv3//5upZunTpzJkz2/69oC0qKyulLsHevXVMeOmg0Xrvb9VVsLlAHPprvfXe3zkYDIb6evxfkpIgCDqd7rbNLPCD7vfffx84cKDp8cCBAy9fvlxZWdmwwfHjx3v16uXh4UFEHh4evXr1On78eHPvdvLkyU6dOgUGBhIRz/MJCQmmxsePHzd/SmJi4okTJ1hTY7AAjmJMOOe4J1HgQtvgTJQtaaTT6S5cuHDz8vj4eHd395KSEj8/P9MSf39/IiopKWnXrp25WWlpqbmBqY1Wq23us5pr3PBT/Pz86urqrl+/3qFDhybfZOvWrX369DE9joqK+uqrr1ryNcGC9Ho9x2EC6Va6uFFZjeuZEn2El1Xi0KqrYFOB63dD6nWOfwaIVdXU1CiVSqWyRbtZsAZRFFvyV9CiNXT+/PnZs2ffvPz777+Piory9vY2GAymJdXV1UTUMAWJqGEDItLr9Y0atKRxo0/hOM7b27u5N0lISHjjjTdMjz08PG7REqyEMYb/7beVHiFklin/GmSVrpX1VsG5CiaQkBDqZY03dybKP0hdiHyJotgwUJrTojXUt2/fQ4cONfevkZGRly9fNj2+fPmyp6dnQEBAcw2IKDc3NzIy8hbvlp+fX19f7+LiYmp877333vwpwcHBbm5uzb2Jv79/3759W/LVACR0XxS/8ITw1+4ONsa4CRfaBudigb/AqVOnrlu37urVq0T0xRdfTJkyhed5Ilq9evXhw4eJKDU11WAwbNmyhYg2b95cW1ubnJzc3Lv169cvICBgzZo1RJSTk3PhwoVx48aZPmX16tVVVVWMsSVLlkydOrXtlQNIK03FnSpnmmqp62iljQVieoSDhTfALVigz56UlDR58uT4+PiOHTsqlUpT4BHRp59+OnHixP79+7u5uS1duvTRRx8NDQ1Vq9XLly93dXUloqysrNTUVFNjjuNGjRq1ZcsWjuMyMjKmTp26aNGiwsLCzz77rH379kR07733btiwITY2tl27dn5+fp999lnbKweQlpuC0sP5H/PEZ+IdJldKDPT7NZaqQo8QnAdnqWMvy8rKKioqYmJizDOTpllK89Pq6uqioqLw8HDT4aO3Vltbm5+fr1KpGk1ylJSU6PX66OjoW7w2IyMjJydn2bJld/pVwAKqqqp8fHykrsIBbMgXPzolZo+1/DSSlVbBp2fEnFK2apjC4u/sfHCwjORMc4ReXreZz7bYGurYsWPHjh0bLjENkJp5enp26dKlhe/m5ubWZOOgoKA7rhDADo0M4x/bI2iqKcRT6lJaZu0V8cW7HKb/CtAS2KABpOSmoNFh/M/5jnF7Qq2BTl7HBWXA2WCDBpDY/VHcuiuOEYTrrojjI3h3DIuCc0EQAkhsdDh//Borvf3JTtJbmytOjsZOA5wNtmkAibkraKQjjI5qDXSqnKXheFFwOghCAOnd38kBRkfX5orjI3k3jIuC00EQAkhvdDh/6Cq7WiN1Hbe09oo4OQp7DHBC2KwBpOeppBGh/C92PDqqqaYz5SwtFOOi4IQQhAB2wc6PHf0hVxwfybtihwHOCNs1gF1Ij+BzStn1WqnraAbGRcGJYcsGsAueSkq119HRYj07W8FSMS4KTgpBCGAv7PbY0bVX2MROGBcFp4VNG8BepEfwe7Ws3P5GRzEuCs4NGzeAvfB2oRQV/2uBfXUKi/TsfAVLxnn04LwQhAB25P4obt0Vy9wZzVLWXWETMC4KTg1bN4AdGR/J7y0RS+zpuqM/5GJcFJwctm8AO+LjQpM68Ssu2Mvo6OVKdrmKpWBcFJwaghDAvjzZlV96ThTtY3z087PijC68C/YT4NSwgQPYl/4BnJ8b7VRLn4QGI626KD7ZFXsJcHLYxAHszqw4fsk56UdHv88VBwRyUT4YFwUnhyAEsDsPd+az1aK6WuJO4ednxdndcNclcH4IQgC74+1Ck6L4lRekDMLj11iJgUaGoTsIzg9BCGCP5nTjl5wTBemi8N+nxTnxvAI5CDKAIASwR3d34ALcKbNYmiSsqKOf8sXHYrF/AFnAhg5gp57sKtkhM1+eF9PD+UAPST4cwNYQhAB26qEYfrdGLNbbulPIiJacE2fHY+cAcoFtHcBOebvQ5Ch+hc0PmcksZh4KSgzE9CDIBYIQwH7Njuczztv6kJnPz4rPoDsIcoLNHcB+9fLngjxoe5HtkrBIz/ZoxIc6Y88AMoLNHcCuze7Gf3RKsNnHZZwTH+rMeylt9oEA0kMQAti1RzrzBTrKssmlR+tFWn6BPYWLi4LMYIsHsGtKnv7Rh3/jiC06hWsui3HtqbsfDpMBeUEQAti7B2P4WoE25Fv3nEKDkf5+VHyrDy4uCrKDIASwdxzRP/sp3jhi3ZsULjwpDA7iBgejOwiygyAEcABjw7kO7vTtZWt1Cov17NMz4jv9sUMAOcJ2D+AYFvRV/L+jYq115gpfPSzO6cZHeqM7CHKEIARwDIODuW6+tOy85TuFx8rYLg17pRdmB0GmEIQADuO9/ooFx4Wqeku+JyOae1B4ux/OHQT5QhACOIye/lxSCP/vM5bsFH53WdTX08O4lAzIGLZ+AEfyz778ov8I12st824GI71+RFycqOAxOQgyhiAEcCSx7bkHovln9lvmQtyLTon9O3JDcMoEyBuCEMDBfDhAcaWKLTzR1gHSEgN9fFp4PwE7AZA7/A0AOBh3Ba1PVfz7jLi58M67hYzo6b3CzC58Jx90B0HuEIQAjifUi1uXqpi+23im4g6z8B9HhdIa9g9cUA0AQQjgoAYGcv8aoLhvh1BR1+rXrr8ifn2J/ZimdEMOAiAIARzXo7H8iDDuwSxjq46cOX6NPb1PWJeqCPKwWmUADgVBCODA/jVAUSfQ/zva0guvZWvY6K3GjCGKvh0xNQjwXwhCAAfmwtMPqcr1V9jEHcKp8tt0DL84K07NMq5OVk6IxB8+wP/g7wHAsXV0pxP3KYeFcMM3Gx/ZJeRWNRGHhXo26zfhk9Pi3nHKFBX6ggB/giAEcHhuCprbg784xSWuPTfwF+OcfcKVKlag5y5Xsh3F7IEsofePxnaudOBeZed2SEGAxnCdXQAn4e1Cf+vNz4nnF54QBm4wunGurkohyIOmxvDLhrj4uEhdH4C9QhACOBV/N1qYoFiYoKiqqvLx8ZG6HAAHgKFRAACQNQQhAADIGoIQAABkzQmDsKamxmAwSF2FrBmNxqKiIqmrkLuioqL6eovezB5a6dq1a5WVlVJXIWt6vb60tPS2zZwwCA8fPnz06FGpq5C1oqKiUaNGSV2F3I0fP/7KlStSVyFrCxYs+Pbbb6WuQtYyMzPnzp1722ZOGIQAAAAthyAEAABZQxACAICscYzd+U2u7dNnn302f/58Ly8vqQuRL0EQbty44e/vL3Uhsnb9+vX27dsrFLjloGSqqqqUSqWHB+53JZm6urpHH3307bffvnUzJwxCIiooKDAajVJXAQAAEgsLC3N1db11G+cMQgAAgBbCHCEAAMgaghAAAGQNQQgAALKGIAQAAFlztvsRnj59et++feankydP9vPzIyLG2HfffXf48OGoqKhZs2a5u7tLV6OTq62t3bFjR05OjiAISUlJI0eONC0/duzYkSNHzM0efvhhT09PiWp0fuvWrdu/f39ERMSsWbNwKpFtlJeX//rrr6dOnfLx8ZkwYcJdd91lWr5ly5bCwkLTYx8fn6lTp0pXo5OrrKz87rvvzE8HDhzYs2dP0+PDhw+vXbvWy8vrsccei4yMbPRCZ+sRZmVlLV68OPcP5osO/+1vf3v77bdjY2M3bdo0YcIEaYt0bt9+++0777zj4uLSrl27GTNm/OMf/zAt37Rp05IlS8yrBue3WM8777zzxhtvdO7cOSsrCxd9tZl58+b99NNPAQEBOp1u4MCBmzdvNi3/+OOPN27caNrszYkI1nD16tW5c+eadzIVFRWm5Xv27Bk+fHhgYGBlZWX//v21Wm3jVzLn8vHHH0+bNq3RwsrKSh8fn5MnTzLGDAaDn5/fkSNHpKhOFgwGg/nxxo0bO3bsaHo8f/78Z555RqKiZMRgMHTo0OHAgQOMsbq6upCQkN27d0tdlCw03PJfe+21e++91/R41KhR3377rURFyculS5f8/PxuXj5mzJiFCxeaHt93333z589v1MDZeoRElJub+/7773/11Vfl5eWmJceOHfPy8jKNVLi7uw8ePHj37t2S1ujMGg4719TU65q34gAABIZJREFUeHt7m5+eO3du4cKF33zzjU6nk6I0WTh16lR9ff2AAQOIyMXFJSkpCVu7bdxiy9+5c+eHH364YcMGURSlKE1G6urqFi9e/Nlnn509e9a8cPfu3cOHDzc9Hj58+M1/Ec4WhH5+fnFxcZWVlatXr+7ateulS5eISKvVBgQEmNsEBQWp1WrpapSLqqqq119/fd68eaangYGBUVFRVVVVy5Yti4+PxyqwEtPWznGc6Sm2dts7derU8uXLX3zxRdPTLl26eHt7l5WVvfLKK2lpaZgUsB6lUjls2LCrV68eOXIkISHhm2++IaLKykq9Xm+OgMDAQI1G0/iVVu+sWtrKlSsVTamurm7Uctq0aY899hhjbN26dfHx8ebl06dPf/31121atNNxc3O7eRV8/PHH5gbV1dUpKSnTp08XRfHml48ZM+b555+3Yb0ysmXLlsjISPPTOXPmzJ07V7pyZKegoCAqKmrJkiU3/1NVVVVoaOjatWttX5UMffPNN4GBgYyx6upqIrpy5Ypp+ffff3/33Xc3aux4PcLp06cbm3LzlW3vueee3NxcIlKpVBqNxjwoUVxcHBISYuu6nUtNTc3Nq+DZZ581/Wttbe19990XEhKyfPlyc9ekoUGDBplWDVicSqUqLS01HyaGrd2WioqKkpOTX3jhhSeffPLmf/X29u7Vqxe2fNsYNGhQaWmpTqfz8PDw8/MrLi42LS8uLlapVI0aO14Q3lpNTY3pgSiKmzZt6tGjBxH169fPw8Nj586dRKTRaPbv3z927Fgpq3RqdXV1U6ZM8fLyMvXdzcsNBoPpgdFo3LJli2nVgMX16NEjJCRk48aNRHTt2rXs7Oxx48ZJXZQslJSUjBgxYtasWX/5y1/MCwVBqKurMzc4dOhQ9+7dJSrQ+Zn3/0T066+/durUyTRTO378+HXr1hGRKIrr16+/+S/C2S66PXToUHd396CgoKNHj/I8v3PnzqCgICJatWrVyy+/PHbs2D179owfP37RokVSV+q0Pvvss2eeeaZnz54uLi6mJXv27PH09OzVq5dKperYseOBAwf8/PwyMzPbt28vbanOau3atXPmzElPTz9w4EBSUtKSJUukrkgWpk2btn79evMvvNjY2DVr1ly9erVHjx733HOPm5vbzp0709PTv/zyyyaHSaDt3nzzzU2bNnXt2rWoqOjkyZOrV682ncd8/vz5pKSkpKSkkpKSqqqqPXv2NDq51tmC8OrVq4cPH75x40ZERMTAgQMb9kguXLhw9OjR6Oho0wF1YCUlJSVFRUUNl/Tu3Zvnea1We+TIkaqqqqioqISEBJ53ttEIu5Kbm5uTkxMRETFo0CCpa5GL3Nxc85HqROTh4REfH09EFy5cOH36tNFo7N69u2kJWInBYDhy5EhRUZG/v39CQoLpaiom5eXlmZmZ3t7eKSkpbm5ujV7obEEIAADQKvhVDgAAsoYgBAAAWUMQAgCArCEIAQBA1hCEAAAgawhCAACQNQQhAADIGoIQAABkDUEIAACyhiAEAABZQxACAICs/X/myX2Sv6u1bwAAAABJRU5ErkJggg==", "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" ], "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" ] }, "metadata": {}, "execution_count": 6 } ], "cell_type": "code", "source": [ "using Plots\n", "using LinearAlgebra\n", "nucleus = ElementGaussian()\n", "plot(r -> DFTK.local_potential_real(nucleus, norm(r)), xlims=(-50, 50))" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "With this element at hand we can easily construct a setting\n", "where two potentials of this form are located at positions\n", "``20`` and ``80`` inside the lattice ``[0, 100]``:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=1}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd2BUVdYA8HPfzKSHhNSZFAKEFCCUEDqkQRJAmlIEBLGtAkp0hbUri20VxW8tu66ugiArKE0UkJJAKp2AhA4hvUwq6ZNk5r37/TEaMQRImZnXzu+vyeMmc+fxzpx3y7uXUEoBIYQQkiuG7woghBBCfMJEiBBCSNYwESKEEJI1TIQIIYRkDRMhQgghWcNEiBBCSNYwESKEEJI1TIQIIYRkjZ9EqNfr09PT09PT9Xr9ncpcvHgxMTGxqqrq1oONjY1JSUmnTp1qsw7AtWvXEhMTS0tLbz3Y3Nycmpp67NgxlmVNW3+EEEKSobT8W1ZXV0dGRtra2hJCGhoaUlJSevbs2abMs88+u2vXrsGDB584cWL79u2RkZEAkJWVFR0dHRwcXFJS4uvru3v3bqVSCQBvvvnmF198MXz48OPHj3/99dczZ84EAK1WGx4ertFoGhoarK2tExMT7ezsLP9hEUIICR21uHfeeWfSpEkcx3EcN3Xq1LfeeqtNgczMTGdnZ61WSyn98ssvhw8fbjy+ePHi+Ph4SmlTU9PAgQN/+OEHSmlBQYGtre2NGzcopTt27OjTpw/LspTSFStWzJ8/n1JqMBjGjRv3r3/96071yc/PLy0t7UjNOY4z/nFkPgaDge8qSByeYXMzfrnxXQuJM+1lzEPX6M6dOxctWkQIIYQsXLhw586dtxeIjY319PQEgAULFpw9ezYvL894/OGHHwYAa2vruXPnGn/x559/HjVqVN++fQFgxowZVVVVZ86cAYAdO3YYCysUigULFtz+Lq0++OCDrVu3dqTmLMs2Nzd37VOjDmpsbOS7ChKHZ9jc9Hr9XQZ9kEmY9jLmoWu0oKDAz8/P+NrPz6+wsLBNgcLCwtYCjo6OPXv2LCws7NGjR319fevxXr16JSQktCmsVCq9vLwKCwvDwsKKi4vv/i6tWJY9c+ZMay709fUdNWpUuyW533XhU6MOwjNsbniGzc14evEkm1XHL2OGuXd7j4dE2NzcrFKpjK+trKx0Ol2bAk1NTcbBPyNra2udTmdsirX+ovHgnQobDAaDwXD3d2llTIQ1NTXGH8eMGTNkyJB2SxoMBr1eTwjpzMdFnXPr5YHMAc+wubW0tAAmQjPr+GVsY2Nzz1zIQyJUq9UVFRXG1xUVFRqNpk0BjUZTWVlpfM1xXGVlpUajcXd3VygUFRUVxpk15eXlxl9Uq9U5OTmtv1teXu7l5aVSqVxdXSsqKgIDA+/0Lq2srKz+8pe/LF++/J41NyZCW1vbzn1g1Bksy+K0JrPCM2xuxltzKysrvisiZaa9jHkYIxwzZkxKSorxdXJy8pgxY9oUGD16dGpqKqUUAE6cOOHo6NivXz+FQjFixIjWX0xJSTH+4pgxY44cOWIwGADg+vXrN2/eNLbnRo8efXthhBBCqA0eWoTPPfdcZGSkn58fwzD/+c9/kpKSjMddXFy2bds2ceLEGTNmvPbaa0899VRMTMw777zz7LPPWltbA8ALL7ywbNkyKyurnJyc1NTUzz//HADCw8P9/f0XLVo0e/bstWvXPvnkk87OzgCwcuXK+++/383Nrba29ocffjh16pTlPylCCCHh46FFGBoaevDgwQsXLmRmZh44cGDYsGHG4y+88EKfPn0AQKlUpqSkeHh47N27d8WKFa+99pqxwKxZs9avX5+cnHzz5s0jR46o1Wrj8X379gUHB+/evfuRRx5Zu3at8WBUVNSuXbtOnjyZk5OTnJwcEBBg8Q+KEEJIBAj98xItMhQfHx8UFIRjhAJRV1fn6OjIdy2kDM+wuRkny+AYoVmZ9jLGtUYRQgjJGg9jhAghC6hqhowKmllFr9XQ/Hpa1Ag3m6HBQCm1JkSvJOBiTdxtwdee9HWEkJ5kqCsJcsZng5AcYSJESDqaWDhUTPcVcIeLaVEDDXMjg13JUFdyvx/jbQ89rcFBSerr6x0cHPQc3Gyh5TrIq6dZtbAzl752mqtuoRFqJtabTPcjvvaYE5FcYCJESPQoQFIx3Xid+zmPG+JKpvoy30Uzg12Ior1cpmyhjtYAAB62JMgJxsMfhbQ6SC7m9hXSv59hA3qQRf2Yh/oxzjjUhaQOEyFCItZogG+ucZ9e5GwV8HgQ88FIlWc35nKpbWG+PzPfHwycIrGYbrzOvZ6hn9uHeT6ECXbGBiKSLEyECImSzgCfX+bWZrJjPZn1EYpxnqZMVEoGJvuQyT6KMp3iiytc1F5DlIZZPQzTIZImnDWKkMhQgM03uODthmNlNPE+5Y4YE2fBW3nYwqpQ5sY81TA3ErnXsOwIW4n7ryDJwUSIkJhcraHRew3/PM9tiVZsn6gY2NMSTTR7Jbw4mLkyR2WtgIHb9euvcXJ/+hhJCyZChMSBo/DReS58t2FOH+bETOVYs7UC76SnNXw8WnFgivKLy9zkfYaiBsyGSCIwESIkAlodxO0z/JzHnZypXD6AYfgbqhviQo5OV0ZomLBdhp/ycKchJAWYCBESunQtHb7LEK5mDk9V9nbkf7qKkoHXhjI/xSqfP869dJJlsWWIRA4TIUKC9vVVbs4hw/oIxd+HMe0+F8iXUR7k9P3Ks5V02gFDTQvftUGoGzARIiRQFODVU+zaTC59ujLOW0g58Hcu1rBvsjLAiYzfbSjAIUMkWpgIERIilsITqWyKlh6ZoezXQ4hZ0EhB4NMxiscCmfG72as1mAuRKOED9QgJjoGDhclsTQtNmKK0E0OMrhjEuNnAhL3s/smKQS7CTdsItUsMQYaQnBg4WJDENrH0p1iltYLv2nTY4gDGRgGT9hsOTlGGWOTpRoRMBRMhQgLCUXg8lW0w0B9jxJQFjR7sy3AUJu9nk6cqhNydi1AbmAgREpCVJ9i8erp/sviyoNF8f6bBAJP2senTlRo7vmuDUMfgZBmEhOKj89yhYvpznNJWzDeoTwQxfwlmph4w1Ov5rgpCHYOJECFB2JXHfXKB+2WSwkn8+/+9MoQZ4U4WJBnwWXskCpgIEeJfZhVdks7+GKvwkcq+8P8eq2hi4ZVTLN8VQejeMBEixLObzTArkf1ktCLMTSJZEACUDPwwQbkzl/6QjeuRIqHDRIgQnyjA4hTDjF5kvr/UgtHFGnbEKOKPspeqsYcUCZrUYg8hcfkwk7vZDGtGinOS6L0McSFrRioePMQ2GviuCkJ3hokQId6cLKf/d57dMkGhkm4gPhbIDHMlzx/HwUIkXNKNP4SErV4PC5PYz8cpfKUyQeZO/j1OcaiY7sLNC5FQYSJEiB8rT7ARajKrt/Rj0FEFm6IUTx9hS3V8VwWh9kg/CBESoAOF9GAR/ecYaQ4N3m6MB3kskFl2BDtIkRBhIkTI0ur0sCSd/Wq8ooeK76pY0KpQxbUa+v0N7CBFgoOJECFLe+UUG+tNYgS51675WCtgfYRixQm2spnvqiD0Z5gIEbKo42X0x1z64Si5dIreaqQ7mdeXefEEdpAiYcFEiJDlGDhYms5+NIpxFv+Col3zVpgioYimafEReyQgmAgRspzPL3PutiC9RWQ6zlEF/xzNPHOUNeBYIRIM+QYkQhZWpoN3zrKfyWam6J3M7sOobeHzy5gJkVBgIkTIQl7PYB8OYIKd5TVHpl0fj1G8+ytb0cR3PRACAEyECFnGuSq6O49bFSr35qDRAGcyry+z+gzOmkGCgIkQIUt44QT7RqgUNt01lb8PU2zN5q7W4KwZxD9MhAiZ3cEiWtAATwVjuP3B1RpeHKJ4+SSOFCL+YWQiZF4chZdOsu8OZ5QYbX+2fABzppIeK8NGIeIZhiZC5rU1m7NWwAMyWFy7s2wUsHoY8+opHClEPMPgRMiMDBz8/Qz37nAFThVt1+IARquDg0XYKER8wkSIkBn9L4vztoOJXpgH26cg8OYwZlUGNgoRnzARImQuBg7ePsu9FYaPTNzNnD5MowH2FWCjEPEGEyFC5vK/LK63I4xXY3PwbhgCq0KZN89ioxDxBhMhQmbBUvjHOXyCvkNm9WbqWnCkEPEGEyFCZrEtm1PbQqQGm4P3xhB4LZR5FxuFiCeYCBEyPQrw3jnulSHYHOyoeX2ZwgY4UoqNQsQDTIQImd7+AsoQmOyLzcGOUhB4YTCz5hwuNIN4gIkQIdP7IJN9cTCDabBTHg1kTpVzl6qxUYgsDRMhQiZ2qpzm1sPcPhhcnWOjgOUDFR9lYqMQWRrGKkIm9n8XuOcG4sqiXbG0P7MrjyvV8V0PJDMYrAiZUkEDTSjkngjCyOoKV2uY15f59yWcPoosCsMVIVP610XukUDGUcV3PUTruRDmv1e4ZkyFyIIwESJkMg0GWH+NWz4Aw6rrgpxImBvZfANHCpHlYMQiZDLfZXHjPJk+jjhdtFviByr+dQkTIbIcTIQImcy/L2Fz0AQm+ZA6PRzFh+uRpWDQImQaaVqq52CiNzYHu4sAPNOf+Tc2CpGlYCJEyDT+c5lb2h8fojeNRwKZfYVcGT5HgSwCEyFCJlCmg/2F3OIADCjTcLaCWb2Zb65hoxBZAj9x+/333/fu3dvR0fH++++vqqq6vcDVq1fHjBnj4OAwePDg48ePtx7/8MMPPT09e/bs+dRTT+n1egAoLS1dsGCBt7e3nZ1deHj4qVOnjCUPHTrkf4uUlBTLfDQkTxuucw/4Mc5WfNdDQpYEM19d5TgcKETmx0MizMvLe/LJJ7/77rvy8nIbG5uXXnrp9jKLFi2aPHlydXX1X//619mzZxtzXlJS0kcffZSenp6bm5uZmfnJJ58AQH19/YgRI44fP37z5s3Y2Nhp06a1tLQAQENDg4uLS8LvRowYYeGPieSDAnx1hVvSH5uDpjTCnfRQwaFizITI7HgI3W+//TYmJmbcuHE2Njavv/765s2bdbo/DQVkZmZeunTpxRdfVCqVjz/+uI2Nzf79+wFg/fr1jzzySEBAgJOT0wsvvLB+/XoA8Pf3X7Fiha+vr7W19YoVK8rKygoLC41/x8bGpu/v7OzsLP9JkUwkFVMHFYx0x/FBE3sqmPnqCvaOIrPjIRFeu3Zt0KBBxtcDBgxobm4uKipqU6Bfv362trbGHwcNGnTt2jUAuHr1ausvhoSEXL9+neP+FCQHDhzQaDS9evUy/piZmdmnT5+wsLA1a9YYDAazfigkZ19fxTXVzOKhfkxiMVfRxHc9kNQpzfFHKysrt27devvxBx54QK1W37x5szWfMQxjb29fWVnZr1+/1mI3b950cHBo/bFHjx6VlZXG446Ojq0HDQZDbW2ts7Oz8Uh2dvYzzzzz3//+V6lUAkBISMjevXv79u17+fLlp556qqWl5Y033mi3tk1NTfHx8fHx8cYf586du27dunZLGgwGvV6POdWs6uvr+a5C59xsIfsKrNYMbqmrE0cnnojOMAG4z0v19UXDM4FiCjrj6IyVFY4Ym1HHL2M7OzuF4h5bZJslETY1NV26dOn245MmTQIAV1fXuro64xGWZevr693c3G4tdmsBAKiurh46dKjxeG1trfFgTU2NSqXq0aOH8ceCgoKYmJjVq1fPmDHDeMTYIwoAXl5eb7311nvvvXenRGhjY/PZZ58tX778np/LmAhbm6rITFpvd0Rh4yXuPl/q6+pw76KCIaIzvGQgXXaEfTlMTEGHidAyTHgZmyURent7f/bZZ3f618DAwJMnTxpfX7hwwc7Oztvbu02BrKysxsZG48De+fPnlyxZAgBBQUGZmZnGMpmZmYGBgQzDAEBpaWlcXNzSpUuXLl3a7jsqlco2nagImco317g1I+5xv4m6bLyaNLNwspziECwyHx4GNh555JHDhw8nJCTU1dW9+eabCxcutLGxAYD333//+++/B4CQkJDBgwe/8847Op3u888/Z1nW2JR84oknNm7ceOHChfLy8g8++OCJJ54AgIqKioiIiLCwsIkTJ2ZkZGRkZBibzHv27Ll48WJdXd2JEydWr149a9Ysy39SJHmZVbSiCSZ44Xe0uRCARwOZDfhAITIns7QI787Hx2fTpk3PPvtsZWVlbGzsmjVrjMdLS0udnJyMrzdv3rxkyRIfH5+AgIBdu3YZh/0iIiJWrVo1ffr0pqamefPmGUf1cnNzHR0dr1y5Ymw1AsC6deuGDBly7dq1lStXlpeXq9Xq+fPnv/LKK5b/pEjyNl7nFgcQXE7GrBYHkNCd7D9HK6yx4Y3Mg1AqjhF+84mPjw8KCsIxQoGoq6sTywiWgQPfLfrUacoAJzFlQhGd4VYxvxiW9mfm9BHH1FwcI7QA017G4riwEBKghCLax5GIKwuK1OIAZlOW3G/ZkflgIkSoizZlcQ/j4qIWMas3k6bFBwqRuWAYI9QVdXrYV8A92BcjyBIcVHCfL/NDNk6ZQWaBYYxQV/yYy0VqGFdrvushG4v6Md9lYSJEZoGJEKGu+C6LW+CPo4OWE+NFsutotkiW70HigokQoU4r1cGpCjq9F4aP5SgZmNuH2XIDEyEyPYxkhDptWw43zZex4+EpXFlb4M9suYG9o8j0MBEi1Gnf3+Dm+2PsWNoYT1Kvhws3sVGITAyDGaHOya+nV2torDcOEFoaAZjXl3yPjUJkapgIEeqcbTn0gd6MCkOHD/P6Mj9kY4sQmRhGM0KdszWbe1Aka31JzzA3QgicqcBciEwJ4xmhTsito7n1NBq3m+DP3D5kWw72jiJTwkSIUCdsy6EP+DEKzIP8mduH2ZaDLUJkSpgIEeqE7TmcWPZAkKqhroQA9o4iU8KQ7pyPLzNfXObOVMh+8ypZyqunOXU0SoPtQZ7N7kO2Y++oLJXpYFsO9+6vJv7fx0TYOSoCZyrpQ0ls4FbDvy9xLRiMcrIzl87wY5QYNHyb04fZkYv3ovJyspzOTGCDt+u/y6I6A+VM+v+Pa2N0zjPBnK2tAgBOlNG3zrL/vsRtiFSMdMcmgizszOVeHYK7pPMvzI00s3DhJg3piaEnffV6WHGC3VdAXw9lvo9W2SoBAOrqTPkWeHPbRaM8yN5JyrfCmBkHDV9cxoah9Gl1cOkmnYjP0QsAAZjVm+zERqEM3KilI38ysBxcmqNcEszYmqfthomwW+b0YY5MV358gVt9huW7Lsi8fsrjpvgyVhgxwjCrN7MThwml7lwVjdjDPhfCrItQOKrM+EYY1t3l34OkTVfuzKXvnMWwlLIfc7kH/LA5KBRjPUmpjubgrkzSdaWaTtlv+GQMsyTY7HkKE6EJuNtAwhTlxuvchmuYC6WppgWOl9HJvhgvQsEQmO7H7MrDRChNWh1MOcCuGamwzNNKGNim4WkLeyYpXjrFHi/DyJSgvQVcpIaxx7llQvKAH7MrF289JUjPwZxEw6MBzMP9LJShMBGaTJATWReunHeYrWrmuyrI1H7KozOxX1RgJniRzCpa0cR3PZCpvXKKdbEmq4ZZLj1hIjSlab3InD7kqXScOCMpzSwkFHHTsF9UYKwVEOvN7MnHRqGkJBTRrdl0Q6RF1zHE2DaxfwxXXKuhm3HLNAlJLqEDexIPW77rgW4zw4/8hMOEElLTAn9JYzdEKlysLfq+mAhNzFoB6yMUK4+z2GMjGT/nc9N7YaQI0VRfJqmE0xn4rgcykZdPsVN8yASL7+6C4W16w93IAn/mpZPYQSoFFGA3DhAKVU9rCHUlh4qxUSgFx8vo7ny6ZiQPizdhIjSLN8MUB4voCZxBKn5nK6iNEoKcMBEK1PRezG4cJhQ/jkL8UfaDkYyTFQ/vjonQLBxV8O5w5q/HWcyEYrengE7vhVlQuKb3InvycTcY0fv2OmetgAX+/KQkTITm8nAAY+BgWzbeq4rbnnxuGg4QCliAE3FU4faE4qYzwKoM7v9G87bjNUa4uRCA90cqXs/gDJgKRUurg6xaOt4TW4SCNq0X2ZOPiVDEPrvEjfYgPG7jg4nQjCZ6ET8H2HgdM6FY/VLAxXkzKowSYZvWi9lbgFEmVrV6+Og8+/ZwPsMMQ9y83gpTvPsrp8cgFae9+XQqDhAK3jhPklVLtTq+64G65NML3GQfht/5aJgIzWuMBwnoAZuyMBOKTwsHh4u5KT4YI0KnYiDWm9mHjUIRqtPDZ5fY14fyHGUY5Gb3eqji/XMczh8VnTQtDXYmbjZ81wN1wH2+5JcCjDHx+eIyF+PFBPD9eBImQrMLVxMPG9iJy+SLzb4C7j5cX1QkJvswh4pxDEJkmln45CL30hD+o4z/GsjBi4OZDzMxRkXmlwI6xRcHCMXB0xb69SBHS7FRKCb/y+KGuMBgF/6jDBOhJUzrxdS2QKoWo1Q0cuvozWYa5sZ/iKIOmuJD9hXi7aZoUIB/XuBWDuJhQbXbYSK0BIbAX0OYjy9glIrGvkI6yYfBNCgiU3yZfThMKB4JRVRJwPLra7cLE6GFLA5g0rRcbh0Gqjjsw35RsRnhToobaVEDhpg4fHqBfS5EKAlIKPWQPDslPBLAfHEFG4Ui0MJBqpaL9cboEBMFgVhv5kARJkIRuFFLT1XQ+X2FEmJCqYccLOvPfHONa8LdmQQvTUsHOBMLbw2Kum+yD9mPvaNi8MVl7tEAxlbJdz1+h4nQcvx7kDA3sj0HG4VCd6CQm4TP0YtQnA9zqBhX9xW6JhY2XueW9hdQiAmoKnKwJJj5EntHBe9AIZ3sgwOE4qO2BT8HcrIcG4WCtj2HG+ZG+jgKKMQwEVrUVF8muxYuVWOgCldJIxQ10OH8LYSPuiPOhxwswntNQfv6KrckWFipR1i1kTwlA48FkvVXMVCF62ARN9Gb4W1jNNQ9k3yYA4V4oylc12vo1WoqtD0+hVUbOXgskPlfFteCqVCoDhTSSdgvKlrjPMnlanqzme96oDtYf41bHCC4rc0EVh0Z8O9Bgp3J3nzMhELEUThUzMV5YyIUKysGwtUksRjjS4hYCt9ep48GCi7vCK5CcvBoALPxOvbeCNGvldTNmvjYYyIUsVhvJgGfJhSkhCLq6wD9nQUXX5gIeTCnD5Oq5cqb+K4Hus3BIhqL/aIiF+dNMBEK08br3CMBQkw6QqyT5DmoYJov8/0N7L0RnIQiLsYLg0Lcgp0Jy8G1GsyFwlKrh/0F3DzBrCZzKyHWSQ4W9mO+w0QoMI0GOFVOozTYIhS9GGwUCs/OHC7aixHmgk2YCPkR403y6ynetApKqpaGuhIHFd/1QN0W600SMREKzKYsbqG/QO8yMRHyQ0Hgwb7MlhsYqwKSWIQLbUvERC8muQTXWhOQ4kb6ayWdKrDHB1sJtFpysKAv8302RqqAJBbRGHxwQhI8bMHPgZyqwBtNodiaTWf6MTaC2IW3HZgIeTPSg7SwcLYSY1UQynSQV0+H45b0UhHjTQ5h76hgfJ/NzfcXbroRbs0kjwDM9yc/YKNQGA4Vc5EaRokBIRUx3gw+Vi8QuXU0p45OEPA0NIx7Pj3Yl9mWTfGuVQgOFWO/qKSEq8mZCtpg4LseCGBrDp3VW9B3mQKumgwMcSEqBjJwJEMADhXTiV6YCKXDXgnD3EiaFoOLf1uzuQcF+fhgK352CE5JSVm/fj0APPbYY1FRUbcXKCoq+vDDD0tKSiIjI5csWaJQ/DbGun379p07dzo4OMTHxw8aNAgA6urq3n333dZfjI2NnThxovH1/v37v/vuO5VKtWTJklGjRpn7Q3XNnD5kWzY33E2og8jycKOWtrBCXPkJdcdEL+ZwMTfZB4OLTzl1tLCBRqgFHVw8ZOnTp09Pnz59/PjxERERM2fOPHnyZJsCer0+KiqKZdkFCxZ89dVXb7/9tvH4Dz/8EB8fP3PmTH9//8jIyOLiYgBoaGhYu3Ztz9/Z2NgYCycmJj700ENxcXFhYWFxcXFXrlyx5GfsuDl9mO05eNPKs8PYHJSiCV7kUDEGF8+25dAHegt+XzNqcYsWLXrppZeMr1999dWHHnqoTYGtW7f279/f+PrEiROurq5NTU2U0pEjR65bt854fPbs2W+++SaltKSkxMrK6vZ3mTJlyvvvv298vXTp0uXLl9+pPsuXL//ss886UnO9Xt/Y2NiRkp3S7wf9mQrO5H9WpGpray3/pvMOGb65ylr+fXnByxnmRQtLnTa2VDRZ+n2bm5ubm5st/a5CNXKXPqHQ9N9vpr2MeWgRnjx5Mjw83Pg6IiLixIkTdykwYsSIxsbGrKwsg8GQkZERERFx+y+yLLty5cqXXnrp4MGDt/6RdgsL0Ow+ZEcOTm/jDQVIKuEm4kwZyVExMM6TJOPcUf4UNtDsOhEsW2iWMcKGhobq6uo2BwkhXl5eAFBSUuLq6mo86ObmptVq25TUarW+vr6tv+Xi4lJSUuLi4sKy7K2/WFJSAgAqleqJJ57w9/evrKxctGhRfHz8G2+80dLSUlVVdfd3adXS0vLpp5/u3LnT+GNkZOTKlSvbLWkwGPR6PcuyHT4THTLZg1l2QvlysM60f1akGhoaCLFo2FysJg4KVU/aUF9vybfljeXPMI/GuSoO5pNJ7hadPNrS0gIAVlZWlnxTYfr+mmKShjQ1mn6j5I5fxnZ2dgxzjyafWRLh1q1bV69e3fadlMobN24Yq9Xc/Nt5aWpqsrOza1PS1ta2tYCxjL29vbHYrb9ob28PAK6url9++aXx4Lhx46ZOnfrKK6+oVCqVSnX3d7m1YpMmTZo5c6bxR7Va7eDg0G5JYyK0tbW91wnonEgHqD9qKGLtg5zk8vV0F5TSO51/Mzmey8X4WPpNeWT5M8yjKX3oQ0msg4ONJd8UE2GrX0oMzw9izHH+TXsZmyURPvbYY4899tid/tXX1zcvL8/4Oi8vr7Xxd2uBc+fOGV/X1tZWV1f7+Pg4OTk5Ojrm5eWp1eo7/WJISEhTU9PNmzfd3d29vb3z8vKMM0vbLdyKYZigoKCYmJjOf1DTIAAz/ciPufTlIZgIeXC4hD4k1LWAUTcNcSHlOlrcSL3s8L/Y0iqb4WwlFcX6vTxUcc6cOTFc3qcAACAASURBVJs2beI4jlL67bffzpkzx3h88+bNxt7OOXPmJCQkGCeFbtq0afjw4cY0NmfOnI0bNwKATqfbunWr8RfLyspa+yrXrVvXu3dvd3d3AJg7d66xsMFg+O6771rfRZju92N+ysORDB6wFNK0XLRGBLGKuoAhEOXFJOHcUT7syedivIW7vuiteHiOcNmyZTt27AgLCyOEKBSKp59+2nj8qaee+vHHHzUaTf/+/ZcsWTJ8+PABAwacP3/+xx9/NBZYtWpVdHR0RERESUlJSEjIjBkzAGDz5s3vv/9+UFBQWVlZQ0PDd999Zyy8cuXKCRMmjBo1qqGhwd3dffHixZb/pB0XqSHXamhJI2ju2IOLzOJsJfW2Ix4m7u1GAhKtIUkldGE/vushP7ty6ew+4miIE8rHCl8cx509exYAQkNDW4cxtVptz549ra1/27cxNzc3Pz8/NDTU0dGx9RdbWloyMjJ69OgxcODA1oN5eXn5+fnOzs7BwcEq1R+7yRkMhjNnzlhZWQ0ZMuQuw6rx8fFBQUHLly+/Z7XNNEZotCiZjVCTp4Ll3jSpq6u79X/c3D7M5Aoa6KdjxHDXaiIWPsO8u1RNpx9gb8yz3E0/jhECgM4Ams36nHmqnubZide0lzE/K8swDBMWFtbmoHHwr1Xv3r179+7dpoyVldWYMWPaHPTz8/Pz87v9XZRK5ciRI7tbV0uZ0YtsvM5hIrSwpBLuySA851LW35k0GmhePfVzEEfrRBoOFdNhrsRMWdDk8CtAKCb5MOlaXCPYogwcHC2lEThAKGkEIMqLSS7BYUKL+jmfm+EnmsgSTUUlz8kKRnmQhCKcMmM5pytoH0fiKpKbVtRlURqCidCSKMCefG56L9E0wTERCsi0XszuPAxXy0kqodGCX/MCdV+UhuDEUUs6XU6drYh/D9EEFyZCAZnWi/xSwHEYsJaSXMwJf/En1H1BTkTPQXYdhpaF7MnnpomnOQiYCAWlryNxsSa4PaFl6Dk4XkbD1RgCsoC9o5a0t4BO6yWmyBJTXeVgai+ytwCHCS3hdAXt5ySaWW2om6I0JAUToUUUN9LcOjrWA1uEqKvu82V+KcBwtYTkEhEsio9MJRJbhJayr4DG+jBKUeUWUVVWBsZ5kqxaqsWNKMwvuZiLFPau2ciEAp2IAYcJLeKXAnqfr8gi626JsLKy8sqVK1lZWfUy2Z9GAFQMxHgx+7F31MxwgFCGsHfUAlo4OFzMTfERWWS1rS6l9PDhw48++qivr6+bm1v//v0DAgIcHR0HDhz4/PPPX7hwgZdaysp9vmRfIYareWVUUP8eOEAoL5GYCM3vSCkNciZuFt32ygT+tMTa3r17X3rppYsXL/bv3z8uLi44OLhnz54Gg6GysvLcuXNbt279+OOP4+LiPvjggyFDhvBVY8mb7MusPKE3cApxdbKLS3IJjcQBQpmJ1JD3zmFfi3ntKxBfcxBuTYTbt29//PHHly1btnXr1gEDBtxelFKanp7+1VdfjRw58urVq7cvBIpMQm0LvR3J8TI6HkewzCalhFvaX3zhirojyIk0s7joqHntK6DrIsQXWX/UODQ0NDs7e82aNe1mQQAghISHh3/77bfnz593dna2VA3laLIP2V+It67mYuDgGA4QylKEmsHeUfMpbKClOjrcTXz3GX98F/j7+7u5uXXkdwIDAzERmtVkH+YADhOazdlK6udAXHCAUH5wmNCs9hfSWG+GEV8evMOsUeNmgbfbsWOHOSuDfjPag9yoo2X4EIV5pGhpBHY7y1KEmqRqMRGay4FCOllsD04YtZ8Ip02btnr1apZlW480Nzc/99xzc+bMsVTFZE3FQLSGwZ0ozCQVZ8rI1YCepKaFFjVgLjQ9AweHi7k4b1GOOLRf6aVLl77zzjuTJ08uLS0FgKtXr44ePfqrr7768ssvLVs9+ZrkQ7B31Bw4CkdKORwglCcCEKFhsFFoDifLqZ8D8bTlux5d0v7XwRtvvHHw4MGLFy8OHTr0zTffHD58eHNz84kTJ5566ikL10+2JvmQhCLciML0zt+knrZiDVfUfeGe2DtqFgcKuUk+Yu1oueN98YQJE44ePdrQ0LB69Wp/f/9Tp04NGjTIkjWTOT8H4mRFzlVixJpYSgkOEMpahIak4nwZMzhYRONE+ASh0R3rfenSpWnTplFKZ86cee7cuWXLluFCaxYW50MOFmHEmliqlkbgAKGMDXEhJTpa3sR3PaTlZjNcrqbjPMUaWe0nwg0bNowYMUKhUJw6dWrXrl1btmz56aefhg8ffu7cOQvXT85ivclBfJrQpChAmpbDFqGcMQTGeZI0LUaWKR0q5sZ7EiuxNgjvkAhffvnlRYsWHT9+PDg4GADmz5+fkZFhb28/evRoy1ZP1qI0zMlyqjPwXQ8JuVpNHZTExx4ToayFq3G+jIklFNFYcc4XNWq/6t98882XX35pa/vHjIJ+/fodPXr0ySeftFTFEDiqINSVpGDEmk6qloZjc1D2ItQkDcPKpA4W0VjRzpSBOyXCKVOm3H7Q2tr6008/NXN90J/EeuPThKaEA4QIAMLcSFYNrWnhux5Scb2G6jkY4CziyPojEVZVVXFch75zGxoadDpc9cQSYr1JIs6XMZ00XFMGAagYGOFOjpRiZJlGYjGN8RJ3WP2RCPfu3Tt48ODNmzc3Nd1xQlVlZeXatWv79u1bUFBgkerJ3XB3UtBAS/GuwxTy6qmeo/16iDtikUmEq5lUnC9jIglFNNZb3GH1xzZMc+bMKSgoWLp06dNPPz116tTRo0cHBga6uLgYDIaqqqpz584dPXo0ISFBrVb/3//9X0BAAI+Vlg8FgSgNk1jELewn4oFogUjV4o4T6DfharIqAxOhCbAUkku4z8ep+K5It/yRCG1tbV999dWlS5du2LBh48aNW7ZsofSPrgOVSjV27Nh169bNnTvX2hrX7becGC9yqJgu7Md3PcQP+0VRq9Ee5FwV1RnAVnnvwuguMiqojz1Ri3ypprZXgYuLy4oVK1asWFFRUXHhwoWysjKVSuXh4REaGmpnZ8dLFWUuxpusycRbVxNI09LlA7BFiAAA7JQw2IWcKKdROHmqexKL6ESRDxBCm0TIsqxCoTC+dnNzi4qK4qFG6M8CnQilcK2GBjqJ/mrjUZkOSnU0pCeeQ/Sb8Z4kVYuJsLsSi7iVgxV816K7/nSDHBYWtmTJkjNnzvBVG9SuiTh3tNvSS7mxHkSMW4YiMwlXM0dwvkz36AxwqoJGin/E4U+JsH///uvWrQsLCxs2bNhnn31WVVXFV7XQrYzDhHzXQtzScKYM+rNxnuREOTVgKuyG9FI61JU4iHuiDECbRLhly5b8/PyPP/7YYDA8++yzXl5e06dP37Ztm8GAy3zxaYIXk1LC4Z5M3ZGGa8qgP+tpDX4O5NcqjKuuO1zMTZBE33Lbe2QvL6/nnnsuMzPz9OnTS5YsOXr06IMPPujn5/fyyy9fv36dlyoijR142pKzuCVTV9Xp4WoNDXOTQsQiEwpX45ZM3XKomE4U8xKjre74GcLCwj755JOCgoJNmzb179//ww8/DAoKioyMtGTlUKuJ3uQw9o521bEyGuZGrEU/oo9MbLwnScf1ZbqqugWuVtPRHlK4v7xHMrezs1u0aNG2bdv+9re/AUBqaqpFaoXamuBFDhXjaEYXpZZw2C+KbjdeTdK1OObQRSkl3Bgxb710q7t9CJZl9+3bN2/ePC8vrw8++GDAgAEfffSRxWqGbhWpZo6V0hZMhV2SXkrHe0oiXpFJ+dgTBxW5Uo2psCsOF9MJXhIJq/aXVbh69eqWLVs2bNiQl5fXo0ePuXPnLl68OCYmxsKVQ616WkOgEzlZRsdjy6aTWjjIqKBjRLt3NjKrcDVJ09L+Yt45gS+HiumGSCkmwpqamp9++mnTpk2HDh0ihIwZM+bVV19duHChvb09X/VDrSZ4kcMlmAg77XQ5DXQiPcQ/wxuZw3hPkqalTwXzXQ+xKdVBSSMNdZXI19GfEmFYWNiNGzf8/PxWrVr16KOP9u7dm6daoXZEezHvn2NXhUrkFsxi0rR0PDYH0R2MV5P3zuGQQ6clFXMRakYhlcD6UyK8//77J0+ePGHCBIbBb1vBCVeTuYdwmeBOSy/lHgnA6xm1L9iZ1OtpYQP1sZfKl7pFJJXQaPEvMdrqT18Qa9eujYmJwSwoTPZKGOJCjpbhwH4ncBSOltLxuKYMugMCMM6TSddiWHXO4WIaLYlH6Y3wC0JMor1IEj5E0RmXqqmrjej3iEFmFa7Gpwk7p6CB1uppiAsmQsSHaA2ThAthdAYOEKJ7GudJsEXYKcklNFItqRXsMRGKyVhPcr6K1uv5rod4pGtxni26h2FuJLuOVrfwXQ/xSCqW1AAhYCIUFxsFhLlhN04npJdiixDdg4qBEe7kKIZVhyWVSGqAEDARik6UhsFhwg7Kr6dNLA3ADY3RvYSryZFSDKsOyamjzSwNltYSBJgIRSbaiyTjMGHHGFdWk1S8IvPAiaMdl1xCozVSSxxS+zySN8qdXK6mtThM2AE4QIg6aLQHOVNJm1m+6yEGySU0Slr9ooCJUHSsFTDSnaTh3WsHpJfiZryoQxxVEORETldgWN1bUjEmQiQAURomGYcJ7+VmM+TV0aESetQJmZVx9W2+ayF0N2opBZDeuDsmQvHBYcKOOFpGR3oQJV7gqGPGeeJ8mXuTZL8oYCIUoxHu5GoNrcHHnu4qXcuNwwcnUIeNVzNHSynu0nt3ySU0EhMhEgIrBkZ5YDfOPeBmvKhT1LbgYk0u4Sa9d5UsuScIjfCbQpSiNExyCXbj3FEzC2cr6GgPCUYsMp9wNa61djdZtZQA+PeQYFhhIhSlKA0OE97NqXLa35k44Ga8qDPGeeKyTXcj1QFCwEQoUiPcyTUcJrwzfHACdcF4bBHeVYpEBwgBE6FI4TDh3eFMGdQFgU5Ex9KCBgyr9km4RcjPZucsyx4+fLikpCQqKqpXr17tlsnMzPz1118DAwNHjx7delCn0x08eLCxsTEuLs7V1RUAGhsbjx49eusvBgcH+/j4lJeXnzt3rvVgaGiosbxkRGmYlBJuWi8F3xURHI7C0TK6LgJv8lDnEIDxnky6li7wl+bXfXdk1VKQ6AAh8JIIKaXTp08vKysbPHjwihUrtm7dOmHChDZlPv/887fffnvGjBmrV6+eN2/ee++9BwD19fVjx451d3dXq9V//etf09LSAgMDq6qq1qxZY/wtlmWTkpJ+/vlnHx+fY8eOPfLII8OHDzf+0z/+8Q/JJULy12M4X6YdF29SdxviiZvxos4zDhMu8Oe7HsKTIt3mIPCSCBMTEy9evHj58mU7O7svvvji9ddfb9Ok0+l0b7zxxu7du8eOHZufnx8UFBQfH+/l5bVhwwYnJ6eEhASGYZ5//vn33nvvm2++8fHxSUhIMP7i3r17L126NHnyZOOPISEhrf8kPa1PEzpZ8V0VgUkvpdgvirpmvJpsvI73l+2Q8AAh8DJGuHv37qlTp9rZ2QHAnDlzjh07VlFRcWuBo0eP2tjYjB07FgB69eoVFha2f/9+4y/Onj2bYRjjL+7evbvNX16/fv3ixYtVqt8mCzY2NiYmJmZkZOj1Elyj2oqBke7kCE5yuw2utY26LNSV5OAmve1JLqGR0g0rHlqERUVFoaGhxtdubm42NjZFRUVubm63FvD29m790dvbu7CwsM1xHx+fyspKnU5na/tbF1hFRcXevXvPnj3b+os6ne7jjz++evWqSqXas2dP3759260Py7JpaWkKxW+Dbf7+/tHR0e2WNBgMer1eqeRnYPV24R5wuMgQq+a7Hial1+u7eeOSpoXXBoNej1sJtK/7Z1jahrtBapF+ik/X/4Lx9BIinbSRUw8shd52BuFcOB2/jJVK5T3/L8zynX7t2rX58+fffnz9+vVDhw5lWbY16wCAQqEwGAy3FmNZ1tjs+62KSqWxwK2/aHzBsn982W3cuDEsLKx///7GH6dOnTpjxgwA4Dhu8eLFf/vb33bu3NlubTmOKygoyMjIMP5ICImIiGi3JPu7u398ixnvQV7OICwrqZ6cbp7hggZoYpk+9pxg/pcER1DXsACNdSfpWojTdL2vxXh6pXSSk4tJuIewvmo6fhkrFAp+EmGvXr02bNhw+/F+/foBgEajKSsrMx6pr69vaGjw8vK6tZhGoykvL2/9sbS0NDw8vM0vlpaW9ujRw8HBobXY+vXrV65c2fpja8pkGGbBggXLly+/U21VKtVDDz10lwKtDAaDQqGwsbG5Z0nLGO8NV5P0eoWNo4SeHNfr9d05w6eKuHA1tbXBgdM76uYZlrxIb/r2WdbGpuvfjcb7eCsr6VyERyrYCT7ExkZAM7FNexmb5YPZ2NgMbo9xXDAyMjIxMZHjOAA4ePBgUFCQWq0GgIaGhqamJgAYOXJkSUlJVlYWANTW1h4/fjwyMtL4iwcPHjS+RUJCgvGg0fHjx/Py8ubOndtufc6ePevj042eDqGyVsBwN3wE+E9wgBB102gPchY36f0zaU8ZBV7GCGfNmvX2228vXLhw1KhRH3zwwfvvv29st86bN2/o0KHvvPOOi4vLsmXLZs+e/Ze//GX79u2TJ082dnguWbJkyJAhzz33nFqtXrNmzd69e1v/5rp16+bNm+fo6Nh65Pnnn6eU9urV69KlS99///2uXbss/0ktIFLDpGq5Kb74NOFv0rT08SAB3bci0TFu0ptRQcfi3GMAAMipoy0cDZTcHoS34uErw8rK6siRIyNGjCgtLf3uu+8WL15sPL58+fKZM2caX69du/bVV18tLCx8+OGHt2zZYjyoVqtPnz7t6enZ2NiYnJw8bty41r85YcKEV1999dZ3efTRR729vUtLSwcMGHD+/PmYmBiLfDhLi9KQJFx09HdVzZBfj5vxou4KV5NU7Gj5XaqWRmkkfnNJKJX7/3d8fHxQUFAHxwj1en3rPFUhaGLB43/64odUkllguq6u7taWfafsyaefXGQTpghlWq8wdecMy8T2HG7jdW53XBcvpJaWFpDQGOHjqexId7K0v7ByoWkvY2F9NtRZNgoY5kaOlsn9bsboSCmHexCi7sNNem8l1c14b4XfGqIXqSbJxQKa1syjNJwpg0xBbQuuNuTiTcyEkF9PGw002FniYYWJUPSivJgUHM8AaGLhXBVuxotMYzzuTQgAAClaGqlhJB9UmAhFb7Q7yayijYZ7l5S2U+V0gDOxx/FBZArhatzmDMC4xKgMelkwEYqerRJCXclR2d+9Yr8oMqHxnviELoDU19puhYlQCiI1JEUr92HCNC03Hh/8QiYS4ERaOJpXL+tcWNhAa/V0QE/phxUmQimIVDPJ8n6akKVwvIyGq/F6RiYTrmZk3ihMKaERaukPEAImQmkY60l+raQ6GQ8Tnq+iXnbEDVfQRKaDw4QpWln0iwImQmmwU8IQF3JMxk8T4gAhMjmcOCqTmTKAiVAyIjUkpUS+w4TppTRcHhGLLGaIKylqoJXNfNeDJ8WNtKqZhshjwUJMhBIRqZH1MGGalsNEiExLQWC0B0mX6zS0lBIaLo8BQsBEKBnjPMmZStoky71jrtdQJSF+DjKJWWQ54WpGtsOE8ukXBUyEkmGvhEE9yXFZDhNivygyEznPl0nR0igvuYQVJkLpiNSQFFn2jqZpMREisxjhTi5V0wb5zcfW6qBcRwfJ4AlCI0yE0hGlYZJlOV8GEyEyExsFhLrKsaMlpYQLVzNyGSHERCgl4zxJRoXshgmLG+nNZlksfoF4Ea4mqfK7v5TD1ku3wkQoHQ4qGNiTnJDZ3Wu6lo6Xzdw2ZHnhakaGu9WnlNAoTIRIpCI1RG5bMqVpaQT2iyKzGetJMipoi5zahKU6KNXRwfJ4gtAIE6GkRGkYuW3Sm6qlEXK6dUUW1kMFgU7kdLmM7i/lNkAImAglZpwnOS2nYcKqZsirp0PldOuKLC9CZg9RyG2AEDARSoyjCgb0JCdlc/d6pJQb7UGUeBUjc4rQkFQ5rS+TLLMBQsBEKD1RGiKftdZSS3DrJWR24z2Zo6WUlUdUGQcIh8islwW/RKQmSsMkyWaYMBVnyiDzc7MBb3tyrlIWmTBZfgOEgIlQesbL5mnCej1crqYj3WUWsogPEWoik4cokktotMz6RQETofTI52nCo2V0mBuxVvBdDyQDkRoZJUK5zZQBTISSFKUhSTJYCyO1hJPP6viIX+FqkqblJJ8JSxqhXGZPEBphIpSgaC8mqVjyMQspWhqhwQsYWYKXHelpTS7elHhYJZdwERrZDRACJkJJGudJzlZSnaSXzG80wK+VdIyH/EIW8SRSLf3dXWT44IQRJkIJslfCYBdyVNLDhMfL6BAXYqfkux5INuSwzVmSLGfKACZCqYrSEGlvyZSq5fDBCWRJkRqSKulhwsIGWt1MQ+Q3QAiYCKUq2os5LOlhwuQSGokDhMiCfO2JnZJcrZZsWCWV0CiNDMcHATARStVYD3K+itbr+a6HeTSxkFFBx3nKM2YRbyIlvWxTUjGN9pJpTGEilCZbJQxzI0dKpRm0J8poSE/ioOK7HkhmoiS9zZlsZ8oAJkIJi9YwhyW61po8n/lFvItUkxSJDr3n1NEmlgY7yzSsMBFK1gQvkiTRbpzkEi4KBwiRxfV2JNYKcrVGgmF1uJhGy3WAEDARStgoD3K1mla38F0PUzMOEI7HAULEB6nu7pJUQifIdYAQMBFKmBUDYzwl2JNzHAcIEX+iNESSyzbJeaYMYCKUtmiNBB+iSC7hZDukj3gXpSEpJVJ7mvByNbVSQF9H+YYVJkIpm+BFpJcIk4pplBdet4gffg7EVkkuS+tpwsPFdKKMm4OAiVDahrmRokZaquO7HqbTaIAzlfgEIeJTtIYkS+v+8nCxTFdWa4WJUMoUBCLVknqI4lgZHepK7HGJUcSfaC9yWELzZVgKySXcBHn3ssj6w8vBBC9ySEJ3r0nFOECIeBatISklnGTGCc9WUo0d0djxXQ9eYSKUuInekkqEh4upzG9dEe+87YmLNTkvlb0JDxXJfYAQMBFK3gBn0sLCjVopBG2dHi7cxD0IEf+kNA3tUDGHiRATofRNlErvaKqWjnQnNgq+64Fkb4IXkcbQexMLx8twIxdMhDIw0ZskFkkhER4u5qKxXxQJQLSGSddSg/hT4dFSGtKTOFnxXQ++4deK9MV4kSRJjO3j005IINxswM+BnK4QfVAlFnETvTGmMBHKgLc9cbchZyvFHbTlTZBbR4e7YdAiQZDGNLTEYhqDvSyYCGUiRvy9o0nFXLiaUeIFi4RhoheTWCTuvtGqZrhaTcfg8hSYCGUi1pskiDxoD2G/KBKSCDXJqKCNBr7r0Q2Hi7nxamKFSQAToUxEaZiT5VQn5qBNLKIxOJiBBMNBBUNcSXqpiDtaEoporDemAABMhDLhqIJQV5KqFWvQZtdRHUsH9MREiAQkRuS9owlFNBZvLgEAE6F8xHoz4u0dTSiiMV7y3T4bCVOsmIfes2qpnoOBeHMJAJgI5SPOhxwUbdBivygSoJHuJLuOljfxXY8uOViIzcE/YCKUizA3UtJIixvFlwtZCoeLOUyESGiUDERqmEPi7Gg5WETjMKZ+h4lQLhQEJnoxCSJsFGZUUG974mWHQYsEZ5I3EWNM6TlIKeFicKbM7/BEyEicDzlQKL6gxSF9JFix3qIccThWRgOciJsN3/UQDEyEMhLnTRKKxLfW2sFCLg5vXZEgBTgRJQMXxbYl0/4CbpIP3lz+gbfvF51Op9Pp7lLAYDDU1NR08K+xLFtdXX378cbGxqYmcY5lm4GPPdHYiWyBxDo9/FpJI9QYtEigxNg7eqCITvbBm8s/8HAuDAbDY489plarNRrNo48+ajC085j3N9984+npGRAQEBoamp2dbTy4ZcuWiRMnenp6PvPMM7cW/vHHHzUaTVBQUHBw8Pnz540Hm5qa5syZ4+3t7enpGR8fT6nIrlQzmeRN9ouqd/RwMTfag9gq+a4HQncwyYccKBTTfJlSHeTW0VHueHP5Bx4S4YYNG86cOVNUVFRcXJyZmbl+/fo2BYqLi+Pj4xMSEsrKymJjY+Pj443HHR0dn3nmmenTpzc2NrYWrq2tfeSRR3744YfS0tLHH3/8iSeeMB7/9NNPtVqtVqvNy8s7cODAjh07LPPpBG6yL7O/QExBe6CQTsJbVyRgE7yYo6W0ieW7Hh22v5Cb6IXL9v4JDyfj22+/Xbp0qYODg52d3dKlS7/99ts2BTZv3hweHj5s2DAAWLFixYEDB8rKygBg2rRps2bNcnd3v7Xwjz/+GBgYGB0dDQDPPPPM+fPnr1y5YnyX+Ph4a2trZ2fnxx9//PZ3kadwNblUTSub+a5Hh+0vpDiYgYTMyQqGuJKUEtF0tOwvpJN9Mab+hIdEeOPGjeDgYOPr4ODg1p7PVtnZ2a0F1Gq1o6Njbm7unf7arYXt7e19fX1v3LjR5ni779KKUlpRUZH9O2PSlSorBqI0zEGR9ORcraEsxcUvkNBN9mH2iySmWAoJhdxkvLn8M7OMvZw4cWLfvn1tDjIMs2rVKgCora21s7MzHnRwcLh9kktNTU3v3r1bf3R0dGx3Ikxr4da/1lq4ublZp9Pd/V1atbS0/Otf/9q4caPxx/vuu+/9999vt6TBYNDr9e0OaopItLvi52xmmoee74q0r76+vvX1rizlRE9SV3e3SVWos249w8gkIlzI48es3hr425BNS0sLAFhZCXHf9+MVjLedypGtr6vjuyrd0/HL2M7OTqFQ3L2MWRIhwzAqler2g8YX7u7urdNBq6urPTw82pS8tcCdytxaODMzs01ha2trJyenu79LK2tr69WrVy9fvvyen8uYCG1tbe9ZUshmBdB3LxjsHGwUQr0pdHR0L/fELgAAFORJREFUNL5IKjcsCWYcHfFxJxNrPcPIJMY5Qn26vhwc+joSEHYiTL7KTveTyAVgwk9hlkQ4YsSIESNG3OlfQ0JCMjIy4uLiAOD06dMhISG3F1i3bp3x9dWrVw0Gg7+//53+2sCBA7/55hvj6/Ly8oKCggEDBhiPZ2RkDB8+/E7vIlu+9kRjR06W0zEeQs2EAADQYICjpXTbRBzTR0JHACb5MPsK6DMDBB1TALA3n34+7h7NIxni4Vtm6dKln3766YkTJ06dOvXJJ58sXbrUeDwuLu706dMAMG/evGvXrv33v//Nycl58cUXFy5caMz8+fn5iYmJubm5xcXFiYmJWVlZADB16tSWlpY1a9bk5eWtXLnyvvvu8/b2BoBly5atWbPm3LlzaWlpX3/99ZIlSyz/SQXrPl/yi+Dnjh4u5ka6E8e2PQsICZEoYqqwgRY10lHCvgPmBQ+J8L777lu9evXTTz+9ZMmSN954Y9q0acbjPXr0UCqVAODo6Lhv375t27ZNmTLFy8vrn//8p7HAqVOn1qxZU1FRwXHcmjVrkpOTAUClUu3bty81NXXSpEkMw7Q2JRctWvTMM88sXrx45cqVH3/88bhx4yz/SQVrqi+zJ1/ok9x+KaBTfLE5iMQhzptJ1wp97+s9+XSyDyPYMREeEXzSPD4+PigoSD5jhADAUlB/pz/zgNLXXnAxUVdX5+joSAH8thgOTlEEOwuuhmJnPMN810KCovca/jZYMdWXCHaMcPpBw6J+zLy+Uri/NO1lLIUzgjpLQWCKj6AbhZlV1FoBmAWRiEztxezJF27vaKMBUktweYr24UmRqel+RMhBuzuPTuuFWRCJyTRfsidfuD1sh4q54e7EWXDNVEHARChTk32YdC2tF+jDhLC3gJvWCy9OJCbBzsRWCb9WCjQV/pxHZ2BM3QGeF5lyVMEYT3JQkJtra3VwrQZ3nEDiM70X2S3IEQeOwp58boYfxlT7MBHK10w/5qc8IQbtnnxukg+jwmsTic1MP+bnPCHeXB4vox62pI8jJsL24ZeNfM3oRfbmcwbhhe3PeXQGDhAiERrnSfLraVHjvUta2K487n5sDt4ZJkL58rYn/ZxIilZYjcIGA0nVcvgEIRIjBYH7fJndBYJLObvy6P29MabuCE+NrD3gx/yYK6wmYUIJGe1BnHBuGxKn+3uTXQIbcbhwk+o5CHUVXHoWDkyEsvZAb7Irj3JCCts9RYoH8NYViVacN3OmilQJacvPnbn0AewXvSv8xpG1QCfS0wqOlwklE7ZwcLCEud8PL0skVnZKmKCBvYV81+MWO3K4WXhzeVd4duRudh+yUzC9o4lFtL8T9RT9GnZI1mb60h/zhNICu15Dy5voWE+h1EeYMBHK3Zw+zPYcoSyHsSOHu99HKFkZoa65zwfSSqFWGKtVbM+ls3szDObBu8JEKHchPYm9Ek4KoHdUz8HP+dwMH2Ev4I/QvThZwXhPKpAlDLdlc3Mlscq2WeEJQvBgX2ZrDv9Be6iYBjoRbzu+64FQt83pDduy+b+5vF5DS3UwHvtF7wUTIYIH+5Kt2fzPHd2azT3YBy9IJAXTfCCphOO9d/T7bDqnD8F+0XvC7x0E/Z2JqzWkl/KZCZtZ+DmPm9sXQxZJgZMVRGmYXXxPQ/shm5vvj1/y94bnCAEAzPdnttzgM2j3F3KDXIiXHSZCJBHz+xJ+Y+pcFW0wwGgPjKl7w0SIAAAW+JMdOZyev7DdcoMuwFtXJCEz/JjjZbRMx1sFNmdxD/kTTIMdgV89CADAz4H0dyb7CvjJhHV62F/IzcEBQiQhdkqY1ou3aWgchc036KJ+GFMdgqcJ/ebhAGZTFj/DhDtyuGgN42LNy5sjZC6L+jH/y+InESaXUA9b6O+MDcIOwUSIfjOnD5NYxN3kY43ETVncwwEYsUhqYrxJQT1creHh/nLjdW4xNgc7DM8U+o2zFUz25WHKTF49zayiU3HfJSQ5CgIP+ZNvr1s6pur08HMe9xAmwg7DM4X+8GgAs8HiQbvxOp3fl7FWWPhtEbKERwKZb69T1rJtwm05XJSGcbex6JuKGiZC9IdYb1Kqg8wqy0UtR2HDNe6xQLwOkTSF9CTe9pBQZNFM+PUV7okgjKlOwJOF/sAQeCyQfH3Vco3Cw8W0pzUMc8MBQiRZjwcyloypizdpfgNM8cWY6gRMhOhPHg9kNmdxjZZa+Pq/V7i/4K0rkrQF/szhYk5rqQcKv7zCPR5IFJgHOwO/g9Cf9HIgYz2Z77MtcQOr1UFiMbcQh/SRpDmqYE4fZr1FGoWNBtichTeXnYbnC7X19ADm80uWCNqvr3AP9mF6qCzwVgjx6en+zH+vcBaYMrP5BjdezfRywPZg52AiRG3FeZNaPRw18xrceg6+vMI9PQCvQCR9Q12Jtz38nGf2+8vPLnLPYEx1Hp4y1BZDYPkA5tOL5g3anblcvx4w2AVvXZEsPDuQ+czMMZVUQlkKMd4YU52GiRC147FAJrGIy6s3Y6Pw4wvccyF4+SG5mNWbyaqFs5VmjKn/O88+NxA3H+wK/CZC7XBUwRNBzMcXzHUDm66llc0woxdefkguVAw8G8J8dN5cMXW5mp4upw8HYEx1BZ411L7nQphvr3OV5ll6dE0m+7dBDN67Ill5Kpg5UMjl1pmlUfhBJrd8oMIGV2jqEkyEqH1edmRuH+bjC6zJ//K5Knq2Eh7BW1ckMz1U8FQw80Gm6RuFuXV0Tz5Ok+k6PHHojl4awnxx2fT7Ubx1hls5CBcXRXL01xDF1myusMHEjcL3znFL+zPOVqb9qzKCiRDdUR9HMqs3s/a8KRuFZyvp8TK6JBgvPCRH7jbwZDDz7q+mbBTeqKU7c7nnQ/DWsuvw+wjdzRuhzJeXuZJGk/3BV0+xrw1l7JQm+4MIicsLgxU7crisWpM1CldlcM8NVOC+1t2BiRDdjY89eSKIWZVhmkZhYhHNqoUnsTmIZMzFGlYMUrx8yjSNwtMVNEVLnx+EMdUtePrQPbw6VLEnn+v+808GDp4/zn44ilHhRYfk7bmBTEYFTS7pbkxRgOeOsW+HMfbYxdI9+J2E7sHJCv4xQvH0EZbrXtj+6xKnsYP7/fCSQ3Jnq4SPRjHxR1l995qFG65xLMUJ2CaAZxDd26OBjIqBzy93PWpz6+i7v7L/Hovj+QgBAMzqzfg5QHcepSjVwSun2C/GKfB53O7DRIjujQB8Fa546wzbtRF+jsJf0tgXBisCnDBkEfrNf8YrPr3IZlZ1sadlSTr7RBAz1BVjygQwEaIOCXIifx+meCiJben8LeyHmVwLBytxPB+hW/jak7WjFA8lsQ2d3wf7P5e5ggb692HYxWIa+N2EOurpAUwvBxJ/tHMzSA8V008uspujFbhlNkJtPNyPGeFOnkzr3E6Fx8vom2fYHyYorPD720TwRKKOIgDfRCiOl9EPOzywcfEmXZhk2BKt9LHHNIhQOz4fq8iuo3/v8BNKWbV0diK7LkLZrwfGlMlgIkSd4KiCXyYp/nOZ+6QDG1NcqqaT9rMfj1ZEajBiEWqfrRJ2xym35dB/dGC5maxaGvML+1YYM9UXY8qUMBGizvG2J0lTFf+5zK08cbfJ3wcK6YS9hg9HMvP98RpD6G7cbeDwfcofsrllR9jmO7cMU0poxB7DG6HME0EYUyaGJxR1mp8DOTZDmVULo382pNz2UHBJIyxJZ59MY7dPVC7ALIhQB2jsIG26srIJhu8yJBS1jakyHTx7jF2YzH4bqcQsaA64IAHqip7W8FOsYssN7sl01pqBWG/iaUt0LD1VTo+V0ccDmfOzlU64Fj5CHdZDBVsnKnbmcn89xlKAOB+isSXNHJwup+ml3KJ+zLlZSldcUNQ8MBGirlvgz8z3Z46X0TQtrWqm1gw8HshsmcD0UPFdM4TEaVZv5oHezKlymqal5U3UioGHA8imKBXeVpoVJkLULQRgjAcZ44FD9wiZBgEY6U5GumNMWQ52NyOEEJI1TIQIIYRkDRMhQgghWcNE2An/+9///va3v/FdCykrKysLCwvjuxYSN23atF9//ZXvWkjZ2rVr165dy3ctpOzXX3+dNm2aCf8gTpbphJaWlsbGRr5rIWUcx9XW1vJdC4mrq6szGDq/zDPqMJ1Ox3cVJM5gMNTV1ZnwD2KLECGEkKxhIkQIISRrhNIubgspGW+99db27ds9PT3vWbKsrKy+vr5v374WqJU86fX68+fPDxs2jO+KSNn58+f79u1rb2/Pd0Ukq7CwEAB8fHz4rohkNTQ0ZGdnDxo0qCOFv/jiC39//7uXwUQIer0+JSWF71oghBAyvVGjRjk6Ot69DCZChBBCsoZjhAghhGQNEyFCCCFZw0SIEEJI1jARIoQQkjVcWaajdDrd/v37Gxsb4+Li3N3d+a6ORNTU1Bw5cqSiomLQoEGhoaHGg9XV1adPn24tExISolareaqg6F26dKm4uNj4mmGYCRMmtP5TWlpaVlbWiBEjQkJCeKqdFGRlZeXm5t56ZOLEiYSQc+fOlZeXG49YW1uHh4fzULn/b+/uQ5r6/jiAn9nUZkRilm238g9RW0Hmw8qlwcxCwa/Lh5UhaJr+EVYkUeE/BQlCUP1R2SOpmE6I0BSdukVFsqUlPqCFZfmAejVlRqgrte7u948Ll5Av/GbF77bd9+uvs8Pgvrkczmf34Zw5ueHh4YGBgdDQUB8fH75zdnbWaDT++PEjPj7e29ub7+/p6eno6AgKCoqKilrugfDWqEPm5uZ2797t5+e3YcMGo9HY0tKyZcsWoUM5PZqmlUqlWq1WKBRGo1Gr1d65c4cQYrFY4uLi1Go197WCgoLY2FhBkzqxo0ePms1mf39/QoiHh4fBYOD6T548aTQa9+7dW1tbW1RUlJOTI2hMJ3b//v2HDx9y7ZGRkZmZmfHxcYlEkpyc3N/fr1AoCCE+Pj78d8BBFEXZbLavX782Njbu27eP67RarZGRkVu3bl25cqXFYmlra9u0aRMh5Pbt24WFhQcOHHj69KlWq7169eryDsaCA4qLi6OjoxmGYVn29OnTmZmZQidyBTabbXx8nGsPDw+7ubn19/ezLGs2m5VKpaDRXEd2dvaVK1eWdA4ODspksomJCZZlnz175ufnt7CwIEQ6V5OcnHzu3DmunZSUVFJSImwep/bx40e73b5x48YnT57wnRcvXtRqtVw7KysrPz+fZdn5+XlfX9+WlhaWZcfGxmQy2ejo6LKOhWeEDmloaEhNTXVzcyOE6HS6hoYGoRO5Ai8vL7lczrXXrVsnlUoXFxe5jwsLCyaT6dWrV/Pz88IFdBGjo6NNTU0fPnzgexobG9VqNXfDWaPRMAzT3t4uXEAXMTU1ZTAYsrKy+J7BwcHm5uahoSHhQjmxgIAAiUSypJObirk2PxW3tbVJpdLo6GhCCEVRKpWqsbFxWcdCIXQITdMURXFtiqI+f/6MDeb/rEuXLkVERCiVSu6jVCotLi7OyclRKpW9vb3CZnNq7u7uHR0dN2/e3LVrV3p6ut1uJ4TQNM1vACaRSORyOU3TgsZ0BQ8ePFCpVPwY9vT0NJvN169fDw0NPXbsmLDZXMaSqZgbt1wnXzUpiuKfizsIL8s4hGEY7nKQELJixQpCCP7I5g+qqqoqLS198eIFd5IjIyPfv39PCGFZNj8//8SJE9gD75fdunWLG7FWqzUsLEyv12dkZDAM8/NvbalUivH8+8rKyn7+v1K9Xs+deZqmd+zY8c8///zZv9ATpyVTMfe46vfHM64IHSKXy6emprj25OTk6tWr/+fmdeCg6urqM2fOGI1GfmNcbvoghEgkksOHD3d1dQmXzunxJ9PX13f//v3cyfx5PBNCJicnuXc64Je1traOjIzodDq+hz/zFEXt2bOns7NToGguZclULJfLuVsaP4/nT58+8c9cHIRC6BCNRmM0Grm2yWTSaDSCxnEdTU1Nx48fr6+v37Zt239+obOzk3srDH6T3W7v7u7mTqZGo7FYLDabjRDS29s7OzsbHh4udEDnVlpampaW9p+/jxcXF9+8ebN58+b/fyrXExMTYzKZuDY/FatUqunp6b6+PkLI3Nxca2vrcqdoLJ9wyOTkZEhIiE6nk8vlly9fNhgMv7BUBZYYHh4ODg6OiIjg17Hl5eWFhIScP3/earUGBAQMDAxUVlaWl5enpKQIG9V5RUVFxcbGenl5NTc3j4+Pt7e3r1mzhhCSmJhos9m0Wu29e/d0Ol1hYaHQSZ2YzWZTKBTNzc38mp+ZmZmEhITY2FhPT8+6urqFhYWXL1/KZDJhczqXoqKikZERvV6v0Wgoirpw4QJFUUNDQ+Hh4bm5uTKZ7Nq1ay0tLdu3byeEFBQUNDQ05Obm1tbWrl27trq6elnHQiF0FE3TFRUV3759S0lJCQkJETqOK5ienl4yXuPi4vz9/d+9e2c0GicmJtavX5+QkBAcHCxUQhdQX1/f1dX1/fv3gICAtLQ0fi5eXFwsLy8fHBzcuXNncnKysCGdHU3TJpMpOzub72EYpq6urqenx263BwcHHzx40MPDQ8CEzqimpsZqtfIfdTodt6x+aGhIr9czDJOWlsav52ZZtqampr29PTAwMDMz093dfVnHQiEEAABRwzNCAAAQNRRCAAAQNRRCAAAQNRRCAAAQNRRCAAAQNRRCAAAQNRRCAAAQNRRCAAAQNRRCAAAQNRRCAAAQNRRCANfX19cnl8vPnj3L9zx+/NjHx6eyslLAVAB/Cew1CiAKN27cOHXqVE1NTVJS0ujoaGhoaExMzKNHj4TOBSA8FEIAsUhNTX3+/Pnr16+PHDkyNjbW1dXFbecPIHIohABi8eXLl7CwsOnp6fn5ebPZrFKphE4E8FfAM0IAsfD29s7IyJiZmUlMTEQVBODhihBALLq7u9VqdWBg4Nu3bw0GQ3x8vNCJAP4KKIQAojA3N6dSqVatWmWxWA4dOtTa2trd3a1QKITOBSA83BoFEIW8vLyxsbGqqipPT8+ysjKZTJaens4wjNC5AISHQgjg+kpKSioqKu7evRsUFEQI4VYQms3moqIioaMBCA+3RgEAQNRwRQgAAKKGQggAAKKGQggAAKKGQggAAKKGQggAAKKGQggAAKKGQggAAKL2L1ABDHwSF48tAAAAAElFTkSuQmCC", "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" ], "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" ] }, "metadata": {}, "execution_count": 7 } ], "cell_type": "code", "source": [ "using LinearAlgebra\n", "\n", "# Define the 1D lattice [0, 100]\n", "lattice = diagm([100., 0, 0])\n", "\n", "# Place them at 20 and 80 in *fractional coordinates*,\n", "# that is 0.2 and 0.8, since the lattice is 100 wide.\n", "nucleus = ElementGaussian()\n", "atoms = [nucleus => [[0.2, 0, 0], [0.8, 0, 0]]]\n", "\n", "# Assemble the model, discretize and build the Hamiltonian\n", "model = Model(lattice; atoms=atoms, terms=[Kinetic(), AtomicLocal()])\n", "basis = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1));\n", "ham = Hamiltonian(basis)\n", "\n", "# Extract the total potential term of the Hamiltonian and plot it\n", "potential = DFTK.total_local_potential(ham)[:, 1, 1]\n", "rvecs = collect(r_vectors_cart(basis))[:, 1, 1] # slice along the x axis\n", "x = [r[1] for r in rvecs] # only keep the x coordinate\n", "plot(x, potential, label=\"\", xlabel=\"x\", ylabel=\"V(x)\")" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "This potential is the sum of two \"atomic\" potentials (the two \"Gaussian\" elements).\n", "Due to the periodic setting we are considering interactions naturally also occur\n", "across the unit cell boundary (i.e. wrapping from `100` over to `0`).\n", "The required periodization of the atomic potential is automatically taken care,\n", "such that the potential is smooth across the cell boundary at `100`/`0`.\n", "\n", "With this setup, let's look at the bands:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computing bands along kpath:\n", " -½ -> ½\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=12}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVxV1f7/8XWYQeAwDwKCCQiIiJkDSjigWX4th3Li5nDtZtat1LzaozLz17dbWnYNuzaYFdhAikOW1VeQq4Lj11kQRFRUQBAOo6AM57B/f+zv48SF43BTzkH26/nX3uusdfZn89DH+7HOXntvlSRJAgAApTIzdQEAAJgSQQgAUDSCEACgaAQhAEDRCEIAgKIRhAAARSMIAQCKRhACABSNIAQAKBpBCABQNAUFoVarPXHihKmrAAB0LCrlPGtUo9GEhoaWlZWZuhAAQAeioBkhAABtEYQAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaMoKwqoq31deWWrqKgAAHYiybqj38jpjZvanxsZLpq4FANBRKGtGKMQrkyaNNXUNAIAORFkzQh6xBgBoRVkzwubmZlOXAADoWJQVhFVVg0NDh5u6CgBAB6KsIBTi6fz8i6auAQDQgSgrCK2sFm/e/E9TVwEA6EBYLAMAUDRlzQgBAGiFIAQAKBpBCABQNIIQAKBoBCEAQNEIQgCAopkmCFesWNGjR48ePXosX77cYIddu3YNGDDA29s7Li6usrJSbmxqapo/f76fn1+vXr2+/fZbfee6urpXXnnlgQce6Nat25w5c4xxAgCAzsLC+IfcunXrmjVrUlJSVCrV6NGjg4KCnnzyyZYdqqqqJkyY8Nlnn40aNeqvf/3rvHnz1q9fL4T48MMP9+/ff+jQofz8/LFjx0ZERERERAgh4uLihBDbt29Xq9WZmZnGPyMAwP3LBDfUP/bYYyNGjFi0aJEQYtWqVf/zP/+zY8eOlh0+++yz7777LiMjQwhx/vz5Xr16Xb16Va1W9+jRY9WqVU888YQQYu7cuVZWVqtXrz58+PAjjzxy6dIlR0fHWx+XG+oBAG2Z4KfRnJycvn37ytuRkZE5OTmtOmRnZ+s79OjRw8rK6vz58w0NDfn5+W0HHj16tF+/fv/4xz8GDhw4ceLEU6dOGes8AACdgQl+GtVoNPrZm1qt1mg0rTqUl5cHBATod52cnDQajUajkSSp5UB5bldUVJSenj548OANGzZs3rx5xIgRZ8+edXFxMXjouro6Z2dn/e6ePXu6d+9+784MANCx2NnZmZub37qPCYLQxcXl2rVr8nZNTY2rq2urDs7OzrW1tfrd6upqV1dXOduuXbumVqvlgW5ubnJnJyenZcuWmZmZLVy4cO3atbt27Wp10VHPzs4uLy9Pv+vo6HjbPxAAoHMzwU+jgYGB2dnZ8nZ2dnaPHj1u0aGoqOj69esBAQG2trZdu3bVt+fk5MgDg4KCbGxszMz+70Ts7OwaGhpudmiVSuXcAikIADBBEM6cOfOTTz4pLy+vrKz85JNPZs6cKbe//PLL8mW/adOmHTx4MCMjQ6fTvffee2PHjpVnjbNmzVq5cmV9ff2FCxe+//57eeCjjz7a3Ny8efNmIcS//vWvc+fOPfzww8Y/KQDA/UoyOp1ON2/ePLVarVarX375ZZ1OJ7eHh4fv3btX3k5OTu7atau9vf2wYcMKCwvlxtra2smTJzs4OLi6un7wwQf6Lzxw4ECvXr2cnJx69uy5ffv2mx23rKzMzc2t3U4LAHBf6tDvI5QkSaVS3UnjneD2CQBAWx36EWsGA++PpSAAAAZ16CAEAKC9EYQAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaAQhAEDRCEIAgKIRhAAARSMIAQCKRhACABSNIAQAKBpBCABQNIIQAKBoBCEAQNEIQgCAohGEAABFIwgBAIpGEAIAFI0gBAAoGkEIAFA0ghAAoGgEIQBA0QhCAICiEYQAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaAQhAEDRCEIAgKIRhAAARSMIAQCKRhACABSNIAQAKBpBCABQNIIQAKBoBCEAQNEIQgCAohGEAABFIwgBAIpGEAIAFI0gBAAoGkEIAFA0ghAAoGgEIQBA0QhCAICiEYQAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaBYmOWpJSckXX3yh0Wgef/zxkSNHtu2g1Wq/+uqrzMzMXr16PfPMM5aWlnJ7bm5uYmJiY2NjXFzcgw8+KDf+8MMPp06dqq2tDQ0NnT59ur29vfHOBABwnzPBjLC2tnbQoEEFBQU9e/Z8+umnN2zY0LbPnDlzEhISIiIikpKSZs2aJTdeuHBh0KBB5ubmnp6ew4cPP3z4sBBCkqTk5GRnZ+eePXtu2rQpJiZGq9Ua83QAAPc1lSRJRj7k559/npiYuH//fiHE999//8EHHxw/frxlh8LCwsDAwIsXL3p5eVVUVPj4+GRnZ3fv3v1vf/tbeXn5119/LYR46623cnNzf/jhh5YD6+vrHRwcTp48GRYW1va4Go0mNDS0rKysPU8OAHCfMcGMMCMjY9SoUfJ2bGzsiRMnampqWnY4cOBASEiIl5eXEMLFxSUyMnLv3r1CiPT0dP3vqCNHjkxPT2/7zWq12s/Pr93PAQDQWZjgGmFJScnAgQPlbTc3NzMzsytXrjg6Ouo7FBcXu7u763c9PDyKi4tbtbu7u5eWlmq1WgsLCyHEo48+euDAAUmStmzZ4uDgcLND19fXP/PMM/rdN99808PD456eHACgA7GxsTEzu82UzwRBaGlpqb+Mp9PpJEmytrZu2cHKyqrldb6mpiYrK6tWA5uamszNzc3NzeXd5OTkurq6H3/8cerUqSdPnvTx8TF4aHNz8/79++t3HR0dWx0aANCZqFSq2/YxQRB27dq1qKhI3i4qKlKpVN7e3q06FBYW6ncLCwvlYPPx8Wk5sGvXrvozdHBwcHBwmDt37rp161JTU/Xra1qxtLScO3fuvT4hAMB9zATXCMeNG7dt27b6+nohxIYNGx555BEbGxshxJEjR/Lz84UQw4cPLy0tPXbsmBAiOzv7woUL8jXFcePGbdy4UV7ds2HDhnHjxgkhbty40dzcLH9zaWlpfn5+QECA8U8KAHCfMsGqUZ1ON2bMmLKysuDg4LS0tF9//VX+uXLYsGGjR49+7bXXhBDx8fHLly9/5JFHdu7cOX/+/EWLFgkhqqqqHn74YTc3N3lp6N69e/38/P71r3/Nnj1bvqcwPT39ySef/Pzzzw0el1WjAIC2TBCEQojm5ub09PTKysro6Gj9+pe8vDxHR0dPT09598yZM5mZmWFhYb169dIPbGho2LVrl1arHTZsmP7G+XPnzmVnZ6tUqrCwsB49etzsoAQhAKAt0wShSRCEAIC2eNYoAEDRCEIAgKIRhAAARbvVfYR1dXUajcbCwsLNzY0bzwEAnZKBGeH+/fufe+65oKAge3v7gIAAX1/fLl269OvX77XXXsvLyzN+iQAAtJ9/WzWalpa2ePHiY8eOBQQEDB48ODg42MXFRavVlpeXnzp16sCBA+Xl5WPHjn3//fdDQkJMWPQfw6pRAEBbv/80umXLlunTp//lL39Zt25d375923bV6XRpaWnr1q3r06fP2bNn/f39jVgnAADt4vcZYV5enoODg/zyo1vLysry8fFxdnZu59ruMWaEAIC2fp8RBgUF3eGY8PDw9ikGAABju/3bJyRJysrKunr1qryrfzUuAACdwO2DcPz48RqNRv+GP4IQANCZ3CYIq6urtVrtvn37jFMNAABGdpsnyzg6Ojo5ORmnFAAAjO9WM8L4+Pj6+vqamppp06ZFRkbKja+++qpRCgMAwBhuf40wOjraCHUAAGASvI8QAKBo/3aNcP369UePHtXvZmZmVldX63ePHTs2b94845UGAED7+7cgfPPNN1NSUuTt5ubmiIiI3377Tf/pmTNnVq9ebdTqAABoZ7yPEACgaAQhAEDRCEIAgKIRhAAARWt9H2FKSkptba1+d8OGDZmZmfL26dOnjVcXAABG8W/3Efr7+1++fPnWA+7f+w65jxAA0Na/zQiPHj2q0+lMVQoAAMb3b0Ho5uZmqjoAADAJFssAABTt9yBMTk6eN29ecXHxrQecOnXqySefzMvLa+fCAAAwht+DsF+/fqdOnfL393/88ccTEhJycnKam5vljxobGw8fPvzRRx9FRUX17dvX3t7e29vbRAUDAHAvtX77xC+//PLxxx+npqbKKejk5KTVauUbKmxsbCZNmrRgwYK+ffuapti7w6pRAEBbhl/DVFxcvHv37szMzKtXr1pbW3t4ePTv3//hhx92dHQ0fon3CkEIAGiL9xECABSNVaMAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaAQhAEDRCEIAgKIRhAAARSMIAQCKRhACABTNcBDu27dv//798vaNGzcWL148dOjQv/3tbzdu3DBibQAAtDvDQTh9+vQjR47I28uWLVu5cmVzc/MXX3wxZ84cI9YGAEC7MxCEtbW1+fn50dHRQgidTvf111/Pnz8/IyMjOTn5hx9+qK6uNnqRAAC0FwNBWFNTI4RwdXUVQhw/frysrGzSpElCiJiYGK1We/HiReNWCABAOzIQhO7u7mZmZufOnRNCbNq0ydHR8aGHHhJCyO+pNzc3N3KJAAC0H4u2TZaWlo8++ujcuXMnTZr02WefTZw40dLSUgiRmZlpZmbm5+dn9CIBAGgvhhfLrF27tmfPngkJCdHR0cuXL5cbExIS+vTpo1arjVgeAADtSyVJkqlrMBKNRhMaGlpWVmbqQgAAHcitbqjXarXnz5/X31AIAEDnYzgIdTrdm2++6eTkFBgYOHnyZLlx3rx5zz33nBFrAwCg3RkOwrfeeuuDDz546aWX9BcIhRCjRo36/vvvGxoajFUbAADtzkAQNjU1rV69evny5e+9996gQYP07ZGRkbW1tYWFhUYsDwCA9mUgCMvKyq5du/bII4+0apfXi1ZUVNz9UZubm48ePXro0CGtVnuzPhcvXty9e3dVVVXLxsbGxv379584caLVGp+ioqLdu3ezEAYA8J8ycB+hg4ODmZlZcXFxWFhYy/bMzEwhhLe3910esq6ubvTo0TU1NTY2NvX19bt27ZKfYtPS0qVL165d26dPn2PHjm3YsGHEiBFCiKKiomHDhnl4eFRUVHTr1u3nn3+2srISQsTHx7/zzjv9+vU7cuTI2rVrJ06ceJcVAgAURDJk6NCh0dHRNTU1e/bs8fHxkSSpsrIyOjo6MjLSYP//SHx8/JAhQ7RarSRJ48ePf+ONN1p1uHDhQpcuXfLz8yVJ+uqrr3r37i23//Wvf50xY4YkSQ0NDX369ElMTJQkSaPR2NnZnTx5UpKkn3/+2dfXt6mpyeBxy8rK3Nzc7r5+AEBnYnixTHx8fGZmZkhIyH//93/X1NTMnDkzJCTkyJEj//znP+8+ejdv3jx9+nT5UW0zZ85MTk5u1WHLli0xMTEBAQFCiKlTp549ezY3N1cIsWnTplmzZgkhrKys4uLiNm3aJITYvn17WFhYRESEEGLMmDENDQ0HDx68+yIBAAph4KdRIUSfPn2OHDny9ttvp6amXrt2bdu2bcOGDVu2bFlkZOTdH/Ly5cv+/v7ydkBAQEFBQasOBQUFcgoKIWxtbT09PQsKCrp3715aWtpyYFJSUqvOZmZm3bp1a/uFelqtNjU1Vb8bFRVlZ2d392cEAOiYzMxu//55w0EohAgMDFy/fv09ref/1NfXy9f2hBDW1taNjY06na7ls7zr6+ttbGz0u9bW1tevX79x44YkSS0HXr9+vdW3tWw3qK7OcvHixc7OzvLup59+6uPjc+/ODADQsdjZ2d02C28ahO3Hy8urvLxc3tZoNO7u7q3eaOHl5XX+/Hn9bnl5ube3t1qttrW1LS8v9/X1lQfKy3Y8PT317xBu2W6QTvfhmTNv3bhx/N6eEQDg/mU4CBctWiS/lbCtzz///C4POWDAgIyMDPkdh+np6QMHDmzbISEhQZIklUqVlZWl1WpDQ0P1A/v06dNy4MCBA998882GhgZra+vi4uJLly49+OCDNz/4j76+Xe+yfgBAZ2I4CNPT0zUajX63pqZGXpzp5eV194d88cUXhwwZEhISYmdnt3Llyh9//FFuDwwMXL169ZgxY0aPHm1vb//ss8+OHTv23XffnTNnjr29vRBi4cKFs2fPdnZ2vnr16rZt244fPy6EGDBgQERExIwZM+Li4lavXj1t2rRbFOnktDsvr/zuTwEA0Gnc6dsnsrOz4+LiFi9eHBcXd/dHPXjw4Lp167Ra7fTp02NjY+XG9957b/z48fLkr7S0dOXKlZcvX46Ojp47d66Fxf8F9vbt2zdu3Ghra/v888/rV+5UVVV9+OGHZ8+efeihh15++WVra2uDB+XtEwCAtv6D1zClp6ePHTu2pKTkPl1pSRACANq6/bpSvZ49e167dk2+pQ8AgM7hPwjCn3/+WdyLR6wBANBx3NGqUa1We+7cuYyMjFGjRt2T9TIAAHQQd7Rq1MLCwtfXd/ny5S+++KKxCgMAwBgMB+GhQ4eMXAcAACbxH1wjBACg8/l9RlhWVnbu3LnbDoiKimrPegAAMKrfg3D79u2zZ8++7YA7v+8QAICO7/cgHDNmzJ49e0xYCgAAxvd7EHp6enp6epqwFAAAjI/FMgAARbvp+wgPHTq0cePGCxcu1NbWtmxv+YZ3AADud4aDMCkp6emnn/bz82tsbLSxsVGr1dnZ2VZWVoMGDTJyfQAAtCvDP40uWbJk0qRJ58+fHz16dFxc3PHjx8+cORMcHBwTE2Pk+gAAaFcGgrCuri4/P3/BggXm5uZCiMbGRiFE9+7d165d+/e///1mb64HAOB+ZCAIb9y4IUmSWq0WQri6uuofOhoWFtbQ0HAnN90DAHC/MBCErq6u9vb2ly5dEkIEBwfv3LlTXi+TlpYmhHBzczNyiQAAtB8DQahSqWJjY7du3SqEmDZt2vXr13v16jVq1Kinnnpq6NChfn5+Ri/ynqmpcVuzZq2pqwAAdCAqg49Mq6iouH79uq+vrxAiKysrPj4+Pz+/b9++b7zxhpOTk9GLvDc0Go2X12Fz8+cbGi6auhYAQEdhOAg7JY1G4+X11ODBrunpm01dCwCgozB8+8Trr79+4MABI5diBM7Op0lBAEBLhoMwOTl58ODBISEhy5cvLyoqMnJNAAAYjeEgPH369E8//RQREbF06dJu3bpFR0evXbu2rq7OyMUBANDebnON8OrVq99//31iYuLJkyednJymTJny2WefGa24e0uj0YSGhpaVlZm6EABAB3Kni2UyMjKmT59+6dKl+3dxDUEIAGjrpm+fkDU3N6elpSUmJm7duvX69etDhgwxTlkAABjHTYMwNzc3KSlp/fr1+fn5Pj4+L7300uzZs4ODg41ZHAAA7c1wEI4ePTolJcXW1nbixImff/55bGysmRmv8AUAdEKGg7BLly5ffPHF5MmTHR0djVwQAADGZDgIt2zZYuQ6AAAwiVstlikvLy8qKmpqamrZ2K9fv3YuCQAA4zEchNnZ2XPnzs3IyGj70f17+wQAAG0ZDsKpU6eWlpauWrWqZ8+elpaWRq4JAACjMRCENTU1WVlZmzdvnjBhgvELAgDAmG56U8R9/QJeAADukIEgdHR0fOyxx37++WfjVwMAgJEZvkY4f/78Z555pqqq6rHHHnN3d2/5EatGAQCdieGHbnt5eV29etXggPt31SgP3QYAtGV4Rrhx48bGxkYjlwIAgPEZDsKYmBgj1wEAgEnc6sky+fn5p0+frq6u/tOf/iSEqKmpsbS0tLW1NVZtAAC0O8O3T9TV1T311FMPPPDA448//uqrr8qNCxYs4M5CAEAnYzgIX3zxxd27d69fvz4pKUnf+PTTT+/atev69evGqg0AgHZnIAjr6+uTkpI++uij6dOne3t769vDwsIaGxsLCgqMWB4AAO3LQBCWl5c3NDQ89NBDrdqtra2FENeuXTNGXQAAGIWBIHRxcbG0tDxz5kyr9oMHD6pUKn9/f6MUBgCAMRgIQltb27Fjx7722mv5+fkqlUpuPHv27IIFC4YPH97qQTMAANzXDN8+8fHHH8fExISEhAQEBFRUVAwcOPD48eOurq4//fSTkesDAKBdGV416uPjc+zYsWXLlnXt2tXPz8/c3HzhwoUnTpwICgoycn0AALQrw88a7ZR41igAoK2bvo8QAAAlMHyNcPLkyZWVlW3bHRwcunfvPmHChOjo6HYuDAAAY7jpjPDw4cO7d+8uLS2VJKmoqCgtLS0rK0uj0XzzzTcxMTHvv/++MasEAKCdGA7CqKiogICA3NzckydP7ty5Mzs7+/jx49bW1vPmzSsoKJgzZ86bb755sxcWAgBwHzGwWKapqcnd3X3btm1Dhw5t2Z6QkPDhhx9mZmbeuHHD2dk5KSnp/noGN4tlAABtGZgRajSa6upqNze3Vu3u7u55eXlCCFtbWz8/P561BgDoBAw/Yq1Lly4JCQktGyVJSkhI0D9fTaPRtE3KO5ecnBwWFubl5fWXv/zlxo0bbTucPXs2NjbWzc1t+PDh2dnZ+vZVq1Y98MADfn5+S5YskeeykiRNmDDhgQce8PT0HDp0aHp6+h+uCgCgRJIhb731lkqleuyxx9asWbNx48aPPvpo8ODBQoivv/5akqTdu3erVKqSkhKDY2/r3Llz9vb2O3fuLCsrGzly5OLFi9v26dOnz9KlS6urq999993Q0NDm5mZJknbs2OHl5ZWZmXnx4sWQkJAvv/xS7rxx48ZLly6VlpbGx8fb29tXVVUZPG5ZWZmbm9sfqxkA0FkZDsLm5uZVq1a1fAdTYGDgd999J39aWVmZn5//hw/5xhtvTJkyRd5OT093c3PT6XQtOxw6dEitVjc0NEiSpNVq3d3d09PTJUmaNGnS0qVL5T7r1q2Liopq9c1ardbCwiIrK8vgcQlCAEBbhleNqlSq+fPnX7lypaysLDMzs7KyMi8vLy4uTv7UyckpICDgD89Bc3JyIiMj5e3IyEiNRtNqAUtOTk5YWJiVlZUQwtzcvHfv3jk5OUKI7Oxs/cA+ffrIjbLDhw//9ttvzz///IgRI0JDQ/9wbQAApTF8Q72em5vbH7gWqNFoduzY0bZ9zJgxzs7O5eXljo6OcouDg4O5ublGo/H09NR3q6iocHBw0O+q1WqNRiOEaDlQrVZXVVU1NTVZWloKIb799tsTJ07k5eW99dZbZmY3vTmytrZW/z4NIcSxY8cCAwP/07MDANwv7OzszM3Nb93n9yA8f/78vn37HnzwwfDw8OTkZINrWIQQM2bMuO2BKysrU1JS2rY//PDDzs7Ozs7O+hWntbW1Op3O1dW1ZTdnZ+fa2lr9bnV1tdzBxcVFP7CmpkatVsspKISIj48XQhQWFvbu3TsiIiIqKspgYfb29jc7LwCAMv0ehOnp6bNnz3733XfDw8NffPHF0tJSgwPuJAiDgoISExNv8enp06fl7dOnT6vV6lbvOAwKCsrNzZUv+EmSlJOTI8/b5IHjx4+XB7adzPn6+gYGBubl5d0sCAEAaOX3IJw6deqjjz4q//aYmZmp0+na6ZCzZs2Kioo6efJkSEjIihUrpk+fLs9b33///fDw8DFjxgwePNjNzW3NmjUvvfTSunXrrK2thw0bJg+cP3/+rFmzunTpEh8f/+c//1kIkZ+ff+XKlYceekgIsWnTppycHHmBKwAAd+L3ILS1tbW1tZW3PTw82u+QYWFhK1eufPTRR+vq6kaNGvXOO+/I7ZmZmWq1WgihUql++OGH2bNnv/baa0FBQcnJyXJSjh8//ujRoxERETqdLi4ubu7cuUKI69evz58/X55i9u7dOzk5mct+AIA7d9P3EUqStH///tOnT2u12hdeeEEIUVRUZG1tfTf30ZsWj1gDALRlOAg1Gs24ceP2798vhPDx8SksLBRCvPTSS6dOndqzZ4+xa7xHCEIAQFuG7zR47rnnCgoKUlNTWy7+nDZt2r59+3jEKACgMzEQhNevX//pp5/i4+NHjhwp39UuCw4O1ul0BQUFRiwPAID2ZSAIq6qqtFptSEiIwQHchwcA6EwMBKGrq6utre2JEydate/evdvMzOyBBx4wSmEAABiDgUesWVtbP/XUU6+++mqPHj30DyTbu3fvggULHn/8cWdnZ+NWCABAOzK8arSiomL06NFHjhyRn2rm4eFRVFQUHBy8a9eurl27Gr/Ke4JVowCAtgw/dNvFxWXfvn1JSUkpKSmlpaWOjo5Dhw595plnunTpYuT6AABoVze9ob7zYUYIAGjrpm8sAgBACQhCAICiEYQAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaAQhAEDRCEIAgKIRhAAARSMIAQCKRhACABSNIAQAKBpBCABQNIIQAKBoBCEAQNEIQgCAohGEAABFIwgBAIpGEAIAFI0gBAAoGkEIAFA0ghAAoGgEIQBA0QhCAICiEYQAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaAQhAEDRCEIAgKIRhAAARSMIAQCKRhACABSNIAQAKBpBCABQNIIQAKBoBCEAQNEIQgCAohGEAABFIwgBAIpGEAIAFI0gBAAoGkEIAFA0C5Mc9dSpU+vXrxdCTJ8+vU+fPm07VFZWfvrpp0VFRcOGDZs0aZK+fffu3Vu3brW3t58zZ46/v7++PT09fdu2bVqtNjo6umV/AABuzQQzwqysrOjoaBcXFzc3t5iYmMzMzFYddDrdsGHDsrKy+vbt+/rrr69cuVJu/+WXXyZOnBgcHNzQ0DBgwACNRiO3r1q1aurUqa6uroGBgQcPHjTqyQAA7nMqSZKMfMhnn33Wzs4uPj5eCLFw4cKqqqovv/yyZYeff/55wYIFZ8+eNTMzS09PnzJlyuXLly0tLWNiYqZOnfrCCy8IIcaOHRsTE7N48eLCwsLg4OCjR4+Ghobe+rgajSY0NLSsrClicjYAAA3ISURBVKz9Tg0AcN8xwYxw7969I0eOlLdjY2MzMjLadhg2bJiZmZkQYsiQIZWVlefOndPpdAcOHGg7cNeuXZGRkXV1de+999769esbGxuNeCoAgPteu1wj1Ol0Wq22VaNKpbKyshJCFBcXu7m5yY0eHh7FxcWtehYXF/v4+Mjb5ubmrq6uV65ccXFx0Wq1bQfm5+cXFBQsXLhwwoQJ33zzzdq1a/fs2WNubm6wsPr6+ieffFK/+8EHH3h5ed3t2QIAOiobGxt5WnUL7RKEn3/++aJFi1ofycKiurpaCGFlZaWPycbGRmtr61Y9W3YQQjQ1NVlbW8sh2naghYVFRUVFVlaWWq2eO3dut27ddu3apZ84tmJubj558mT9rouLS9ujAwA6DZVKdds+7RKEL7zwgnwlzyAfH5/CwkJ5u7CwUD/5a9khNzdX3q6rq6usrPTx8XF2drazsyssLPTw8Gg50NfX18vLS61WCyFsbGy6d+9eUFBws0NbWlpOmTLl7k4OANCpmOAa4bhx43744Qd5e8OGDePGjZO3U1JSysvL5Q6pqakVFRVCiE2bNvXq1at79+5CiCeeeEIe2NTUtGXLFnngmDFjysrKLl++LIQoLy/Pycnp3bu38U8KAHCfMsGqUY1GEx0d7ePjY2Zmdvny5b1797q7uwsh7O3tt27dOmrUKCHE7Nmz9+7d++CDD+7cufO7774bPXq0ECInJ2f48OFRUVGFhYV2dnYpKSnyD5sffPDBxx9/PGLEiIyMjLFjx8rrUQ0el1WjAIBWTBCEQoj6+vo9e/YIIYYOHWpjYyM3nj592t/f397eXt49fPhwUVHRwIEDvb299QOrq6szMjIcHByio6Nbrog5e/ZsdnZ2cHBwWFjYzQ5KEAIA2jJNEJoEQQgAaItnjQIAFI0gBAAoGkEIAFA0ghAAoGgEIQBA0QhCAICiEYQAAEUjCAEAikYQAgAUjSAEACgaQQgAUDSCEACgaAQhAEDRCEIAgKIRhAAARSMIAQCKRhACABSNIAQAKJqyglCr1dbX15u6CgBAB6KsILx27U9+fv1NXQUAoANRVhBKUg9mhACAlixMXYBRdeny/44c+V9TVwEA6ECUNSO0trYMDg42dRUAgA5EWUEIAEArBCEAQNEIQgCAohGEAABFIwgBAIpGEAIAFI0gBAAoGkEIAFA0ghAAoGgEIQBA0ZQVhE1NTbW1taauAgDQgSgrCGtrZwcEDDJ1FQCADkRZQShJXo2NjaauAgDQgSjtNUx/P3nyhKmrAAB0IMqaEVpbW3Xv3t3UVQAAOhBlBSEAAK0QhAAARVNQEFZUVFRWVpq6CgBAx6KgIBRCSJJk6hIAAB2LsoIQAIBWCEIAgKIp6D5Ca2trGxubUaNGmboQAICRfPbZZz169Lh1HwUFob+//7FjxwoKCkxdCADASDw8PG7bR8X6EQCAknGNEACgaAQhAEDRCEIAgKIRhAAARVNEENbX1//v//5vSkpKy8bi4uJWLQAABVJEEJ47dy49PX3RokUtG+fPn29ubm6qkgAAHYQigjA8PHzBggUtW3bu3GlmZhYbG2uqkgAAHYSCbqjXa2xsfP3117du3WrqQgAApqeIGWFDQ8OFCxcaGxsvXLgghPjHP/4xefLks2fPLlu27IsvvmhoaDB1gQAAk+m0T5bZtWtXRUWFECIqKkoIsXr1arn9hRdeeOqpp/bt27d7925/f/+1a9c6OjouXbrUlLUCAEyn0/40mpGRkZ+fL4QICAjo16/f8uXL5fZJkyatWLHC0tJSfvp2dXX1gAEDTFkoAMCkOm0QGpzkpaamWllZDR8+XN5dsWJFt27dJk+ebNzSAAAdSKf9abSthoaGIUOG/PTTT127dhVCvPXWW0ePHp0/f76Xl1d4eLipqwMAmEannRG2lZubu2TJEjkFhRC9e/f28fG5cOECdxMCgJIpaEYIAEBbirh9AgCAmyEIAQCKRhACABSNIAQAKBpBCABQNIIQAKBoBCEAQNEIQgB/RFVV1dq1ay9evGjqQoC7RRAC+CNKSkqee+6548ePm7oQ4G4RhICJabXa0tLS5ubmPza2pKSkrq7u7stobGwsKSm5ceOGwU+bmppu8el/WpVOpyspKeFVoOggCELgHps9e/aIESP0u6NHj3ZxccnOzpZ3ExMT3dzcKisrhRC//vrrgAEDrK2tPT09u3Tp8thjj126dEnutmLFCk9PT7mb3tSpUwcOHChvNzY2Llq0yN3d3dvb29HRMTY2Vn7vdFvPPvtseHh4q6Dt379/XFycvF1ZWTlr1iwnJydvb2+1Wj158mT5XZ6yqqqqZ5991sXFxdvb287OLjw8/NixY4cOHZIrmTFjhouLi4uLy7Zt24QQOp1uyZIl+qpGjhx59uxZ/VcNGTLkueeeW7lypdzhu++++wN/XuDekwDcU5988om5uXlFRYUkSTU1NZaWllZWVvHx8fKn06ZNi4yMlLe//PLLFStW7N27NycnZ+PGjd27d+/Tp09zc7MkSXl5eSqV6tNPP9V/7dWrVy0tLd955x15d9KkSY6Ojh9//HFmZuaOHTsiIiK6d+9+7dq1tvX88ssvQoi0tDR9y8GDB4UQGzdulCSpoaGhf//+Pj4+33777enTpzdv3uzn5xcdHS2X0dDQMGDAAHt7+1WrVp04ceLAgQPvvvtuenp6ZWXlV199JYRYtmxZampqampqSUmJJEmvvPKKSqV6/fXXjx8/vmXLFh8fHx8fH/lPIUlSeHi4h4dHRETEpk2b9u7de+bMmXv8pwf+EIIQuMdyc3OFEFu3bpUkafv27TY2NrNmzXriiSckSWpubvby8lq4cKHBgampqUKIzMxMeTc6OnrQoEH6T1etWmVmZnbp0iVJkvbv3y+EWL9+vf5T+T0qX3/9dduv1Wq1vr6+M2bM0Lc8//zzLi4u9fX1kiQlJCQIIQ4cONCqjH379kmStH79eiFEUlJS26/NyckRQmzZskXfUl5ebmVl1fJA+/btE0KsWLFC3g0PD7e1tS0qKrrJXw4wDQW9hgkwjuDgYF9f37S0tPHjx6elpT388MP/9V//NXv27KampjNnzpSUlMTGxuo7nzlz5scff7xy5UpDQ4N8Ue3cuXPyCzJnzpz57LPPnjlzJiQkRAiRmJgYGxvbrVs3IURKSooQwsHBYefOnfqvcnd3z8rKaluPubl5XFzcmjVr/vnPfzo4ODQ2Nm7cuHHatGnW1tbyVzk5OdXW1uq/qrGxUaVSZWVlDR48ODU11d7efsqUKXdy4jk5OY2NjS3fdD148GA/P7/09PTFixfLLVFRUfpXoQEdBEEI3HsjRoyQcyUtLS0uLm7EiBF1dXVHjx49ePCghYVFdHS03O29995bsmTJoEGDIiIinJ2dLS0thRDV1dXyp5MnT543b94333zz97//PSsr68SJE99++6380dWrV1Uq1ezZs1sdVz+2lZkzZ77//vubN2+eNWvWtm3bysvLZ86cqf+q2traluklhHBycpIvT2o0Gh8fH5VKdSdnLV/gbJVzXbt21Wg0+l1PT887+SrAmAhC4N6LjY1dv379iRMnsrKyRo4c6eLiEhkZmZaWdvDgwUGDBjk4OAghtFrt22+//fLLL69atUoedeLEiTVr1ui/xNHRceLEiYmJiW+//fZXX33l6Og4YcIE+SO1Wq1SqS5duiR/1W2FhYUNGDAgMTFx1qxZiYmJYWFh/fv313+Vp6dnYWGhwYFOTk5Xr169w7Pu0qWLEKKsrKxlY1lZWWBgoH7XzIwFeuhw+EcJ3HsjR44UQrzxxhtqtToyMlIIERsbu2PHjvT0dP3voiUlJfX19f369dOP+vXXX1t9z8yZM4uKilJTU5OSkqZMmWJnZye3Dx06tLm5edOmTXde0syZM/fs2XPw4MGUlJRZs2bp22NiYoqKiuSLjm3FxMRUVVXJv8S2Ym9vL4Sor6/Xt/Tv39/c3FxemyPLysrKz8+Pioq68zoBEzD1RUqgc5Iv7E2cOFHe3bFjh/w/bs+ePXKLTqfz8PAYPHhwfn5+TU1NQkKC/LNhQkKC/kt0Op2/v39AQIAQIiMjQ9/e3Nw8fPhwtVq9du3akpKSmpqa48ePL126ND09/Wb1VFRU2NjYBAQEWFhYXLlyRd9+7dq1oKAgPz+/zZs3l5eXl5eXHzhw4OWXX75w4YIkSXV1dcHBwd7e3lu3bq2qqiopKUlKSjp69KgkSVqt1tXV9ZFHHtm9e/eRI0fkpaEzZsywtrb+4osvKisrjx07FhERoVar9atjwsPD//SnP92rvzBwrxCEQLt44YUXhBCffPKJvFtXV2dtbW1nZyev1ZT99ttv7u7uckD6+/tv3bq1VRBKkrRkyRIhRGBgoHw/g15NTc2f//xn+bKiLDIy8tixY7coadKkSUKIMWPGtGovLCwcO3as/kdLMzOz6Ojo4uJi+dPLly+3XN3j6uq6d+9e+aOffvqpV69e8qIb+WaMurq66dOn678qKCho//79+gMRhOiYVJIkteN8E8At1dfX5+bmWlhYhIaG/oHrZ7W1tbm5uebm5r6+vm5ubndTSXl5+fnz5+3s7Pz8/NRqdatPr1y5UlhY6ODgEBQUZGFxm7UFGo3mwoULjo6OPXv2vMOFNoAJEYQAAEVjsQwAQNEIQgCAohGEAABFIwgBAIpGEAIAFI0gBAAo2v8H5/WWbk6GnHQAAAAASUVORK5CYII=", "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" ], "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" ] }, "metadata": {}, "execution_count": 8 } ], "cell_type": "code", "source": [ "using Unitful\n", "using UnitfulAtomic\n", "\n", "plot_bandstructure(basis; n_bands=6, kline_density=15)" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "The bands are noticeably different.\n", " - The bands no longer overlap, meaning that the spectrum of $H$ is no longer continuous\n", " but has gaps.\n", "\n", " - The two lowest bands are almost flat. This is because they represent\n", " two tightly bound and localized electrons inside the two Gaussians.\n", "\n", " - The higher the bands are in energy, the more free-electron-like they are.\n", " In other words the higher the kinetic energy of the electrons, the less they feel\n", " the effect of the two Gaussian potentials. As it turns out the curvature of the bands,\n", " (the degree to which they are free-electron-like) is highly related to the delocalization\n", " of electrons in these bands: The more curved the more delocalized. In some sense\n", " \"free electrons\" correspond to perfect delocalization." ], "metadata": {} } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.3" }, "kernelspec": { "name": "julia-1.6", "display_name": "Julia 1.6.3", "language": "julia" } }, "nbformat": 4 }