{ "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\n", "--- --------------- --------- --------- ----\n", " 1 -0.143607104651 -0.42 7.0\n", " 2 -0.156040143100 -1.91 -1.10 1.0\n", " 3 -0.156772188827 -3.14 -1.57 1.0\n", " 4 -0.157037422201 -3.58 -2.24 1.0\n", " 5 -0.157005723168 + -4.50 -2.11 2.0\n", " 6 -0.157056388062 -4.30 -3.68 1.0\n", " 7 -0.157056405877 -7.75 -4.25 1.0\n", " 8 -0.157056406906 -8.99 -5.17 1.0\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380288 \n AtomicLocal -0.3163454\n LocalNonlinearity 0.1212602 \n\n total -0.157056406907" }, "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.05567723029220664, 0.0, 0.0]\n [0.055675847250452506, 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+gvaeTAAAgAElEQVR4nOzddXxUV/438HNl3OMe4kYI7hqcpkALlEILlCpLDUq93Zb9dbtbY8uWbt3dgUJwhwLBk0A8xH3c5crzR/JQCgRCMskd+b5f/SOkl5nDmbn3c+9RjGVZBAAAAPgrnOsCAAAAAFyCIAQAAODXIAgBAAD4NQhCAAAAfg2CEAAAgF+DIAQAAODXIAgBAAD4NQhCAAAAfg2CEAAAgF+DIAQAAODXPC4IX3jhBZvN1sWDKYqCJeI45HK5uC6CX4P65xbUP7fcWP8eF4RffvmlVqvt4sEul4thmF4tD7gOu93OdRH8GtQ/t6D+ueXG+ve4IAQAAAD6EgQhAAAAvwZBCAAAwK9BEAIAAPBrEIQAAAD8GgQhAAAAvwZBCAAAwK9BEAIAAPBrJNcFAL3LxaBSA1tmYKtNyMUgmkVSHgoVoQQ5lhWA8eBGCADu0Cwq0bOlBrbRinQOxMcRhqEYCUpWYKlKTAyX574CNe2bmm3oh0pmVwPzRzMbLsZSlFiCDPEJRGCoxYaOtqBiPVNpZAcGYnNi8flxWJwM47rIAPiLNjvaWM38WsUcb2XDxViGCosQI5UAWShEM+hUGyozMJUmdlAgNjUSvzMeS1LA6dm7IAh9ze4G9p0L9B8t7G2x+P0p+LcTcZXg2kdaKHS0hf2lihn5Oz0gAFvdn5gZjcEJB0DvOd7Kvn2e2d3AzIzCV6ThP07GlfxrH2ml0OFmdkc9M24rnaTAHk7DF8TjBJyfvQPztEWro6Ki8vLyIiMju3KwzWbj8/kEQfR2qbzC9jp27Rna7EJPDcAXxOOSLt/kuBj000XmP+cZF4PeGE7MiOrq2WYymWQyWTeLC3oM6p9bN1X/Z9TsUyfoahN6vD++PBmX8br6Li4Gbatj1hUyzTb0fBa+NAnHIQ4RQm79/kMQ+oISPbv6OF1tRq8OxefGdv882VzDPHuSiZGg/40hEuU3fhW4EHML6p9bXax/tR2tyaP3NLAvD8bvS+n+U92hZvaFU7SVQv8dSYwNgzB05/cfBkt4NxeD/u8sM34rNT0KL7idvL1fj+4W58TihbeTM6PxUb9T6woZ2rPukQDwPj9UMgN+c4UIUekC8sHUHrVtjg/DDuWQT2Xii/fTDx6hjbAHlPtAEHqx8zp2+GYqr5U5exu5qj/uliGgJI5W9cfz5pC5tUx2LtVggTAEoDtMLrR4P/3qOWbzVPLNEYS0y22h14EhdGcCfn4+iSE04FdqbyOcnu7hniD85ptvkpKSQkJC7rvvvutsq7tt27ahQ4du3brVLW/q5z4pZbJzqccy8NzpZKTEze0k8TJszyxyaiQ+dBO1sx5ONgBuTr6WHbqJkvLQiTnksGA3n55yHvpwLPHxOGLZQfrvp2louek5NwRhcXHxww8//MUXX5SWltbU1LzyyivXPMxgMDzzzDNarVaj0fT8Tf2ZlUJ3H6A3XGAO5ZDLk3vrmR7H0IuD8B+yyXsP0esKYfdjALrqt2pm2nZq7WD8o7GEqNcG5k+NxM7MJU+0spO3US2dPn2ALnHDZfSzzz6bO3fumDFjVCrViy+++Omnn15zAM6aNWtWrVoVHBzc83f0Z7VmdtxWisDQ8dlkqrLXO8wnhGPH5xDfVDDLD9FOSEMAbuSfZ5nVx5kdM8hFCb3e8RQiQttnkBPDsRGbqdNqeDDsPjd8VCUlJVlZWe0/Z2Vltba26nS6K47Zt29fZWXlvffe2/O382fHWtmRv1N3JeBfTujFO80rREuwI7eSegfK2UmZoH8egE7QLHrwCP17LZM3hxwU2EejOnEMrR1MvD0Sn7WT+uki3Kt2kxuuplqt9tIYVrlcjhBSq9UBAQGXDrBYLI899timTZuwLkzXbmtri4qKuvTHnJyc77//vrOD/Wr6xKY6/InT5AfDqWkRDrO5r9/9i5HoiVPkhC3Mr+NdwcKOe0+LxdKVzxT0Eqh/bl1e/zYKLT/GczBoywRKQtv7+AydGoQ2T8DuOMwv0zpXpVJ9+t7c6eL3XygUkuQNks4NQRgQEGA0Gtt/bv8hMDDw8gPWrl07bNgwg8Fw+vRpi8VSXV1dU1MTGxt7zVcLDg7u+jxCgiD8JAjfPs+8XcjsmkkMDOxkIYre9+kktPYMPesAvmcm0T48h2VZqVTKVXkA1D+3LtW/2YUWHqQiJNjn4wke3slKTr1spBQdm83m7MIaHbx3RxP+MOnejd9/NzSNJiUlXbhwof3n8+fPBwYGXv44iBBiGKawsPChhx566KGHqqurv/jii6+//rrn7+snWISeO0l/WsocnU0M7Kv2ls6sHUzcl4JPyKVrzNAhAQBCCBmcaPoOKlGBfTWB4HYV+0gJdiiHLDeyd+6jHTSXJfE6blhZprCwcNy4cQcOHEhOTl6wYEF6evqbb76JEHrllVeGDh06c+bMyw8eMWLEypUrly1b1tmrwcoyl6NZ9NARukjHbp1OBnBzo3kN7xYx6wqZ/bOIQGSGlU04BCvLcMtkMjEC2bTt1MgQbP0oT1kH1MmgJQdojZ3dNJV0y+RFj+VZK8tkZma++eabOTk5ERERcrn85Zdfbv99aWlpS0vLFQdHRkbCqdtFTgYt2kfXmdndszwoBRFCj6TjT2Xi2dvoequHnPsAcMBEYTN2UCNDsP96TAoihPg4+n4SkSjHpmyndA6uS+MlYK1RD2Wj0Py9FB/HfsgmBB7579twgXm7kDp0Ky/K3dP5QRfBEyGHLBSastUxOIR8d7QHpeAlLEJP59G7G9idM8lQEdel6R2e9UQI3M5CoVt3USoB9vNkD01BhNCjGfj9ifTU7XQrTOYFfsZOo7m7qRQ565kpiBDCEHpzBHF7HD4xl2q0etbTjgeCIPQ4FgrN3kWFibEvJxCkZ38+j6ZQC+OxqdspLbTAAL9Bs2jJAVrOw9YPdXlmCl7y0iD8vhR87Ba6ygRZeD2efaH1P0YXmradSlZgX00gPPRW86/WDiYmR2Czd1FWf5m8BPwai9A9B2k7zf6Q7R1n6JOZ+OMZ+ORtdDVkYecgCD2IwYmmb6cyVdh7Y7xpGtC6kUSqEpuzm4I12IDPezKPrjSyP2ST3M6UuCmP98efy8In5tKVRsjCa/OeD9PX6Z1o6nZqRAj2/livuNH8E4bQB2MIEYHddwjWwQe+7NVzzN4GdtsMUtJXCxy6ywOp+HMD8cnbIAuvDYLQI+gcaOo2alwYtn6kl6VgOxJHP2YT1Wb2qTyYxwt806elzGelzI6ZpJKzxZ165KFU/IWBePY2ugKy8CoQhNzTOtCU7dTEcGzdCE8dIdoFIhL9PpXcUc/Cnk3A92ypZf5+mt4xgwjz5qkID6TiLw3Cs3PpcgNk4V942xO+z9E40NRt1NRI7PXhXpyC7VQCtGMGMWYLHS5Gi3t/DxoA+saxVvb+w/TWaWSSwhvba/7ivhQcw9DkbfSeWUSy9/9z3AWCkEs6B5qxnZriEynYLkqCbZ9BZOdSQUJsWiScZsDrVRjZ+XvoLye4f6N5rtybjIsINCmX3jWTyFD5yD+qh+C2nTOtNjQxl5oehb3hKynYLl2J/TSZXHKAytdC8wvwbk1WNG07/epQfEaUTwXGogT8teH49B10kR5OUoQgCLnSbEOTcqnb+mH/HOpTKdhufBj2v9FEzk66FjapAF7L5EK37KTuS8HvSfbB6+SSRPyN4fjUbXQh3LBC0ygnGizs5G300iT8+YE+eIK1mx+HN1jQzB30kVtJlSetGA5AV7gYtGAvNTwYe8F3T9LFCTiJoek7qNzp5CCut3jjls9+xh6r2sROyKXvT/XlFGz3eH98ZjQ2dzdlhykVwKuwCD1wmObj2P/G+GCDzeXuiMf/N5qYuYPKa/Xr50IfvxZ7mnIDOzGXXt0ffzLTL2r+jeFEhARbcoBm/PosA17m76foEoPXLKLWQ7f1wz8bT87eTR1u9t+z1C8uxx6iUMtO2ka/NBh/ON1fqh3H0BfjCbWdfQIm2gMv8UEx83MVu3UaKfabjqNZ0dj3k8j5e6kd9X6ahf5yReZcXis7dTv1nxH4vb7Y8X4dAgJtnErua2Tfgon2wONtrmFeOctsn0EECbkuSt/KjsA2TyXvOUj9UuWP56l/XZS5sruBnb2b+nw8eUe8P1a4ko+2TyfevcB8V+mP5xjwFkdb2AeP0FumEfEyP2gSvcrIEGzXTHLVcebjEr87T/3xutzHNlYzSw9Qv04mZ0b749nVLlKCbZtBPHGc3umvbS/Aw1UY2QV76a8mkIOD/Pc8HRCAHc4h3ixknj3pX30ZEIS967/nmcePMbtnkWPD/PfsapeuxH6dQi49SJ1SQxYCz1JvYadso18fjk/3rYnz3RAnww7lkDvr2UeP0v6zmwwEYW9hEVp7hv6ohDl8K9Ef1jFCCCE0JhT7dBx5606qFNb8BR7D4EQ5O+lHM/C7E+F6iBBCYSJ0OIesMLLz9tB+suE2fPC9wkahhXvp/Y3skVvJWCmk4J9yYrBXhxEzd9CNVshCwD0rhW7ZSU2NxNb4x4ymLpLy0OZppJyHJm+jWmxcl6b3wWfvfs02NDGXEhBo10xYVOUa7k3GH0rFp2+ndQ6uiwL8W/vyMYly7A1v3gGtl/Bx9OVEYnoUNup36oLOx29bIQjd7KyGHbmZuiUG/2oiIYCTqxPPZOEzorBbdlIW/2h4AR6IYdE9B2kejn0yzh/mzXcHhtDawcQrQ/DsbVRunS9nIQShO/10kZmxg3prBP7SIBxOret7YwSRpsRu3005/Gt4GvAUjxylG63sD9kECVfB67orEd88lVxxhH4932dXiIKvgHtQDHrmBP3sSWb3THJ+HNTqjWEIfTSOkPOxxftpyu+mLQGOPXuSPqVmf59GCqHZpgtGhmB5c4jfqpk799FmF9el6QVwyXaDFhuatp0q0LIn55IDAuBRsKsIDH07ibDR7L2HYDFS0Hf+dY7JrWW3zyBlPK6L4j0ixNihHFLFR8M3U8U+t4shBGFP7W1kh2yiJoTjudPJQBgac5P4OPp1MllnYVf84T9zlgCX3j7PfFnO7JoJZ+tNExDog7HEM1n4xFzqizKfasaBIOw+F4NeOEUvO0h/NYF4eTAOvYLdIyLRlmlkkY59/Bj0FoLe9V4R8+4FZu8sIlzMdVG81rIkfN8scl0hc/cB2ugrzaQQhN1UomdH/07la9izt5HZEZCBPSLlodzpZF4r+8RxyELQWz4qYd4oYPbOIqIkcML2SIYKOzGHlPHQwN+oQz6xeRME4U2jWfT2eWb8VuqBVHzrdDLYz1ap7yUKPto5kzzczK6BDZtAL/i4hHn1HLN3FtHPLxfUdjsRid4fQ7w7mli8n34yj7Z5+TwoCMKbc0HHjt1Cbalhjs0mH0yF2nMnJR/tnkUeamLX5EF/IXCnj0uYf55j9s0iEuSQgu40Kxo7dxvZbEMDfqMONnnxWQuX8q6yUuj5k/SkXGp5Mr73FhLOqN6g5KNdM8kjzexjRyELgXu8W8T8Kx9SsLcECdE3E4m3RxJLD9LLD9Ftdq4L1C0QhF3ySxWT8StVa0EF83gPpsKwmF6kEqDds8izGvahIzCnAvTUW4XM+vPMgVsgBXtXTgx2YR4ZKED9f3X9r4jxupnBEIQ3cLKNnbCVevUc8/l44puJRJiI6wL5ATkP7ZhBVhjYJQdol7edUcBzvHya/rSUOXgLAQvf9wEpD701gtg3i9xcw2T9Rm3zqiXZIAg7VaRn5+2hb99DL03CT88lJ4bDudR3pDy0bQZpcqF5e2g7jJ4BN4lF6PFj9JZa9lAOGQljRPtQhgrbNZN8fTjx9Al6/FbqiJeMKYUgvIYCLXvnPjo7lxodipUtIO9LgTmCHBAS6NcphIyHZuyg9E6uSwO8h4tBSw/QZzTs/ltgUDc3cmKwgtvJ+1PwZQfpyduo/R4/jgaC8C/2NLCzdlIzd9DDgrHKhbw1mbiI5LpMfoyHo68nEgMDsQlbqSYr16UB3sBCodm7KKML7ZpBKvhcl8aP4RhamoSXLiCXJeF/O0KP/J366aLn9h1CECKEkN6JNlxg+v9KrT5Oz4/DLy4k12TiEohAD4BjaP1IYlECPmaLD65wCNyr1Yayc6lICfbrZAJuYT0BiaOlSXjRfPK5LPx/RUz8j9T/nWU8cFNuv/6yuBi0u4H9poLZXsfMiMbfHU1MCMegEdQDPZuFR0nQpFzqx2xyAnTWgmspNbC37KSXJuF/h03QPAyOoTmx+JxYPF/LfljMZP5KjwzB7k7E58TiYs+IIM8oRd+yUGhPA7Opht1Sw6QosbsT8XdH8wJgBV7PdnciHinBFu6j3hxBLEmElgzwFwea2EX7qNeGE8uS4LvhubICsPfGEG+NIDbVMF9XMCuO0FOj8NtisRnROLdroPtLEDpodLKNPdjMHmhi8lrZESHYrTH4K0NIWHXQi0wKx/bfQt66iy7Wsf8cSsBtP2j3aSnz4in6u2xyErQWeAMxiRYn4IsTcK0D/V7D/FLFPnzUlabEJkdg48Px0SGYtM+3x8JY1rOaa6OiovLy8iIjI7tysM1m4/P5BHHl3po0ixosbIURFenZ81r2tJot0rMZKmx8GDYxHJ8UgUH/n1uYTCaZTNbHb6pxoNt3UwEC7KuJhJ/vJ8dJ/XsUikFPnaC31bFbpxFJir5OQah/d3Ey6I8W9kAjc6CJPaNhY6XY0CCsfwDWX4UlylGsFONd6znfjfXv3YGwvxk72MZiGG1xIaMLGZ2o1c42WVGjlQ0WYolylKbEBgRgy5LxQYEYbEXtGwIFaPcs8rGj9KjfqU1TiURYMcRfaRxo4V6Kh6O8OaQSBoh6Mz6OJoVjk8IJhBDFoPM69rSavaBjd9UzlUbUYGWDhViEGIWKMAUfyXjo4XQ8Q+XOE9+7g1CAIyUP4TgWKUZyPlLwUYgQDxejCDEmgNjzXXwcfTCW+LCEGbuF+mQcmRMDWeh3zqjZBXvp+XHYv4YRBHz+PoTE0cBAbGDgnx8qxaBmG9tgQa12Vu9EFhdye9updwfh6BB2YhRGENA97o8eSsWzArCF++i8NmztYLga+pFPS5nnT9HvjSbmxcG57/tIHEVJsCgJQqi3TnL4GgEvNjIEOzWXPN7KTt1GeeDkJOB2ZhdaeoBef545lENCCgJ3gW8S8G7BQrRzBpkdgQ/dRO2ohyz0Zec07NBNFJ9AJ+aQKX0+NAb4MAhC4PVwDL04CP8xm1xxhF59HBbp9kEsQv8pZKbvoNYOxj8ZB6vGADeDIAQ+YlwYdvY2stGKhm+m8rXwaOg76izs9O3UxhombzZ5ZwJcsoD7wbcK+A6VAP2YTTw1AJ+2nXr1nOeu8Au67styZugmamI4fuAWsp8MmkNBr4AmBuBrliTik8Kx+w/Tm6qZT8cTAwLg6umV6i3siiN0gxXtnknChwh6FTwRAh8UJcG2zyBXpuNTt1N/Pw29hl6GYdH7xczgjdTIEPzEHEhB0OsgCIFvwhBanoyfu41XokeZv1J7GqDX0Dvka9kxW6jvKpkDOeSLg/Brrq0FgHtB0yjwZeFi9PNkIreOffAIPSwYWzcCh2XWPZbBiV46Tf9wkXl1KHFfCqypDvoO3G4B33dLNHZhHpmqQIM2Uv86x9gorgsE/oph0aelTNovLjuNiubz7ocUBH0LghD4BRGJ/jGEODGHPKdh036hvqtkoKnUQ+xtZIdsor4qZ7ZMIz8cS3C7Lx3wT9A0CvxInAz7aTJxpJldk0evK2ReH0ZMiYRnD86c07DPnqQrjei1YTislwY4BEEI/M7YMOz4HPLXKubho3SkGP1zKDE6FOKwT5Xo2bVnmINNzIuDiAdTYUQM4Bh8AYE/whCaH4dfmEfenYTfdYCeuYM62gJtpX2hRM8uPUBPyKUGBmIVC3kPp0MKAu7BdxD4LxJH9ybjpQvI2/vhdx+gJ2+jdsMsi15zRs3euY+emEulKrHyO3jPZuESaJACngG+icDf8XH0QCq+PBn/rpJ54jjNJ9CTmfj8OHhScQ8WoV317LpCukSPVmfin4zjuX1XVQB6CIIQAIQQInG0NAlfkoTn1rJvn6efPsE8nI7fl4IHC7kumdeyUOibCuad8wyfQKv744sS4N4CeCgIQgD+hCGUE4PlxJD5WnbDBSblZ1dONP5gKj42DEbT3IQLOvajEubbCmZCOP7uGGJSONQe8GgQhABcQ1YA9sk44s3hxBflzINHaAyhe1PwuxLxMBHXJfNgRhf66SLzeRlTY0b3JmNnbyejYR0f4A0gCAHolEqAVvfHV/fHjzSzn5Ux6b+4RodgixPxObEw0ONPLgbtamC/q2C21THZEfizWfisaJyABATeA85mAG5sbBg2NozYQBEbq5lvK5iH/6CnR+EL4rCZ0bjYX88hikH7m9ifq5hN1UyyAlucgL8zmgfrwgBv5K8nMQA3T0KiuxPxuxNxjQNtrGY+KmHuO0xPCsdnx2KzovFQ/2g1NTjRrgbm9xp2ex2TpMAWxOGn5pIxUngABF4MghCAmxYoQPen4Pen4DoH2lbH/F7LrslzJciw6VHY1Eh8dCjmY8MjGRad1bC7Gtjd9cwpNTs2DJsdg78+nIwQQ/4BX+DOIHQ4HALB9VpGXC4XjwdziIDvUAnQXYn4XYmIYog/Wtid9cyTeXSZgR0Vik0Mx8eGYsOCMQHBdSm7hWZRvoY90sIeaGIPNjGhImxqJLZmADEpHPPb1mDgq9zzja6vr1+0aFFBQQGfz3/77bfvvvvuKw54+eWXP//887a2NolE8sgjj6xdu9Yt7wuAhyBxNCEcmxBO/GsY0jnQoWbmQBO7+jhTrGezArERwdiwYGxIEJYoxzx5h6FaM3tazZ5Ss8db2VNtbJQEGxuGzY/D/jeaFy7munAA9Br3BOGqVav69+9/8ODB06dPZ2dnT5o0KTIy8vIDgoODd+/enZKSUlJSMn78+IyMjAULFrjlrQHwNCoBmhOLz4lFCCELhU61sXlt7K/V7AunGI2dzQzAMgOwTBWWqsTSlBiHAaN1oFIDe0HHXtCxhVr2nIbl4WhIEDYsGH8yEx8RggXAyBfgHzCW7eniinq9Pjg4uLy8vF+/fgihnJyc8ePHP/30050dP3/+/IyMjH/84x/X/L9RUVF5eXlX5GhnbDYbn88nCO9se/J+JpNJJpNxXQpvonOgAi1boGUv6NgSA1ukY200SpJjcTIsToZipFiMFEWIsQgxChbeuE21K/VPMajVzrbYUIMF1VnYWjNbbUZVJrbcwNIsSlFgGSosTYllBWKZKi5T2RvB959bbqx/NzwR1tTU8Pn89hRECKWnp1dWVnZ2sFarPXLkyIoVKzo7gGVZg8EgFneckUKhUCTyj9F4wA+oBO0tqH82jxqcqNzIVpnYahOqMLL7GlGTlWm0olYbKyBQkBALECAlH8n5mJhEEhKJSCT8/wHpdJJ8Pt3+s4tBZheyUchGI72TNTiR3onUdtbkQsFCFCrCoiQoSoJFSbBbY1CcDE+UY7B6HADt3BCEer1eIpFc+qNMJquqqrrmkS6Xa+nSpdOnT58yZUpnr6bRaEaNGoXjHaPupk6d+tlnn3V2MDwRcstisWCYB3d5eQMCoVQhShUiFHzl/zK6MI0D6Z1I78RMFLJRyEpjNgo5mI46J1mnAP3/nwm2nxAJcFZEIiUPyflIwWMDBCiA30mTD4XM5t76R/kJ+P5zq4v1LxQKSfIGSeeGIAwKCjIYDJf+qNPpQkNDrz6MpuklS5YghD7++OPrv1rXm0YJgoAg5BDLslKplOtS+CwpQhHXPcBkcslk8FjHGfj+c8uN9e+G6U79+vUjSfLChQvtfzx79mxaWtoVxzAMs3z5cq1W+8svv/D5/J6/KQAAAOAWbghCiUSyePHiF154oamp6fvvvz979uyiRYsQQgUFBTk5Oe3HrFy5cu/evY8++uiRI0f27NlTXFzc8/cFAAAAes490yfWrVu3evXq0aNHh4eHb968OSAgACHEsixFUe0HmM3m9PT0d955p/2Ps2bNuvqpEQAAAOh7bpg+4V4wfcKLwPBxbkH9cwvqn1turH/fWhIRAAAAuEmwaKDP0tn159uKz7eV1BjrmswtGpvOxVB2yi7nyxRCeYg4KEEVl6SKGxI2MFCk4rqwAPgXs9NypqWgVFNeqa9uNLcY7Eaj00jiPAHBVwjk4dLQaHlEelBKZnB6mCSE68L6PghCX6O2anZXHzxcd7zWWN8/ODUjKHV20oxwSWigOICP84Sk0OAwGhymZktLpa76j/oTG05/EiIOGh89ekZ8dqjkqrlsAAD3MTpMu6sP7q85fFFfMyAkPS0wJSdxepQsQiGQKQRyF0M5aafObmg2t1Qb647U5b135nMFXzYuetSUfuNjFdFcF99nQR+h7zjdnL+xLDe/5cL4mFETY8YMCh1A4jeuGYZlijVlu6sO7qs5nBKYuCjt9sFhA7r4jtBHwi2of27dVP1X6Kq+L/rteOOpURHDpsVPHBSSySNuvBUPi9hiddmhumO7qvZHySLmJs+aGDMGx6BLCyG3fv8hCH3B0YYTX5//yeqyLUidPbnfBBHZnUnWLtq1p+bQD0W/CUnhA1lLhoYPvOFfgQsxt6D+udXF+i/TVn6S/02lvnpB6uxbE6dLeN1Z0ZVi6D/q834p3aK16e7OmD8tfhKB+ft1D4KwAwRhkbr0vTOf2Sj7ssw7x0aNxHu84BOL2EO1xz7J/zpEEvzw4PvilbHXORguxNyC+ufWDeu/1ar+8OyXZ1sKlmXeOSthKg93Q1dUfuuFLwt/aLNqHhq0bGzUiJ6/oPeCIOzgz0GoteneO/N5fuv5+7LunhY3qecReDmapbdW7Pq84LtpcZOWD1jc2SMmXIi5BfXPrevUP83SP5f8/t2FX29LvuXO9Nu610hzHScaz7x/9nOlULF62IoYeZR7X9xbQNOzem4AACAASURBVBB28M8gZFh2U9m2Lwt/yEmatiTjDiHZW7vG6R2GD858caal8JmRjw4Jy7r6ALgQcwvqn1ud1X+lvvrfR9erRMrVw1ZESMN66d0ZltlYtu2rwh9zEqcuy7yTT/jd0pUQhB38MAhrjQ1vHH8Hx/AnRzwSI+9SLfXQ6eb8149vGBkxZOXg5cK/3tjChZhbUP/curr+GZb59sIvv5ZuWTFo+Yz47D4og9am23D6kwrdxadHPpYZ7F/LdUEQdnBXENIOxmlwuSw07WBoB40Qwnk4wccFKp5AxcNwj9hphUXsLyVbvjn/8z0D7pyTNMu9baHXZ3FZN5z6+IK65KWxTyWp4i/9Hi7E3IL659YV9d9iaXvlj3UCkv/cyMeDxIF9WZLDdcfXn/pwcuy4+7Pu9pxHQ4feZdc4GRfbflEl+DjOx3lSki8lSYkbnl4gCDt0IwhdZsra4rC1Oe1qh63Nadc4HToXS7N8BcmTkjgfJ4UEQoh2MoyTsWudLhMljhCqkqWqNJk8jrMNvFut6n8fW++iqedHr+q9xpbr21tzeMOpj5b0v2Neyq3tv4ELMbeg/rl1ef0frju27sR7C9NuW5h2W1/epF5idJjWnXivxlj/4ujViZfdrfYxa5NdW2TSlZpNtTZSgAuD+IQAJwQdF1XawbjMlMtEMTQrVPGEgXxhEF8ULBAF80XBAoGSh26m5iAIO1wnCBmKdRpcDp3LrnPaNS672mlTO+ytTkQgUbBAHCIQBfOFQQJhEF+o4pHiTqOUoVhTjVVfZtYUGDECixgXGDJUiRF9+kU/UPvH+pMfLEidsyj9dm6nEDWZW14+/Hq4NPTpkY9KeGK4EHML6p9b7fVPMfSHZ784VHfsH+OeSQ1M4rZIu6oOvHfm08UZ8xekzsZuKlV6iEWa88bGwxq72hmUJVcmS+XxEkLQ6cWKdjAOrdOmcdrVTlub09bmsLU6KDsjCuKLgvnCwI4rs0DF4yt5BP/arwNB2KH5nMZUZscwjGVY2s7QToay0ZSVdpkoxsXw5bz2ts2O+44gvihI0P1Hchbpy80NB9QOnSvxjsi+eTq0U44Npz8+13L+72PWcH6OtXPRrg2nPznTUvDK+OeCcBVciDkEQcgtk8nkJKmXD78u40ueG71KzveIz6LJ3PLKH+skPPHzo1ephMo+eEdbq6Pi50bayURlBwVmyrvdl0Q7GFubw6522tQdbXUOvcupdyGEeDKSJyFIMUEICEKAR4wLlEQKIQg7aCr1jmYKx3GEYaQIx3k4KSZ4YoInI0lRb42g0RQYL25qChwgj5sd1qvdh1X6mrVH3kgJTFo9bIXbh1/3UPuN50P9l81MmcJ1WfwXBCG3TtSceePMu3OSZtzdf0GfPn7dCM3SXxT8sO3inudHrbrmeG83ajykqdvTFjMtOHxMYC/VAe1gnEaKstKUjaYdNO1gVClSvoIHQdiBq1GjlJ0u/aoOYVjq0ujrPP73xJaKnZ+c+2bl4OXT+2TsWTdU6KpeOPBqdr9xDwxcAms+cQKCkEO/l+/4NP/bF0avHh4xmOuyXNvZlsJXj749PX7SvQMW98YyNCzDXtzYZKyypt8XK1DdeLk4t4Mg7MDh9AmWYSt/azLVWDP/FnedLsZusLis6/L+V2Osf3ns030zQaLbGjVNb517j8TJl8Y8KeVLuC6O34Eg5ATF0O+c+qig9cLzQ1clhyVyXZzr0TsMrx592045XhqzJlgc5MZXZmm2+ItalmZTl8X00sPADcF+hNzDcCxxfoQyWVr0WS1Due1mokxb+cD21VK+9P3pb3p4CiKEZHzpm9lrY+SRK3Y+WWts4Lo4APQ6vcPwxN4XNTbte9PfDJeEcl2cG1AKFG9MenlkxJAHd6w51nDKba/LovKfGhBC6ffHcpWC7uUL/wYOxeWECQN4JV/WsowbsvC30tyn9699IGvJE8P/5jmTga6PwIhHhtx/V/q8x3Y/d6LxDNfFAaAXVeqrV+x4Miuk/yvjnxfzRFwXp0swhN2VMf//xj379sn33z/zOcXQPX/Nmu0t1hZHypJoD5lj3XMQhD2DocQ7Imk7U7enrScvY3KaXzz07x1Ve9+b/uak2LHuKl2fmZkw5Z/jn3/9+Ds/FW/iuiwA9IrDdcfW7P37gwOX3pd1FyczBXsiMzjtk5nra40Nj+5+psnc0pOXUucb284ZMu6P7WxWgzfynX8JV3ASS10W3XxUa6qxde8VCtuK79u2KlwS8r9pb3A1Wb7n+genvj/jrd3VB/99bL2LdnFdHADchkXsl4U/bjj9yRuT1mbHjuO6ON0kF8j+NfGFybET/rbzyf01R7r3Ik4jVflbY8rd0TypT23qDkHoBjwpmTAvovSbOtrB3NRfZFjmy8IfXjr82uphKx4ecp9bdmnhUIg4aMPU1xy08/E9L2hsOq6LA4Ab2Cn72sNvnGg6/cH0t5IDErguTo9gCJufeusbk9Z+WvDt68ffsVP2m/v7LCr/sSFiTIAsxjuahbsOgtA9AjPl8jhxde5NtDk0W1of2/18fuuFT2auHxU5tPfK1peEpODlsU+Nihy6YseaEk15770RQ7HWFoe+zNz+n7XZcbN3IcAbMRRra3UYKiztn7ulwd6rn3uzpfXhXc+IeaL1k18NEKl67436UnJAwscz32ZZ9oHtq0s1FV3/i03HtJSVjpoS3Htl44p3P4J4lPjbwk+/Vh42SiUJv/Hk950X97135vNFGbffkTrX6/obrg9D2JL+d8Qr+z174P/+NsjN8yBtrQ71OYO60GhrcQhU/I7FCVnkNLrsOhcpImSxInk/sTJF2pVPAXgFu9apLzEbq6ymWqtD5xKoeHxFx1L4LhNlUzsFSjKwvzwwS+HeJ5UzzQWvHF13d8aCeSk5bnxZTyAihc+Oenx/zZFnDvxjfsrsxRnzbjgVmLLStTtaM1fG+cwAmcvBPEJ3ajys0RaZ+j/U7zrH6B2GdXnv1ZkaXxz9RKIqrq+K1iuuP4+nxlD3wqF/DQ8f/PCQe3s+n9eucdZsbzWUm4MGKYIGyGX9xFefkHat01RtM1ZbdMVmlmYDM+VBgxTyWLEnLfrhTr49j9Da7Gg7q9cUGF1WWpUqVSRIZLFicYjgyk+TRZZGu7rQ2HZaLwzk98sJlUa5IQ5/Lvn9u6JfXx7z5MDQzM6O8YH6b7Oq/33sv3bK8fzoVVGyiOsceXFTE0OxifOvd0wfgwn1HTwtCFmaPfNmRfycMFXatT+eg7VH/3vqw+lx2fcOWMwjOFiLwb1u+EU0Oy3/PLrO6rL9Y9wzPVn2sOGgum5PW8S4wMgJQV2ct9T+7Nh2zsA42ZBhypChSmGgd8xI6TofuBBfzWWmWk/rW0/qKRsdNFARnKWQRou6civDMmzLcV3trtagLEXc7LBur4xvpxxv5r1ba6z/5/jnQiUh1znSN+qfRezG0m1fFH6/LHPhbck512ygsrU5C96pHPxMkkeNkYEg7OBpQYgQ0l4wVec2D3oy8YrnFZ1dv/7kh1WG2mdGPpoRlMpV8dyrK19EhmW/LPwht3L32rFP9w++6X847WQqfmywaZxp98QIlN25dbA02FtO6trOGCQRwrBRqsD+8j7ePKT3+MaFuAOL9OXm5mNafZkloL8sdLhKES/pxqM8bWfKvqt3WenUZdF82U1ftRtMTX8/9O+kgPgnhq8U3Ggury/Vf72p8fXjG1iWeXrkY1cv5VH8ea0sVhSV7Vm9gxCEHTwwCBFCBRsuRowLDBqoaP8ji9htlXs+PvfVLQlTl2Xe6S0z5bui61/EYw2nXj/+zl0Z8+en3tr15YlpB3P+gypxqDBhfgRO9ii9GIrVnjc2H9NaWxyhw1VhowI4WR3RvXzjQkxZ6JaTuqajWlKAh40KCB6sJIQ9G8THoro9bS15usxH4m7q5ql9T8HlAxbPSZrZleN9o/4vYVh2c/m2zwu+n5dy6+KMeZcGsVsa7Rc+rhn6QnIPz0G3gyDs4JlBqC0y1e5oHfhEAkKoylD7nxPvU4xrzfCHvb1H8Go39UVstrS+dPi1MEnIMyMfk/BuvIkVQ7FFn9QIA3iJCyLd2Mlna3U0HdW2ntbL+4nDxwSqUqTe24Po7RdiU62t6Q+N9rwpoL8sfHSgLNadQ10aDqpbjusyH4nndWHnNRdDfXT2y5vdU9Db6/+aWixt/z31Yb2p6YlhK9r7R8u+rReHC6Oy3blUqVtAEHbwzCBELDrzZnlETsCvtk27qw8sH7D41sQZPjY0tN3NfhFdtOu9s58fbzj10tgn0wKTr3coi0q+qkUIS1kS1Ruj1Bgn03bW0PSHlrLRYaMDQoerunK59DReeiGmnUzbGUPzUS1lp8NHBYQOV3V/l9DrqtnWois1Zz4cd/01UJrMLf848magSPXMqMduak9BL63/rjhcd/zd059kBKU8EL+05n3tsBdSevqY3gsgCDt4ZhAyLLM79w/1GUPT9JoHBi5RCORcl6i3dO+LeLju2H9OfnBH6pyFabd1dn/QeEijPmfIfDiut/vzTLW25qNaTaFRlSYNGxXQvX4prnjdhdjSaG8+pm07a1DES8JGB/TB43jZd/U4D09c0OlYx701hzec+vjujAXzUnNudk9Br6v/m+Kgnd9e+MW005kQFDtlySihh22JiiAIL/HAIDxSn/fJua+DBIHz/rgj84EEaaTHfXvcqNtfxFar+tU//oPj+POjVl29O4yl0X7+g+qsVfHCgD7qT6VsdOspffMxLUuj0BGqkGHKboyz6HveciGm7LT6rKElT+c0UqEjVKEjVN0b99QNtIM5u64i7tawwMwr70ctLuv6kx+Waiv+PmZNkiq+Gy/uLfXfbZSFPvGv0iPZ+08YT9+TuXBmwpTe2New2yAIO3hUEJ5uzv80/1sH7Xhg4JKREUPr96ttrY6khZ6+lVJP9OSLyLDsd0W//lKyeeXg+6bFTfzz9xR77j+VUZODQoZ0f7pFt5lqrM3HdZoCozxOHDJMGZAh97QBApfz8Asxy7CGckvLKb2uyKRMloSO6ItHwKuZamzFn9UMXJPIl/95c3O2pfC1Y/8dETFk5eB7haSgm6/s2fXfcw0H1NYmR9KiyBJN+Ufnvmq1qpdnLpoUO85DOnogCDt4SBCebDr7ZeEPBofpnsw7L31LXGbq9L/Lh72U4hv7dV1Tz7+I5bqLrx59O1oW8cTwv7VPNKzZ0WprdaQujXZTGbuDcTLqAmPrab25zhaYKQ8epFAkSjxwQQ0PvRCzyFRrU5/Tt501CFS8kCHK4EHKXuoF7KLana3WJnvqPTEIITvl+Dj/q4O1x54e8UgPN5f30Pp3nzOvlycujJT36xjadro5/7OCb01Oy9L+d0yKHcv50yEEYQdug5Bm6YO1R78v+s3FUHdnzM+OHXfFMkUlX9QqU2VhI31kicKrueWL6KJdXxR+n1u5Z+XgeyeoxpxbVzFwTWKfNZ1dn9NItZ3Vq88Z7BpXYKY8KEuuSJB4zjREz7oQs8hUa1UXGDX5RpyHBQ1UBA9WioI9YrIQQ7FnXitPWhRZKal8M+/dzJD0R4bcd1PjYq7Js+rf3YxV1oqfGgY/c+UY2lNN574+/1OrVb0wbe6M+Mndfp7uOQjCDlwFodFhyq3cvbEsN1wSujD9tlGRQ6/Zza4rNtXuas163LtXrL8ON34Ry7SVrx9/Z0rJ9AFJqRm3dqfDplfZtU5NvlFdYLSpHaoUWWB/mTJFSop854642xgnY6i0aM6btEUmUoQHDVAEDpBLIjyua7z2ZHPRjqrPMz99YvjfRkQMcctrekL9957y7xvEEYLICdeeNVGkLv2u6LfC1qJZCVPmJs+8/hI8vQSCsEMfByGL2PzWC1vLdx1rPDkuetS8lJzr97GzDHvqn2UZD8SKfXQBaPdeCPTVlvxPy9/J/O9t6bMWps31zCXonEZKW2TSXjAaKi2ScKEqRapMlkpjRJw0nHJ2IWaRpdmuLzPrS83GKqs0WhSQIQvMkAuDPOL57wosYndU7v3o7NcrylZkZCdEjXTbJduHg5C2MydfKR3y3A3WVGsyt/xWlrvj4t70wJTZSTNGRAwh8b67O4Qg7NBnQVhvatxTfXDnxf0CUpCTMG16/CQZX9qVv1izo5W20/Fzw3u7hJxw74Wg8L2qkGFKNt317ulPqvS1jw69f2SE5+5OxVCs8aJFV2LWl5kdOpc8TqxIkMjjxNJoUZ+1nfblhZhlWGuzw3jRarhoMVRaCAGuTJaqkqXKZKkHzjC7pFRTsf7UhxjCVg17KNwUUfJl7dAXkt31AflwEDYf0+rLzKnLYrpysIN27q85srViV72pcWq/CVPiJqQEJPZ2CREE4SW9HYS1xoYj9cf31xzR2LQTY8ZOj590sx+wXevM/+/F4S+neOBQi55z4xfRVGMr/bpuyPNJ7RV1ounMhlMfh0lCVw5eHqeMdctb9B6XmTJctBorLYaLFrvaKYkQymLFshiRNEokDOT33jjJ3r4QO/QuS4PdVGs11djMtTaenGzPe0WixEM6ca+jzar+OP+bU01nHxy4bHr8pPbOiwsfVQcNVIQOd0+3vQ8HYcGGi9FTgjvbPKAzDaamnVX79lQfwhA2MXbM+KhRyYEJNzs7s+sgCDv0RhA6aWdhW3Fe4+ljDadslH1s1IgJMaOzQjJuuF9XZ/LXV8beEqpM6tITpHdx4xex+ItaZaIkfGzgpd9QDL25fNvX538aEzVieeaiIHHgdf6656AdjKnWZq61mmpt5nobbWMkkUJJhFAcLhSHCsShAlLstq+rey/EtIOxtTqszQ5rs93caLc02hGLpNEiWbRIGiOSxYq9ZfEdi8v6fdFvm8u3z06acVf6fDHvz5Xb9OWWi781Dn46yS0XZ18NQqfBdfatiuFrU7v96FyqrThYe/Rw3XE7ZR8RMWRExJAhYVmXfxBuAUHYwV1BaKccJZqy/Nai/NbzxZqyeGVs+4eXHOCG25mG/WqbxulR+3i5i7u+iLY2R8GGqmEvJuNXLYVldlq+K/p1S8XOWxKm3pl+m1Kg6Pnb9SWXhbY02iyNdmuzw9rssLU4EIFEwQJREF8YyBcG8gUqnkDF48t53Ziw2L36ZxnWZaLsWpdD57JrnXaN065x2lodlJ0RhwhEoQJJmKA9uT3/se8Kdsq+sWzbj8UbR0UOWz5gcchVazUghM69XRkzLSQgww3fW18NwsZDGkuT3S1zoOuMDccaT+U1ni5Sl/ZTxAwOzcwMSe8flCblS3r+4hCEHbodhHbKflFfW6mvKtNWlmjKa431Ccq4zOC0rND+A0P6u/fOxYdbR931Raz4uZEvJ2OmdzqKQW3TfnP+5701h25NnL4gdY5K6GVxeDmXibK1OWwap13jcmiddq3LoXM6jRQpIvgyktf+n5gg2/8TEoQAJ4Q4ISQwHLWPU8UIrH1yqtlslkqlCCHGxTIuBiFE2xmWYSkbTTsZ2sHQNoay0S4LRVlpl4V2GV1OE+Wy0DwJIQjgC5Q8YQBfGMgTBvFFQQKBkudFy8tdwU7ZN5fv+KF4Y1ZIxvLMRbGKTuehqvMNjYc0Ax51w8hkXw3Cgg0Xo6eGqFLd2YjlpJ0X1KVnWwoL24pKNOXB4qC0wKTkgIQEVVyiMq57uejG+veCdaR6yOAwtlnVrVZ1o6mlwdxUb2qsNdTrHYZYeXSCql+SKmFm/JQkVVzvjVEUBvCFKp6h0qpMcsNNkO9xmSl1vmHIc9dbgztIFLBq2EOLM+Z9e+GXpVtWTo2buDBtbqjEs3ZH66L2qJPHX/llcJkop5lymSiniaIsNGWl7WonZaNpB0M7GNpOswyibDRCiKVZ2sEghFiWxTAMIYTzMJyHI4QIIY5hGCkmcB5GCAlSiJMiQqDiSSNFPCnBk/P4MpInJXzpnszoNG0s3baxbOvA0Mz/ZP/fDXuUAzPl1VtbTDVWWeyNt0DxQ06Dy9bqcPvFik/wB4VmDgrNRAgxLFNlqC3RlJdqKvbVHK7S14pIYawiOlIWHiENi5CGhUqCQ8RBSqGyz5aw8e4nwlP150r05TiOW102mqHttN1OOawum8lpNjhMBodRb9eLeKJgUWCoJDhMGhIhDY+SRcQqosIkoX25SlD9frXdF1tH3XJH1nBQbW1yJN3Z1XYYrU33U8nm3Mrdw8IHLUyb2zfj0zyTrz6RdFGjufnnks17qg6NjR6xKH3e1dvJdsZdyx/6ZP27sV2061qt6lpDfb2pscHc3GRuabW0tVrVJqdJIZArBHK5QC7nS8U8kYAQSHhiDMNmJUyJkkXAE2EHJ+00Oy0Yhol5IoJHBBGBQlIg4YmlPIlCKFcK5Eqh8tL2khwKypLnr7+YcHu4L92Ju0tLni5xwU2cdQEi1YpB9yzpf8fWil1/P/RakChgXkrO+JjRnvBBgz7AsOyp5rO/lW4tVpfnJE79IufdQNHNjQINHao8/Vp5/NxwH17+sNvU+YboqX09Oz5EHBQiDhoaPvDyX9Isrbcb9Q6D0WEyOIw2yu6gHFbK1hsF8O5rx6CQzBFRQzhfa/SGhAF8gYpnrLIqEqB19C9MNVaWQZcWM+w6CU+8MG3ugtTZf9Sf2FiWu+H0JzPis29JmBrd5ccC4HU0Nt32i3tzK3bJ+NK5ybP+Me5ZAdGdKfw8GalIlKjPGUJH+Ozyh93jNFHWFve3i3YPgRGBItXN3uV0j3cHoRcJSJfpik0QhFdoydOFjlB1e4wGjuHjokeOix7ZYGraUrHzsT3PR0rDZsRPnhgzxi3D0oAncNLOow0nd1zce76tZGLMmLVjn04J7Gl7eOgIVf3eNgjCK+iKTMpkqeespttnIAj7iCpVVv5jQ78crsvhSWgnoy4wDn76ylV9uyFSFr5i0D0PDFyS13h6x8V97535bHBY1pR+40dGDOVwUWDQExRDn2nJ31d9+Eh9XnJAwoz47LVjn3HXp6lKlVb83GhtcYhD4evxJ12JKSDdZzcSvw4Iwj4iixG5zJRD5xKovGxuVu9R5xsU8ZLLd4nrIQIjRkcOHx053OKyHqo9mlu5+43jG4aFDxoXPXJkxFB4RvQKDtp5uvnc4brjf9SfiJZHTIwZ+8DApW5vH8NwLGSosvWErt+tYe59Ze/F0qy+zJJwu6+N6esKCMK+giFVqlRXbAobHcB1UTxF2xlDL+1RJeGJZyZMmZkwxegw/VGft6/m8H9OvJ8SmDgqctjIiKFdH14I+kybVd2+nNPZlsKUwMSxUSM6mxHvLiGDFRc+rumXE+a9syfdy1hlFQbzeTJ/DAV//DdzRZUqazurhyBsR9loc41NtbxLq/p2m1wga09EO+U43Zx/rOHkLyW/YwgbFj5oSFjW4LABCoE/NgR5CBtlP9dy/kxz/snmc1qbblj4oOzYcc+MeqznOwV2hThcSAhwU61NFuvmpb+8lK7EFHCTi4v6DAjCvqNKk1b+0shQbDcW0/I9mgKjMkVCXLWmWi8RkoIxUcPHRA1HCNUY6k42ndtVdeDNvHdDJcEDQ/sPCM7IDE7zluVMvZrRYTqvLi5oLcpvvVClr0kLSh4SmvXsyMeSAxL7cmpvu8AshbrAAEHYTltk7vp0Xh8DQdh3SBEhDhcYKy3KFB9cgPtmqQuMocOUnLx1rCI6VhE9P/VWhmXKtJX5rRd2Vx9Yf/JDASnoH5SSFpSSGpiUpIqHUTZuQTH0RX11iaa8SF1apClTWzVpQckDgtNXDLonLTCJ3635D+4SNEBe/FltHLSOIuTQu1xmShbtp/cEEIR9SpUm0xabIAgpG22qsqYu7XRByL6BY3hqYFJqYNLCtLkIoTpjQ7GmrEhdtq/m0EV9bbg0NEkVl6SKT1DFxSv7efUCp33J4rJe1FdX6qrLdRfLtRdrjPWR0rCUwMSM4NT5qbPjlbHd3sjF7SQRQozAzPU2qb8GwCW6YpMqVeq3NwQQhH1KmSSt+LmB61JwT3PeqEiSeNq6HtHyyGh55LS4SQghiqGrDDXl2ovluot/NJy8qKvGMTxOGROriI6VR8fII6PkESHi4L5vzfM0apu23tRYZ2yoNdRXGWprDHVml6WfIiZB2S85IGFWwpQEZZwnP1sHZcnV+UYIQn25xW87CBEEYR+TRgsdOpfLTPGkfl3zmgJj8ECPfsAicSJJFZ+k+nOPAo1NV22orTHU1xjrjtQfrzM26B3GcGlohDQ0QhoeLg0Nk4SESoJDJEFet1dUV1hc1lZLW5Oltdnc2mxpaTS3NJmbG0xNQlIQJYuMVURFyyOHhGf1U8SESoJ7by9WtwvMUpR8UdsvJ5TrgnCKRYYKS5wfzyTx68tx38NwTB4nNlRagrJ88FrZRYyLMVRaku+K4rogN6d9tachYVmXfuOgnY2mpkZzS6O5ucnccraloMWibrW22SlHsDgwSBQQLA4KECoDRKoAoVIpVCgFCqVQoRDIhKSQw3/INTlpp9Fh0juMOrte7zDo7AatTaex6TQ2rdqmbbW04RgeIg4KlYSESoLDpaHpQSntewVIeN69h4M0UohY5Ocz663NdkKI+/MUZwjCvqZIkhgq/DoIDRUWaZSIFHr6CrE3JCD4ccrYq/f9cdDONqtaY9O1WdVau15j1VYZanV2vc6mNziMRqeJYVkZXyLjS2V8qYQnkfBEYp5YypeISCGf4Ev5Eh7OE5ICPsETEAIcwy9tkCnj/6V32WK1mDDLX37jsjIsgxCyUXaaoV2My045nLTLQTtsLpuDdlpdNovLanXZLC6L2WU1O80mp8XkNNEMoxDIFAK5UqgIEKpUQkWASBWvjA0UBQSKAkLEQW7fXtxzKFOlumKTPwehvsKiTPLrgQsQhH1NkShtPlbHdSm4pC02q1J9uTdCQPCjZBFRsk5X6HDQTpPTbHKazU6z2Wm1uqxWymZymu2UXW83NJianIzLQTmctNNBOxmWzYalmgAAIABJREFUsbo6Vtw3Oc2Xvw7DMDj+l35WCU/cPhRFSApInLw8UEWkkE/yZQJpmDRERIqkfLGEJ5H+/zwWed5Dap8JSJM1HtJETuzFyfsezlBhCfLsroreBkHY16QRQspCOQ0uvsJPGyJ0Jaa0Xp5H7+EEBF8gCggS9XRpBZ/cD6/vKZIkpd/W0Q7G00Zv9REWGS9aEnxut9Sb4pcfPLcwJE+QGCotNz7SF9laHSzFSsL89/kDeBqCj8tixPpy840P9UXmBhtPSvL9cmW1SyAIOaBMlBgq/DQIdSVmVZrMewYVAr+gSpPqiv00CA3lFoV/dxAiCEJOKJKk+nI/DUJt+7xdADxJQJpMV2LiuhTcMFRaFIn+vjELBCEHxCEC2sE49C6uC9LXaCdjqrEqkyEIgWcRhQgwHLM2O7guSJ9jkbHaqoj37jkwPQdByAUMyWJFphor1+Xoa8ZKizRK5KdDEoBnU6VIdaV+91BobbbzJISfr++BIAi5Io+TGKv8LggNFRal3zfCAM/UPsGX61L0NWOVVR4HpyQEIUfkcWI/DEJ9BXTLAw+lSJAYL1pZhuW6IH3KWGWV9fP3dlEEQcgVabTI1uqgHQzXBek7lJ22tTpgdWPgmXhSkq/gWertXBekTxmrrfI4CEIIQo7gJCYJF5rrbFwXpO8YK62yWDFsSgw8ljJJoven1lGXiaKstDjEf9eWu8RtQWi1Wmtraxmm00ccmqZramrsdv+64boOmZ+1jurLzcok6I0AnkuRKDFU+NFsQkOVVR4nhkm9yF1B+M4770RGRmZnZ6emphYXF199wKlTp+Lj46dOnRoREfHVV1+55U29nbyf2FjlR7efhgqYrgQ8miJBYqy2srS/dBOaqizQLtrODUFYU1Pz/PPPHz16tKKiYvHixY8//vjVxzz00ENr1qwpKyvbtWvXypUrNRpNz9/X28njJaZqm590zlNW2qF1SaOggxB4LlJMiAL5/tNhASNlLnFDEP7www+TJk1KS0tDCK1cuXLv3r0tLS2XH1BcXFxUVHT//fcjhIYOHZqVlbVx48aev6+340kInoywtvjFHF5DhUUWJ8YIaIUBHk2RKPWTbkLGxVibHTIYvIYQcsvuE1VVVcnJye0/h4SEyOXympqa0NDQyw+IjIwUiztuPZKSkqqrqzt7NZqmCwoKmpub2/8YHh4eEdHpsuisw+pqraVxbx3yIwmmdGdreJS3lp+2Wp3iLt1Ras9R0iDMWVfe20XyK12vf9BFEiXTfJYOTdF35WCvrn9TIysMYKmWSq4LctMwDCfDYzHCnYsAuOG1TCZTUNCfW3lJJBKDwXDFASKR6DoHXM5gMDz11FM8XscWRUOGDHn77bc7O9j8x3Z04Wj3i845a39tVRBx/gDX5egmhmFsXbsL0WvuCJYd0pQ393aR/ErX6x90EcMKzK1L1N//D0M37rPw6vrXWwYQtELz/WGuC9Idwml3kQkDzOYujWwSCoWXAqUzbgjCkJAQnU536Y86ne7yx8H2A/R6/eUHpKend/ZqAQEBO3fujIyM7Mpbk9nz+DMWEYS37nUuq7WV/9wQvuYOrgvSTV3cD492MtUvlcQ8+3eYO+FesB9hb2h5vVxx1zpJ5I13CvPq+jd+UxeWKgsZuojrgvSIu+rfDbczWVlZJ06caP/5/PnzOI4nJCRcfkB6enpra2tDQ0P7H0+cOJGVldXz9/UBkkihvc1JO318Wr2pxiqJFEIKAq/gJ6s+mWpsshjoIOzghiBcsGBBbW3tW2+9de7cuVWrVi1fvlwikSCEnn322ddeew0hFBoaOn/+/EceeSQ/P/+ll15iWXbmzJk9f18fgBGYOExgafDxuZWmKli9AngNWT+xsdrHg5Cy0pSVFgXDVPoObghCiUSyZ8+e48ePr1ixYtiwYW+88Ub776OiosLDw9t/fv/992NjYx988MGKioqdO3eSpL8vdn6JNEZsqvXxs85YDaO0gdeQ9xObfD0ITbU2aZQIptJf4p5AyszM/OWXX6745SOPPHLpZ7lcvn79ere8l4+RxYh8fEdQFplqbcmLIQiBdxAFC2gn49C7BMobjLDwXuZaqxTaRS/jrUOefIY0WmSq9eUJvJZmO09CwoZnwGv4wXahpjobzCC8HAQhx8QhAspCUxaa64L0FlO1VQ7tosCr+Px4GXOdDZ4ILwdByDUMSaOEJt9d1clYZZXBSBngVeQ+PV7GoXMhFvlww283QBByTxot9uHlDY3wRAi8jTRaZGt2MD46r8kEj4NXgSDknjRG5KsDRykLTZlpcSiM0gbeBOfh4jCB2UfnNZmhg/AqEITck0YJfXUqoanOJo2GUdrA+0hjRL7aTmNuPyvBZSAIuSdU8Wkn4zJRXBfE/UwwSht4J1m02FeHc1sabJJIOCv/AoLQA2BIGin0yXYYcy0s4wS8kq92WDh0LoRjfDlMZ/oLCEKPIIkUWRp88PbTXA+NMMAr+eq8JnO9DfbHvhoEoUeQRgrN9b72RAijtIEXw5AkUmiu97XbU0uDXRp14401/A0EoUeQRInMPvdEaKq1yWJh4gTwVrIYH+wmNDfYpNBBeBUIQo8gDhG4TBRl96l2GHMdjJQBXswnuwnN9XYJPBFeBYLQM2BIEu5rkyhMMFIGeDNZjMjsW0+ELjPFOBmhis91QTwOBKGnkPjYbEIWWert0AgDvJdAyUMYcuhdXBfEbcz1dmmUEOb1Xg2C0FNII0W+1DNvbXXwZAQpIbguCADdJ432qYdCS4NNAkNGrwWC0FNIokS+NHDUXAejtIHXk0aLfGlBfHODXRoBHYTXAEHoKSRhArvWybh8ZJ1fcz3cewKvJ43yqQm+8ETYGQhCT4ERmCiIb212cF0Q94DpSsAHSCOFPrPiKO1knAZKFAwjZa4BgtCDSCKElkafaB1lkaUBRsoAr8dX8DAc843xMtZGuyhMgOEwVOYaIAg9iM8EoU3tIMUEKYaRMsDrSaJEvjGc29IIHYSdgiD0IJIIoaXJJ065Bpi0C3yE1FcWWrM02sXhcFZeGwShB/GZJ0IzzCAEvkLiK/OaLI12CTwRdgKC0IPwpCRO+EKHhLneBiNlgG/wkX2zWWRphiDsFAShZ5FE+sJDIez8CXyGMMAX9s22a52kiCBF0G1/bRCEnkUS7vVBCDt/Ap/iE/tmWxrtEugg7BwEoWfxgW5C2PkT+Bgf2DcbOgivD4LQs4i9PwhhKj3wMT6wbzYE4fVBEHoWcYjAoXd59UJrZuggBL7FB/bNtjRBEF4PBKFnwQhMFOzdC61ZGu3SSDjlgO8QBfNdRoq2e+vtKe1gXEZKGASLq3UKgtDjePV4GcpGU1ZaGACnHPAdGI6JwgTWZm89K61NdlEoLK52PRCEHkcSLvTeU65jcBqcccC3ePXtKbSL3hAEoccRhwstTd7aNAp98sAnSSKEZq8NQmuzQxwm4LoUHg2C0OOIwwRWr11x1NL4/9q788AmyrQB4G/us/eRpCdQWrqA0APKUZajCFKxHJWzwuLCClK8YGFXVz9AcQVFBQRZ5CqLB1VQXI6i4nIqp5WjpdDWUmjapkeapE0yuTPfH6Oxm0KbNsdMkuf312T6ZuZJmsmTmfeZ99ULoIMQ+Byvvq9JK4ObCLsAiZByOMEsqwU3abxyJAttHZwRAh8kiOJiMj3CyY6jR7AGAwy33TlIhFTEF3O9sXAUt+K6JgNfDIcc8DVMHoPJZ+hbjGQH0m1GtRnhODsARnrqDCRCKhKIOd44H5Ou2cgOYjI48KECPshLuwkxGcy+1DX4zqIivsQrzwihUgb4MEE0D/PGRNhggA7CLkEipCKBxCvrZbT1ekEUjCkDfJOX1stoZXooGe0SJEIq4ku4Wi/smYczQuDDhN6ZCOHSqCMgEVIRk8dgcBleN0Ovtk4HiRD4Km4Y26Q1m3UWsgPpDhxhjQa+CM4IuwCJkKIEEi+rlzFrLVYTzglmkR0IAO5BQ3wJ17v6LPQKI5MP8/F2DRIhRXndIaclrsDA4GrAdwkkXK1XVbHBrfQOgkRIUQKxNx5ycAUG+DKBhONdhaNwK72DIBFSFN/bCkexBriVHvg4vpjrXR0WWAOUjDoEEiFF8SM5umYjbvWaylFtPVTKAB8niOJiMoMXlXNjMoMAfp46ABIhRdHZdHYg02uGdCKK0+C3J/BpTD6DzqZ5Szk3bsV1LUZeJEwO2jVIhNTFF3O8ZXwZvcLI5EFxGvB9giivuTqqlxs5QUw6C77kuwbvEXV50dDbUJwG/IQXlXNrZdBt7yhIhNQlEHO8Zap6TAbFacAvCLxn3myolHEcJELq8rYzQjjkgO/zojNCKOR2HCRC6uKJOLoW7ygchfEMgZ/gE0elxRuOSjgjdBgkQuqiM2mcIJaumeqFo1Yzrlea+JFwyAHfR2fSOMEsXRPVL9XgFlyvMPEi4Kh0CCRCSvOKwlFdo4EbyqYxYHQ14Be8optQ12zkhLDoTDgqHQKJkNL43lAvo22ADkLgR/gSLzgqsQa9AK6LOgwSIaV5Rb0MlIwCvyLwhoHWtDDKaHdAIqQ0rzgjhD554Ff4Ei/osICjslsgEVIaP5KjV5goXqIGd9MDv8INY5vUZovBSnYgnYF7J7oFEiGl0Rg0biilS9QsBqsZs3BDYTxD4C9odBovkoM1UveotJpxg9LEC4ej0lGQCKmOL6L0IYfJ9DwRB+bjBX6FL6b0LGm6JgM3jAWF3I6DREh1PDGXyolQ2wDzvAC/Q/EqNqzRwBfBUdkNkAipji+idM889MkDPySQcLBG6p4Rwpxo3QWJkOooXjgKffLAD/HFlL6nHn6edhckQqqjeOEoBsNtA//DCWZZjVaLjqKFo1gDXBrtHkiEVEdj0DghLJ2ciiOOmjGr1Yyzg1hkBwKAZ9EQX8TRN1HxqMQtuEFp4kVAyWg3QCL0ApS9OqpvNMLoFcA/8SVcXZOJ7CgeAGsycEOhZLR7IBF6AcrWyxiaTTCeIfBPfDHH0EzFRKiDSpnug0ToBfgiit5BoWsyQaUM8E98MVffSMVECPVrPQCJ0AtQ+tIo/PYEfkkg4ego2UeINer5Ijgqu8dliVCn0/300081NTUPayCXy4uLi6VSqav26D94kRx9CxULR/XNJkiEwD+xhEwajWZsM5MdiD2swQCJsLtckwiLi4sTEhJeeumloUOHvvLKKx0bzJkzJykpKT8/Py0tbfr06UYjFX9JUdavk2JTrHDU2GqiM2gsIZPsQAAgBzeSRbVLNcTE9FyYmL6bXJMIV65c+eKLL/7www/Xrl3bsWNHWVmZXYOFCxc2NDRcvnz57t27paWlBQUFLtmv/6DgVPVamYErghJt4L84ESyqHZW6JgM3FCam7zYXJMLGxsazZ88uWrQIIRQVFfXYY4998cUXdm0mTpzIZrMRQgEBAQMGDGhsbHR+v36Fgt2EWKOeEw6ng8B/8UQsqs3QizXCddGecMEXmVQqDQgICA8PJx727t27k57CX3755ezZs2vXrn1YA5PJdP78edvWoqKi+vfv73yQ3o4v4raUtpEdxf/AZAauGM4Igf/iRrJbSyh2VELJaI84mgjXr19fUlJit3LIkCErVqzAMIw42yNwuVytVvvAjSgUiieffHLlypUpKSkP25FWq926dSuH8+uPmsGDB69bt+5hjXU6HZvNZjAYDr4KLxZo0dTrNBoN2XH8Tl2nDU3kUSokf6PVamk0uAhGGqvQpK3XadQa6kxD1lanCR4g8JOj0sHPP5fLZTK7yHSOJsLMzMyEhAS7ldHR0QghkUikUqmsViudTkcItbS0iMXijltobW2dNGnSo48++uqrr3ayo+Dg4C+++ILYcpcYDIafJEJ+L7xS2SDgCagyYASODHJzcHygUCgkOxT/heM4vP8kwnGcyWtjmTmcEKqMMmiUy0Ljg/hCvzgpdOHn39FEOHr06If9qXfv3sHBwZcuXRo5ciRC6Mcff1y+fLldG61WO2XKlJSUlHfffbfHsfozW+EoRToA9Aojk89gcOA+VODXBBIO1qCnSCKEktEec8EXGZvNzs/Pf+6557755ptXXnmlubl5xowZCKELFy7ExcURbWbMmFFTU5Oenr5r166dO3eeO3fO+f36G0oVjkJXBAAIIb6IQvMxQcloj7mm6m/16tVhYWHbt2+Pjo4+e/Ysl8tFCIlEotmzZxMNxowZM3jw4OrqauKhrRYGOI4vIuYCDSQ7EIQQwhpg9iUAEF/Maa16cEmE50HJaI/RcJxa45XExMRcvnzZwT5CPyqWQaj559aW0rbkP8WSHQhCCJV/WhvST8jrxwgICCA7Fv+lVqvh/SeRWq1GSmbVofqUFfb1E6So+aYJIRQ3KZLsQDzEhZ9/6OPxGtS6NCqDKbABQHwxB2sy4FZKnE5gjXBU9hAkQq/Bi+ToW4xUGHEUt+K6ZiMvEg454O8YbDpbyDQoKDENBYwy2mOQCL0GnUmVqer1ciM7iMlgw4cHAMSXcKgwvgxRMgo/T3sGvsu8CUVm6NXKDAKYmB4AhBBCfDGXCkelDiamdwIkQm/yW+EoybAGPR8SIQAIIYQEYkqcEWIwMb0TIBF6E4r89sRkegEccgAghBDiS7hUGBAfazDwRfDztIcgEXoTisxBoYW76QH4DTFvttVMchUblIw6AxKhN6HCVPVWM25QmngRMO8EAAghRGfSuKEsXTPJl2q0Mvh52nOQCL3Jr4WjzWQWjmINel4EG/rkAbDhS7gYqd2E8PPUSZAIvQwxyC+JAcAoowDYEZA92IWuycANg5+nPQeJ0MvwxVwtqYccjDIKgB2+hEtu4Sgmg6PSKZAIvQxfzCH3Igx0RQBghy/mYqTOQQH1a06CROhl+GKSa7UxmZ4Pvz0BaIcXzjZpzRaDlawAsAYoGXUKJEIvw4tgG1rNVhM5h5xZbzHrLNwQ6JMHoB0a4kWQ2U2IwXUa50Ai9DI0Oo0XzsYayTnkMJmBL+Yg6JIH4H/xyatisxitJo2ZG8YiZe++ARKh9yFxfBmsQQ+jjALQkUBMWr0M1mDgRXJodPh92nOQCL2PQEJavQxUygDwQHwJafUy0EHoPEiE3ocv5pJ1aVRbrxdEQSIEwJ4giqut15Gya+ggdB4kQu/DJ2+0e5iYHoAHYgcyEY1mbDN7ftdYI9xE6Cwm2QG4xsqVK2tqasiOwnNUldqgG3wP9wpYTbi6BguqFNjWTJ069amnnvJkDABQFtFnwQ4Ueni/0GHhPB9JhJ988smbb74ZFBREdiB+pKio6MKFC5AIASDwJVxtvT64n0cToVlnsRgsnGAoGXWKjyRChFBOTo5IJCI7Cj9SX19fUVFBdhQAUIVAwm2t0np4p9p6vUDChTuanAR9hAAA4AKCKBLuoPg1EQLnQCIEAAAX4Is5umajh6cLxWR6PhRyOw0SIQAAuACdRecEe3q6UK0MzghdABIhAAC4hkDC1dZ78OooTswPCvdOOAsSIQAAuIZA4tF7fPUtRpaAyeQxPLZHXwWJEAAAXIMfxfXk8IfaepgTzTUgEVLL5cuXx44da3v4+uuvf/nll520Lygo2LRpk9vDAgA4wMOXRrUyGPLQNSARUovRaJTL5cRybW3tv//975ycnE7az5w58/3331cqlR6JDgDQGW4o26y3mDGLZ3aHQaWMi0AidIuSkpIzZ85cuXJlzZo1JSUlCKHi4uJ33nlny5YtdXV1RBuVSvXZZ5+tWbPmgw8+qK+v77iR3bt3T58+nc1mI4S+/vprqVSq1+t3796NECovL//2228RQkKhcMKECfv37/fcawMAPAzNoyeFmno9HxKhK0AidIuTJ0/+5S9/WbNmTXBwsNls3rlz57x58zgcjlKpHD58eHV1NULoxIkTxcXFUVFR9fX1aWlpMpnMbiPHjh3LysoiltevX19aWqrRaJYuXYoQunDhwvbt24k/ZWVlHT9+3IMvDgDwUMIYnqbOE9NQWIxWU5uZF8H2wL58nu8MsWbn2R8sVWoP3dkazKYVZjEY/zvKEZ1OP3bsGIPB0Gq1Y8eOvX79eu/evRFCDAZj8+bNW7ZsmTt37ty5c4nGra2tBw4cWLFiRfst3Lp1KzExscu99+vXjzjpBACQThDloYHWMBnMx+syPpsIl/yB3uKpTmseEzE6fBrT09MZDAZCqKKiQqfTLV68mFjf2NgokUgQQrdv337++eelUqnZbFYqlQsWLGj/dKPRqNfrBQKB/XY7EAgEbW1tLnkhAAAnCaK5dWdbPLAjuJXehXw2EaaGkfxDicvl2hY4HM6OHTtoNFr7Py1dunTWrFn5+fkIoRUrVlgs/9PBzmazg4KClEpldHQ0QojD4RgMv0/Gq9frOZxfy6aVSmVERIT7XxAAoGsCCVevMFrNOJ3p3q8gbZ1OEA2J0DWgj9Dt+vbtGxkZefny5T6/IRJhQ0NDcnIyQkitVn/99dcdnzhkyBDbNc/BgwefOHGCWLZarSdOnBg8eDDx8MaNGxkZGZ54JQCArtAYNF442wN3E2rq9JAIXQUSoduxWKzCwsI1a9aMHTt22rRpycnJn376KUJo2bJleXl5s2bNGjVq1B/+8IeOT5wxY8Y333xDLP/f//1fSUlJZmamxWJJTk62WCwvvfQS8acTJ07MmDHDYy8HANA5YTRPU+fmRIgjTAaJ0GV89tIoufLz89tf6hw6dGhZWVllZaVKpUpKSgoPD0cIPf/881OmTJFKpYMGDWKz2VarFSE0fPjwM2fOEM+aN2/exo0bW1pawsLCIiMjL1y4cP369SFDhpw5cyYqKopoI5VKy8vLp02b5ulXCAB4CEE0V1unQyjEfbvAmgysACaTC4OruQYkQrewdRDasFis/v37262Mj4+Pj4+3a0akSYSQUCjcsGHDxYsXn3jiCWJNTEwMjUazZUGE0MWLFzdv3kzcawgAoAJBNFd+vdWtu9DW6YTRPLfuwq9AIqS0mTNntn/I4/FWrVrVfs2sWbM8GxEAoAvCaJ5Wpkc4ct/E8Zo6vTAGrou6DPQRehOBQPDWW2+RHQUAoDMMLp0lZOqaDV037SltrU4AZ4SuA4nQLT788EPivgh32Lp16zvvvNNJg6KiomeffdZNewcAdEkYzXVrvYymHiplXAkSoVuMHj169uzZ7tiyRqN5++23n3nmmU7aZGdnnz9/vqyszB0BAAC6JIzlaWrdNdCaXmGkM2jsAOjYchlIhG5hsVjMZjNCSC6Xf/bZZ5WVla+99tpbb72lUqlaWlo2bNiwZs2a+/fvE41ra2s//PDDVatWvffee01NTbaN3Lp1a+3atevXr5dKpcRY2wihwsLCzMzMkJAQhND3339PZLuCggKtVltfX3/o0CGEEI1Gy8vL27Fjh4dfNQCAIIzlaWrclQi1dXphDFwXdSVIhG5x6tSpTz75BCEklUqXLVv20ksvxcfHX758+cknn5w/f35gYKBSqRw3bpzJZEIIHT9+XKVSpaamNjU1paenE+OllZaWjh07lsvlhoaGzp49mxhrGyFUVFQ0btw4Ynnnzp1nz55FCL3wwgsqlaq8vHzdunXEn8aNG1dUVOT5Fw4AQAgJ43iaOh1udctwx5o6vQASoUv57Mm1snCzWdHomX3RWOzwRWsQ/cG/KjQazb59+yIiIqZMmSIWi//73/9mZWXhOP7VV1+VlpampqYuWbKEaJmXl1dRUVFUVDRnzpxNmzYtXrz45ZdfRgiJRCJb+WhpaakjvY+JiYlVVVU6nY7HgwMGAE9jchnsAKau2cgXuX4GeW2tTjTMjTcp+iGfTYTCMdOsag9NV0vjCh6WBRFCYrGYGAtUJBLRaLSBAwcihGg0mkgkamlpQQidO3fuhRde0Ov1QqHw/v37o0aNQgiVl5dnZ2cTW0hLS7NtTavV8vn8LkMiRutWq9WQCAEghTCWr6nRuSMRamp1CU9Gdd0OOMxnEyFL0gtJepEdBUIIEXNQdHxIo9FwHEcIPf300x999NGECRMQQjNnziSGpAkODm5t/fWeXJVKZXt6RESEQqEglrlcrl7/e2Va+/M/uVzOZDJDQ0Pd9KIAAJ0LiOOppVjk0GDXbtagMiEccUJYrt2sn4M+QvLJ5XJiigmpVPrdd98RKx9//PGdO3dqNBocx7ds2WJrnJGRcePGDWI5LS2tqKiI6GhECH399depqanE8s2bN9PT05lMn/2hAwDFualeRn1fFxDf9TUh0C2QCMm3cuXKRx99dPLkyVOnTh02bBix8plnnklJSUlISOjbt29sbKztVG/mzJm2KpilS5eGhIT07dtXp9NlZmbevXv39ddfJ/50/PhxGHQGABIJo7naBoPV7OJ6GU0NJoyD/g5XwykmOjq6trbWwcYYhpnNZhzHRSJRQ0ODO+PqHoPBgGEYjuNms1mlUtnWKxQKq9VKLLe2thqNRmK5pqbm+vXrRqNRq9XqdDpbe6PRaLVajxw5kpaWRqyxWq0pKSklJSW2Ni0tLXw+/8aNG7Y1ra2tffr0USgUbnt9OI7jmzdvzs/Pd+suQOfa2trIDsGvdfn+/7yxUi3FXLvTm9vuKivUrt2ml3Lh5x8unbmFbRRsBoMRFBRkW0/c/0cIDAy0LcfGxsbGxiKEWKxfL/0bjcZFixaNGjVKoVBs27Zt06ZNxHoajbZ9+/Zbt24RRTcIodDQUDqdHhYWZttaSUnJhg0b2u8LAOB5xNVRF97zh1txTa0rNwgIkAgpisViTZ8+vaysjM/nFxUV2abhRQiNGDFixIgR7RsvX748ICDA9jAzM9NzgQIAHiIgjqeu0Yldt0GswcAJZjF5MPuSi0EipCgajZabm5ubm+tI4zfeeMPd8QAAuksYx68/1+LCDaprdNBB6A5QLAMAAG4hkHCMrWaz1tJ1U8eo72NQMuoOkAgBAMAtaHSaMJ7Xdg9z1QY1NboAOCN0A7g06i4OMjmwAAAVD0lEQVQGg2H37t3l5eVJSUlLliyxVcEAAPxHYG9+WzUWOiCg66ZdsRis+hajQAKzL7kenBG6BY7jkydPvnLlypgxY06cOJGXl0d2RAAAEgT2FrTd1bpkU+oanSCaS2O4bdp7P+azZ4Rv/LCxVi3zzL44TM6WR/9Jp/3+q+KHH36oqan57rvv6HT6xIkTY2JiKioqkpKSPBMPAIAiAuN5WpneasbpTGcTWFuVNihB4JKogB2fTYR/SZmvNmo8sy8Og9M+CyKEbt26lZqaSqfTEUIBAQFJSUl37tyBRAiAv6Gz6XwRR1ODBfZxNoe1VmljH41wSVTAjs8mwiihC+/e6TYMwzSa39OwRqNpf58fAMB/BPYWtN51NhFazbimVhfQC0pG3QL6CN3l9OnT1dXVCKEff/yRmHGX7IgAACQg6mWc3Ij6PsYXcxkc+MZ2C589IyTdyJEjZ86cyePxbt++vXPnzvYDqgEA/EdgH35lYR1uxWn0nncTtkIHoTtBInSX2NjYb775pqqqKj4+nsuFimcA/BRLyGQFMjGZQRDd8++BtiosemxY1+1Aj8CJthsxmcx+/fpBFgTAz4X0EyrLe167h1twdQ0W0Bs6CN0FEqFb/PGPf5w5cybZUQAAKCEkWagqV/f46eoaHS+Sw+TCWNvuApdG3WLo0KHOPN1qtVosFpcMRmMwGDgcjvPbAQD0WGCC4M7HUovB2rNql9YqbZDTd1+ATsAZobtcunSpqqqq8zanT5+ur6/vuP7QoUOPP/648zFYLBYul9vW1tbjLSiVyqKiIucjAcCfMdj0gFh+a1UPh5hRlqlDkoWuDQm055pEaLVaL168ePz4cZVK1Ukzg8FQXFzc0uLKeUkoa9OmTSdPnuy8zerVq3/66SfPxNMz1dXVy5YtIzsKALxecLJQeacn3YRmrQVrMARCyag7ueDSqMViycnJkUqlffr0WbRo0cmTJx955JEHtly9evW77767d+/eBQsWOL9fKjt//vzNmzdbW1vv3buXkZGRm5ur1Wq3b99++/btfv36LVu2TCgUFhUV3bt3b//+/RcuXJgwYcL48eMfuCm1Wr1t27bKysoBAwbk5+fzeL+OPX/kyJGTJ0/q9fqsrKy5c+dqNJr9+/dfu3aNxWJNnjx58uTJnYT3yiuvzJs3b8+ePUajccGCBbYLuVevXi0sLNTpdFOmTJk0aRJCaPv27Uql8uWXX0YIrV69ms+H7noAeiKkn/DOv2sQknT3iYoydVCSwPkR2kAnXHBGeOzYsV9++eXKlSv/+c9/Fi9evGbNmgc2u3Llyvnz51NTU53fI/WJRKKQkJC4uLj09PT4+Hir1Tp+/PgbN27k5ubeunVr7NixFoslJiZGKBQmJCSkp6dLJA8+PMxm86hRo3755Zfc3NwrV65MnDgRx3GE0OrVq//xj39kZmbm5ORUVlYihOrr65uamqZNmzZmzJgVK1YUFhZ2Et4777yzaNGikSNHDhw4MDs7++bNmwihc+fOZWdn9+vXb/To0UuWLNm7dy9CKDk5mc1mp6enp6enM5nQowxADwkkXGL6iO4+UVHW5pLJK0AnXPDVdvjw4enTpxNnKnl5eY888ojZbLb70jQajc8+++zevXuXLFni/B4dUbbnPtZg8My+mDzG4Jf6tL9bNikpKTY2Ni0tjagd/e677xoaGn788UcGg5GdnZ2YmPjtt98+/vjj4eHhmZmZU6ZMediWjx49SkznRKPRJk6c2KtXr7Nnz6akpLz99ttlZWUJCQkIIeLpSUlJa9eu1el0jY2Ny5YtKywsnDNnTicxv/rqqzk5OQihmpqaDz74YPfu3Rs3bly1atXixYsRQlwud+XKlQsXLszKyvrwww+hAhYAZ9FQSL8A5R2NJDPU8SfhFlxVoU3IjXJfXAC5JBHW1tYOGTKEWI6LizObzQ0NDTExMe3bvPnmmzk5OSkpKV1uzWAwHDlyJDT0189KTEzM8OHDH9bYYrFYLA+e/TlxVrTFaHX0NTiHxqB1PmZEZWXl4MGDGQwGQojBYKSmplZUVDhSDlNZWZmWlkaj0RBCbDZ70KBB5eXlAoEgODiYyII2Mpls1qxZzc3NsbGxCoWC2FcnbP+L1NTUDz74ACFUUVHx4osvEiuHDBlSXV1tNHbx65Wobu3yVQA36eTzDzygu+9/cH+B7AdF5PAgx5/SWqnlRrDpfBr8ozty8P2n0+nEV2gnHEqE169fX7FiRcf1BQUF8fHxRqPRVuhPLBgM/3MqdvPmzcOHD1+9etWRfen1+qNHj9p6wtLS0jq5mmowGHAcZzAYxAXD9lgBTOrMhBsUFNS+jEipVAYHBzvyxODgYLsnhoSEhIaGtrW12Z12v/322xkZGe+99x5C6MCBA++//37nW7aVkqpUKiKY9vtSqVQCgYDNZnf+AbJarXb/a+BJRqMR3n8Sdff95/dhYQf16iYtO8jRMxB5SWtgPx78lx/IwfefzWZ32a3j0P+jd+/ea9eu7bg+IiICISQWi+VyObGmubkZIWTX4/XPf/4zOjqa2IJUKj148GBAQEBubu4D9xUUFLRr167o6GhHAqPRaGw2m8FgdJnwPS88PFwqlRLLY8eOXbp06c8//5yWlnb9+vUrV64UFBTYtXmgrKysv/71r2VlZf3797906dLt27f/+Mc/ikSi5OTk99577+9//ztCSC6Xh4eHa7Xa8PBwhJBOp9uxY0eX4e3YsWPr1q0Gg6GgoICYN3jSpEk7d+6cOnUqk8nctm0bUSwTFhbW0tKCYdgDy2SYTCaUz5DIYrHA+0+iHrz/YY8Eae8Yg8c5NvIwjtrKa/+wKI7Ph9GpHsCFn3+HimWCgoJGPwgRRGZm5qlTp4iWp0+fTk1NtQtu4cKFubm5ffr06dOnD4fDiYiIIL6yfdvixYuPHj0qFouXL18eExOza9euJ554IiUlJTs7+6OPPoqPj0cIvfjii1u3bhWJRO+8807759LpdOInTN++fbdt25aVlZWSkpKbm1tQUCCRSOh0+oEDBw4dOtS3b9+UlJQZM2YQm/r444+HDBmSmpo6ePBgYjs0Go3JZD7wVwKDwRg0aFBCQkJcXNwzzzyDEFq1ahVxxTUpKamsrGzTpk0Iobi4uHnz5vXr1y80NFSpVLr5PQPAx0WkBTdfa3WwcWuVlsGhCySQBd0Pd5pCoZBIJMuXL9+1a5dYLC4sLCTWjxw5cuPGjXaNMzIy9u3b18nWoqOja2trHdw1hmFmsxnHcZFI1NDQ0P3YPcpisdTX11sslp490Wq12q1XKBTtX7XJZLp//z5xubhzdDq9ublZpVIpFAq7P6nV6o4rH2jz5s35+fmOtARu0tbWRnYIfq0n778Vv7z2jrZB70jbOx/X1J2Td3sXfsOFn38X3D4REhJy6dIlLpd77dq1goKC2bNnE+vz8/PHjBlj13jZsmV+OzMfnU4nzud69sSOJ3YhISEikcj2kMlkxsXFsdlsBzcbFBQUEhJit1IoFHZcCQBwDRqKSAmUO3BSaNZalLc1kUMcKiYATnLNnWFxcXFvvfWW3cqnnnqqY8s//elPLtkjcMb69esFAhioAgASRKQH39knjZ0Y0XmpeVOxKnRAAJMHA217Aow16o/+9re/2epyAQCeJIzhcUJZzT93cVLYcEkhHtGNOw6BMyARAgCAR8U9Fik92YRb7W/6slFVahGOAmECQk+BRAgAAB4VlCBgBz30pBC34tX/kcU9FunhqPwZJEIAAPC0+GxRzXcPPimU/aBgCZnhKd0YgAY4CRIhAAB4WmBvPi+cff94o916Y5tZ+n1zQm63J6kAzvCR+QQ4HM6wYcO6HGDTb+EW3GrBGSw6+q1ODbfiVhPOYP++prtaW1vnz5/vqggB8Df95sXe3HqXFciMHvPrACNmraXi01rRsBBeJIfc2PyNjyTC4uJiZ+Zh9weNV1WyH1riJkawA1lamUH2Y0vi7GhBVM8HraDRaEFBcPUGgB5i8hkDlvS6ufWuGbME9RUgK6r8oi48JSh+EvQOepqPJMLw8HB/GLbNGX36IE2mrrKwDmfSIqJDBvw1ObCXszVparXaJbEB4J84wayBS3rJflRIv2s2tJoSnowK7Q9TD5LARxIhcIQwhpe6si/ZUQAAfseL5PSZDj2CJINiGQAAAH4NEiEAAAC/BokQAACAX/PuRHjkyJGysjKyo/BfBQUFTU1NZEfhvzZt2mQymciOwk+Zzeb333+f7Cj8l1wu37Nnj6u25t2J8Pjx41evXiU7Cv914MCByspKsqPwX1u3blWpVGRH4afUavXmzZvJjsJ/VVVVffrpp67amncnQgAAAMBJkAgBAAD4NUiEAAAA/BoNxx86JxYphEJheHi4g6OGNjc383g8oVDo7qjAA8lkstDQUA4HxkUkR01NTUxMDJ0OP2dJYLVapVJpfHw82YH4KaPRKJfLo6KiumyZl5e3bt26zttQLhHW19fr9XoHG5tMJgaDAV8EZDEYDJAFSQTvP7ng/SeXg++/RCLh8Xidt6FcIgQAAAA8Cc6lAAAA+DVIhAAAAPwaJEIAAAB+DRIhAAAAv+Y18xE2Nzfv2bOnubl58uTJWVlZHRtYLJb9+/eXlJQkJyf/+c9/ZrFYng/SVzU3Nx89evT27dvBwcGzZs1KTEy0a2C1Wnfv3m17OHDgwJEjR3o2Rl925syZiooKYpnJZC5cuLBjm/r6+r1796pUqmnTpo0aNcqzAfq4Xbt2tS8q7PjxrqioOHPmjO3htGnTIiNhlnmnNDU1FRcXS6XSsWPHJiUl2dbX1tYWFBS0tbU9+eSTw4cP7/hEo9G4Z8+eysrKlJSUefPmOXhPgXecEWIYNmLEiPLy8vj4+Ly8vM8//7xjm2XLlm3fvj0xMfGTTz55+umnPR6jL3vuuee+/fZbiUTS1NSUkpJy4cIFuwZms3nJkiV37ty5e/fu3bt3W1paSInTV+3bt+/zzz8n3tvq6uqODVQqVUZGRl1dXUxMzNSpU48fP+75IH1YdXX13d8sX768pKTErsHFixc3btxoa2MwGEiJ05eMHz/+jTfeeOWVVy5evGhbKZfLhw4d2tzcLJFIsrOzT5482fGJc+fO/fzzzxMTE7ds2bJixQpH94d7g7179w4dOtRqteI4/tlnnw0aNMiugUwm43A4UqkUx3GFQsHlcquqqkgI1EfpdDrb8pIlSxYuXGjXgDjy1Wq1Z+PyFwsWLHj//fc7abB58+asrCxieceOHZmZmR6Jy+/89NNPPB5PqVTard+3b19OTg4pIfkqi8WC4/iwYcP27dtnW/n2229nZ2cTy1u2bLF95m1u377N5/NbW1txHK+urubxeHK53JHdeccZ4blz5yZMmECj0RBCEyZMuHnzplKpbN/g4sWLffv2jYmJQQiFhISkpaWdP3+enFh9EZfLtS3r9fqHDeXz0UcfffDBBz///LOn4vIjly5d2rhx48GDBx8479K5c+cmTpxILE+YMOHixYswPZM77NmzZ+bMmcHBwR3/VFtbu3Hjxr179zY3N3s+MN/zwEuadp/z8+fPW61WuwYZGRmBgYEIoV69esXGxl6+fNmh3TkdsCfIZLKIiAhiOSwsjMlkymSyhzVACIlEovr6eo+G6B8uXrx4+PDh559/vuOfxo8f39LScvv27aysrI0bN3o+Nh8WFxcXERGhUCjWr1+fkZGBYZhdg/af/8jISKvV2tDQ4PEwfZxOpyssLHxgB21QUNDAgQNbW1sPHz6cnJzc8dopcAm7z7nJZJLL5e0bNDQ0tE8EkZGRDiYC7yiWYTKZZrOZWLZYLBaLhc1m2zWwWCy2hyaTya4BcF55efnMmTN37drVt29fuz+x2ezvv/+eWJ47d+6jjz6an58vEAg8HqNveuONN2wLaWlpe/fufe6559o3aH+AEAvw+Xe5L7/8MiQkZPTo0R3/NG3atGnTphHLS5cuff311w8dOuTZ6PxCl5/zHicC7zgjjI6OtiV2YkEikdg1qKursz2sq6tzZDBW4LiKiorx48dv2LBh1qxZnbccMWKE2Wyura31TGB+hcViZWRkdKyXaX+A1NXVsVis8PBwj0fn4/bu3bto0SKig6YTI0eOvHv3rmdC8jd2n3M+n293mbrHicA7EmFOTs7Ro0eJwbgPHTqUlZVFnG1cv379/v37CKExY8bI5fLi4mKEUGVl5Z07d2yXkoHz7t+//9hjj61evXrevHnt11+9epX42Ol0OtvKY8eO8fn8Xr16eThIH2Z7e9Vq9ZkzZwYMGIAQMhqNp06dIsqUcnJyDh8+TPQLHjx4cPLkyQ7O3wIcVF1dff78+fnz59vWYBh26tQp4rzE9g/Ccfz48eMDBw4kJ0pfl5OT89VXXxHnfAcPHszJySHWX716lUiQjz322I0bN4hfipcuXcIwLDMz06FNu6TCx93MZvPEiRPT09Pnz58fFhZ24cIFYn1WVta6deuI5c2bN0skkoULF8bGxtpWApeYNGmSQCBI/82zzz5LrE9NTd22bRuO47t27RowYMBTTz01adKkwMDAjz/+mNR4fU1oaGhOTk5eXp5EInniiSdMJhOO48SRf+/ePRzH9Xr9qFGjRowYkZeXFx4efu3aNbJD9jWvvvrq5MmT26+5c+cOQqilpQXH8UmTJo0fP37evHmPPPJIYmLi/fv3SQrTdyxfvjw9PV0gEPTq1Ss9PZ34zscwbNiwYaNGjZozZ05kZGRpaSnReNCgQTt27CCWX3vttfj4+IULF4rF4n/9618O7s5rZp+wWCxnzpyRy+VjxowRi8XEyvLy8oCAANvJb2lpKXFDfWpqKnmR+qCKigq1Wm17GBAQQNzieuvWrYiICKLXuri4uLq6OigoaMiQIXA3sWvdu3fv+vXrRqMxKSkpJSWFWGkyma5du5aSkkL0gphMptOnT6tUqnHjxrWvFwAuUVFRERgYaPvmQQjp9fqbN2+mp6czGAyFQnHlyhWlUhkdHT1ixAgYzcN5VVVVKpXK9jAxMZGoBSUuhLS1tY0fPz4sLIz4a2lpqUgksn3si4uLKyoqBg8e3L9/fwd35zWJEAAAAHAH7+gjBAAAANwEEiEAAAC/BokQAACAX4NECAAAwK9BIgQAAODXIBECAADwa5AIAQAA+DVIhAAAAPwaJEIAAAB+DRIhAAAAvwaJEAAAgF/7f0l8cT44QMeQAAAAAElFTkSuQmCC", "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.4" }, "kernelspec": { "name": "julia-1.8", "display_name": "Julia 1.8.4", "language": "julia" } }, "nbformat": 4 }