{ "cells": [ { "cell_type": "markdown", "source": [ "# Custom potential" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "We solve the 1D Gross-Pitaevskii equation with a custom potential.\n", "This is similar to Gross-Pitaevskii equation in 1D example\n", "and we show how to define local potentials attached to atoms, which allows for\n", "instance to compute forces.\n", "The custom potential is actually already defined as `ElementGaussian` in DFTK, and could\n", "be used as is." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using DFTK\n", "using LinearAlgebra" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "First, we define a new element which represents a nucleus generating\n", "a Gaussian potential." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "struct CustomPotential <: DFTK.Element\n", " α # Prefactor\n", " L # Width of the Gaussian nucleus\n", "end" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Some default values" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "CustomPotential() = CustomPotential(1.0, 0.5);" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "We extend the two methods providing access to the real and Fourier\n", "representation of the potential to DFTK." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function DFTK.local_potential_real(el::CustomPotential, r::Real)\n", " -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)\n", "end\n", "function DFTK.local_potential_fourier(el::CustomPotential, q::Real)\n", " # = ∫ V(r) exp(-ix⋅q) dx\n", " -el.α * exp(- (q * el.L)^2 / 2)\n", "end" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "We set up the lattice. For a 1D case we supply two zero lattice vectors" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "a = 10\n", "lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "In this example, we want to generate two Gaussian potentials generated by\n", "two \"nuclei\" localized at positions $x_1$ and $x_2$, that are expressed in\n", "$[0,1)$ in fractional coordinates. $|x_1 - x_2|$ should be different from\n", "$0.5$ to break symmetry and get nonzero forces." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "x1 = 0.2\n", "x2 = 0.8\n", "positions = [[x1, 0, 0], [x2, 0, 0]]\n", "gauss = CustomPotential()\n", "atoms = [gauss, gauss];" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "We setup a Gross-Pitaevskii model" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "C = 1.0\n", "α = 2;\n", "n_electrons = 1 # Increase this for fun\n", "terms = [Kinetic(),\n", " AtomicLocal(),\n", " LocalNonlinearity(ρ -> C * ρ^α)]\n", "model = Model(lattice, atoms, positions; n_electrons, terms,\n", " spin_polarization=:spinless); # use \"spinless electrons\"" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "We discretize using a moderate Ecut and run a SCF algorithm to compute forces\n", "afterwards. As there is no ionic charge associated to `gauss` we have to specify\n", "a starting density and we choose to start from a zero density." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n Energy log10(ΔE) log10(Δρ) Diag Δtime\n", "--- --------------- --------- --------- ---- ------\n", " 1 -0.143697315381 -0.42 7.0 \n", " 2 -0.156054236401 -1.91 -1.10 1.0 2.05ms\n", " 3 -0.156771127149 -3.14 -1.57 1.0 1.75ms\n", " 4 -0.157013765528 -3.62 -2.11 1.0 1.63ms\n", " 5 -0.156930211636 + -4.08 -1.91 2.0 1.93ms\n", " 6 -0.157056376920 -3.90 -3.61 1.0 1.46ms\n", " 7 -0.157056404273 -7.56 -4.03 1.0 1.48ms\n", " 8 -0.157056406868 -8.59 -4.80 1.0 1.53ms\n", " 9 -0.157056406914 -10.33 -5.34 1.0 1.74ms\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380293 \n AtomicLocal -0.3163463\n LocalNonlinearity 0.1212606 \n\n total -0.157056406914" }, "metadata": {}, "execution_count": 8 } ], "cell_type": "code", "source": [ "basis = PlaneWaveBasis(model; Ecut=500, kgrid=(1, 1, 1))\n", "ρ = zeros(eltype(basis), basis.fft_size..., 1)\n", "scfres = self_consistent_field(basis; tol=1e-5, ρ)\n", "scfres.energies" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "Computing the forces can then be done as usual:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "2-element Vector{StaticArraysCore.SVector{3, Float64}}:\n [-0.055676487596999225, 0.0, 0.0]\n [0.05568269960062737, 0.0, 0.0]" }, "metadata": {}, "execution_count": 9 } ], "cell_type": "code", "source": [ "compute_forces(scfres)" ], "metadata": {}, "execution_count": 9 }, { "cell_type": "markdown", "source": [ "Extract the converged total local potential" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1" ], "metadata": {}, "execution_count": 10 }, { "cell_type": "markdown", "source": [ "Extract other quantities before plotting them" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=4}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzddXwUd/4/8M/Iusc9xBNCCO7u0BTaAqVQoKVy9FulpVcX7tre1bgKlau7C7TF3QoET5A4cc+6787M5/dH8uMoEAjJZmfl/Xz0cY+EG3Y/fHZnXjMfJTDGCAAAAAhWJN8FAAAAAPgEQQgAACCoQRACAAAIahCEAAAAghoEIQAAgKAGQQgAACCoQRACAAAIahCEAAAAghoEIQAAgKAGQQgAACCo+VwQPv3003a7vYsHMwwDS8TxyO12812EoAb1zy+of355sP59Lgi/+OILnU7XxYPdbjfHcb1aHnAFDoeD7yIENah/fkH988uD9e9zQQgAAAB4EwQhAACAoAZBCAAAIKhBEAIAAAhqEIQAAACCGgQhAACAoAZBCAAAIKhBEAIAAAhqNN8FAL3LzaESIy414kozYjjk4pBKiCLEKE1F9A8hBHAjBAB/OIyKDLjEiBttyOhCFIFIAiXIUZqSyFQTUrg8ewvUdGBqsqPvK7it9dyfTThaSmSqiWQFElJISKIKEzrQjIoNXLkJDwglbuhDzutD9FEQfBcZgGDRYkdrq7lfK7lDLThKSvTTENFSpBIiFiOWQ4dbUZmRqzDjgaHE1FjylmQiTQWnZ++CIAw02+rx22fYP5vxjYnkXRnkNxNIjejyR9oYdKAZ/1LFDf+d7R9CPNyPmhlPwAkHQO851ILfOM1treNmxJN/yyS/n3Sl03NfE95cx41dz6apiPuyyPnJJAXnZ+8gfG3R6ri4uPz8/NjY2K4cbLfbhUIhRVG9XSq/sKkWrzrOWtzo7/3J+cmkrMs3OW4O/XiO+89pzs2hV4dRM+K6eraZzWaFQtHN4oIeg/rn1zXV//E2/PfDbJUZPdSPXJZOKgRdfRc3hzbWcqtPcU129FQuuTSNJCEOEUIe/f5DEAaCYgN++BBbZUEvDSFvSOz+efJbNffEES5Bht4dTaUqr/4qcCHmF9Q/v7pY/20OtDKf3V6Pnx9E3pnR/ae6vU346aOsjUFvjaDGREEYevL7D4Ml/JubQ/88wY1bz0yPIwtvom/q06O7xTmJ5Kmb6Jnx5MjfmdWnONa37pEA8D/fV3D9f3VHiFHJfPpvmT1q2xwXRezNo/+eQy7axf5tP2uCPaA8B4LQj53W42G/Mfkt3Ikb6RX9SI8MAaVJtKIfmT+H3lDDTdrA1FshDAHoDrMbLdrFvniS+20q/dpwSt7lttArIBC6JYU8PY8mEOr/C7OjAU5Pz/BMEH799ddpaWkRERF33nnnFbbV3bhx45AhQ9avX++RNw1yH5dwkzYwD2aTG6bTsTIPt5MkK4jts+hpceSQdcyWOjjZALg2BTo8ZB0jF6Ajc+ih4R4+PZUC9MEY6qOx1G172GePsdBy03MeCMKioqL77rvv888/Lykpqa6ufuGFFy57mNFofPzxx3U6nVar7fmbBjMbgxbvZtec4fbm0cvSe+uZniTQ0wPIHybTd+5jV5+C3Y8B6KpfKrmpG5nnB5EfjqEkvTYwf2oscfwG+nALnryRae706QN0iQcuo59++ukNN9wwevRojUbzzDPPfPLJJ5cdgLNy5coVK1aEh4f3/B2DWY0Fj13PUAQ6NJvOVPd6h/m4KOLgbOqbcu6OvawL0hCAq3nxBPfwIW7zDHpRSq93PEVI0KYZ9IRoYvhvzLE2eDDsPg98VMXFxbm5ue0/5+bmtrS06PX6i47ZuXNnRUXFHXfc0fO3C2YHmvGI35nFqeQX43vxTvMi8TJi3/W03onytjBm6J8HoBMsRn/bz/5ew+XPoQeFeWlUJ0mgVYOoN0aQs7YwP56De9Vu8sDVVKfTnR/DqlQqEUJtbW0hISHnD7BarQ8++OC6deuILkzXbm1tjYuLO/9rXl7ed99919nBQTV9Yl0t+fBR+sMRzNRop8Xi7Xf/fAR65Cg9/g/ul3HucHHHvafVau3KZwp6CdQ/vy6sfzuDlh0UODn0x3hGxjm8fIZODUO/jSdu3ics1blWZDJefW/+dPH7LxaLafoqSeeBIAwJCTGZTO0/t/8QGhp64QGrVq0aOnSo0Wg8duyY1Wqtqqqqrq5OTEy87KuFh4d3fR4hRVFBEoRvnObeOMVtv47KDRHyVYZPJqJVx9mZu8gds6j24TkYY7lczld5ANQ/v87Xv8WNFuxhYmTEZ+MoAdnJUjG9bIQcHZyN87YSDU7BO6OoYJh078HvvweaRtPS0s6cOdP+8+nTp0NDQy98HEQIcRx36tSp5cuXL1++vKqq6vPPP//qq696/r5BAiP05BH2kxLuwGwqN4Tnb/eqQdRdmeT4DWyVGTokAEAIIaMLzdjMpKqIL8dT/K5iHysj9ubRZSZ8y07WyfJZEr/jgZVlTp06NXbs2N27d6enp8+fP79v376vvfYaQuiFF14YMmTIzJkzLzx4+PDh995772233dbZq8HKMhdiMVq+nz2rx+un0yH83GhexrtnudcKuV3XUWHIAiub8AhWluGX2WzmRIppm5jhEcRbI31lHVAXh5bsZrUOvG4q7ZHJiz7Lt1aWycnJee211/Ly8mJiYpRK5fPPP9/+5yUlJc3NzRcdHBsbC6duF7k4tHAnW2vB22b5UAoihO7rSz6eS07eyNbafOTcB4AHFoaYuZkZ4UspiBASkui7iVSqkpiyidE7+S6Nn4C1Rn2UnUHzdjBCkvh+EiXyyX/fmjPcG6eYvdcL4jw9nR90ETwR8sjKoCnrnYMi6HdG+VAKnocReiyf3VaPt8ykIyV8l6Z3+NYTIfA4K4PytjIaEfHTZB9NQYTQA9nkXans1E1sC0zmBUHGwaIbtjHpSuybKYgQIhB6bTh1UxI5YQPTYPOtpx0fBEHoc6wMun4Lk6QgvhxP0b79+TyQwSxIJqZuYnTQAgOCBovRkt2sUkC8NcTtmyl43nMDyTszyDF/sJUwuu2KfPtCG3xMbjR1I5OpJj4a6x8DoFcNoibHELO3MrZgmbwEghpG6PY9rIPF30+ifPRh8K8ezSEfyiYnb4SR3lcCQehDjC40fRPTP4R4d7RfnGIdVo+gMtXEnG0MrMEGAt6j+WyFCX8/ieZ3psQ1eagf+WQuOWEDW2GCLLw8//kwA53eiaZuYkZEEO+P8acURAgRCP13NCWhiDv3wjr4IJC9dJLbUY83zqBl3lrg0FPuziSfHEBO3ghZeHkQhD6hPQXHRhFvjPCzFGxHk+iHSVSVBf89H+bxgsD0SQn3aQm3eSat5m1xpx5Znkk+PYCctJEthyy8BAQh/3RONGUTMzGaWD3cV0eIdoGERr9PpTfXYdizCQSeP2q4Z4+xm2dQUf48FeHuTPK5geSkDWyZEbLwL/ztCT/gaJ1o6kZmaizxyjA/TsF2GhHaPIMa/QcbLUVe2IMGAO842ILv2seun0anqfyxveYv7swgCQJN3shun0Wl+/8/x1MgCPmkd6IZm5gpAZGC7eJkxKYZ1KQNTJiYmBYLpxnwe+UmPG87+8V4z280z5c70kkJhSZuYLfOpLI1AfKP6iG4bedNix1N2MBMjyNeDZQUbNdXTfw4mV6ymynQQfML8G+NNjRtE/vSEHJGXEAFxsIU8uVh5PTN7FkDnKQIQRDypcmOJm5gbuxDvDgkoFKw3bgo4t1RVN4WtsYCpxnwV2Y3um4Lc2cGeXt6AF4nl6SSrw4jp25kT8ENKzSN8qLeiidvZJemkU8NCMATrN28JLLeimZuZvdfT2t8acVwALrCzaH5O5hh4cTTgXuSLkohaQJN38xsmE4PDA2oR95rFbCfsc+qMuPxG9i7MgM5Bds91I+cGU/csI1xwJQK4FcwQnfvY4Uk8e7oAGywudDNyeS7o6iZm5n8lqB+Lgzwa7GvKTPiCRvYh/uRj+YERc2/OoyKkRFLdrNcUJ9lwM88e5QtNvrNImo9dGMf8tNx9OxtzL6m4D1Lg+Jy7CNO6fDEjexzg8j7+gZLtZME+nwc1ebAj8BEe+An/lvE/ViJ/5hGS4Om42hWPPHdRHreDmZzXZBmYbBckXmX34KnbmL+M5y8IxA73q9ARKG1U+mdDfh1mGgPfN7v1dwLJ7jNM6hwMd9F8a5JMcRvU+nb9zA/VwbjeRpcF2W+bKvHs7cxn42jb04OxgpXC9Gm6dQ7Z7hvK4LxHAP+4kAzvns/+/s0KlkRBE2ilxgRQWydSa84xH1UHHTnaTBel71sbRW3dDfzy2R6Znwwnl3tYmXExhnUI4fYLcHa9gJ8XLkJz9/BfjmeHhwWvOdp/xBiXx712inuiSPB1ZcBQdi73jrNPXSQ2zaLHhMVvGdXu75q4pcp9NI9zNE2yELgW+qseMpG9pVh5PTAmjjfDUkKYm8evaUOP3CADZ7dZCAIewtGaNVx9sNibt/1VD9YxwghhNDoSOKTsfT1W5gSWPMX+AyjC+VtYR/IJhenwvUQIYSiJGhfHl1uwnO3s0Gy4TZ88L3CzqAFO9hdDXj/9XSiHFLwf/ISiJeGUjM3sw02yELAPxuDrtvCTI0lVgbHjKYukgvQb9NopQBN3sg02/kuTe+Dz97zmuxowgZGRKGtM2FRlcu4I51cnklO38TqnXwXBQS39uVjUpXEq/68A1ovEZLoiwnU9Dhi5O/MGX2A37ZCEHrYCS0e8RtzXQL55QRKBCdXJx7PJWfEEddtYazB0fACfBCH0e17WJokPh4bDPPmu4NAaNUg6oXB5KSNzIbaQM5CCEJP+vEcN2Mz8/pw8rmBJJxaV/bqcCpLTdy0jXEG1/A04CseOMjW2/APkygaroJXdGsq+dtU+p797CsFAbtCFHwFPIPh0OOH2SeOcNtm0vOSoFavjkDow7GUUkgs2sUyQTdtCfDsySPskVb8+zRaDM02XTAigsifQ/1axd2yk7W4+S5NL4BLtgc029G0TUyhDh+5ge4fAo+CXUUR6JuJlJ3Fd+yFxUiB9/zrJLe+Bm+aQSsFfBfFf8RIib15tEaIhv3GFAXcLoYQhD21owEPXseMjyY3TKdDYWjMNRKS6JfJdK0V3/Nn8MxZAnx64zT3RRm3dSacrddMRKH/jqEezyUnbGA+Lw2oZhwIwu5zc+jpo+xte9gvx1PPDyKhV7B7JDT6Yxp9Vo8fOgi9haB3vXeWe+cMt2MWFS3luyh+67Y0cucsevUpbvFu1hQozaQQhN1UbMCjfmcKtPjEjfSkGMjAHpEL0IbpdH4LfuQQZCHoLR8Wc68WcjtmUXEyOGF7JFtDHJ5DKwRowK/M3oDYvAmC8JqxGL1xmhu3nrk7k1w/nQ62Vep7iUqItsyk9zXhlbBhE+gFHxVzL53kdsyi+gTlgtoeJ6HR+6Opd0ZRi3axj+azdj+fBwVBeG3O6PGYP5g/qrmDs+m/ZULteZJaiLbNovc24pX50F8IPOmjYu7Fk9zOWVSKElLQk2bFEydvpJvsqP+vzJ5GPz5r4VLeVTYGPXWEnbiBWZZO7riOhjOqN6iFaOtMen8TfvAAZCHwjHfOcv8qgBTsLWFi9PUE6o0R1NI97LK9bKuD7wJ1CwRhl/xcyWX/wtRYUeFcwd8yYVhML9KI0LZZ9AktXr4f5lSAnnr9FPfmaW73dZCCvSsvgTgzlw4VoX6/uN89y/ndzGAIwqs40orHr2deOsl9No76egIVJeG7QEFAKUCbZ9DlRrxkN+v2tzMK+I7nj7GflHB7rqNg4XsvkAvQ68OpnbPo36q53F+ZjX61JBsEYafOGvDc7exN29mlaeSxG+gJ0XAueY9cgDbOoM1uNHc764DRM+AaYYRWHGL/qMF78+hYGCPqRdkaYutM+pVh1GOH2XHrmf1+MqYUgvAyCnX4lp3spA3MqEiidD59ZwbMEeSBmEK/TKEUAjRjM2Nw8V0a4D/cHLptN3usDe+8DgZ18yMvgSi8ib4rg7xtDzt5I7PL58fRQBD+xfZ6PGsLM3MzOzScqFggWJlDSmi+yxTEBCT6agI1IJQYv55ptPFdGuAPrAyas40xuNDWGbRayHdpghhJoKVpZMl8+rY08v/2syN+Z34857t9hxCECCFkcKE1Z7h+vzAPH2LnJZHnFtArc0gZRKAPIAn05ghqYQo5+o8AXOEQeFaLHU3eyERLiF+nUHAL6wtoEi1NI8/Oo5/MJd89yyX/wPzzBOeDm3IH9ZfFzaFt9fjrcm5TLTcjnnxnFDU+moBGUB/0RC4ZJ0MTNzA/TKLHQ2ctuJwSI75uC7sklXxuEHRl+BaSQHMSyTmJZIEOf1DE5fzCjoggFqeScxJJqW9EkG+UwrusDNpez62rxn9UcxlqYnEq+c4oQQiswOvbFqeSsTJiwU7mteHUklRoyQB/sbsRL9zJvDyMui0Nvhu+KzeEeG809fpwal0191U5d89+dmoceWMiMSOe5HcN9GAJQieLjrTiPU14dyOX34KHRxDXJ5AvDKZh1UE/MjGa2HUdff1WtkiPXxxCwW0/aPdJCffMUfbbSfREaC3wB1IaLUohF6WQOif6vZr7uRLfd8CdpSYmxxDjosmREYTC69tjERj7VnNtXFxcfn5+bGxsVw622+1CoZCiLt5bk8Wo3orLTeiMHp/W42NtuMiAszXEuChiQjQ5MYaA/j+PMJvNCoXCy2+qdaKbtjEhIuLLCZT3Txifwkv9+xSGQ38/zG6sxeunUWkqb6cg1L+nuDj0ZzPe3cDtbsTHtThRTgwOI/ppiJwQIlWJEuWE4HLP+R6sf/8OhF1NxJ5WTBCsxY1MLmRyoxY7brKjBhsOFxOpSpSlJgaEEsvSyQGhBGxFHRhCRWjbLPrBA+zI35l1U6lUWDEkWGmdaMEORkCi/DkwQNS/CUk0MZqYGE0hhBgOnTHgY234jB5vO8VVmFC9DYeLiRgpipAgpYBQCtF9fcl+Gk+e+P4dhCISqQWIJIk4GVIKkEqIwsVktBTFSAkRxF7gEpLov2OoD4q5MX8wH4+l8xIgC4PO8TY8fwc7L4n411CKgs8/gNAkyg0hckP+96EyHGqy4wYbarZjkwtZGKT0dFOQfwfhqAg8IY6gKOgeD0bLM8ncEGLBTja/lVg1CK6GQeSTEu6po+x7o6i5SXDuBz6aRHEyIk6GEOqtkxy+RsCPjYggjt5AH2rBUzcyPjg5CXicxY2W7mbfPM3tzaMhBYGnwDcJ+LdwMdoyg54UQw5Zx2yugywMZCe1eMg6Rkihw3PoDK8PjQEBDIIQ+D2SQM8MJH+YRN+zn334ECzSHYAwQv85xU3fzKwaRH48FlaNAR4GQQgCxNgo4sSNdIMNDfuNKdDBo2HgqLXi6ZuYtdVc/mz6lhS4ZAHPg28VCBwaEfphEvX3/uS0TcxLJ313hV/QdV+UcUPWMROiyd3X0X0U0BwKegU0MYBAsySVnBhN3LWPXVfFfTKO6h8CV0+/VGfF9+xn621o20waPkTQq+CJEASgOBmxaQZ9b19y6ibm2WPQa+hnOIzeL+IGrWVGRJCH50AKgl4HQQgCE4HQsnTy5I2CYgPK+YXZXg+9hv6hQIdH/8F8W8HtzqOfGUhedm0tADwLmkZBIIuWop8mUxtq8d/2s0PDidXDSVhm3WcZXei5Y+z357iXhlB3ZsCa6sB74HYLBL7r4okzc+lMFRq4lvnXSc7O8F0g8FccRp+UcFk/ux0sOjtPcBekIPAuCEIQFCQ0+sdg6vAc+qQWZ/3MfFvBQVOpj9jRgAevY74s4/6YRn8whuJ3XzoQnKBpFASRJAXx42RqfxNemc+uPsW9MpSaEgvPHrw5qcVPHGErTOjloSSslwZ4BEEIgs6YKOLQHPqXSu6+A2ysFL04hBoVCXHoVcUGvOo4t6eRe2Yg9bdMGBEDeAZfQBCMCITmJZFn5tKL08hbd7MzNzMHmqGt1BuKDXjpbnb8BmZAKFG+QHBfX0hBwD/4DoLgRZPojnSyZD59Ux9y8W528kZmG8yy6DXH2/AtO9kJG5hMNVF2s+CJXFIGDVLAN8A3EQQ7IYnuziSXpZPfVnCPHGKFFHo0h5yXBE8qnoER2lqHV59iiw3o4Rzy47ECuad3VQWghyAIAUAIIZpES9PIJWnkhhr8xmn2scPcfX3JOzPIcDHfJfNbVgZ9Xc69fZoTUujhfuTCFLi3AD4KghCA/yEQyksg8hLoAh1ec4bL+MmdF08uzyJHw2iaa3FGjz8s5r4p58ZHk++MpiZGQ+0BnwZBCMBl5IYQH4+lXhtGfVbK3b2PJRC6I4O8NZWMkvBdMh9mcqMfz3GflXLVFnRHOnHiJjoe1vEB/gCCEIBOaUTokRzykRxyXxP+rJTr+7N7VASxKJWckwgDPf7HzaGt9fjbcm5jLTcphnwil5wVT1KQgMB/wNkMwNWNjSLGRlFrGGptFfdNOXffn+z0OHJ+EjEznpQG6znEcGhXI/6pkltXxaWriEUp5NujBLAuDPBHwXoSA3DtZDRanEouTiW1TrS2ivuwmLtzHzsxmpydSMyKJyODo9XU6EJb67nfq/GmWi5NRcxPIo/eQCfI4QEQ+DEIQgCuWagI3ZVB3pVB6p1oYy33ew1eme9OURAz4olpseTICIIOrOGRHEYntXhbPd5Sxx1tw2OiiNkJ5CvD6Bgp5B8IBJ4MQqfTKRJdqWXE7XYLBDCHCAQOjQjdmkremooYjvqzGW+t5x45xJYa8ahIYnw0OTaKGBJGiCi+S9ktLEYFWry/Ge9pxHsauQgJMS2OeLQ/NSGaCNrWYBCoPPONrqurW7hwYWFhoVAofOONNxYvXnzRAc8///xnn33W2toqk8nuv//+VatWeeR9AfARNInGRxPjo6mXhiC9E+1t4nY34hUHuWIjzg0hhoUTw8KJwWFEqorw5WeoGgs+1oaPtuFDLfhoK46TEWOiiLlJxDujBNFSvgsHQK/xTBCuWLGiX79+e/bsOXbs2KRJkyZOnBgbG3vhAeHh4du2bcvIyCguLh43blx2dvb8+fM98tYA+BqNCM1JJOckIoSQxY2OtuH8FvxTJX7yKKdz4P6hRI6GyAkhMtVElprgcT6GzolKjPiMHp/R41M6fFKLBSQaHEYMDScfzSGHRxAhMPIFBAcC454urmgwGMLDw8vKyvr06YMQysvLGzdu3GOPPdbZ8fPmzcvOzv7HP/5x2f83Li4uPz//ohztjN1uFwqFFOWfbU/+z2w2KxQKvkvhT3ROVKDDp3T4rB4XGXCRATtYlKokkhREkgIlyIkEOYqREjFSFC6+eptqV+qf4VCLAzfbUb0V1VpxjQVXWVClGZcZMYtRhorI1hBZaiI3lMjREPDYd03g+88vD9a/B54Iq6urhUJhewoihPr27VtRUdHZwTqdbv/+/ffcc09nB2CMjUajVNpxRorFYokkOEbjgSAQIkITo4kLV1oxuFC5CVeacZUZlZvwzgbUaOMabKjFjkUUChMTISKkFiKVkJDQSEYjCY3E/z8gXS5aKGTbf3ZzyOJGDhbZGaR3YaMLGVyozYHNbhQuRpESIk6G4mREnIy4PgElKchUJQGrxwHQzgNBaDAYZDLZ+V8VCkVlZeVlj3S73UuXLp0+ffqUKVM6ezWtVjty5EiS7Bh1N3Xq1E8//bSzg+GJkF9Wq5Xw6T4vP0AjlClGmWKEwi/+v0xuQutEBhcyuAiTG9lZZGcJO4OcXEed09glQv//Zwr3ESMhiaU0UgmQSohUAhwiQiHCTpp8GGSx9NY/KkjA959fXax/sVhM01dJOg8EYVhYmNFoPP+rXq+PjIy89DCWZZcsWYIQ+uijj678al1vGqUoCoKQRxhjuVzOdykClhyhmCseYDa7FQp4rOMNfP/55cH698B0pz59+tA0febMmfZfT5w4kZWVddExHMctW7ZMp9P9/PPPQqGw528KAAAAeIQHglAmky1atOjpp59ubGz87rvvTpw4sXDhQoRQYWFhXl5e+zH33nvvjh07Hnjggf3792/fvr2oqKjn7wsAAAD0nGemT6xevfrhhx8eNWpUdHT0b7/9FhISghDCGDMM036AxWLp27fv22+/3f7rrFmzLn1qBAAAALzPA9MnPAumT/gRGD7OL6h/fkH988uD9R9YSyICAAAA1wgWDQxYeofhdGvR6dbialNto6VZa9e7OcbBOJRChUqsjJCGpWiS0jRJg6MGhEo0fBcWgOBicVmPNxeWaMsqDFUNlmajw2RymWhSIKKEKpEyWh4Zr4zpG5aRE943ShbBd2EDHwRhoGmzabdV7dlXe6jGVNcvPDM7LHN22oxoWWSoNERICsS02Og0GZ3mJmtzhb7qz7rDa459HCENGxc/akbypEjZJXPZAACeY3Kat1Xt2VW975yhun9E36zQjLzU6XGKGJVIoRIp3RzjYl16h7HJ0lxlqt1fm//e8c9UIuW4+BGTE8clquL5Ln7Agj7CwHGsqWBt6YaC5jPjEkZOSBg9MLI/TV69ZjjMFWlLt1Xu2Vm9LyM0dWHWTYOi+nfxHaGPhF9Q//y6pvov11d+d/bXQw1HR8YMnZY8YWBEjoC6+lY8HMbF2tI9NQe2Ve2OV8TOSZ85IWE0SUCXFkIe/f5DEAaCA/WHvzr9o81tn585e3Kf8RK6O5Os3ax7e/Xe78/+KqbFd+UuHho98Kp/BS7E/IL651cX679Mf+6jk19VGKrmZ86+PnW6TNCdFV0Zjt1fd+iXkvU6u35x9rxpyRMpItivexCEHSAIz7aVvHf8UzvjuC3nljFxI8geL/iEEd5bc/Djgq/CpWH3Db4zRd3nCgfDhZhfUP/8umr9t9raPjjx5fHmgttybpmVMlVAeqArqqDlzBenvm+1aZcPvG1M3PCev6D/giDsEMxBqLPr3zv+WUHL6TtzF09LmtjzCLwQi9n15Vs/K/xuatKEO/ov6uwREy7E/IL659cV6p/F7M/Ff3xz5ucb0mct7HtT9xppruBww/H3T3ymFqseHnpPgjLOsy/uLyAIOwRnEHIYryvd+BQzJQcAACAASURBVMWp7/PSpi3JvllM99aucUan6f3jnx1vPvX4iAcGR+VeegBciPkF9c+vzuq/wlD17wNvasTqFUOXxyqie+ndOcytLd345akf8lKn3pZzi5AKuqUrIQg7BGEQ1pjqXz30NkmQjw6/P0HZpVrqoWNNBa8cWjMiZvC9g5aJ/3pjCxdifkH98+vS+ucw982Zn38p+eOegctmJE/yQhl0dv2aYx+X6889NuLBnPDgWq4LgrCDp4KQdXIuo9ttZVknxzpZhBApICkhKdIIRBoBQfrETisY4Z+L//j69E+3979lTtosz7aFXpnVbVtz9KMzbcXPjn40PSTl/J/DhZhfUP/8uqj+m62tLx74j5ASPDnioTBpqDdLsq/20JtHP5icOPau3MW+82joNLgdWhfnxu0XVUpIkkJSIKeFcpqWeeDpBYKwQzeC0G1hbM1Oe6vL0ea0t7ocWpdT78YsFqpogZwmhSQtphBCrIvjXJxD53KbGWmMWJMu12QplEm8beDdYmv798E33Szz1KgVMfIoXsqwo3rfmqMfLs6+eW5mHoEIBBdivkH98+vC+t9Xe2j14XcXZN24IOtGb96knmdymlcffq/aVPfMqIdTNcneL0A7W6NDd9asL7GYa+y0iBSHCSkRSYk6Lqqsk3NbGLeZ4Vgs1gjEoUJxmFASLpKECyXhIpFagK6l5iAIO1whCDkGu4xup97t0LscWrejzWVvczpaXIhCknCRNEIkCReKw0TiMKFYI6ClnUYpx2Bztc1QatEWmgiKiBkbGjFETVBe/aLvrvnzzSP/nZ85Z2Hfm/idQtRoaX5+/yvRssjHRjwgE0jhQswvqH9+tdc/w7EfnPxib82BVWMfywpN57dIWyt3v3f8k0XZ8+ZnziauKVV6CCPtaVPDPq2jzRWWq1Sny5XJMkrU6cWKdXJOncuudTnaXPZWl73VaW9xMg5OEiaUhAvFoR1XZpFGIFQLKOHlXweCsEPTSa251EEQBOYw6+BYF8fYWcbGus0M5+aESkF722bHfUeYUBIm6v4jOUaGMkv97jan3p16c6x3ng4djHPNsY9ONp9+dvTKzNA0L7zjVblZ95pjHx9vLnxh3JNhpAYuxDyCIOSX2Wx208xz+16RC6VPjXpYKfSJz6LR0vzCn6tlAulTo1ZoxGovvKO9xVn+UwPr4uImhYXmKLvdl8Q6OXur09Hmsrd1tNU5DW6XwY0QEihogYyipRQloigRGTM2VBYrhiDsoK0wOJsYkiQRQdASkhSQtJQSSCmBgqYlvTWCRltoOreuMbS/Mml2VK92H1YaqlftfzUjNO3hofd4fPh1D7XfeC7vd9vMjCl8lyV4QRDy63D18VePvzM7bfqSfjd79fHraljMfl74/cZz258aueKy4709qGGvtnZ7a8K08OjRob1UB6yTc5kYxsYydpZ1sqyT02TIhSoBBGEHvkaNMg625MtaRBCZS+Ov8PjfE3+Ub/n45Nf3Dlo23Stjz7qhXF/59O6XJvUZe/eAJbDmEy8gCHnUfoY+NWrF8JjBfJfl8k40n3rpwBvTkyfe0X9RbyxDgzl8bm2jqdLW985Ekebqy8V5HARhBx6nT2AOV/zaaK625fxf0hW6GLvB6ratzn+32lT3/JjHvDNBotsatI2vn3yPJunnRj8qF8r4Lk7QgSDkBcOxbx/9sKDl9NNDHk6PSuW7OFdicBpfOvCGg3E+N3pluDTMg6+MWVz0eQ1mceZtCb30MHBVsB8h/wiSSJ0Xo06Xn/20hmM8djNRqqu4e9PDcqH8/emv+XgKIoQUQvlrk1YlKGPv2fJojame7+IA0OsMTuMjO57R2nXvT389WhbJd3GuQi1SvTrx+RExg/+2eeXB+qMee12Myn6sRwj1vSuRrxT0rED4N/AoKS9KHCIo/qIGcx7Iwl9LNjy2a9XduUseGfZ/vjMZ6Moogrp/8F239p374LYnDzcc57s4APSiCkPVPZsf7R+R/cK4p6QCCd/F6RICEbdmz/vn2CfeOPL++8c/Yzi2569ZvanZ1uzMWBLvI3Osew6CsGcIlHpzLOvgare39uRlzC7LM3v/vblyx3vTX5uYOMZTpfOamSlTXhz31CuH3v6xaB3fZQGgV+yrPbhyx7N/G7D0rtzFvMwU7Imc8KyPZ75ZY6p/YNvjjZbmnrxUW4Gp9aQx+67EzmY1+KPA+ZfwhaSJzNvimw7ozNX27r3CqdaiOzeuiJZFvDvtVb4my/dcv/DM92e8vq1qz78Pvulm3XwXBwCPwQh/ceqHNcc+fnXiqkmJY/kuTjcpRYp/TXh6cuL4/9vy6K7q/d17EZeJqfi1IWNxvEAeUJu6QxB6gEBOp8yNKfm6lnVy1/QXOcx9cer75/a9/PDQe+4bfKdHdmnhUYQ0bM3Ul52s68HtT2nter6LA4AHOBjHqn2vHm489t/pr1+4vqA/IhAxL/P6Vyeu+qTwm1cOve1gHNf29zEq+6E+ZnSIIsE/moW7DoLQM0JzlMokadWGa2hzaLK2PLjtqYKWMx/PfHNk7JDeK5s3iWnR82P+Pip26D2bVxZry3rvjTgG25qdhlJL+3+2Jue13oUAf8Qx2N7iNJZb2z93a72DdfXi595kbblv6+NSgeTNyS+FSDS990belB6S8tHMNzDGd296uERb3vW/2HhQx9jYuCnhvVc2vvj3I4hPSb4x+tjLZVEjNbLoq09+33Ju53vHP1uYfdPNmTf4XX/DlRGIWNLv5hRNnyd2/9Pja/DbW5xtBSZtodHW7BRphB2LE2LkMrkdejctoRSJEmUfqTpD3pVPAfgFh85lKLaYKm3mGptT7xZpBEJVx1L4bjNjb3WKNMLQHEVYf5Xco08qJ5pP/fPP1xdnz5ubcb0HX9YXSGjxEyMf2lW9//Hd/5iXMXtR9tyrTgVmbGzN5pace5MCZoDMhWAeoSc17NPqzpr7Le9zhWMMTuPq/PdqzQ3PjHokVZPkraL1iivP46k21j6991/DogfdN/iOns/ndehc1ZtajKWWsIGq9ufvS09Ih85lrrKbqqz6IgtmcWiOMmygSpko9aVFPzwpsOcR2pqcrScM2kKT28ZqMuWqFJkiUSqNEF38aWJkaXBoC42tx43iUGGfvEh5nAfi8Kfi3789+8tzox8dGJnT2TEBUP+ttrZ/H3zLwTifGrUiThFzhSPPrWvkGJw670rHeBlMqO/ga0GIWXz8tfLkOVGarMt/PHtqDrx19IPpSZPu6L9IQPGwFoNnXfWLaHFZXzyw2ua2/2Ps4z1Z9rB+T1vt9taYsaGx48O6OG/J3uJsO2lsPWnkXDhiqDpiiFoc6h8zUrouAC7El3JbmJZjhpYjBsbOhg1Qheeq5PGSrtzKYBY35+trtraE5aqSZkd1e2V8B+N8Lf+dGlPdi+OejJRFXOHIwKh/jPDako2fn/rutpwFN6bnXbaByt7qKny7YtDjaT41RgaCsIOvBSFCSHfGXLWhaeCjqRc9r+gdhjePfFBprHl8xAPZYZl8Fc+zuvJF5DD+4tT3Gyq2rRrzWL/wa/6Hsy6u/Id6u9aVdXuCSN2dWwdrvaP5iL71uFEWI44aqQntp/Ty5iG9JzAuxB0wMpRZmg7qDKXWkH6KyGEaVbKsG4/yrIMr/bbObWMzb4sXKq75ql1vbnx277/TQpIfGXav6GpzeQOp/uvMDa8cWoMx99iIBy9dyqPosxpFoiRukm/1DkIQdvDBIEQIFa45FzM2NGyAqv1XjPDGiu0fnfzyupSpt+Xc4i8z5bui61/Eg/VHXzn09q3Z8+ZlXt/15YlZJ3f6v5XSSHHKvBiS7lF6cQzWnTY1HdTZmp2RwzRRI0N4WR3RswLjQsxY2eYj+sYDOlpERo0MCR+kpsQ9G8SHUe321uZ8fc79Sdd087Sv9uDqw+8t679oTtrMrhwfGPV/Hofxb2UbPyv8bm7G9Yuy554fxG5tcJz5qHrI0+k9PAc9DoKwg28Goe6suWZzy4BHUhBClcaa/xx+n+HcK4fd5+89gpe6pi9ik7XluX0vR8kiHh/xoExw9U2sOAaf/bhaHCJInR/rwU4+e4uz8YCu5ZhB2UcaPTpUkyH33x5Ef78Qm2vsjX9qdafNIf0U0aNCFYmeHOpSv6et+ZA+5/5kQRd2XnNzzIcnvthbe/AfYx/v+n5n/l7/l9VsbX3r6Ad15sZHht4zIDIHIVT6TZ00Whw3yZNLlXoEBGEH3wxChNHx18pi8kJ+sa/bVrV7Wf9F16fOCLChoe2u9YvoZt3vnfjsUP3R58Y8epUtTDEq/rIGISJjSVxvjFLjXFzrCWPjnzrGzkaNCokcpunK5dLX+OmFmHVxrceNTQd0jIONHhkSOUzT/V1Cr6h6Y7O+xJJzX9KV10BptDT/Y/9roRLN4yMfvKY9Bf20/rtiX+2hd459nB2WcXfy0ur3dUOfzujpY3ovgCDs4JtByGFu24Y/244bG6dX3z1giUqk5LtEvaV7X8R9tQf/c+S/N2fOWZB1Y2f3Bw17tW0njTn3JfV2f565xt50QKc9ZdJkyaNGhnSvX4ovfnchtjY4mg7qWk8YVcmyqFEhXngcL/22jhSQqfM7Heu4o3rfmqMfLc6ePzcz71r3FPS7+r8mTtb1zZmfzVtcKWGJU5aMFPvYlqgIgvA8HwzC/XX5H5/8KkwUOvfPm3PuTpHH+ty3x4O6/UVssbW99Od/SJJ8auSKS3eHsTY4Tv+3KndFsjjES/2pjJ1tOWpoOqjDLIocrokYqu7GOAvv85cLMeNg204Ym/P1LhMTOVwTOVzTvXFP3cA6uROry5OujwrNufh+1Oq2vXnkgxJd+bOjV6Zpkrvx4v5S/93GWNnD/yrZP2nXYdOx23MWzEyZ0hv7GnYbBGEHnwrCY00FnxR842Sddw9YMiJmSN2uNnuLM22Br2+l1BM9+SJyGH979pefi3+7d9Cd05Im/O/PGXzyPxVxk8MiBnd/ukW3mattTYf02kKTMkkaMVQdkq30tQECF/LxCzHmsLHM2nzUoD9rVqfLIod74xHwUuZqe9Gn1QNWpgqV/7u5OdF86uWDbw2PGXzvoDvEtKibr+zb9d9z9bvbbI3OtIWxxdqyD09+2WJrW5azcGLiWB/p6IEg7OAjQXik8cQXp743Os2359xy/lvitjDH/l029LmMwNiv67J6/kUs05976cAb8YqYR4b9X/tEw+rNLfYWZ+bSeA+VsTtYF6ctNLUcNVjq7KE5yvCBKlWqzAcX1PDRCzFG5hp76wlD20mjSC0IH6yOGKTupV7ALqrZ0mJrdGTenoAQcjDOjwq+3FNz8LHh9w+LGdSTl/XR+vec46+UpS6IVfbpGNp2rKng08JvzC7r0n43T0wcw/vTIQRhB36DkMXsnpoD35391c0xi7PnTUoce9EyRcWf16gzFVEjAmSJwkt55IvoZt2fn/puQ8X2ewfdMV4z+uTq8gErU73WdHZlLhPTfkF3aN2hOcqwXKUqReY70xB960KMkbnG1lZo0haYSAERNkAVPkglCe/mw5ZncQw+/nJZ2sLYc/Jzrx5a0y8864Ehd13TuJjL8q369zRTpa38x/pBj188hvZo48mvTv/YYmtbkHXDjOTJ3X6e7jkIwg58BaHJad5QsW1t6YZoWeSCvjeOjB1y2W52fZG5ZmtL7kP+vWL9FXjwi1iqq3jl0NtTiqf3T8vMvr47HTa9yqFzaQtMbYUme5tTk6kIzVaoM+S0JHDuiLuNc3HGCqv2tFl31kxLyLD+qtD+SlmMz3WN1x5pOrO58rOcTx4eds+IGM+sce8L9d97yr6rl8aIYsdfftbE2baSb8/+eqrl7KyUKTekz7zyEjy9BIKwg5eDECNc0HJmfdnWgw1HxsaPnJuRd+U+dszhoy+WZt+dKA3QBaA9eyEwVFkLPil7O+etG/vOWpB1g28uQecyMbqzZt0Zk7HCKosWazLk6nS5PEHCS8MpbxdijKxNDkOpxVBiMVXa5PGSkGxFaLZSHOaLi0VghDef2/nRia+WlyzPnpwSN9xjl+wADkLWwR15oWTwk1dZU63R0vxr6YbN53b0Dc2YnTZjeMxgmvTe3SEEYQevBWGduWF71Z4t53aJaFFeyrTpyRMVQnlX/mL15hbWwSbfEN3bJeSFZy8Ep96rjBiqxn3d7xz7uNJQ88CQuzx1594bOAabKqz6Eouh1OLUu5VJUlWKTJkklcdLvNZ26s0LMeawrclpOmcznrMaK6yUiFSnyzXpcnW63AdnmJ1Xoi1/8+gHCKEVQ5fHmGOLv6gZ8nS6pz6gAA7CpoM6Q6kl87aErhzsZF27qvevL99aZ26Y2mf8lKTxGSGpvV1CBEF4Xm8HYY2pfn/doV3V+7V23YSEMdOTJ17rB+zQuQreOjfs+QwfHGrRcx78Ipqr7SVf1Q5+Kq29og43Hl9z9KMoWeS9g5YlqRM98ha9x21hjOdspgqr8ZzV0eaSxYgViVJFgkQeLxGHCHtvnGRvX4hdRrelzmGusZlr7JZqu0BJt+e9KlXmI524V9Bqa/uo4OujjSfuHrB0RvKk9s6LMx9WhQ1QRQ7zTLd9AAdh4Zpz8VPCO9s8oDP15sYtlTu3V+0lEDEhcfS4uJHpoSnXOjuz6yAIO/RGELpY16nWovyGYwfrj9oZx5i44eMTRuVGZF91v67OFLxZkXhdpDqtS0+Q/sWDX8Siz2vUqbLoMaHn/4Th2N/KNn51+sfRccOX5SwMk4Ze4a/7DtbJmWvslvbwqLWzTk4WI5bFiKXRYmmkSBol8mDPomcvxKyTs7c4bc1OW6PD0uCwNjgQRvI4sTxBqkiQKBKl/rL4jtVt++7sr7+VbZqdNuPWvvOkgv+t3GYos577tWHQY2keuTgHahC6jO4Tr5cPW5XZ7UfnEl35npoD+2oPORjH8JjBw2MGD47KvfCD8AgIwg6eCkIH4yzWlha0nC1oOV2kLU1WJ7Z/eOkhHridqd/VZte6fGofL0/x1BfR3uosXFM59Jl08pKlsCwu67dnf/mjfMt1KVNv6XujWqTq+dt5k9vKWhvs1gaHrclpa3LYm10EhSThInGYUBwqFIcKRRqBSCMQKgXdmLDYvfrHHHabGYfO7dS7nXqXQ+uyt7nsrS7GxkojRZJIkSxK1J7cvv/YdxEH41hbuvGHorUjY4cu678o4pK1GhBCJ9+oSJgWEZLtge9toAZhw16ttdHhkTnQtab6gw1H8xuOnW0r6aNKGBSZkxPRt19Yllwo6/mLQxB26HYQOhjHOUNNhaGyVFdRrC2rMdWlqJNywrNyI/sNiOjn2TuXAG4d9dQXsfynBqGSTpje6SiGNrvu69M/7ajee33q9PmZczRiP4vDC7nMjL3F6dC6HFqXU+926NxOvctlYmgJJVTQgvb/pBQto2gJRUsoSkRSIpISUwSJ2p8mCYpon5xqsVjkcjlCiHNjzs0hhFgHhzFmbCznxqyDZRwcY2fdFoaxsW4r6za5XWbGbWUFMkoUIhSpBeIQoThUIA4VSsJFIrXAj5aXu4iDcfxWtvn7orW5EdnLchYmqjqdh9p20tiwT9v/AQ+MTA7UICxccy5+aoQm05ONWC7Wdaat5ETzqVOtZ4u1ZeHSsKzQtPSQlBRNUqo6qXu56MH694N1pHrI6DS12tpabG0N5uZ6S2OduaHGWGdwGhOV8SmaPmmalJnJU9I0Sb03RlEcIhRrBMYKmzrNAzdBgcdtYdoKjIOfvNIa3GGSkBVDly/KnvvNmZ+X/nHv1KQJC7JuiJT51u5oXSRU0EIFrUq5+MvgNjMuC+M2My4zw1hZt421t7pYB8s6OdbJMQ4WcYixswghzGLWySGEMMYEQSCESAFBCkiEECUmCYKgpRQpICgRRYtJWkqJNAJ5rEQgpwRKgVBBC+RUIN2TmVzmtSUb15auHxCZ859J/7xqj3Jof2XVhmZztd2zm10EDJfRbW9xevxiJaSEAyNzBkbmIIQ4zFUaa4q1ZSXa8p3V+yoNNRJanKiKj1VEx8ijYuRRkbLwCGmYWqz22hI2/v1EeLTuZLGhjCRJm9vOcqyDdTgYp81tN7ssRqfZ6DQZHAaJQBIuCY2UhUfJI2Lk0XGKmERVXJQs0purBNXtanMEYuuoR+7I6ve02Rqdabd0tR1GZ9f/WPzbhoptQ6MHLsi6wTvj03xToD6RdFGDpemn4t+2V+4dEz98Yd+5l24n25m6XW2OVmfqzT1t+gvI+vdgu2jXtdjaaox1deaGektTo6W5xdraYmszu8wqkVIlUipFSqVQLhVIRJRIJpASBDErZUqcIgaeCDu4WJfFZSUIQiqQUAIqjAoV0yKZQCoXyFRipVqkVIvV57eX5FFYrrLgzXMpN0UH0p24pzTn61PnX8NZFyLR3DPw9iX9bl5fvvXZvS+HSULmZuSNSxjlCx808AIO46ONJ34tXV/UVpaXOvXzvHdCJdc2CjRyiPrYy2VJc6IDePnDbmsrMMZP9fbs+AhpWIQ0bEj0gAv/kMWswWEyOI0mp9noNNkZh5Nx2hh7bxTAv68dAyNyhscN5n2t0asShwhFGoGp0nZpg1iQM1fbMIfOL2bYdTKBdEHWDfMzZ/9Zd3ht6YY1xz6ekTwpL3VanCLQHrvBeVq7ftO5HRvKt8qFshvTr/vH2CdEVHem8AsUtCpV1nbSGDk8YJc/7B6XmbE1e75dtHsoggqVaK71Lqd7/DsI/UhIX4W+yAxBeJHmfH3kcE23x2iQBDk2fsTY+BF15ob15Vsf2PZkrDxqRvLkCQmjPTIsDfgCF+s6UH9k87kdp1uLxyeMWjXmsYzQnraHRw7X1O1ohSC8iP6sWZ0u953VdL0GgtBLNJmKsh/q++TxXQ5fwrq4tkLToMcuXtW3G+IUMfcMvP3uAUvyG45tPrfzveOfDorKndJn3IiYITwuCgx6guHY480FO6v27a/LTw9JmZE8adWYxz31aWoy5eU/1tuandJI+Hr8j77YHNI3YDcSvwIIQi9RJEjcFsapd4s0fjY3q/e0FRhVybILd4nrIYqgRsUOGxU7zOq27a05sKFi26uH1gyNHjgufuSI2CEywTU3wALvc7Guo00n99Ue+rPucLwyZkLCmLsHLPV4+xhBEhFDNS2H9X2uj/LsK/svzGJDqTXlpmDsXIAg9BYCaTLl+iJz1KgQvoviK1qPG3tpjyqZQDozZcrMlCkmp/nPuvztVXtXH34vMzRtROyQETFDuj68EHhNq62tfTmnE82n0kNSxsaP6GxGvKdEDFKd+ai6T16U/86e9CxTpU0cLhQogjEUgvHfzBdNpqL1hAGCsB1jZy3Vds2yLq3q221KkaI9ER2M41hT4cH6Iz8X/04gYmj0wMFRuYOi+qtEwdgQ5CPsjONk8+njTQVHmk7q7Pqh0QMnJo55fOSDPd8psCuk0WJKRJprYEJhB32xOeQaFxcNGBCE3qPJklf83MAxuBuLaQUebaFJnSGjLllTrZeIafHouGGj44YhhKqNtUcaT26t3P1a/juRsvABkf36h2fnhGf5y3Kmfs3kNJ9uKypsOVvQcqbSUJ0Vlj44MveJEQ+mh6R6c2pvu9BcVVuhEYKwne6spevTeQMMBKH30BJKGi0yVVjVGQG4APe1ais0RQ5V8/LWiar4RFX8vMzrOcyV6ioKWs5sq9r95pEPRLSoX1hGVlhGVmh6qiYJRtl4BMOx5wxVxdqys20lZ7WlbTZtVlh6//C+9wy8PSs0Tdit+Q+eEtZfWfRpTRK0jiLkNLjdFkYRH6T3BBCEXqXJUuiKzBCEjJ01V9oyl3a6IKR3kASZGZqWGZq2IOsGhFCtqb5IW3q2rXRH1d5KY020PDJNk5SmSU7RJCWr+/j1AqfeZHXbzhmqKvRVZfpzZbpz1aa6WHlURmhqdnjmvMzZyerEbm/k4nGyGDFBEZY6uzxYA+A8fZFZkykP2hsCCEKvUqfJy3+q57sU/NOeNqnSZL62rke8MjZeGTstaSJCiOHYSmN1me5cmf7cn/VHzumrSIJM1iQmKOMSlfEJytg4ZUyENNz7rXm+ps2uqzM31Jrqa0z1VYaaKmONxW3to0pIUfdJD0mZlTIlRe3Tz9Zhucq2AhMEoaHMGrQdhAiC0Mvk8WKn3u22MAJ5UNe8ttAUPsCnH7BokkrTJKdp/rdHgdaurzLWVBlrq4w1++sO1ZrqDU5TtDwyVh4Vo4iKkkVGySIiZeERsjC/2yuqK6xuW4u1tdHa0mxtabK0NFiaGixN9eZGMS2KU8QmquLilbGDIvsnqRMiZeG9txerx4Xmqoo/r+mTF8l3QXiFkbHcmhTEM0mC+nLsfQRJKJOkxgprWG4AXiu7iHNzxgpr+q1xfBfk2rSv9jQ4Kvf8nzhZV4O5sd7S1GBparQ0n2gubLa2tdhaHYwzQhoWKtGES8NCxOoQiSZErFaLVWqRSi1WqUQKMS3m8R9yWS7WZXKaDU6T3mEwOI16h1Fn12vteq1d12bXtVhbSYKMkIZFyiIiZeHR8sissPQYeVSsItrfZ2fKY8UIoyCfWW9rclBiMpinOEMQepsqTWYsD+ogNJZb5XESWuzrK8RelYgSJqkTL933x8m6WqytWoe+1damdxi1Nl2lsUbvMOjtBqPTZHKZOYyVQrlCKJcL5TKBVCaQSAVSuVAmocVCSigXyoSkQESLhJRQRAkpgpT8/w0yFcK/9C5bbVYzYb3wT2xuG4s5hJCDcTIcw3CMnXG4ObeDcdrddifrsrntVrft//+vzeyymFwWs8vMcpxKpFCJlGqxKkSs0YhVIRJNsjoxVBISKgmJkIZ5fHtx36HOlOuLzMEchIZyqzotqAcuQBB6mzpVXnywlu9S8ElXZNFkBnJvhIgStnc3dnaAk3WZnWazy2J2Waxuu81tszF2s8viYBwGh7He3OhkXS7W5WSdLtbNsazVtQAAIABJREFUYs7utiOEMMIW119ij+M4kvxLP6tUIKUIEiEkokUCkqZISkpLBKRATIsktFhICxUieZQ8QkJLZAKpTCiVC2QKoVwpUkh87yHVa0KyFA17tbETenHyvo8zllvDfLurordBEHqbLEbMWBmX0S1UBWlDhL7YnNXL8+h9nIgSiqShPZ+2GJD74XmfKk1W8k0t6+R8bfSWl2BkOmdNCbjdUq9JUH7w/CKQMkVmrLBe/chAZG9xYgbLooL3+QP4GkpIKhKkhjIL3wXhh6XeLpDTwqBcWe08CEIeqFNlxvIgDUJ9sUWTpfCfQYUgKGiy5PqiIA1CY5lVFdwdhAiCkBeqNLmhLEiDUNc+bxcAXxKSpdAXm/kuBT+MFVZVarBv3glByANphIh1ck6Dm++CeBvr4szVNnU6BCHwLZIIEUEStiYn3wXxOoxMVTZVsn/Pgek5CEI+EEiRKDFX2/guh7eZKqzyOEmQDkkAvk2TIdeXBN1Doa3JIZBRQb6+B4Ig5IsySWaqDLogNJZb1UHfCAN8U/sEX75L4W2mSpsyCU5JCEKeKJOkQRiEhnLolgc+SpUiM52zYQ7zXRCvMlXaFH2CvV0UQRDyRR4vsbc4WSfHd0G8h3Gw9hYnrG4MfJNATgtVAmudg++CeJWpyqZMgiCEIOQJSROyaLGl1s53QbzHVGFTJEphU2Lgs9RpMkMwtY66zQxjY6URwbu23HkeC0KbzVZTU8NxnT7isCxbXV3tcATXDdcVKIKsddRQZlGnQW8E8F2qVJmxPIhmExorbcokKUzqRZ4Kwrfffjs2NnbSpEmZmZlFRUWXHnD06NHk5OSpU6fGxMR8+eWXHnlTf6fsIzVVBtHtp7EcpisBn6ZKkZmqbJgNlm5Cc6UV2kXbeSAIq6urn3rqqQMHDpSXly9atOihhx669Jjly5evXLmytLR069at9957r1ar7fn7+jtlssxcZQ+SznnGxjp1bnkcdBAC30VLKUmoMHg6LGCkzHkeCMLvv/9+4sSJWVlZCKF77713586dzc3NFx5QVFR09uzZu+66CyE0ZMiQ3NzctWvX9vx9/Z1ARgkUlK05KObwGsutiiQpQUErDPBpqlR5kHQTcm7O1uRUwOA1hJBHdp+orKxMT09v/zkiIkKhUFRXV0dGRl54QGxsrFTaceuRlpZWVVXV2auxLFtYWNjU1NT+a3R0dExMp8uiY6fN3VLDkv465EcWzuhPVAsYfy0/a7O5pF26o9SdZORhhKu2rLeLFFS6Xv+gi2RqrukEG5lh6MrBfl3/5gYsDsFMcwXfBblmBEHS0YkE5clFADzwWmazOSzsf1t5yWQyo9F40QESieQKB1zIaDT+/e9/Fwg6tigaPHjwG2+80dnBlj83oTMHul903tn66SrDqNO7+S5HN3EcZ+/aXYhBe3O4Yq+2rKm3ixRUul7/oIs4LLK0LGn77l0CXb3Pwq/r32DtT7Eq7Xf7+C5Id4in3Uqn9LdYujSySSwWnw+UznggCCMiIvR6/flf9Xr9hY+D7QcYDIYLD+jbt29nrxYSErJly5bY2E43Nb0QPWmucMZCivLXvc4VNfayn+qjV97Md0G6qYv74bEuruq54oQnnoW5E54F+xH2huZXylS3rpbFXn2nML+uf9PXtVGZioghC/kuSI94qv49cDuTm5t7+PDh9p9Pnz5NkmRKSsqFB/Tt27elpaW+vr7918OHD+fm5vb8fQOALFbsaHWxrgCfVm+utslixZCCwC8EyapP5mq7IgE6CDt4IAjnz59fU1Pz+uuvFxQUrFixYtmyZTKZDCH0xBNPvPzyywihyMjIefPm3X///YWFhc8//zzGeObMmT1/3wBAUIQ0SmStD/C5leZKWL0C+A1FH6mpKsCDkLGxjI2VhMNU+g4eCEKZTLZ9+/YDBw4sX758yJAhr776avufx8XFRUdHt//8/vvvJyQk3H333aWlpVu2bKHpYF/s/Dx5gtRcE+BnnakKRmkDv6HsIzUHehCaa+zyOAlMpT/PM4GUk5Pz66+/XvSH999///mflUrlW2+95ZH3CjCKBEmA7wiKkbnGnr4IghD4B0m4iHVxToNbpL7KCAv/ZamxyaFd9AL+OuQpYMjjJeaaQJ7Aa21yCGQ0bHgG/EYQbBdqrrXDDMILQRDyTBohYqwsY2X5LkhvMVfZlNAuCvxKwI+XsdTa4YnwQhCEfCOQPE5sDtxVnUyVNgWMlAF+RRnQ42WcejfCKIAbfrsBgpB/8nhpAC9vaIInQuBv5PESe5OTC9B5TWZ4HLwEBCH/5AmSQB04ylhZxsJKI2GUNvAnpICURoksATqvyQIdhJeAIOSfPE4cqFMJzbV2eTyM0gb+R54gCdR2Gkv7WQkuAEHIP7FGyLo4t5nhuyCeZ4ZR2sA/KeKlgTqc21pvl8XCWfkXEIQ+gEDyWHFAtsNYamAZJ+CXArXDwql3I5IQKmE6019AEPoEWazEWh+At5+WOmiEAX4pUOc1WerssD/2pSAIfYI8VmypC7QnQhilDfwYgWSxYktdoN2eWusd8rirb6wRbCAIfYIsTmIJuCdCc41dkQgTJ4C/UiQEYDehpd4uhw7CS0AQ+gRphMhtZhhHQLXDWGphpAzwYwHZTWipc8jgifASEIS+gUCy6ECbRGGGkTLAnykSJJbAeiJ0WxjOxYk1Qr4L4nMgCH2FLMBmE2JkrXNAIwzwXyK1ABHIaXDzXRCPsdQ55HFimNd7KQhCXyGPlQRSz7ytxSlQULSM4rsgAHSfPD6gHgqt9XYZDBm9HAhCXyGLkwTSwFFLLYzSBn5PHi8JpAXxLfUOeQx0EF4GBKGvkEWJHDoX5w6QdX4tdXDvCfyePC6gJvjCE2FnIAh9BUERkjChrcnJd0E8A6YrgQAgjxUHzIqjrItzGRlJOIyUuQwIQh8iixFbGwKidRQjaz2MlAF+T6gSECQRGONlbA0OSZSIIGGozGVAEPqQgAlCe5uTllK0FEbKAL8ni5MExnBuawN0EHYKgtCHyGLE1saAOOXqYdIuCBDyQFlozdrgkEbDWXl5EIQ+JGCeCC0wgxAEClmgzGuyNjhk8ETYCQhCHyKQ0yQVCB0Sljo7jJQBgSFA9s3GyNoEQdgpCELfIosNhIdC2PkTBAxxSCDsm+3QuWgJRUug2/7yIAh9iyza74MQdv4EASUg9s22Njhk0EHYOQhC3xIA3YSw8ycIMAGwbzZ0EF4ZBKFvkfp/EMJUehBgAmDfbAjCK4Mg9C3SCJHT4PbrhdYs0EEIAksA7JttbYQgvBIIQt9CUIQk3L8XWrM2OOSxcMqBwCEJF7pNDOvw19tT1sm5TYw4DBZX6xQEoc/x6/EyjJ1lbKw4BE45EDgIkpBEiWxN/npW2hodkkhYXO1KIAh9jixa7L+nXMfgNDjjQGDx69tTaBe9KghCnyON/n/t3XlgE2XaAPA3dzJJ7yNJT6C0dAGhB5SjLEcRpGIRKofUsriwgoAXLOzK6gcorlwqIMgiR2E9q6C4IkXF5VTuytFSaGspNG3SI03S5j7n+2PcWNvSps0xk+T5/TWZvpl5kmbyZOZ95n25Wpm3XhqFPnngk/hRXI3XJkJdvRETcciOgtIgEVIOJuLovHbEUa3UwIcOQuBzvPq+Jq0MbiLsBiRCyuEEs2xW3KzxypEstHVwRgh8ED+Kq5MZEE52HL2iqzfCcNtdg0RIRZiI642Fo7gN1zcaMREccsDXMHkMJsYwNJvIDqTHTGoLwnF2AIz01BVIhFTEF3G8cT4mfZOJHcRkcOBDBXyQl3YT6mQw+1L34DuLijCxV54RQqUM8GH8aJ7OGxNhvRE6CLsFiZCK+GKvrJfRSg38KBhTBvgmL62X0coMUDLaLUiEVISJuVov7JmHM0LgwwTemQjh0qgjIBFSEZPHYHAZXjdDr7ZOD4kQ+CpuGNustVj0VrID6Qkc6RqMmBDOCLsBiZCi+GIvq5exaK02M84JZpEdCADuQUOYmOtdfRYGhYmJwXy83YNESFFed8hpiSswMLga8F18MVfrVVVscCu9gyARUhRf5I2HHFyBAb6ML+Z4V+Eo3ErvIEiEFIV5W+Gorh5upQc+DhNxvavDQlcPJaMOgURIUVgkR99kwm1eUzmqlUKlDPBx/CiuTmb0onJunczIh5+nDoBESFF0Np0dyPSaIZ2I4jT47Ql8GhNj0Nk0bynnxm24vtnEi4TJQbsHiZC6MBHHW8aXMShMTB4UpwHfx4/ymqujBrmJE8Sks+BLvnvwHlGXFw29DcVpwE94UTm3Vgbd9o6CREhdfBHHW6aq18mgOA34Bb73zJsNlTKOg0RIXd52RgiHHPB9XnRGCIXcjoNESF08IUff7B2FozCeIfATGHFUWr3hqIQzQodBIqQuOpPGCWLpm6heOGqz4AalGYuEQw74PjqTxglm6RupfqkGt+IGhZkXAUelQyARUppXFI7qG4zcUDaNAaOrAb/gFd2E+iYTJ4RFZ8JR6RBIhJSGeUO9jLYeOgiBH8HEXnBU6uoNfLgu6jBIhJTmFfUyUDIK/ArfGwZa08Iooz0BiZDSvOKMEPrkgV/BxF7QYQFHZY9AIqQ0LJJjUJgpXqIGd9MDv8INY5vVFqvRRnYgXYF7J3oEEiGl0Rg0biilS9SsRptFZ+WGwniGwF/Q6DReJEfXQN2j0mbBjUozLxyOSkdBIqQ6TEjpQ04nM/CEHJiPF/gVTETpWdL0jUZuGAsKuR0HiZDqeCIulROhth7meQF+h+JVbLoGIyaEo7IHIBFSHSakdM889MkDP8QXc3QN1D0jhDnRegoSIdVRvHAU+uSBH8JElL6nHn6e9hQkQqqjeOGoDobbBv6HE8yymWxWPUULR3X1cGm0ZyARUh2NQeOEsPRyKo44atHZbBacHcQiOxAAPIuGMCHH0EjFoxK34kalmRcBJaM9AInQC1D26qihwQSjVwD/hIm5+kYz2VF0Qtdo5IZCyWjPQCL0ApStlzE2mWE8Q+CfMBHH2ETFRKiHSpmeg0ToBTAhRe+g0DeaoVIG+CdMxDU0UDERQv1aL0Ai9AKUvjQKvz2BX+KLOXpK9hHqGgyYEI7KnnFZItTr9VevXq2pqXlQA7lcXlxcXFtb66o9+g9eJMfQTMXCUUOTGRIh8E8sAZNGo5nUFrIDaU9Xb4RE2FOuSYTFxcUJCQkvvfTS8OHDV69e3bHBk08+mZSUtHTp0tTU1BkzZphMVPwlRVm/TopNscJRU4uZzqCxBEyyAwGAHNxIFtUGWiMmpufCxPQ95JpEuHLlyhdffPHHH3+8du3a7t27y8rK2jVYsGBBfX39pUuX7t69W1paeuDAAZfs139QcKp6rczIFUKJNvBfnAgW1Y5KfaORGwoT0/eYCxJhQ0PDmTNnFi5ciBCKiop65JFHPv/883ZtJk+ezGazEUIBAQGDBg1qaGhwfr9+hYLdhLoGAyccTgeB/+IJWVSboVfXANdFe8MFX2QSiSQgICA8PJx42Ldv3y56Cn/55ZczZ86sW7fuQQ3MZvO5c+fsW4uKiho4cKDzQXo7TMhtLm0lO4rf0cmMXBGcEQL/xY1kt5RQ7KiEktFecTQRbtiwoaSkpN3KYcOGrVixQqfTEWd7BC6Xq9VqO92IQqF44oknVq5cmZKS8qAdabXaHTt2cDi//qgZOnTo+vXrH9RYr9ez2WwGg+Hgq/BigVaNVK/RaMiO4zdqqS40kUupkPyNVqul0eAiGGlsArNWqteoNdSZhqy1ThM8iO8nR6WDn38ul8tkdpPpHE2EmZmZCQkJ7VZGR0cjhIRCoUqlstlsdDodIdTc3CwSiTpuoaWlZcqUKQ8//PArr7zSxY6Cg4M///xzYsvdYjAYfpIIsT54pbKez+NTZcAIHBmbzMHxkQKBgOxQ/BeO4/D+kwjHcSavlWXhcEKoMsqgSS4LjQ/CBH5xUujCz7+jiXDs2LEP+lPfvn2Dg4MvXrw4evRohNBPP/20fPnydm20Wu20adNSUlLeeuutXsfqz+yFoxTpADAoTEyMweDAfajAr/HFHF29gSKJEEpGe80FX2RsNnvp0qXPPffct99+u3r16qamppkzZyKEzp8/HxcXR7SZOXNmTU1Nenr63r179+zZc/bsWef3628oVTgKXREAIIQwIYXmY4KS0V5zTdXfmjVrwsLCdu3aFR0dfebMGS6XixASCoVz5swhGowbN27o0KHV1dXEQ3stDHAcJiTmAg0kOxCEENLVw+xLACBMxGmp6rwkwvOgZLTXaDhOrfFKYmJiLl265GAfoR8VyyDU9HNLc2lr8p9iyQ4EIYTKP64NGSDgDWAEBASQHYv/UqvV8P6TSK1WIyWz6rA0ZUX7+glS1HzbiBCKmxJJdiAe4sLPP/TxeA1qXRqVwRTYACBMxNE1GnEbJU4ndA1wVPYSJEKvwYvkGJpNVBhxFLfh+iYTLxIOOeDvGGw6W8A0KigxDQWMMtprkAi9Bp1JlanqDXITO4jJYMOHBwCEiTlUGF+GKBmFn6e9A99l3oQiM/RqZUY+TEwPAEIIIUzEpcJRqYeJ6Z0AidCb/K9wlGS6egMGiRAAhBBCfBElzgihZNQZkAi9CUV+e+pkBj70yQOAEEIIE3OpMCA+3NrrDEiE3oQic1Bo4ZAD4H8oMm82lIw6AxKhN6HCIWez4EalmRcB804AgBBCdCaNG8rSNZJ8qUYrg0ujvQeJ0Jv8WjjaRGbhqK7ewItgQ588AHaYmEvuVPW//jyFktHegkToZYhBfkkMALoiAGiH9MEu9I1Gbhj8PO09SIReBhNxtaQecjDKKADt8MVccgtHYaQnJ0Ei9DKYiEPuRRitDM4IAfgd0su5tfVGKOR2BiRCL4OJSK7V1skMGJwRAtAGL5xt1lisRhtZAcCtvU6CROhleBFsY4vFZibnkLMYrBa9lRsCJaMAtEFDvAgyuwmh595JkAi9DI1O44WzdQ3kHHI6mRETcRB0yQPwexh5VWxWk82ktnDDWKTs3TdAIvQ+JHZI6OoNMMooAB3xyeuz0NUbsUgOjQ6/T3sPEqH34YtJq5eBShkAOoWJuVopaT9PoWTUSZAIvQ8m4pJ1aVQrNfCjIBEC0B4/iquV6knZNXQQOg8SoffByBrtHofblQDoHDuQiWg0k9ri+V3Drb3OY5IdgGusXLmypqaG7Cg8R1WpDbqBebhXwGbG1TW6oEq+fc3jjz/+1FNPeTIGACiLL+bopAb2AIGH9wsdFs7zkUT40UcfvfHGG0FBQWQH4keKiorOnz8PiRAAAibmaqWGYM8mQoveajVaOcFQMuoUH0mECKGcnByhUEh2FH5EKpVWVFSQHQUAVMEXc1uqtB7eqVZq4Iu5cEeTk6CPEAAAXICUEUd/TYTAOZAIAQDABTAxR99k8vB0oTqZAYNCbqdBIgQAABegs+icYE9PF6qVwRmhC0AiBAAA1+CLuVqpB6+O4sRNhHDvhLMgEQIAgGvwxR69x9fQbGLxmUwew2N79FWQCAEAwDUwMdeTwx9qpTAnmmtAIqSWS5cujR8/3v7wtdde++KLL7pof+DAga1bt7o9LACAA/hRHi0c1cpgyEPXgERILSaTSS6XE8u1tbX//ve/c3Jyumg/a9asd955R6lUeiQ6AEBXuKFsi95q0Vs9szsdVMq4CCRCtygpKTl9+vTly5fXrl1bUlKCECouLt68efP27dvr6uqINiqV6pNPPlm7du27774rlUo7bmTfvn0zZsxgs9kIoa+++koikRgMhn379iGEysvLv/vuO4SQQCCYNGnSBx984LnXBgB4EBrii7naOg+dFGqkMDG9a0AidIsTJ0785S9/Wbt2bXBwsMVi2bNnT35+PofDUSqVI0eOrK6uRggdP368uLg4KipKKpWmpaXJZLJ2G/nmm2+ysrKI5Q0bNpSWlmo0miVLliCEzp8/v2vXLuJPWVlZx44d8+CLAwA8kCCGp6nzxDQUVpPN3GrhRbA9sC+f5ztDrLXz7I/WKrWH7mwNZtMKsxiM349yRKfTv/nmGwaDodVqx48ff/369b59+yKEGAzGtm3btm/fPnfu3Llz5xKNW1paPv300xUrVrTdwq1btxITE7vd+4ABA4iTTgAA6fhRHhpoTScz8mA+Xhfx2US4+A/0Zk91WvOYiNHh05iens5gMBBCFRUVer1+0aJFxPqGhgaxWIwQun379vPPPy+RSCwWi1KpnD9/ftunm0wmg8HA5/Pbb7cDPp/f2trqkhcCAHASP5orPdvsgR3BrfQu5LOJMDWM5B9KXC7XvsDhcHbv3k2j0dr+acmSJbNnz166dClCaMWKFVbr7zrY2Wx2UFCQUqmMjo5GCHE4HKPxt8l4DQYDh/Nr2bRSqYyIiHD/CwIAdI8v5uqbTTYLTme69ytIW6fnR0MidA3oI3S7/v37R0ZGXrp0qd//EImwvr4+OTkZIaRWq7/66quOTxw2bJj9mufQoUOPHz9OLNtstuPHjw8dOpR4eOPGjYyMDE+8EgBAd2gMGi+c7YG7CTV1BkiErgKJ0O1YLFZhYeHatWvHjx8/ffr05OTkjz/+GCG0bNmyvLy82bNnjxkz5g9/+EPHJ86cOfPbb78llv/v//6vpKQkMzPTarUmJydbrdaXXnqJ+NPx48dnzpzpsZcDAOiaIJqncXfhKI50MkiELuOzl0bJtXTp0raXOocPH15WVlZZWalSqZKSksLDwxFCzz///LRp0yQSyZAhQ9hsts1mQwiNHDny9OnTxLPy8/O3bNnS3NwcFhYWGRl5/vz569evDxs27PTp01FRUUQbiURSXl4+ffp0T79CAMAD8KO52jo9QiHu24Wu0cgKYDK5MLiaa0AidAt7B6Edi8UaOHBgu5Xx8fHx8fHtmhFpEiEkEAg2btx44cKFxx57jFgTExNDo9HsWRAhdOHChW3bthH3GgIAqIAfzZVfb3HrLrR1ekE0z6278CuQCClt1qxZbR/yeLxVq1a1XTN79mzPRgQA6IYgmqeVGRCO3DdxvKbOIIiB66IuA32E3oTP57/55ptkRwEA6AqDS2cJmPomY/dNe0tbq+fDGaHrQCJ0i/fee4+4L8IdduzYsXnz5i4aFBUVPfvss27aOwCgW4JorlvrZTRSqJRxJUiEbjF27Ng5c+a4Y8sajWbTpk3PPPNMF22ys7PPnTtXVlbmjgAAAN0SxPI0te4aaM2gMNEZNHYAdGy5DCRCt7BarRaLBSEkl8s/+eSTysrKV1999c0331SpVM3NzRs3bly7du39+/eJxrW1te+9996qVavefvvtxsZG+0Zu3bq1bt26DRs2SCQSYqxthFBhYWFmZmZISAhC6IcffiCy3YEDB7RarVQqPXz4MEKIRqPl5eXt3r3bw68aAEAQxPI0Ne5KhNo6gyAGrou6EiRCtzh58uRHH32EEJJIJMuWLXvppZfi4+MvXbr0xBNPzJs3LzAwUKlUTpgwwWw2I4SOHTumUqlSU1MbGxvT09OJ8dJKS0vHjx/P5XJDQ0PnzJlDjLWNECoqKpowYQKxvGfPnjNnziCEXnjhBZVKVV5evn79euJPEyZMKCoq8vwLBwAghARxPE2dHre5ZbhjTZ2BD4nQpXz25FpZuM2iaPDMvmgsdvjCtYje+a8KjUZz8ODBiIiIadOmiUSi//73v1lZWTiOf/nll6WlpampqYsXLyZa5uXlVVRUFBUVPfnkk1u3bl20aNHLL7+MEBIKhfby0dLSUkd6HxMTE6uqqvR6PY8HBwwAnsbkMtgBTH2TCRO6fgZ5ba1eOMKNNyn6IZ9NhIJx021qD01XS+PyH5QFEUIikYgYC1QoFNJotMGDByOEaDSaUChsbm5GCJ09e/aFF14wGAwCgeD+/ftjxoxBCJWXl2dnZxNbSEtLs29Nq9ViGNZtSMRo3Wq1GhIhAKQQxGKaGr07EqGmVp/wRFT37YDDfDYRssR9kLgP2VEghBAxB0XHhzQaDcdxhNDTTz/9/vvvT5o0CSE0a9YsYkia4ODglpZf78lVqVT2p0dERCgUCmKZy+UaDL9VprU9/5PL5UwmMzQ01E0vCgDQtYA4nlqiixwe7NrNGlVmhCNOCMu1m/Vz0EdIPrlcTkwxIZFIvv/+e2Llo48+umfPHo1Gg+P49u3b7Y0zMjJu3LhBLKelpRUVFREdjQihr776KjU1lVi+efNmeno6k+mzP3QAoDg31cuo7+sD4ru/JgR6BBIh+VauXPnwww9PnTr18ccfHzFiBLHymWeeSUlJSUhI6N+/f2xsrP1Ub9asWfYqmCVLloSEhPTv31+v12dmZt69e/e1114j/nTs2DEYdAYAEgmiudp6o83i4noZTY1OEAf9Ha6GU0x0dHRtba2DjXU6ncViwXFcKBTW19e7M66eMRqNOp0Ox3GLxaJSqezrFQqFzWYjlltaWkwmE7FcU1Nz/fp1k8mk1Wr1er29vclkstlsX3/9dVpaGrHGZrOlpKSUlJTY2zQ3N2MYduPGDfualpaWfv36KRQKt70+HMfxbdu2LV261K27AF1rbW0lOwS/1u37//OWSrVE59qd3tx5V1mhdu02vZQLP/9w6cwt7KNgMxiMoKAg+3ri/j9CYGCgfTk2NjY2NhYhxGL9eunfZDItXLhwzJgxCoVi586dW7duJdbTaLRdu3bdunWLKLpBCIWGhtLp9LCwMPvWSkpKNm7c2HZfAADPI66OuvCeP9yGa2pduUFAgERIUSwWa8aMGWVlZRiGFRUV2afhRQiNGjVq1KhRbRsvX748ICDA/jAzM9NzgQIAHiAgjqeu0Ytct0FdvZETzGLyYPYlF4NESFE0Gi03Nzc3N9eRxq+//rq74wEA9JQgDpOea3bhBtU1euggdAcolgEAALfgizlGldmitXbf1DHq+zooGXUHSIQAAOAWNDotIA5rvadz1QY1NfoAOCOzwR7WAAAVHUlEQVR0A7g06i5Go3Hfvn3l5eVJSUmLFy+2V8EAAPxHYD+stVoXOiig+6bdsRpthmYTXwyzL7kenBG6BY7jU6dOvXz58rhx444fP56Xl0d2RAAAEgT25bfe1bpkU+oaPT+aS2O4bdp7P+azZ4Sv/7ilVi3zzL44TM72h/9Jp/32q+LHH3+sqan5/vvv6XT65MmTY2JiKioqkpKSPBMPAIAiAuN5WpnBZsHpTGcTWGuVNiiB75KoQDs+mwj/kjJPbdJ4Zl8cBqdtFkQI3bp1KzU1lU6nI4QCAgKSkpLu3LkDiRAAf0Nn0zEhR1OjC+znbA5rqdLGPhzhkqhAOz6bCKMELrx7p8d0Op1G81sa1mg0be/zAwD4j8C+/NZqZxOhzYJravUBfaBk1C2gj9BdTp06VV1djRD66aefiBl3yY4IAECCwL5Yy11nC0fV93WYiMPgwDe2W/jsGSHpRo8ePWvWLB6Pd/v27T179rQdUA0A4D8C+2GVhXW4DafRe99N2AIdhO4EidBdYmNjv/3226qqqvj4eC4XKp4B8FMsAZMVyNTJjPzo3n8PtFbposeHdd8O9AqcaLsRk8kcMGAAZEEA/FzIAIGyvPe1e7gVV9foAvpCB6G7QCJ0iz/+8Y+zZs0iOwoAACWEJAtU5epeP11do+dFcphcGGvbXeDSqFsMHz7cmafbbDar1eqSwWiMRiOHw3F+OwCAXgtM4N/5UGI12npX7dJSpQ1y+u4L0AU4I3SXixcvVlVVdd3m1KlTUqm04/rDhw8/+uijzsdgtVq5XG5ra2uvt6BUKouKipyPBAB/xmDTA2KxlqpeDjGjLFMHDxC4NiTQlmsSoc1mu3DhwrFjx1QqVRfNjEZjcXFxc7Mr5yWhrK1bt544caLrNmvWrLl69apn4umd6urqZcuWkR0FAF4vOFmgvNObbkKL1qqrNwb1hzNCN3LBpVGr1ZqTkyORSPr167dw4cITJ0489NBDnbZcs2bNW2+9VVBQMH/+fOf3S2Xnzp27efNmS0vLvXv3MjIycnNztVrtrl27bt++PWDAgGXLlgkEgqKionv37n3wwQfnz5+fNGnSxIkTO92UWq3euXNnZWXloEGDli5dyuP9Ovb8119/feLECYPBkJWVNXfuXI1G88EHH1y7do3FYk2dOnXq1KldhLd69er8/Pz9+/ebTKb58+fbL+ReuXKlsLBQr9dPmzZtypQpCKFdu3YplcqXX34ZIbRmzRoMg+56AHojZIDgzr9rEBL39ImK2+qgJL7zI7SBLrjgjPCbb7755ZdfLl++/J///GfRokVr167ttNnly5fPnTuXmprq/B6pTygUhoSExMXFpaenx8fH22y2iRMn3rhxIzc399atW+PHj7darTExMQKBICEhIT09XSzu/PCwWCxjxoz55ZdfcnNzL1++PHnyZBzHEUJr1qz5xz/+kZmZmZOTU1lZiRCSSqWNjY3Tp08fN27cihUrCgsLuwhv8+bNCxcuHD169ODBg7Ozs2/evIkQOnv2bHZ29oABA8aOHbt48eKCggKEUHJyMpvNTk9PT09PZzKhRxmAXuKLucT0ET19ouJWq0smrwBdcMFX25EjR2bMmEGcqeTl5T300EMWi6Xdl6bJZHr22WcLCgoWL17s/B4dUbb/vq7e6Jl9MXmMoS/1a3u3bFJSUmxsbFpaGlE7+v3339fX1//0008MBiM7OzsxMfG777579NFHw8PDMzMzp02b9qAtHz16lJjOiUajTZ48uU+fPmfOnElJSdm0aVNZWVlCQgJCiHh6UlLSunXr9Hp9Q0PDsmXLCgsLn3zyyS5ifuWVV3JychBCNTU177777r59+7Zs2bJq1apFixYhhLhc7sqVKxcsWJCVlfXee+9BBSwAzqKhkAEByjsacWao40/CrbiqQpuQG+W+uABySSKsra0dNmwYsRwXF2exWOrr62NiYtq2eeONN3JyclJSUrrdmtFo/Prrr0NDf/2sxMTEjBw58kGNrVar1dr57M+Js6OtJpujr8E5NAat6zEjKisrhw4dymAwEEIMBiM1NbWiosKRcpjKysq0tDQajYYQYrPZQ4YMKS8v5/P5wcHBRBa0k8lks2fPbmpqio2NVSgUxL66YP9fpKamvvvuuwihioqKF198kVg5bNiw6upqk6mbX69EdWu3rwK4SReff+ABPX3/gwfyZT8qIkcGOf6UlkotN4JNx2jwj+7IwfefTqcTX6FdcCgR3rhxY8WKFR3XFxQUxMfHm0wme6E/sWA0/u5U7ObNm0eOHLly5Yoj+zIYDEePHrX3hKWlpXVxNdVoNOI4zmAwiAuGbbECmNSZCTcoKKhtGZFSqQwODnbkicHBwe2eGBISEhoa2tra2u60e9OmTRkZGW+//TZC6NNPP33nnXe63rK9lFSlUhHBtN2XSqXi8/lsNrvrD5DNZmv3vwaeZDKZ4P0nUU/ff6wfS3fIoG7UsoMcPQORl7QEDuDBf7lTDr7/bDa7224dh/4fffr06bTnLyIiAiEkEonkcjmxpqmpCSHUrsfrn//8Z3R09Lp16xBCEonk0KFDAQEBubm5ne4rKCho79690dHRjgRGo9HYbDaDweg24XteeHi4RCIhlsePH79kyZKff/45LS3t+vXrly9fPnDgQLs2ncrKyvrrX/9aVlY2cODAixcv3r59+49//KNQKExOTn777bf//ve/I4Tkcnl4eLhWqw0PD0cI6fX63bt3dxve7t27d+zYYTQaDxw4QMwbPGXKlD179jz++ONMJnPnzp1EsUxYWFhzc7NOp+u0TIbJZEL5DImsViu8/yTqxfsf9lCQ9o4peIJjIw/jqLW89g8L4zAMRqfqhAs//w4VywQFBY3tDBFEZmbmyZMniZanTp1KTU1tF9yCBQtyc3P79evXr18/DocTERFBfGX7tkWLFh09elQkEi1fvjwmJmbv3r2PPfZYSkpKdnb2+++/Hx8fjxB68cUXd+zYIRQKN2/e3Pa5dDqd+AnTv3//nTt3ZmVlpaSk5ObmHjhwQCwW0+n0Tz/99PDhw/37909JSZk5cyaxqQ8//HDYsGGpqalDhw4ltkOj0ZhMZqe/EhgMxpAhQxISEuLi4p555hmE0KpVq4grrklJSWVlZVu3bkUIxcXF5efnDxgwIDQ0VKlUuvk9A8DHRaQFN11rcbBxS5WWwaHzxZAF3Q93mkKhEIvFK1as2Ldvn0gkKiwsJNaPHj16y5Yt7RpnZGQcPHiwi61FR0fX1tY6uGudTmexWHAcFwqF9fX1PY/do6xWq1QqtVqtvXuizWZrt16hULR91Waz+f79+8Tl4q7R6fSmpiaVSqVQKNr9Sa1Wd1zZqW3bti1dutSRlsBNWltbyQ7Br/Xm/bfhl9bd0dYbHGl750NJ3Vl5j3fhN1z4+XfB7RMhISEXL15ks9lXr14tKCiYM2cOsX7JkiVjx45t13jZsmV+OzMfnU4nzud698SOJ3YhISFCodD+kMlkxsXFsdlsBzcbFBQUEhLSbqVAIOi4EgDgGjQUkRIod+Ck0KK1Km+rI4c5VEwAnOSaO8Pi4uI2bNjQbmV+fn7Hln/6059cskfgjA0bNvD5MFAFACSISA++c1ASOzmi61LzxmJV6KAAJg8G2vYEGGvUH/3tb3+z1+UCADxJEMPjhLKafu7mpLD+okI0Eq7NeAgkQgAA8Ki4RyIlJxpxW/ubvuxaftHiNhTYFy7beAgkQgAA8KigBD476IEnhbgNv/sfWfyUSES5m8J8FiRCAADwtPhsYc33nZ8Uyn5UsPjM8JQeDEADnASJEAAAPC2wL8YLZ98/1tBuvUltkfzQlJDb40kqgDN8ZD4BDoczYsSIbgfY9Fu4FbdZcQaLbr/YgttwmxlnsOm9vvzS0tIyb948V0UIgL8ZkB97c8ddViAzetyvA4xYtNaKjyTCESG8SA65sfkbH0mExcXFzszD7g8arqhkPzbHTY5gB7K0MqPsp+bEOdH8qN4PWkGj0YKC4OoNAL3ExBiDFve5ueOuRW8LSsAQjio/qwsfGhQ/JZLs0PyOjyTC8PBwfxi2zRn9+iFNpr6ysA5n0iKiQwb9NTmwj7PD9KnVapfEBoB/4gSzBi/uI/tJIfm+ydhiTngiKnQgTD1IAh9JhMARghhe6sr+ZEcBAPgNL5LTbwb0CJIMimUAAAD4NUiEAAAA/BokQgAAAH7NuxPh119/XVZWRnYU/uvAgQONjY1kR+G/tm7dajabyY7CT1kslnfeeYfsKPyXXC7fv3+/q7bm3Ynw2LFjV65cITsK//Xpp59WVlaSHYX/2rFjh0qlIjsKP6VWq7dt20Z2FP6rqqrq448/dtXWvDsRAgAAAE6CRAgAAMCvQSIEAADg12g4/sA5sUghEAjCw8MdHDW0qamJx+MJBAJ3RwU6JZPJQkNDORwYF5EcNTU1MTExdDr8nCWBzWaTSCTx8fFkB+KnTCaTXC6PiorqtmVeXt769eu7bkO5RCiVSg0Gg4ONzWYzg8GALwKyGI1GyIIkgvefXPD+k8vB918sFvN4vK7bUC4RAgAAAJ4E51IAAAD8GiRCAAAAfg0SIQAAAL8GiRAAAIBf85r5CJuamvbv39/U1DR16tSsrKyODaxW6wcffFBSUpKcnPznP/+ZxWJ5Pkhf1dTUdPTo0du3bwcHB8+ePTsxMbFdA5vNtm/fPvvDwYMHjx492rMx+rLTp09XVFQQy0wmc8GCBR3bSKXSgoIClUo1ffr0MWPGeDZAH7d37962RYUdP94VFRWnT5+2P5w+fXpkJMwy75TGxsbi4mKJRDJ+/PikpCT7+tra2gMHDrS2tj7xxBMjR47s+ESTybR///7KysqUlJT8/HwH7ynwjjNCnU43atSo8vLy+Pj4vLy8zz77rGObZcuW7dq1KzEx8aOPPnr66ac9HqMve+6557777juxWNzY2JiSknL+/Pl2DSwWy+LFi+/cuXP37t27d+82NzeTEqevOnjw4GeffUa8t9XV1R0bqFSqjIyMurq6mJiYxx9//NixY54P0odVV1ff/Z/ly5eXlJS0a3DhwoUtW7bY2xiNRlLi9CUTJ058/fXXV69efeHCBftKuVw+fPjwpqYmsVicnZ194sSJjk+cO3fuZ599lpiYuH379hUrVji6P9wbFBQUDB8+3Gaz4Tj+ySefDBkypF0DmUzG4XAkEgmO4wqFgsvlVlVVkRCoj9Lr9fblxYsXL1iwoF0D4shXq9WejctfzJ8//5133umiwbZt27Kysojl3bt3Z2ZmeiQuv3P16lUej6dUKtutP3jwYE5ODikh+Sqr1Yrj+IgRIw4ePGhfuWnTpuzsbGJ5+/bt9s+83e3btzEMa2lpwXG8urqax+PJ5XJHducdZ4Rnz56dNGkSjUZDCE2aNOnmzZtKpbJtgwsXLvTv3z8mJgYhFBISkpaWdu7cOXJi9UVcLte+bDAYHjSUz/vvv//uu+/+/PPPnorLj1y8eHHLli2HDh3qdN6ls2fPTp48mVieNGnShQsXYHomd9i/f/+sWbOCg4M7/qm2tnbLli0FBQVNTU2eD8z3dHpJs93n/Ny5czabrV2DjIyMwMBAhFCfPn1iY2MvXbrk0O6cDtgTZDJZREQEsRwWFsZkMmUy2YMaIISEQqFUKvVoiP7hwoULR44cef755zv+aeLEic3Nzbdv387KytqyZYvnY/NhcXFxERERCoViw4YNGRkZOp2uXYO2n//IyEibzVZfX+/xMH2cXq8vLCzstIM2KCho8ODBLS0tR44cSU5O7njtFLhEu8+52WyWy+VtG9TX17dNBJGRkQ4mAu8olmEymRaLhVi2Wq1Wq5XNZrdrYLVa7Q/NZnO7BsB55eXls2bN2rt3b//+/dv9ic1m//DDD8Ty3LlzH3744aVLl/L5fI/H6Jtef/11+0JaWlpBQcFzzz3XtkHbA4RYgM+/y33xxRchISFjx47t+Kfp06dPnz6dWF6yZMlrr712+PBhz0bnF7r9nPc6EXjHGWF0dLQ9sRMLYrG4XYO6ujr7w7q6OkcGYwWOq6iomDhx4saNG2fPnt11y1GjRlksltraWs8E5ldYLFZGRkbHepm2B0hdXR2LxQoPD/d4dD6uoKBg4cKFRAdNF0aPHn337l3PhORv2n3OMQxrd5m614nAOxJhTk7O0aNHicG4Dx8+nJWVRZxtXL9+/f79+wihcePGyeXy4uJihFBlZeWdO3fsl5KB8+7fv//II4+sWbMmPz+/7forV64QHzu9Xm9f+c0332AY1qdPHw8H6cPsb69arT59+vSgQYMQQiaT6eTJk0SZUk5OzpEjR4h+wUOHDk2dOtXB+VuAg6qrq8+dOzdv3jz7Gp1Od/LkSeK8xP4PwnH82LFjgwcPJidKX5eTk/Pll18S53yHDh3Kyckh1l+5coVIkI888siNGzeIX4oXL17U6XSZmZkObdolFT7uZrFYJk+enJ6ePm/evLCwsPPnzxPrs7Ky1q9fTyxv27ZNLBYvWLAgNjbWvhK4xJQpU/h8fvr/PPvss8T61NTUnTt34ji+d+/eQYMGPfXUU1OmTAkMDPzwww9JjdfXhIaG5uTk5OXlicXixx57zGw24zhOHPn37t3DcdxgMIwZM2bUqFF5eXnh4eHXrl0jO2Rf88orr0ydOrXtmjt37iCEmpubcRyfMmXKxIkT8/PzH3roocTExPv375MUpu9Yvnx5eno6n8/v06dPeno68Z2v0+lGjBgxZsyYJ598MjIysrS0lGg8ZMiQ3bt3E8uvvvpqfHz8ggULRCLRv/71Lwd35zWzT1it1tOnT8vl8nHjxolEImJleXl5QECA/eS3tLSUuKE+NTWVvEh9UEVFhVqttj8MCAggbnG9detWREQE0WtdXFxcXV0dFBQ0bNgwuJvYte7du3f9+nWTyZSUlJSSkkKsNJvN165dS0lJIXpBzGbzqVOnVCrVhAkT2tYLAJeoqKgIDAy0f/MghAwGw82bN9PT0xkMhkKhuHz5slKpjI6OHjVqFIzm4byqqiqVSmV/mJiYSNSCEhdCWltbJ06cGBYWRvy1tLRUKBTaP/bFxcUVFRVDhw4dOHCgg7vzmkQIAAAAuIN39BECAAAAbgKJEAAAgF+DRAgAAMCvQSIEAADg1yARAgAA8GuQCAEAAPg1SIQAAAD8GiRCAAAAfg0SIQAAAL8GiRAAAIBfg0QIAADAr/0/dCp1HY3/cvAAAAAASUVORK5CYII=", "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "execution_count": 11 } ], "cell_type": "code", "source": [ "ρ = scfres.ρ[:, 1, 1, 1] # converged density, first spin component\n", "ψ_fourier = scfres.ψ[1][:, 1] # first k-point, all G components, first eigenvector\n", "ψ = ifft(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n", "ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\n", "\n", "using Plots\n", "x = a * vec(first.(DFTK.r_vectors(basis)))\n", "p = plot(x, real.(ψ), label=\"real(ψ)\")\n", "plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n", "plot!(p, x, ρ, label=\"ρ\")\n", "plot!(p, x, tot_local_pot, label=\"tot local pot\")" ], "metadata": {}, "execution_count": 11 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.5" }, "kernelspec": { "name": "julia-1.8", "display_name": "Julia 1.8.5", "language": "julia" } }, "nbformat": 4 }