{ "cells": [ { "cell_type": "markdown", "source": [ "# Gross-Pitaevskii equation in one dimension\n", "In this example we will use DFTK to solve\n", "the Gross-Pitaevskii equation, and use this opportunity to explore a few internals." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## The model\n", "The [Gross-Pitaevskii equation](https://en.wikipedia.org/wiki/Gross%E2%80%93Pitaevskii_equation) (GPE)\n", "is a simple non-linear equation used to model bosonic systems\n", "in a mean-field approach. Denoting by ψ the effective one-particle bosonic\n", "wave function, the time-independent GPE reads in atomic units:\n", "$$\n", " H ψ = \\left(-\\frac12 Δ + V + 2 C |ψ|^2\\right) ψ = μ ψ \\qquad \\|ψ\\|_{L^2} = 1\n", "$$\n", "where C provides the strength of the boson-boson coupling.\n", "It's in particular a favorite model of applied mathematicians because it\n", "has a structure simpler than but similar to that of DFT, and displays\n", "interesting behavior (especially in higher dimensions with magnetic fields, see\n", "Gross-Pitaevskii equation with magnetism)." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "We wish to model this equation in 1D using DFTK.\n", "First 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": 1 }, { "cell_type": "markdown", "source": [ "which is special cased in DFTK to support 1D models.\n", "\n", "For the potential term V we just pick a harmonic\n", "potential. The real-space grid is in [0,1)\n", "in fractional coordinates( see\n", "Lattices and lattice vectors),\n", "therefore:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "pot(x) = (x - a/2)^2;" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "We setup each energy term in sequence: kinetic, potential and nonlinear term.\n", "For the non-linearity we use the PowerNonlinearity(C, α) term of DFTK.\n", "This object introduces an energy term C ∫ ρ(r)^α dr\n", "to the total energy functional, thus a potential term α C ρ^{α-1}.\n", "In our case we thus need the parameters" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "C = 1.0\n", "α = 2;" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "... and with this build the model" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using DFTK\n", "using LinearAlgebra\n", "\n", "n_electrons = 1 # Increase this for fun\n", "terms = [Kinetic(),\n", " ExternalFromReal(r -> pot(r[1])),\n", " PowerNonlinearity(C, α),\n", "]\n", "model = Model(lattice; n_electrons=n_electrons, terms=terms,\n", " spin_polarization=:spinless); # use \"spinless electrons\"" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "We discretize using a moderate Ecut (For 1D values up to 5000 are completely fine)\n", "and run a direct minimization algorithm:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter Function value Gradient norm \n", " 0 1.696984e+02 1.609230e+02\n", " * time: 0.0004928112030029297\n", " 1 1.685692e+02 1.241667e+02\n", " * time: 0.0017278194427490234\n", " 2 1.205450e+02 1.336305e+02\n", " * time: 0.0038568973541259766\n", " 3 8.836168e+01 1.181966e+02\n", " * time: 0.005646944046020508\n", " 4 2.358290e+01 6.401658e+01\n", " * time: 0.007350921630859375\n", " 5 1.849204e+01 5.369975e+01\n", " * time: 0.008742809295654297\n", " 6 4.844388e+00 2.381481e+01\n", " * time: 0.00999903678894043\n", " 7 3.232381e+00 7.725110e+00\n", " * time: 0.010939836502075195\n", " 8 1.929858e+00 1.037768e+01\n", " * time: 0.011874914169311523\n", " 9 1.489733e+00 6.579945e+00\n", " * time: 0.012948036193847656\n", " 10 1.240645e+00 2.854978e+00\n", " * time: 0.014071941375732422\n", " 11 1.181890e+00 7.487593e-01\n", " * time: 0.014979839324951172\n", " 12 1.164747e+00 4.232490e-01\n", " * time: 0.015902996063232422\n", " 13 1.154497e+00 2.919645e-01\n", " * time: 0.01692485809326172\n", " 14 1.150000e+00 3.531799e-01\n", " * time: 0.017912864685058594\n", " 15 1.147268e+00 1.875993e-01\n", " * time: 0.018877029418945312\n", " 16 1.144921e+00 1.009468e-01\n", " * time: 0.01997089385986328\n", " 17 1.144742e+00 1.068119e-01\n", " * time: 0.02079486846923828\n", " 18 1.144390e+00 6.054911e-02\n", " * time: 0.021623849868774414\n", " 19 1.144129e+00 6.319764e-02\n", " * time: 0.022655963897705078\n", " 20 1.144057e+00 1.269546e-02\n", " * time: 0.0234220027923584\n", " 21 1.144053e+00 9.805785e-03\n", " * time: 0.02440786361694336\n", " 22 1.144040e+00 6.332917e-03\n", " * time: 0.02537703514099121\n", " 23 1.144039e+00 5.551739e-03\n", " * time: 0.02636098861694336\n", " 24 1.144038e+00 4.880190e-03\n", " * time: 0.027254819869995117\n", " 25 1.144037e+00 2.694711e-03\n", " * time: 0.028239965438842773\n", " 26 1.144037e+00 1.783312e-03\n", " * time: 0.02936387062072754\n", " 27 1.144037e+00 1.534160e-03\n", " * time: 0.030261993408203125\n", " 28 1.144037e+00 1.421684e-03\n", " * time: 0.031137943267822266\n", " 29 1.144037e+00 6.666474e-04\n", " * time: 0.031848907470703125\n", " 30 1.144037e+00 2.176300e-04\n", " * time: 0.03268003463745117\n", " 31 1.144037e+00 1.642421e-04\n", " * time: 0.034372806549072266\n", " 32 1.144037e+00 1.307222e-04\n", " * time: 0.035455942153930664\n", " 33 1.144037e+00 1.194693e-04\n", " * time: 0.0364689826965332\n", " 34 1.144037e+00 4.667283e-05\n", " * time: 0.03755593299865723\n", " 35 1.144037e+00 4.223087e-05\n", " * time: 0.03855299949645996\n", " 36 1.144037e+00 2.881812e-05\n", " * time: 0.039528846740722656\n", " 37 1.144037e+00 1.428124e-05\n", " * time: 0.040448904037475586\n", " 38 1.144037e+00 7.671627e-06\n", " * time: 0.041380882263183594\n", " 39 1.144037e+00 5.224334e-06\n", " * time: 0.04231882095336914\n", " 40 1.144037e+00 4.781252e-06\n", " * time: 0.0430450439453125\n", " 41 1.144037e+00 5.081461e-06\n", " * time: 0.043782949447631836\n", " 42 1.144037e+00 2.191185e-06\n", " * time: 0.044551849365234375\n", " 43 1.144037e+00 1.239762e-06\n", " * time: 0.04538583755493164\n", " 44 1.144037e+00 9.169111e-07\n", " * time: 0.04632401466369629\n", " 45 1.144037e+00 7.352664e-07\n", " * time: 0.04721999168395996\n", " 46 1.144037e+00 4.971911e-07\n", " * time: 0.04829287528991699\n", " 47 1.144037e+00 3.084620e-07\n", " * time: 0.049362897872924805\n", " 48 1.144037e+00 2.350292e-07\n", " * time: 0.05026698112487793\n", " 49 1.144037e+00 1.129135e-07\n", " * time: 0.051175832748413086\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Energy breakdown:\n Kinetic 0.2682057 \n ExternalFromReal 0.4707475 \n PowerNonlinearity 0.4050836 \n\n total 1.144036852755 \n" }, "metadata": {}, "execution_count": 5 } ], "cell_type": "code", "source": [ "Ecut = 500\n", "basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\n", "scfres = direct_minimization(basis, tol=1e-8) # This is a constrained preconditioned LBFGS\n", "scfres.energies" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "## Internals\n", "We use the opportunity to explore some of DFTK internals.\n", "\n", "Extract the converged density and the obtained wave function:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "ρ = real(scfres.ρ.real)[:, 1, 1] # converged density\n", "ψ_fourier = scfres.ψ[1][:, 1]; # first kpoint, all G components, first eigenvector" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "Transform the wave function to real space and fix the phase:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n", "ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "Check whether ψ is normalised:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "x = a * vec(first.(DFTK.r_vectors(basis)))\n", "N = length(x)\n", "dx = a / N # real-space grid spacing\n", "@assert sum(abs2.(ψ)) * dx ≈ 1.0" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "The density is simply built from ψ:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "1.3250535760055254e-15" }, "metadata": {}, "execution_count": 9 } ], "cell_type": "code", "source": [ "norm(scfres.ρ.real - abs2.(ψ))" ], "metadata": {}, "execution_count": 9 }, { "cell_type": "markdown", "source": [ "We summarize the ground state in a nice plot:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=3}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3zU9f0H8M/ne3fJZe+7y15kT0LCTAh7yBAVK8W6By6stbZqHbW1/qyj2jrQ2iG1dRQnCDJkJgECBJJAEsggIfMuCdk7d/f9/P44ioggAT7fm6/n4/foL2S8v994d3ndZ1PGGAEAAHBUgqVvAAAAwJIQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NCuNghramoGBgbG/v2iKF7lFYEvxhi22bM2eJlYG6PRaOlbgO/h+4hcbRDeddddRUVFY//+wcFB/Nm1Knq9Xq/XW/ou4DuMscHBQUvfBXwPHhFrw/cRQdcoAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NLmlbwAAvjNgILW9rK1XSBRIoKul7wbAMSAIASyvoZ+9X8U+qBa1QyzSnbrJZI2D+l49WRwm3BUnzAmiArX0LQLYLwQhgCX1jJInDhk/qxNXRgtfzJWl+VLG2MDAsLu7e/co+W+t+HSR8Rd68l6ObJoaYWg/vvzyy9bWVkvfhS3x9/dfvny5RMURhAAW800ju6/AeE0orfmJwsvp/K96O5FV8cKqeOGLU+JNO43XhdOXJ8pc8JK1Cy+99FJQUJBKpbL0jdiGjo6OmpoaBCGAvXm7QvxjqfifGbLcwEs09a6PEGYGCvfvNc7bYtgwV+7jbJ4bBGn9+te/njx5sqXvwjYUFxffeeed0tXHrFEAc2OEPHfE+HqZuGvRpVPQxMeZfDxLNjuIZn9taOjH+S0APKFFCGBuzx42bmlkhUvl/srL+ClKyHMZMneFOGezcd+Sy/tZAPgRaBECmNX7VeJHNWzj/CtMssdShOURdNFWw6CB950BOCoEIYD5bG9mTx4yfrNApna58iIvZMmiPOkdeUb0kAJwgSAEMBPtILllt+G/s+RxXle1EIIS8v50WWM/+0uZyOveABwZghDAHBghd+cbHkgc6+yYH6eUkf/Olr1UaiztRLMQ4GohCAHM4fVjYs8o+U06t1dcqBt9ZZJs5U7jEAYLAa4OghBAcuVd7KWjxv/MkMm4bg7zs3FCuh99ssjIsyiA40EQAkiLEbJ6n/G5DFmEB/890t6aKvvvSbG4Ax2kYDGHDx/euHGj6eOPP/74YlvH7d2799ChQ2a8r8uAIASQ1n9qxF49uTdekteajzP5Q6ZsVYFRRBSChezdu3fdunWEkGPHjr322msX2zdOo9HcddddBoM1duUjCAEk1KsnTx4S35rKuVP0XHfGCXJK3q/CDFLgZmRkpKur6+w/e3p6zvuGjo4Oo/H8PvlXX3313nvvpfTCz/Xo6OiwsLANGzbwvVUuEIQAEnr2sHFxKJ2skvDgCErImmmyp4qMXSPSXQQcRXR09LPPPhsTE7Nw4UJCyBtvvBEWFjZ16tSYmJj8/HxCyHvvvRceHj5jxoygoKDf/va3Z39wZGTk888/X7p06eDgoK+v78jIyIYNG66//npCyP33379mzRpCyLJlyz788EML/WY/BlusAUilro99WCMeX66Q+kLpfvTacOGVo8b/y5JJfS2QzuZGZs4ebh9nMuUHb9G6u7tPnjxZW1srk8m2bt36zjvvHDlyxN/fv6CgYOXKldXV1fPnz7/zzjvlcnlXV1dWVtaiRYtMP1hcXOzj46NWqwcGBrq6uhhjo6Oj/f39hJDBwcHh4WFCSEZGxtNPP23GX3GsEIQAUnn2sPhwksw8m4L+NkNI/dzwYKIQ7IZjC23VmuNGc24XlOhNp6gu8M5p9erVcrmcELJu3brJkyeXlJSYPi+K4vHjx5OSkj799NPjx4+Pjo4qlcri4mLTV1taWi55qpRarW5ra9Pr9QqF5O8OLwuCEEASZV3s22ZxzTQzveCDXOkdscL/lYpvT0Wj0FZ9Pc8q/iD7+/ubPmhvbx8ZGdm+fbvpnz/72c/c3d1vvvlmg8GwYsUKX1/fw4cP9/X1OTk5EUKUSuXo6CghxDRGyM5p3DLGTJ8cGRmRy+WmlLUqVndDAPbhqSLxyTSZhxnf+P4mXRb/mf7RZCHaE41C4CA+Pr6jo+OPf/zj2c8wxjZs2HDq1CmNRsMYe+yxx85+KSYmpqGhQRRFV1dXb2/v+vr6s1+qra1dvHgxIeTUqVMxMTEXm01jQZgsA8DfoXZW0sFWJZj19eXjTB5MlD1fjOmjwMejjz76zTffPPvss3l5eevXr1+1apUoigkJCX/6058KCgoeeuihlpaWs98cExPj5+d39OhRQsjdd999xx13FBQUtLW1/frXv25sbJw/fz4hZO/evXPmzLHY73NxCEIA/v5YKj6WIijN3kn5SLKwqVGsx8m9cKV+/vOf+/j4mD7WaDSlpaWCILz77rsbN27Mzs4WBOGLL74YHBx8++23p0yZ8s4770ycODEzM9PU4Lvnnns+/vhjQsjLL7/82GOPHT9+vK2tzc/Pr7i42MvLizH2ySef3HPPPZb89S4CXaMAnFX2sHyd+MEMC0wH8FSQO2OF18vEP0/GSCFciWefffbcf6pUqueee+7cz0RGRr799tsX/NnVq1dnZWU9/vjjvr6+N9xwA2Psvffee/zxx01f3bBhQ3p6enJysjQ3flXQIgTg7OVScXWSzM1CbzIfTZH9p1psH7bM1cGRubu7FxUVeXh4XPCrc+fOff/99818S2OEIATgqXmAfVkvPpBosVeW2oXcECm8XYGduMEC3Nzczi6NWL58+bZt285+ydXV1TS/1AohCAF4er1MvCNW8HO25D38KlVYUyEOWOOejgDWCEEIwM2AgaytEn+eZOGX1ThPmq0R/l2N6aMAY4IgBODm39VibqAQ5m75ZVIPJQpvluNECrhsn3/++cDAAK9qBoNh06ZNF/vqli1bRkasYodcBCEAN2uOiw9ZujloMiuIygSyR4sohMtz7733trW18ar2zjvv7N69+2Jf3bt371/+8hde17oaWD4BwMcuLTOKZEag5ZuDJvcnCG9ViDMCsY4CLkNHRwevUgaD4eWXXzadWXFBDz30UHp6+sMPP6xUmmVD3ouzinevAHbgrXJxdZJgLTFIyC3jhF0tWFwPl2fatGlNTU233Xbb22+/nZ6e7u3tfc899zQ1Nc2ZM8fPz2/x4sW9vb2nT5+eP39+SEhIcHDwypUrOzs7CSGbNm2Kj48PCgp64IEH7rzzzry8vO3btwcFBUVERGzcuPEXv/gFIeTJJ5/84osvCCHZ2dmtra1qtTouLm7z5s0W/p3RIgTgommA7dGK/8q1oj313RXklhjh75Xi8xPQKLQNel0DIeZ740IVznI/zXmfNB0rceLEiZqamu3bt1NKs7Kyjhw58tFHH0VGRi5ZsuTdd99dtWrVCy+8kJmZOTo6+uCDD/7ud7974oknVq5cuXHjxpycnLVr1951113Lli0rLCzMysoihHR0dNTU1BBC6urqIiIiCCHFxcWmHbonTpyYl5d33XXXme23viAEIQAHa6vYTdGCuxXlICGE3BcvzNlsfC6DyKynoQoX1/mv/yPMfHN9FcHjfG/59cW++uCDD5qOocjNzfXz84uLiyOELFq06NChQ15eXiEhIR9++KFOp5PJZPv27du6devkyZNzcnIIIbfffrtpM5rGxsb4+Pgfv4fAwMCCggKOv9SVQRACXC1GyNpqcd0sq2t4xXvTUDeytYldE4oktAHqx9+19C18x9fX1/SBUqk89+Ph4eFDhw5de+21t912W3h4uLe3d3d3d0dHx9nDmwghAQEBhBBBEIxGI/nfwUwXZDQareFUJowRAlytXS3MTU4y/K0xbO6KE/5RiQWFwNPHH398zz33vPjii/fdd19CQgIhJDo6uqyszPTVoaEhU0doZGRkU1MTIUSlUul0urM/3t3dPTIyYgrOpqYmU2epZSEIAa7WPyrFe+Kt9KW0IkrYpRV1Q5a+D7AjwcHBO3bsqK2tzcvLe/XVVwkhixYt0uv1Dz/88KZNm26//Xa5XE4pnTVrVmFhISFk6tSpWq32pZde6urqqq2tve+++xYtWuTi4kIIKSwsnDVrloV/HwQhwFXqGSXfNIoroqz0peSuIMvChQ9r0CiEMbnhhhvc3d1nz56t0ZyZR5ORkWFq9hFCoqKipkyZ8uCDD2ZlZa1cuXLNmjV/+tOfFi5cqFAo8vLyvL29v/rqq3vvvVej0QQEBOTk5IyMjFRUVHh6ehYUFLS0tBQXF2/atGn8+PEfffQRIaSurk6n01nFCYXs6kyfPn337t1j//6+vj5RFK/yosDRyMjIyMiIpe/Chr1dblyx08CxoCiKfX19HAsW6MTET/UcCzqg3t5evgUnTZq0f/9+vjUtq62tzfRBfn6+n5+f6b/YRx99dNddd539nptuuundd989+88HH3zw73//+1iKHzlyJD09/dzP8H1ELD9KCWDT/lUtPp9pddNkzjVNTQ2MHGxnEwOscRQT7MMTTzxRUFCgVCq7urr++c9/mg5jWrFiRWJi4sV+ZNWqVT/yVXNCEAJcuZpeVt/PZlnNbjIXszJa+OikODHAqgMbbNo//vGP/v7+4eHhc6ePUkrT0tLO/vONN944dxOZlJQUs97ixVnpwAaATfh3tfjTaEFu9S+jn42jH58UDRgoBCm5u7ufm4I/pFKpPD09zXY/Y2f1r2AAK/ZJLbt5nA28iKI9aaQH3dGC7dYALsAGXsMA1qmwjTFCMq1y+eAP3RyNuaMAF4YgBLhCH9aIt9hCc9BkRbSwoUHs11v6PgCsj828jAGsikEkn9aJK6NtozlICAlQkqlquqEBjUKA8yEIAa7ELi2L9KDRnjYThISQldHCuloME8IlvPnmm+np6ZmZmW+99Zal78VMsHwC4EqsqxV/Yq27yVzMteHC6n36nlGZl5OlbwUuZNgwbM73KQIVnGXnPxX+8Y9//Otf//r8888HBwevv/56f3//FStWmPGmLANBCHDZDCLZ0CA+Pd7GXj4eCpIbKGxsEG1ipqsDuv6L241mPIYpyT/utdnPn/fJ995777nnnouOjiaE/OpXv1q7di2CEAAuYJeWRXvQcHdb6hc1uTGSflrHbh5n6fuAC/nmJ59Y+hZIfX19ZGSk6eOoqCjT8RF2D28MAS7bp3XijbbWL2qyNFzYoxV7Ri19H2Ct9Hr92fBrbGwMDAy07P2Yh02+mAEsyCCS9fXi9RG21xwkhHgoyPRAYSPmjsLFvfTSS52dnc3Nza+99trKlSstfTvmgCAEuDy22y9qYuodtfRdgPW68cYbZ82alZOTs3z58ttuu83St2MOGCMEuDyf1YnLI234HeTScOGhffo+vcxDYelbAas0d+7c+++/39J3YVY2/HoGMD+Rka8bxGW22S9q4qkg09R0cyN6RwHOQBACXIb9bUylpFEeNhyEhJBrw4X19egdhQv44osvgoKCLH0X5oYgBLgM6+ttuzlocm24sLlJHEWbEH4gNzfX1dXV0ndhbghCgMuwvp5dG27zrxq1C0nwpru1aBQCEIIgBBi78i42YiTpfjbfIiRnekfRJAQgBEEIMHZf1bNlEdQeYpCQ6yPoV6cYmoQABEEIMHbr60U76Bc1GedJvZ3IoXZEIQCCEGBsWgZZbS/LUdtHg5AQQq4Np+gdBSAIQoAx2tjAFoQKcjt6xSwOEzY1oEUIgJ1lAMZmUyNbEWU/zUFCyGQV1Q2x+n5mu9vF2a5JkybNnz9fLsdf4DExGo1nz8SQAh4GgEsbMZI9WvGf0+1qUzKBknnBwjeN7P4EBKG5vfrqq88++6yl78KWSPqmQf7II4+kpqbedtttMpns3C9otdoPPvhAq9VGRUXdcccdHh4e0t0EgJXb2cLSfKmfs6Xvg7dFYfTf1eL9CXbU4WsjFAqFn5+fpe8CzhASEhLee++9Bx544NzPlpeXp6am1tTUREVFVVZWarVaS90fgDXY1CguCrPDtFgQIuTr2KDB0vcBYFHyVatWLV26NDIy8rnnnjt7BuNDDz20evVqtNwBTDY1sk3z7TAIvZxIhj/d2cIWh6F3FByXQAgJDAyMi4vbt2+f6VODg4N79uxZsmTJX//613feeae5udmidwhgYWVdjDGS6G2fUbEoTNiEkyjAsZ0ZflSpVGf7P+vq6iil99xzz5IlS3Q63TPPPLNv377Y2NgL/vzQ0NCLL774/vvvn/vJefPmXX/99Rf8/uHhYZlMRu1kdw57MDo6SggRRfwpvKj1tcI1wWR4eNg8l2OMDQ8Pm20+4VwV+csx2esT9Oa5nI0aHh5WKOxqqpStG/sj4uTkJAiX6M4582LT6/VOTk6mj2UymSiKv/jFL26++WZCyNDQ0J///Oc1a9Zc8OdlMllKSsp5MRkVFXWxW1QoFAqFAkFoPRhjhBC8yH/Et1rxsRRBoZBd+lt5YIyZXibmuVyyP3GSsaoBRZK3eS5ok8z5iMBYjP0RGUvcnAnC5ubms2dQBQUFUUoTExNN/0xMTMzLy7vYzzs5OS1evDg3N3csN0QIkclkaBFaFdNs4fPmDMNZfXpS3CHOCjbffyHGmOllYqbrEbIg1Lithab62eEgKC9mfkTgkvg+IgIhpLi4uK2tbcaMGZ2dnbt37/b09Jw9e3ZhYaHpO/bv35+UlMTregC2ZXuzOEVNXe16we38YLq1CX3j4Ljkd95555YtW55//nl3d/ft27cvW7asv7//hRdeWLJkSWFhoU6na2xs/Nvf/mbp+wSwjK1NbH6InTeVZgcLt+w29uuJOzr/wCHJ582b98tf/tLU5ps4caKpF3TixInl5eW7d+/28fHJzs52dra7hcQAY7Otmf082c6D0E1OsgLobi0WUYCDkq9YseLsPzw9PTMyMkwf+/v7L1++3EJ3BWAVjnczg0gS7HThxLnmhwhbmsTFYRgGA0dk5291Aa7GliZ2Taj9pyAhZEEI/aYRJ1GAg0IQAlzU1iZxgWMEYYovNYikugdZCI4IQQhwYUMGsr+VzQpylNfI/BC6pQlBCI7IUV7kAJcrT8fS/Kinw0yknBtMtzVjEQU4IgQhwIV92yzODXagF8jsYCFfx/SIQnA8DvQ6B7gs3zazucEOMUBo4udMxnnSwjb0joLDQRACXEDrEGkcYFkBDhSEhJB5wfRb9I6C40EQAlzAt83izEBB5lg5SOaGCN82o0UIDgdBCHABjtYvajJNTSu6WOeIpe8DwLwQhADnY4R82yzOC3G4IHQSyFQ13a1F7yg4FgQhwPnKu5izjEZ5OFwQEkLmBqN3FBwOghDgfN82s/mO1y9qMi+EIgjB0dj1MWsAV2RHs3h7rAXeIw4bhqu76qo6TnYOdMarYmN9o9VuAWa+hyQfOmhgp/pYhEM2iMExIQgBvscgkr2tbG2uWYPQIBq/qNr4YdlnGndVrG+0kjpvrNl2oqM6RZXwQMadQe4as90JJWRmoLBTy+5EEILDQBACfM/BdhblSf2V5rtiTVfd7wpe0bir3pz3YphnCGNsYGDA3d1db9SvO7H+vi2P3Ri/9Jbkn5jtfmYF0R3N7M5Ys10QwMIQhADfs6OFzQo0X2PoeEfVk7v/sHrC3bMjpp/3JYVMcXPS8gVRs5/c/XzXcM/qzLspMceNzQmmTxUZGXG0VZTguDBZBuB7draIs821xejRtvIndz//xJSHf5iCZ/m5+Lw+5w+VndWvHnhbZOaYxhLuTt3ktKILU2bAUSAIAb4zaCCHT7NstTnaQrqBtmfzX3p22q8mB2X++He6KVxfnfX7+p7GD8s/NcONEUJmB9OdLQhCcBQIQoDv7G1laX7UXfqjlwyi8fcFr/4s6cYMTepYvt9Frnx++pNfVW8+rCuV+t4IIbMC6Q4EITgMBCHAd3a0iLODzNEcfLf4fS9njxviF4/9R3yU3k9PffT/9r3eOdQl3Y2ZzAoS9mhFI6IQHAOCEOA7O1vYbOmPpD+kLS5oOvDU1Ecvd/LLeHXKNdFzXz34tkQ3dpbKhYS608OnkYTgEBCEAGd0j5LKbjZJJW2LUC8a3ij6288zV7k7uV3Bj9+aclNjb0thSxH3GzvP7CAME4KjQBACnJGnFSerqJPEr4nPTmwI8QiaEnyJCTIXoxDkj2SteqPob3qjnu+NnWdmIHbfBkeBIAQ4Y7eWzZC4X7RruPuTii8fyLjjaopM0KRFeoV9WrmB111dUG6gsK+VjSIKwQEgCAHO2KVlMyVeSv9u8dolMfNDPYOvss5DE+7+pOJLSWfNeDmRGC9a1I7eUbB/CEIAQgjpHCG1vWyCv4RB2NKvO9ByeGXiDVdfKtBdPS9yxifHv7z6Uj9iZiDdpUUQgv1DEAIQQsgerThNQxVSviA+OPbf62IXuypcuFRbmXjD5todXcM9XKpd0MwgAcOE4AgQhACEELJby2YGSvhyaB1o29d86Ia4y1g4+ON8XXzmRExfd/wrXgV/aLqGHmxjI0bprgBgFRCEAIQQskvLZkg5QPhB2bplsddc2ZKJi1mZeMPGmm3dI1I1Cj0UJN6bHsQwIdg7BCEA6Rghjf0sQ7IBwrbB03mN+5fHL+FbNsDVf0b4tM9PbORb9lwzMEwIDgBBCEB2tYjZGirdsUPrq76ZHznT08mDe+WbEpZ9XbNl1DjKvbLJzCBhVwuGCcHOIQgByB4tmyHZAKHeqP+mdse1MQulKB7iERTjG72rvkCK4oSQbDUtOs2GMUwIdg1BCEB2SzlAuKM+P8Yn6urXDl7MdbGLPqv8WqLi7gqS5INhQrBzCEJwdB0jpGmApftJFYRfVm26LnaRRMUJIZODMvtHB453VElUP1dDd2OYEOwaghAc3R6tOE0t1QDh8Y6qruGeSUEZklQnhBAiUHpt7MIvq76RqH5uoLAHqwnBriEIwdHt0bJcyQYI11dtvj5ukUClfaFdEz1nX9PB3pE+KYrnaOihdqwmBHuGIARHt0fLcqUZIBw2DOc3Fc6PnClF8XN5OnlMDsrcfipPiuIeChLnRQ9hmBDsF4IQHFrnCKnrY+OlGSDcVV+Qrkr2UXpLUfw8C6JmbanbIVHxGYF0jw5BCHYLQQgOLV8nTlVLtcXoltqdC6JmSVL6BzI0aV1D3bXd9VIUxzAh2DcEITg06QYIWwfa6noaJgdd4QG8l0ugdF7kjG11u6QonqOhhW04mxDsFoIQHNoenVQDhFtqd86OmK6QKaQofkELomZvrdtlZPyntXg5kXGe9PBp9I6CfUIQguPqGSU1PSxTgi1GGWHm7Bc1CfUM1rgFFGlLpCg+I5DuwWpCsFMIQnBcBa1sokqSAcLy9kpnmVOc7zj+pX/UgqjZ2+p2S1E5R0PzdegbBfuEIATHla8TczSSvAR2NeTPisiRovKPmxE27UDL4REJ9uDODRT2tTIj2oRgjxCE4LjytGy6hn+/qMjY7oZ9M8Kmca98SV7OnnF+4w60HOZe2deZBLvR0g4kIdghBCE4qEEDOdbFJgbwD8Ky9gpvZ88wzxDulcdiZti03dIcRjFdQ/OwmhDsEYIQHNT+NpbuR13l/CvvathrkeagyfTQqYUth4cNw9wr52hoPoIQ7BGCEBxUvk6UqF80r2FfbthU7pXHyNPZI8k/vlCC3tEZgUK+TkQSgv1BEIKDytMyKWbKHG0v91F6W6pf1GRG+DQpjuoNdCWeTvR4N6IQ7A2CEBzRqEgOn2ZT1fxbhHsa9s4Iz+Ze9rLkhE4u0pVI0Ts6XUPzsJoQ7A6CEBzRoXYW60U9ee/6wggraDwwPXQK57qXydPJI94v5pAEK+sxTAh2CUEIjihfx3IkGCCs6jyplCvDPIO5V75c2SGT9jYd4F42NxATR8EOIQjBEeXrxOkSbDG6t+lAdsgk7mWvQE7I5H3Nh7jvOxrlQQkhp/qQhWBXEITgcERG9rexqSr+T/6CxgPZoVYRhP6ufipX//L2E9wrT1OjUQj2BkEIDudoJ9O4UJUL57KtA22dw90JfnGc616p7NBJBRL0jmKYEOwPghAcTr5Okp3V8hsPTAuZKFBJDnW6Atkhk/IbC7mXRRCC/UEQgsORaKZMQVPhNOsYIDQZ5xPFCKvvaeRbNsWXnh5mbUN8qwJYEoIQHM7eVv5B2DfaX9V5coImjW/ZqzQlOIt77yglZIqaFrTiSCawHwhCcCw1vUygJMydcxAe1Banq5OdZU58y16lKcGZUuy1lqMR0DsK9gRBCI4lT8dyJVg4caDl8KSgCdzLXqXxqpSTXXW9o318y2KYEOwMghAcixQDhCJjB1uOWGEQKmSKNFXSYW0p37KZ/rSqh/XwP/0XwDIQhOBYCiQIwurOk17OHho3Fd+yXEwKmsD9nF6FQDL96f42NArBTiAIwYG0DLKuEZbgzTkIC1uKJgdl8q3Jy6SgCQe0RxjhHFrZGpqvw3wZsBMIQnAgBTqWrRG4jxAWWuUAoUmgu9pd4VbdWcu3bLZGKMAwIdgLBCE4kHwdy+bdL9o72lff05gSkMC3LEeTgzMLW4r41pyqokc62DDnrUwBLANBCA6koJVl8z6D8GDLkQxNqkLG+0gnfiZLMEzoriBxXvTwaTQKwR4gCMFR9OpJTS/L8OcchAdajkwMzOBbk6/UgMS67gYsogC4GAQhOIq9OpblT524PuUZYYd1JVmB43kW5U0hU6QEJB7RHeVbNltNCzBfBuwCghAcRUGryH3hRG13vVKuDHRX8y3LXWZgehHvA+unBwp7W5kRbUKwfQhCcBT5Opat4fyEP6QttvLmoElmYPohbTHfmgFKonKh5V1IQrB5CEJwCCNGUtzBJqs4twiLtCWZmnS+NaUQ6RVmZGJzn5Zv2Ww1hgnBHiAIwSEUnWZxXtSD69TOUeNo+ekT6epknkUlk6lJK9Jx7h3N0dCCVgQh2DwEITgEKXZWO9Z+PMo73MPJnW9ZiVuwOvkAACAASURBVGQGph/iPUyIiaNgHxCE4BDydSL3FYSHtMWZGhsYIDSZoEk/ois1iDzXwEd7UsbIqT5kIdg2BCHYP0bI/jY2Vc352V6kLckMtIEBQhMfpVegu6ays5pv2Wlqmo/eUbBxCEKwf2WdzF9JA1151uwe6dEOtCb6x/IsKrEsCeaOZmsoNh0FW4cgBPuXr+O/s1qx7liaKklGZXzLSipDk1qsO8a3ZrYaQQg2D0EI9q+glf9e20daj45Xp/KtKbXUgKTKzpphwzDHmml+tGWQtfMsCWBuCEKwf1JMGT2iOzpBk8a3ptSUcucY3+hj7cc51pRRMklF97dirzWwYQhCsHP1/WxUZOM8eQbh6cGOAf1ApHcYx5rmkaFOLW7l3jsqYDUh2DQEIdi5fB3L4b2zWpGudLw6lRLuR/xKLkOTelhXyrdmNlYTgo1DEIKdK5BgpsyR1qPj1Sl8a5pHon9cY29z32g/x5qTAuixTjZo4FgSwKyEBx544IMPPhDFC3Tx9/f3v/zyy/n5+ea/LQBeJJkpY4MDhCYKQZ7oH1faVs6xpoucpPrSA+1oFIKtEiZMmPDaa6/95je/+eHXnnjiiRdffHHbtm3mvy0ALjpHSGM/S/PlGYSNvc2UkGCPQI41zSlDk8r/bEKsJgRbJtx1110fffTR22+/3df3vQOs9+3bV1FRsWDBAkvdGcDVK9CJk1VUznUE4EirrTYHTTLUqUda+R/SuxcTR8FmCYSQxMREDw+P4uLvtpwYGhpatWrVmjVrBAGDiGDDClr5n0FY3HrMRgcITWJ9o08PdnQN93Csma0R9rcyA6IQbJPc9P9UKlVLS8vZzz799NMrV66Mj4+/5M8PDAw8//zzf/3rX8/95Pz582+88cYLfv/g4CCllFLbm25nr0ZHRwkhBoN9TnXIa5H/Ll0cHOT2F5oRVqw7dnv8TYODg7xqnn8JxgYHByV9D5roG3ew8XBO8GReBZWEBLnID7YMpfvaZwfp0NCQTGZLuwjZvbE/Ikql8pKvpjNBqNfrnZ2dTR8fOHBg165dBw4cGMs1FApFVlZWXFzcuZ9MTEw8W+08pgshCK2H6bFwcnKy9I3wN2QgZd3i1ECFs5xbzVM9ja4Kl1CfEG4Vf4AxZjAYLvYK4mK8JqWiq2pOVC7HmtMD2cEuYVKgfb60R0dHJX1E4HKN/REZS9zICSGiKGq12uDgYNOnPvvss/b29ilTphBC6urqnJ2d29vb33333Qv+vJOT04IFC3Jzx/qKkslkMpkMQWg9TO+q7PLdblEbS/FlHs48f7WyjuPp6mRJ/3MxxkwvE+kukRGY+uK+P/O9RE6guKGePZJih08k8r8/XJa+C/gO30dETgjZunWrm5vbhAkTtFptVVXVE088cf/995u+/Itf/CI0NPSpp57idT0As5FiZ7WS1rLJQRP41jS/GJ+o00OdXcM9PkovXjWz1fRXB4yEIC3A9gg333zzrbfe+sorr8hkst27d992221+fn5R/+Pu7u7j4xMYaKszxcGRFbTyP4z3aFt5ujqZb03zE6iQEpBY2lbGsWaEB1UItKbXPscIwb7Jb7jhhhdeeCEiIoIQsnDhwvT07x00+oc//AE942CLjIwcaGP/nsFzykl9T6NCplC7qTjWtJR0dXJpW9mMsGkca5pWE/Ld1hXADITrr7/elIKEEG9v74SEhHO/HBkZGRQUZIH7Arg6pR0syJUGKHnWLGkrS7flhRPnSlclF7fybBES09mE2H0bbBCWCYJ9kmJntZLWsnSVzfeLmsT4RkmwmhC7b4NNQhCCfcqXYKbM0bby8bY/QGgiUCElIOEo101Hk33o6WGmlWqBJYBUEIRgn/a2inyDsKG3WS7I7WOA0CRNnVzCdb6MQMkUFd2HvdbA1iAIwQ5V9zAZpeHuPIPwaFuZHcwXPVeaKolvi5AQkqMR8jFMCLYGQQh2KL+VTefdL1raVp6qSuJb07JifaO1/a29o32X/tYxy8EwIdggBCHYoXwd/5kypW3lafYVhDIqS/SPK2s/zrFmZgCt6mE9oxxLAkgOQQh2iPtMGd1Am140hHjY21KiVFXi0bYKjgWdBJLhRwvb0CgEW4IgBHujGyKdIyzRm2cQltrRwolzpauS+e4vQwjJ0dACzJcBm4IgBHuTrxOz1YLAtWf0aHuFnfWLmiT4xdZ1NwzqhzjWzNEIGCYE24IgBHsjxQpC+xsgNFHIFHG+4ypOV3KsOU1ND59mI0aOJQGkhSAEe1PAe6ZM51BXz3BvuFcYx5rWI02dVMp1EYW7gsR50aLTaBSCzUAQgl3p1ZOaXpbpzzMIS9rK0tRJgp0eopkakMR9mHC6hhagdxRsB4IQ7EqBjk0MoAquz+ujbRWpAYk8K1qT5ID4qs6To0aeKx6yNTRfh/kyYDMQhGBX8nVijobzs7q0rczOltKfSylXRniFneio5lhzukbY28qMaBOCjUAQgl3J17HpgTz7MPtG+3UDbeN8ojjWtDbcVxP6K0mgKz3WiSQE24AgBPsxZCClnWxSAM8gPNZekeAXKxdkHGtam1RVUmk7501Hp2toHoYJwUYgCMF+HGhnqb7UVc6z5tE2+1xBeK5UVWLF6UqR8RzVw6ajYEMQhGA/8qRZQWjHA4Qmnk4eAS5+J7tPcayZG0jzdCKSEGwCghDsB/eZMiPG0dru+gS/WI41rVMq7yOZQtyoq5xW9SAKwQYgCMFOGERyqJ1NU/NsEVacroz2jlDKnTnWtE7c58sQQqajdxRsBIIQ7MTh0yzKg3o78ax5tK08VWW3KwjPlaZKKmktY4RnbmGYEGwFghDsRH4rBgivXICrv1Lu3NTbwrHmdA3do0UQgg1AEIKd2KMV+Z5Kb2TGEx3VyQHxHGtas1QV501HY72oXiT1/chCsHYIQrAHIiN7WxnfmTLVnbVqtwBPJw+ONa1ZqirxaDvnYcJsrCYEW4AgBHtQ2skCXajKhWfNow7TL2qSEpB4TIL5MnnoHQWrhyAEe5Cn5byzGiHkaPtxO95r+4fCvUIG9UPtg6c51swNpHvQIgSrhyAEe5CnY7lcBwgZYWXtFSmOMWXUhBKaHBB/rP04x5pJPrRnlLUMIgvBqiEIweYxQgpaRb4twsbeZmeZs8rVn2NN65eiSjzGdZiQEjJNLWARBVg5BCHYvIou5qmgQa48g/BoW4WDrCA8V2pA0tE2ni1Cgt23wRYgCMHm5elYLv8BwgqHmiljEusb3dKv7R8d4FgTqwnB+iEIwebt0TK+KwgJIcfaKlIcaaaMiVyQxfmOKz99gmPNND/aMsjahjiWBOAMQQg2L0/HeSl9x1BX/+hAmGcIx5q2IlWVxHc1oYySqSpa0MrzjCcAvhCEYNsqe5hCoBEePIOwtK0sVZUoUM6tTJuQGsB/9+3cQAG9o2DNEIRg23Zr2awg3v2iDrZw4lxJAfHVnSf1Rj3HmjMC6W4EIVgxBCHYtj1azisIiWnKaIDDzZQxcZErwzxDTnTWcKyZ4U8bB9jpYY4lAXhCEIJt26PlPGV0QD/Y3KeN8Y3iWNO2pKoS+R7SK6Nkiorm6TBMCFYKQQg2rLKHySiJ5DpAWNZ+PME/ViHIOda0LSkBnJfVEwwTgnVDEIINk2KA8Gi7Iy6cOFeaOqms/YTIeOYWhgnBmiEIwYZx7xclZwYIHToIvZ29vJVep3oaONacgGFCsGIIQrBhe7RsBtcg1Bv11Z0nE/3jONa0RSkBCRgmBMeBIARbJcUA4YnO6jCvEFcF14MNbVBqAP9DejFMCFYLQQi2ao+WzeQ+QOjw/aImqaqkUq4tQkLIzEC6C0EIVglBCLZqZwubyXuA8JhD7rX9Q8EegYQx3UAbx5oZ/rRpgLVi01GwPghCsEmMkD1akW+LUGSsrP1ESkACx5q2K1mCYcJstYBhQrBCCEKwSeVdzE1Bw915BmFdT7230stH6c2xpu1KVSXy7x0Norta0DsKVgdBCDZpV4sEKwjbyjFAeFaqKulYO+dDemcF0Z0IQrA+CEKwSbu0/AcIj7ZhgPA70d4RHUOdXcM9HGum+tLOEdY0gCwE64IgBNsjMpKnFbkvpT/WfjzVUQ+d+CGBCon+cWVcF1FQQnIDBWwxA9YGQQi2p6SDqV1okCvPIGzp14nMGOSu4VjT1qUFJHE/m3BmIIYJweogCMH27JRgBWFpW3maKplvTVuXqkosbec8X2ZWEN2OIAQrgyAE27OrRZRggLAc/aLnSfCLbehpGtTzXPoX702NjNT2IQvBiiAIwcboRbK3lc0I4vzUPdpWgRbheRQyRYxvdMXpSr5lZwbSHc0IQrAiCEKwMQfa2DhP6ufMs2bXcHfPSG+EVxjPonYhVZV4VILe0Z2YLwPWBEEINmZ7izgnmP8AYUpAgkA5l7UDaQFJpbzny8wLptubRRFRCFYDQQg2ZkcLm827XxRbjF5MUkB8VWeNXjRwrBnsRn2d6bEuJCFYCwQh2JI+PSnpYFPVvFuErZgyemFuCtdQz+DKjhq+ZecE0+0YJgSrgSAEW5KnYxMDqJucZ80B/WBLvy7WN5pnUTuSpkoubSvjW3N2EN3Rgt23wVogCMGW7GgWufeLHm2rSPCPlQsyvmXtRmoA/923ZwUJe3VsxMi3KsAVQhCCLdnewiSYKVOWGoABwotKUyeVnz4hMp4NOG8nEutFD7SjdxSsAoIQbEbrEGkaYBP8JdhTRo0gvChPJ48AV//qrlq+ZecE0x3N6B0Fq4AgBJuxs0XM1Qgyrjk4bBg51dOQ4BfLs6jdSVMlce8dnRMsYK81sBIIQrAZ25rZvBDOzcHy0ydifKKcZU58y9qZNBX/3bez1bS8i/WM8q0KcCUQhGAztjezeRIMEGLhxCWlq5JL28pExrMB5ywjU1R0J+aOghVAEIJtKOtiShmJ9uS/ghB7bV+Sr4uPp5NHfU8D37LzQoRvsZoQrACCEGzDtib+/aJ6o76q62RyQALfsnZJimHCecF0SxOCECwPQQi2YVuzOJd3v2h5R2WEV6iLXMm3rF1KVyeX8F5Wn+RDDSI52YssBAtDEIINGDaS/a1sZiDnp2tpa1m6KoVvTXuVrk4paT3GCOfQmhNMt6JRCJaGIAQbkK9jqX7Ui/fUzpK2snQ1ZsqMicrV30Xu0tDTxLfs3GCKYUKwOAQh2IBvm8W5wZyfq3rRUNlRgwHCsUuToHd0brCwWyvqMXUULApBCDZgWxP/hRMVpyvDPEPcFK58y9qxdFVySSvnIPRXknGetLANjUKwJAQhWDvtIGkaYFkBnIOwpBX9opcnQ5Na3HqU+zDh/BC6pQlNQrAkBCFYu28axXkhnHdWI4SUYoDwMqlc/ZVyZWNvM9+yC0OFzY1oEYIlIQjB2m1uYgtDea8gFA3HO6pSArCU/vKkq5KLW4/xrTlZRRv6WcsgshAsBkEIVs0gkp0t4jzeM2VOdFRhgPAKpKmTS3kPE8oomR0sbMMiCrAcBCFYtX1tLNqDql04ly1uPZaOLUYvX4Y6pbitjPsw4cIQuhlBCJaDIASrtrlR5N4vSgg5ojs6XoOl9JdN7aZyljlxX024IFTY3iwaMGMGLARBCFZtcyNbGMr5WTpqHK3srMEA4ZXJUKceaT3Kt6bGhYS748B6sBgEIVgv7SBpHGATeS+cKGs/EekVjgHCKzNek8J9vgwhZGEo3dyIJiFYBoIQrJdECyeKW49laFI5F3UYGeq04tZjfM8mJIQsDBG+wSIKsBAEIVivrxvYYikGCFtLx6sxQHiF/Fx8vJ29arvr+JadoqaNA6y+H1kIFoAgBCs1YiS7teIC3gOEQ4bh2u56bDF6NTI0qUd0nIcJZZQsCMHKerAMBCFYqR0tLNWX+jlzLnu0rTzOd5yzjPdJFo5kvFqSYcLFYXRjA4YJwQLkq1evdnJyuvvuuxMSvnuP3NjY+OWXX1ZWVnp6eq5YsSItLc2CtwiOaVOjuDiM/xu14tZj49UYILwqGerUlwvfNIhGuSDjWHZBiHBvvnHQQFzlHKsCXJoQHx/v5eU1derUU6dOnf3sk08+WVpampKSIpPJpkyZsmvXLsvdITiobxrZkjBJVhBmYAXh1fF09gh0V1d11vAt6+VEJvjTHS1oFIK5yR988EFCyMmTJ995552XXnrJ9Nm1a9fK5WfelfX09Kxdu3bmzJkWu0dwPKWdTCAkwZtzEPaO9jX1tcT7xfIt64AyNGmHdaWJ/nF8yy4OEzY2sCVhfKsCXMKZrqdZs2bt3bv37GfPpiAhpLu728fHx9z3BY5tYwNbEi5JczAlIFEhoOvtak3QpBbpSrmXvTacbmzgvTID4FLO/EUICAjQarU//PKePXu+/vrrkpKSi/38wMDAb3/724CAgHM/ec0119x0000X/P7BwUFKKaX8/8bBlRkdHSWEGAwGS9/I93x9Sv50qjg4OMK37IGmw6l+iYODg3zL8sUYGxwcFASrnsgW6x5V2VHd0dvpIldyLBsoJ24y+b6mofG+1pWGQ0NDMhnPAVG4SmN/RJRK5SVfTWeCcHR0VKk8/wldUlJy0003ffjhhxERERf7eScnp5ycnMTE7+1WFRcX98NqJgaDQalUIgith+kp4uRkRbMoW4dIdZ9xTpiTE+8sKGkvWx6/5GJPTivBGDMajVZ+k0qijPMdV91XOzEwg2/l6yLFLTphSpB1vQ/Q6/VW/og4mrE/ImN5T3kmCBsbG4ODg8/9Qnl5+TXXXPP2228vWrToR35eoVDMmTMnNzd3LDdkuidBEBCE1sP0LLGq9sf6BvGaUEEp53xLuoG2YeNIlG8EJVb99GOMmV4mlr6RS5gQmHa49ejk4Ey+ZZdF0FUFxuczrevXt4lHxKHwfUQEQojRaFy3bt2yZcsIIbt3725oaKiqqlqwYMErr7xyww038LoSwBitrxevlWCA8JC2OFOTbuUpaEMyNeMPay86aHLFJgXQ08Osuse6ukbBvgkLFy6cNm2aXq+/7bbbCCEPPfTQ1q1bH3744Z6entdffz0zMzMzM/O+++6z9H2Co+jTk/1tbH4I/3ffRdqSTA1WxHIT7zeufbCjY6iLb1mBkqVhwoYGBCGYj/yRRx5RKpXTpk0zzRTdsGGDr6/v4sWLh4aGzn6Tiwvvc1EBLmJTg5ijoR4KzmVFxkpayx6acBfnug5MoEKaOvlIa+nciBl8Ky+LEP6vxPjLFHRFgpnI58+ff+6/o6KiCCHe3t4Wuh9wdF/Vs2Xh/P8CVnee9FZ6Brj6c6/syLIC0w9r+Qfh7CB68y6mGyIavAMHs8B7LrAiI0aytUlcJMHOakW6kszAdO5lHVymJv2QroQRzt2YCoHMCxGw7yiYDYIQrMjOFpbsS6VoBxxsOcJ9oj8EewQ6CYq67gbula8No1+dQhCCmSAIwYp8Wicuj+D/nBzUD1V31aapkrlXhklBEw60HOZednGYkK9jXZw3VAC4MAQhWAu9SDbUi9dF8F/ecFhXkuQfr5TzPtIJCJkYNP6Qtph7WXcFmRUkfI3eUTALBCFYix0tLNaLhrnzD8KD2uKsoPHcywIhJEOderyjasgwzL3y8kj6WR0WUYA5IAjBWnxWJ94YJckT8pC2eHLQBCkqg1KuTPCLLW7lfGA9IWRpuJCnE3tGuRcGOB+CEKyCQSRfN4jXS9AvWt/TaBSNYZ4h3CuDSVbg+IMt/HtHPRQkNxC9o2AOCEKwCju1LNKDhkvQL3pAe2RKcBb3snDWpKAJhS1FUlReHkk/Re8oSA9BCFbh8zpxeaQkz8aDLUeyAjFAKKEo73CjaGzqa+FeeWmYsLtF7NVzLwzwPQhCsDy9SL48Jd4Yyb85OGwYLj99IkOTyr0ynCsrcLwUiyi8nMiMIGF9PXpHQVoIQrC8rU0szluSftFD2pIEv1g3hSv3ynCuKcGZ+5oPSVH5p1H045MIQpAWghAs7+OT4k+lmS9a2FKEAUIzyAocX3G6ckA/yL3y0nBhfytrHbr0dwJcMQQhWNiggXzTKMkAISOssOXwFN4nx8IPKeXKZP+EIgmOJ3SVk0Vhwmd1aBSChBCEYGHr68WpaqqSYH/Ryo4aN4VriEcQ/9LwA1OCs/ZL1DsaLaB3FCSFIAQL+/gk+2m0JM/D/c1FU9Evai7TQiYWthSJjP9qh/nBtKaX1fVhHQVIBUEIltQ1QvJ04lIJDiAkhOxvPjQZ/aLmonYL8FZ6n+io4l5ZLpDrI4RPahGEIBUEIVjSf2vFBSGCJ+/z6AkhHUNd2v7WlIAE/qXhIqYGZ+1vlmRl/cpo4T/V6B0FqSAIwZLWVou3xUjVHMwKTJdRmRTF4YKmBGftaz4oReVsDTUycqgdjUKQBIIQLKaqhzX2k3kh/JcPEkLyGwtzQqdIURkuJsk/vmu4u6VfJ0Xxm8cJ/0KjEKSBIASL+WeVeEsMlUmQg4P6oWPtFRODcCS9WQmUTgnOKmgslKL4rTH0v7XiiFGK2uDoEIRgGSIjH9Wwn42Tpl+0pShNlYQNZcxveuiUPGmCMNydpvpSHEYBUkAQgmVsa2aBriTZR6J+0f3oF7WICZr0+p7GjqEuKYrfFoPeUZAEghAs41+STZPRG/VF2pKpIVhBaAFyQZYVNF6iKTM3RAp7W5kO260BbwhCsICOEbK1SZRoHf0hXUm0T6S3s5cUxeGSpodOyWvcL0VlNzlZHimsrUKjEDhDEIIFrK0Srw0XfJwlKZ7fuH966GRJSsMYTAqaUN5+om+0X4ri98QJfzshilhGAVwhCMEC/lEp3hMnyXPPyIz7mg5lh0ySojiMhYtcma5OlujM+qwA6u1MdrQgCYEnBCGY224tI4RMVUsyTeawrjTYI1DtppKiOIzRzPDsnacKJCpuahRKVBwcE4IQzO1vJ8T7E6R64u2sL5gVniNRcRijacGTStvKekf7pCj+s3HCjhZRy//oQ3BcCEIwq9PDZHOTKNHyQb1o2Nd0cEbYVCmKw9i5KlwyA9MLGg9IUdxdQa6PwDoK4AlBCGb190oJp8kcbDkS6R3u7+onSXW4HLPCc3bW50tU/L4E4d3johEDhcAJghDMxyCSNRXig4nS9YvmzwrPlqg4XJapwVmVHTVdwz1SFJ/gT4PdyPp6NAqBDwQhmM+X9WKEB8n0l2SazKhx9EDL4emh6Be1Ck4yp4lBGXmN+ySqvzpReLMcQQh8IAjBfN4sF1cnSfWU29d8KN4vxkeJdfTWYlZ4zs5TUvWOLo8UanpJaSe6R4EDBCGYSUkHq+0jy6Q5jJ4QsrV219yIXImKwxWYGJRxqqdR298qRXG5QFbFC2+hUQg8IAjBTN4oFx9IEBTSPOO6hnuOtpdPx3xRa6IQ5DPDs789tVui+vfGC5+fEjtGJCoPDgRBCOagGyIb6sV746V6vm0/tSc7ZJKLXClRfbgyC6JmbandyYgkHZgqF3JDhPB2BRqFcLUQhGAOrx0z3hIj+EuWU1trdy6Imi1VdbhS8X4xTjKn8vZKier/KlV4q9w4YJCoPDgKBCFIrldP3q8SH0mW6slW19PQO9qfpkqWqD5cjXmRM7bW7ZSoeKwXzdEI7+M8Crg6CEKQ3JoK8ZpQIdxdklUThJDNJ7fPj5wpUKnqw9WYFzlzd/3eEeOoRPWfSBP+dEw0IArhKiAIQVojRvJmufhYilTPNINo3H5qz7zIGRLVh6vk7+Ib5zeuoLFQovpZATTCnayrQxLClUMQgrT+WSVO8KcpvlI11wqaCkM9Q0I9gyWqD1dv8bh5G6q3SFf/iTTZiyU4pBCuHIIQJDQqkpePis+Ml/BptqF6y9KYBdLVh6uXEzq5qa+lrqdBovrzQ6i3M/kMjUK4UghCkNDfToipvjQrQKrmYEu/rrb71PQQnEdv1WRUtjB6zqaabdJd4pnxst8dQaMQrhCCEKQyYiR/LBWfTpfwOba+avPCqDkKmUK6SwAXS8bN31q3a9gg1er3ecHUT0k+RaMQrgiCEKTyt0pxgr+EzUG9Ub+1bteicXMlqg8cqd0CEv3idjfsle4Sz4yX/R6NQrgiCEKQRL+evFgiPpch4RNsV8PecT6RIR5B0l0COLo2dsH66m+kqz83mPoryb9r0CiEy4YgBEm8esw4J5im+0m4tu/zyq+Xxy+Rrj7wNTkoq2ekr6z9hHSXeGWS7OkicRAbzcBlQhACf21DZE2F+PsJEj67iluPDeoHJwZOkO4SwJdA6fK4petOfCXdJSYG0CkqinMK4XIhCIG/3x4x3hEr4VYyhJB1x9fflLAMu8nYlmui5xxtK2/p10l3iT9OFF4rM+JICrgsCELgrKKbfXFKfCJNJt0lGnubj3dUzY2cKd0lQApKufOi6LmfV34t3SWiPOhPIoXfHzFKdwmwPwhC4Gz1PuMz42U+zhJe4tMTG66NWegsc5LwGiCN6+IWb63d1TvaJ90lfjdBtq5WxOH1MHYIQuDpv7Vi2xC5T7JzBwkhHUNdu+oLlsUulO4SIB1/F99pIRO/qpJw+qivM3luguzBvUYkIYwRghC4GTSQxw+Kb02VyaV8Wn1c8fn8qJk+Sm8JrwFSujX5pi8qNw7oB6W7xD1xgl4kH2IpBYwNghC4ee6IMTeQ5gZKOIGla7h7a+2uFQnXSXcJkFqwR2CmZvyXlZuku4RAyZtTZY8fFLulOv0J7AqCEPgoOs3+XS2+MlHCOTKEkI8qvpgbOcPf1U/Sq4DUbk9d8VnlhkH9kHSXmBhAl0XQxw5g1gxcGoIQODCI5N5842uTZSoXCa/SM9K7pXbHysTrJbwGmEWIR1CGJu0rKTeaIYS8PFG2W8u2NmGsEC4BQQgc/KHEGOxGfhot7dNp7bFP5qE5aC9uT1mx7vhXfaP90l3CTU7ey5bdt9fYp5fuImAPEIRwtY6cZu8cF/+aLW2naEu/bld9/q3JN0l6FTCbMM+Q6aFT/1P2qaRXmRVEgRGteQAAFkZJREFUZwXSX6GDFH4UghCuyoCB3Lzb+PpkWZCrtJu8rDny/k0J13k5e0p6FTCnO1JXbq7dIelGM4SQv0yR7dSydbWYQQoXhSCEq7J6n3GKiq6UuFO0/PSJ6s6TN8QtlvQqYGY+Sq/l8Uv+XvJvSa/iriAfzpCt3m9s6MdgIVwYghCu3H9rxb2t7I0p0naKioy9UfS3e9JvccJWMnbnJ/HLjrYfl/RICkJIVgB9JFl22x4j1tjDBSEI4QqVd7GH9xs/mSVzl/h8+K+qvnGWOc2OmC7tZcASlHLnBzPufPXg2wZR2mG8x1MFhUB+cwiDhXABCEK4Er16cuMO4x+zZOOlPHGQENI51PVB2Se/nPgAJThowj7NDM8OdFNLejwTIUSg5JNZ8s9PsU9OYrAQzocghMsmMrJyp2FeCL0jVvLnz1+K3rs25ppwr1CpLwQW9EjWqv9WfCX1rBlfZ/L5HNnD+41HsR83fB+CEC7bI4XGYSN5VeJNZAgh+Y2FJ7tP3Zx0g9QXAstSuwWsSLzuTwfXMCJtRKX50remypZsMzYNIAvhOwhCuDx/LBV3tbDP5sgl3VmbENI13P3nor8+OeXnmCPjCG5KuM5gNHx2QsKjCk1+EiWsThIWbDF24fBe+B8EIVyGf9eI7x4XtyyQeUucTYywP+5/Y+m4BUn+8dJeCayDQOlT0x79qOLzmq5aqa/1WIowJ5hev90wZJD6UmAbEIQwVp+cFJ84KG5eIAt2k3zeyqfH1w/oB29JvlHqC4H1ULn6P5Bxx/N7/zRilPzMiNcmyYLd6HXbDcOYRgoIQhijz+vERw8YtyyQJXhLnoIlbWUfH//y6WmPChTPT8cyN2JGvF/My4VvSj1YKFDyr1yZv5Jeuw1ZCAhCGIN/VYur9xu3LJCn+Eqegq0Dbc8XvPrM1F9q3FRSXwus0C8nPtDSr/244gupLySjZO10mbczXbrN0I9duR0bghAu4S9l4u+OiDuvkadKn4JDhuHf7HlhZdLyDE2q1NcC6+Qkc3o+58nPKzfuby6S+lpygXw0UxbtSWdsMrRJeDYiWDsEIVyUXiQP7jP+o0osWCKLl75HVG/UP7XnhUT/OOwp6uD8Xf2ez3nypcK/SL31GiFERsk702SLw+i0rw0V3VhT4aAQhHBhbUNkzjeGhn6Wv1gu9ckShBCRiX/Y95qbwvWRrFVSXwusX6J/7DPTHns674WqzpNmuNxzGbJnxgszNhq+PIV9ZxwRghAuIF/HstYbcgPp+rlyL+lX8RlE4x/2vTZiHPlt9q9lVPJ1+mATJmjSHp34wBO7f1/bXW+Gy90aI2xeIP9FofjEIaMeaehgEITwPXqRPF1kXLHT+M402e8nyATpN/gcNow8teeFIf3w73KekAtIQfjO9NApD2fe++iOZ462lZvhchP8adEy+fFuMnWDoaoH3aQOBEEI3yk6zSatN5R0sOLr5NeEmmOT656R3kd3POOj9PpD7pPO2EEGfmBG2LRnpv3y2fw/5jfuN8Pl/JVk/VzZ3fFC9teGV4+JBjQNHQOCEAghpHuUPFJoXLrN8MsUYeN8ucrFHBetOF11z+ZHMwPTHp/yMHpE4WImaNJenvncm4f//vfS/4jMHNG0Kl4ovFa+vVnMWm/Y14qmof1DEDq6UZG8dZzGfaofNpCyGxQ3jzPHU4IR9nnl17/Z84dHsu69M/VmHLEEPy7WN/q9ha9VnK58bOdvTw91muGKUR50ywL5r1OFn+4y3rjDWNuPp6g9QxA6rmEjWVMhJn1Jt7ewndfI382W+Tqb47rNfdpHvn1qZ33+mvkvTw2eaI5Lgu3zdvZ6ddbvUlWJd3/z880nt5vnoj+NFk4sl2f409nbne7IM2Lg0F7JLX0DYAHaQfLeCfGvJ4xZAcJHuSzLnzg5meMN75Bh+JOKL7+s2nRr8k3Xxy0WKN5lw2UQqHB7yk9zQqe8tP+Nbad2P5BxZ4xPlNQXdZGTJ9OEW0JH3693zdloyFYLDyUJMwPx1LUrCEIHohfJlibxX9VsV4t4U5Sw/Rp5ojcdHTXHTosjxtHNJ7f/u2xdhib1bwtfV7sFmOGiYJeivSPeWfDKxpptj+/63QRN+q3JPwn1DJb6ol4K9sx44dEU4d/V4up9RsbI7bHCzeOoGZbYghkgCO3fiJHsaGFfnBI31Ivx3vTWGOGf0xWeCjNdvWOoa2PN1q+qNyf6xb4445lY32gzXRjsl4zKro1ZOC9y5rrj61d/+2RyQPzy+KVpqiSpB5vd5OS+BOG+BKFAxz6oEVM+N6b50usihGURNFT6I1lAOvI77rgjISFh9erVLi7fmym4bt26zZs3q9Xqhx56KCQkxFL3B1fGyEhJB8vXsW3N4l4dS/ej10UIz4yXh7ub6eU6oB880HJ4W92usvYTM8Oz/zLn/8Kkf9sODsVFrrwt5aYViddtqd3x50N/HTGMLIiaPTN8Wpin5H+vsjU0WyN7Y4psW5P4ZT37fbFR40LnhdAZgXSaWjDPWDtwRD/++OO1a9e6urp+8cV3272/++67L7300vPPP3/o0KGvvvrq+PHjrq6uF/z53Nzc3//+97m5uWO8Xn9/v5ubG8XgEG8GkdT0smNdrKidHT7NDrWzUHeao6Fzg+nsIOFHdocZHR0lhDg5cVjDJzJW01Vb0nrskLak/PSJlIDE2RE500OnKuX4w3AZGGMDAwPu7u6WvhEbU9lRs6Vu596mA04ypylBmeM1KWmqZDfFhf9wXa6+vj4PD4+LfVVk5EgH29rE9mjFwjYW7k6zAugEf5rmS5N9qdSnWDumH39ELhdljPX29mo0muLi4ri4OEIIYywmJua1115bunQpIWTy5Mn33nvvnXfeecGfRxCaX6+eNA2whn5S38/q+lh1D6nuZTW9LNiVpvnR8X40059mBdAxvi29miAc0A829Dad6mms666v7DxZ3XkywNVvvDo1Q5OaFTjeRa68gpqAILxKVZ0nD2qLS1qPlZ8+oXL1j/UdF+MTGeEdFu4ZEuAacGVTtMb+Z9cgkrPvR0s6WUUX83KisV4kxpNGe9IIDxLuTkPciMaFmmHbJjvGNwjlhBBPT8/x48cXFhaagrC1tfXkyZMzZ840fcfMmTP37dt3sSCEqzRkIMNGMmwkQ0Y2oCfDRtIzSvoNrE9PekdJ9yjpHmEdI6RjhHQMs9Yhoh1kAiUhbjTEjUR40Ah3uiKaxHgKsV7UleuAb//owIhxdNgw3D860K8f6B3p6xnp7Rnp7Rzubh88fXqwU9vfamCGUI/gcK/QSO+wW5N/Euc7zt3JjedNAFy+WN/oWN/onyUtNzLjqZ7Gyo6ak911B7RHGnqauoZ71G4BarcAf1c/fxdfb2dPL2dPT2dPdyc3dyc3F7nSTeGqlCsVwpW/luQCGe9Hx/vRewghhDBC6vtYdS+p6mF1fexAO2noF5sGWMcw8VfSACXRuBI/Z+qnJN5OxNuJejoRLyfiqaBucuIiJ95OREaJlxOllKBlKZ0zj3dAQIBWqzV9rNPpnJ2dz4ZtQEDA0aNHL/bzAwMDTz31lL+//7mfXLx48cqVKy/4/QVlRV8ce5/DjZudOOYVRD/8xnN/9uzHph0y6Nn/o9/7X+HM/zKBUnfKPAiJokRGiUCJQAgZIGSAkDbSS8hhQg5f6DYGBCO7yD2LhA3TMzt0MEaGBFEkjBAyQIyEsGFqNBDiygQnInNm1J3IXZnMncg9RJknkfsxRSxz8mUKFYv0YHLSQ0hTDyHHCDk2SMjgWP8jwUUxxoxG46AcE9k48CAkk5BMQggRCAkzENY6ONLRPthFuzuEylNU30cN/cQ4QA0DxDhCxCFqPPv8Fwh1JjI5oYQQpSjI/9eUFAhVsouuwFYSQc4u3NaTEzKOkHGEEEKYjIijRBwlxh4iMiIy0kFoOyOMEZEQxggjhP3vz8WZ/z1z9TPoOVODTH80zjPGBie1tf0sVuX8KiEsamhoSCYb03ZUSqVSEC6xYv7Mi21kZESpVJ79Mb1eL4ri/7d3/zFNnH0AwJ8eRRARUUfDz2FFxCGgNBpEbGWACqmIVCe6LCaDdxp/xGiCWbb5Y4sajJMsZguTgBjlh8YqGioatWhAaNWJjjjG8AdICkULqdDyo/WuvfePx/foW+NrmbTH3vt+/rp7nuPumx7Xb++e554H/7HJZLLrR2NrwoQJy5Ytmzt3rm1heHg4szc7wYLgaP8lDp+j8YKHkIdjQ4DxeGjCW5+5BzHyb8rsx7bQGbwJz3ftn+ARE4k33UYtlGUiMWEC3x0h5O3myUPIg3B3hwHP2EPTtO31CMaWr2ObDVrMNKKHrSRFWxBCfcMGvsebOzILbTVZ3zmlPfMnTkJakYV+s8D8qmYKGRYrohz77U7R6J81pGqQIAAnKQevkfdmQcQkQo1Gw3QNDQwMpGlaq9XiEo1GExT0zv5+7u7uiYmJjrcRzhD4bxf+C9oIx48x7CwDxgRN0/TgoBe0EbLKro/N2LZIgQ9HEIQjGc7RvSGEHj582NHRsXz5cr1eX1lZ6ePjk5ycXFZWhhAyGAwKhUImk43V8QAAAIBxhS+TyRoaGvLy8nx8fFQq1Zo1a2iazsvLk0qlDQ0Nra2tEolEIpGwHScAAADgFPxNmzbl5+cLhUKEkEgkevLkCUJowYIFra2td+/eFQgEsbGxbAcJAAAAOAuRmpqKsyBCyNPTc9Ys3KcJ+fr6rlixYsyzoEajGRgYGNt9gg/R09PT09PDdhRgxMDAgEajYTsKMMJiseA7BDB+tLW14f4NY8LV0zDt2LGjvr7exQcF/8Ovv/5aWFjIdhRgRF1d3c6dO9mOAozo6+tLSUlhOwrwX2QyWUdHx1jtDeYjBAAAwGmQCAEAAHDah45eQVFUU1OT49u/evXq0aNH7xrCG7heR0cHQRC1tbVsBwLe+OOPP/R6PZyR8cNgMFAUBWdkXDGZTPfu3dNqte/dcsGCBZMmvWfoRx79rmG4HHPmzJlffvmF7/BwUAMDA56eno5vD5zNZDIhhGAck/GDoiiTyQSDbo8feGaCKVOmsB0IGGEwGLy9vR15p/7EiRNMJ9B3+dBECAAAAPyjQRshAAAAToNECAAAgNMgEQIAAOA0SIQAAAA4zbm9N0mSvHnzpsFg+PTTT+0m70UI3blz5/nz5yKRaPbs2U4NAzAMBoNarTYajdHR0REREbZVjY0j8/v6+fl9/PHHLo+Oi5qbm3HHXYTQ5MmT7a6FZ8+e/fbbbyEhIQkJCWxEx0W2ZwQhNGXKFKbPYVdX14sXL5iq+fPnOzg3LPh7Xr582dnZGRYW5uv7Zh5Jo9FYU1NDEERycrLdSxEWi+XWrVt6vV4sFgcEBIzqQE7sNWo2m5OSkmiaDg0NVSqVNTU1MTExTO2uXbuqqqokEsmVK1d+/PHHjRs3OikMwPjrr7/i4uIWLlzo5+d3/fr1rVu3HjhwgKl1c3MTi8Xu7u4IIalUCqN8uUZkZOSkSZPwdT5v3ryjR48yVefOndu2bZtUKlWpVElJScePH2cvTA758ssvOzs78fKDBw9kMllRURFe/e6778rLy8PDw/HqxYsX4S0X54mKimpvbydJUi6XZ2RkIIS6u7vj4+Ojo6Mpinr69KlKpfLz88MbW63WtLS03t7eyMjIq1evVlVVLV68eBQHo52mtLR03rx5JEnSNL1nz561a9cyVW1tbV5eXl1dXTRNX7t2LTAwEG8GnOrVq1darRYvP3z4kMfj9fT0MLUEQdiuAtf45JNP6uvr3y63WCxCofDixYs0Tet0usmTJ//5558uj47ThoaGfH19b9++zZR8++23ubm5LIbEKS0tLSRJRkREXLp0CZd8/fXXWVlZeDkjI2P//v3MxpcvX545c+bQ0BBN0/n5+fgezHFObCOsqqrKzMzE785/9tlnCoWC/s/dZ3V1dXx8fGBgIEIoJSXFZDLdv3/feZEAzNfXl3liEBISQtO07SMghFBjY6NKperv72cjOu5qaWmpq6uzmwPk0aNHL1++lEqlCCE/P7/ExESFQsFSgBwll8sFAoHdQ+ne3t6amprHjx+zFRV3zJkzx27olaqqqnXr1uFlnFOYKoVCkZ6ePnHiRFx169Yto9Ho+LGcmAi7urqCg4PxckhIiNls7u3tfbuKIIjAwMCuri7nRQLelpeXl5KSwpwFhNC0adOOHDmya9euGTNmXLhwgcXYOMXLy6u0tHTv3r1CofCnn35iyrVarb+/P35SjRAKDg52ZDQpMIZKSkpycnJ4PB5TwuPxmpqaDh8+nJCQsHLlSrPZzGJ4HGSXU2yzhm1VUFAQj8cb1fXixM4yFEUxLcl4gSRJvEqSpO3QOHw+n6kCLlBSUiKXy+3mw9Jqtfhrt6KiIjs7OzU19b0D9IEPp1ar8cd+586dxMREqVSK+8uQJGnbEYPP54/h7Gvgvdra2lQq1ZkzZ2wL9+/ff/DgQYSQwWBYvHjxzz//nJuby1KAXGSXU2yvCNsqHo9HEMSocooT7wj9/f2Zpz06nc7NzU0gEODVgIAA2wdBOp1utJ18wN9WXl6+b9++GzduhISE2JYzNx/r1683mUwwE6lrMB/7okWLwsLCfv/9d7yKrxGmNUGn0+GmBOAaxcXFaWlpdt9LzMny8fFZvXq1bUdr4AK2icPuirBNN3q9nqKoUeUUJyZCiUSiVCrxslKpTEhI4PP5ZrPZZDJJJJL6+nrcQNXc3Gw0GkUikfMiAYzz58/v3r372rVrTDf9169fDw8P226D26iDgoLYCJC79Hq9RqMJCgqiKGpwcDAqKoogCNx2TpJkbW2tRCJhO0ausFgsZWVl2dnZTInRaLRarbbbNDU12bYsABeQSCQ3btzAy0qlcunSpQih4eFhkiRxusE/HJVKZVRU1PTp0x3fsxNfn9Dr9dHR0ZmZmWFhYQcPHiwrK0tLS9uyZcvg4ODp06dXrFjB4/EyMjIKCgpWrVp16NAhJ4UBGM3NzfPnz09NTZ07dy4u2bx5c2lp6e3bt7dt21ZZWSkSifr7+4uKitasWXPs2DF2o+WCpqamH374IT4+nqbpU6dOhYaGVldXnz179ptvvnn+/Dm+anbs2HH16tX+/v66ujq24+WK6urqnJwcjUaDbwFfv37t4eHR2Nj4/fffx8bGfvTRRzdv3mxoaGhsbLR7rALGUGFhYXt7e1FRUWJiYnh4+NatW/v7+5csWZKbm0tR1LFjx+7evRsRESEWi9PT07dv3x4TEyMWi2NjY/Py8vLz8z///HPHj+Xc2Se0Wu3Jkyf7+voyMzPxWx21tbUkSeKeosXFxW1tbXFxcevWrbNtkQZO0t3dbdfzMD09vbOzs7u7Oy4uTi6Xt7e3e3t7i8XilJQUtoLkFKPRKJfLW1pa3N3dRSKRTCYjCOLJkydqtRq/WXvhwoWGhobQ0NCvvvoKZvF0GbVaPTQ0lJycjFetVmtBQUFWVtb9+/fVavXAwMCsWbM2bNgwdepUduP8/6ZQKLq7u5nV1atXCwSC5ubmiooKgiC++OILPCRIZWWlUCiMjY3t6ekpKSnR6XRSqTQpKWlUx4JpmAAAAHAajDUKAACA0yARAgAA4DRIhAAAADjt34QGJgDhj9yzAAAAAElFTkSuQmCC", "text/html": [ "\n", "\n" ], "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "execution_count": 10 } ], "cell_type": "code", "source": [ "using Plots\n", "\n", "p = plot(x, real.(ψ), label=\"real(ψ)\")\n", "plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n", "plot!(p, x, ρ, label=\"ρ\")" ], "metadata": {}, "execution_count": 10 }, { "cell_type": "markdown", "source": [ "The energy_hamiltonian function can be used to get the energy and\n", "effective Hamiltonian (derivative of the energy with respect to the density matrix)\n", "of a particular state (ψ, occupation).\n", "The density ρ associated to this state is precomputed\n", "and passed to the routine as an optimization." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "E, ham = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)\n", "@assert E.total == scfres.energies.total" ], "metadata": {}, "execution_count": 11 }, { "cell_type": "markdown", "source": [ "Now the Hamiltonian contains all the blocks corresponding to kpoints. Here, we just have one kpoint:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "H = ham.blocks[1];" ], "metadata": {}, "execution_count": 12 }, { "cell_type": "markdown", "source": [ "H can be used as a linear operator (efficiently using FFTs), or converted to a dense matrix:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "ψ11 = scfres.ψ[1][:, 1] # first kpoint, first eigenvector\n", "Hmat = Array(H) # This is now just a plain Julia matrix,\n", "# which we can compute and store in this simple 1D example\n", "@assert norm(Hmat * ψ11 - H * ψ11) < 1e-10" ], "metadata": {}, "execution_count": 13 }, { "cell_type": "markdown", "source": [ "Let's check that ψ11 is indeed an eigenstate:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "1.5515748768512322e-7" }, "metadata": {}, "execution_count": 14 } ], "cell_type": "code", "source": [ "norm(H * ψ11 - dot(ψ11, H * ψ11) * ψ11)" ], "metadata": {}, "execution_count": 14 }, { "cell_type": "markdown", "source": [ "Build a finite-differences version of the GPE operator H, as a sanity check:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "0.00022341923778219596" }, "metadata": {}, "execution_count": 15 } ], "cell_type": "code", "source": [ "A = Array(Tridiagonal(-ones(N - 1), 2ones(N), -ones(N - 1)))\n", "A[1, end] = A[end, 1] = -1\n", "K = A / dx^2 / 2\n", "V = Diagonal(pot.(x) + C .* α .* (ρ.^(α-1)))\n", "H_findiff = K + V;\n", "maximum(abs.(H_findiff*ψ - (dot(ψ, H_findiff*ψ) / dot(ψ, ψ)) * ψ))" ], "metadata": {}, "execution_count": 15 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.3" }, "kernelspec": { "name": "julia-1.5", "display_name": "Julia 1.5.3", "language": "julia" } }, "nbformat": 4 }