{
"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 external magnetic field)."
],
"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 `LocalNonlinearity(f)` term of DFTK, with f(ρ) = C ρ^α.\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",
" LocalNonlinearity(ρ -> C * ρ^α),\n",
"]\n",
"model = Model(lattice; n_electrons, terms, spin_polarization=:spinless); # 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.924889e+02 1.998694e+02\n",
" * time: 0.0005788803100585938\n",
" 1 1.753189e+02 1.639821e+02\n",
" * time: 0.0019469261169433594\n",
" 2 1.335265e+02 1.860744e+02\n",
" * time: 0.004890918731689453\n",
" 3 1.301303e+02 1.693085e+02\n",
" * time: 0.008224010467529297\n",
" 4 8.979427e+01 1.684478e+02\n",
" * time: 0.011606931686401367\n",
" 5 4.492616e+01 1.025466e+02\n",
" * time: 0.014878988265991211\n",
" 6 1.785666e+01 5.936727e+01\n",
" * time: 0.0177609920501709\n",
" 7 1.302049e+01 4.164894e+01\n",
" * time: 0.019841909408569336\n",
" 8 1.002173e+01 1.376787e+01\n",
" * time: 0.021733999252319336\n",
" 9 4.911450e+00 3.058122e+01\n",
" * time: 0.02381587028503418\n",
" 10 3.056549e+00 3.091546e+00\n",
" * time: 0.08835101127624512\n",
" 11 2.192968e+00 2.363784e+00\n",
" * time: 0.09001684188842773\n",
" 12 1.479981e+00 3.439475e+00\n",
" * time: 0.09158992767333984\n",
" 13 1.366366e+00 4.329482e+00\n",
" * time: 0.09291481971740723\n",
" 14 1.311442e+00 3.411253e+00\n",
" * time: 0.09414982795715332\n",
" 15 1.190634e+00 8.936627e-01\n",
" * time: 0.09536194801330566\n",
" 16 1.158207e+00 8.559547e-01\n",
" * time: 0.0965878963470459\n",
" 17 1.157925e+00 4.627043e-01\n",
" * time: 0.09747886657714844\n",
" 18 1.150702e+00 2.518736e-01\n",
" * time: 0.09869384765625\n",
" 19 1.145582e+00 1.206243e-01\n",
" * time: 0.09995889663696289\n",
" 20 1.144415e+00 4.100512e-02\n",
" * time: 0.10120487213134766\n",
" 21 1.144236e+00 2.809725e-02\n",
" * time: 0.10332393646240234\n",
" 22 1.144116e+00 2.785552e-02\n",
" * time: 0.10426688194274902\n",
" 23 1.144092e+00 5.356267e-02\n",
" * time: 0.10560488700866699\n",
" 24 1.144050e+00 1.200704e-02\n",
" * time: 0.10647797584533691\n",
" 25 1.144048e+00 1.041670e-02\n",
" * time: 0.1076970100402832\n",
" 26 1.144039e+00 1.188159e-02\n",
" * time: 0.10891890525817871\n",
" 27 1.144038e+00 5.399591e-03\n",
" * time: 0.11018085479736328\n",
" 28 1.144037e+00 3.148366e-03\n",
" * time: 0.1114039421081543\n",
" 29 1.144037e+00 6.557557e-04\n",
" * time: 0.11276102066040039\n",
" 30 1.144037e+00 5.044121e-04\n",
" * time: 0.11402487754821777\n",
" 31 1.144037e+00 3.123932e-04\n",
" * time: 0.11527895927429199\n",
" 32 1.144037e+00 3.733782e-04\n",
" * time: 0.11652779579162598\n",
" 33 1.144037e+00 2.605472e-04\n",
" * time: 0.11780691146850586\n",
" 34 1.144037e+00 1.952548e-04\n",
" * time: 0.11910796165466309\n",
" 35 1.144037e+00 1.004852e-04\n",
" * time: 0.12041401863098145\n",
" 36 1.144037e+00 4.559759e-05\n",
" * time: 0.12174797058105469\n",
" 37 1.144037e+00 4.029865e-05\n",
" * time: 0.12320685386657715\n",
" 38 1.144037e+00 4.045215e-05\n",
" * time: 0.12451982498168945\n",
" 39 1.144037e+00 1.970928e-05\n",
" * time: 0.12584590911865234\n",
" 40 1.144037e+00 1.326918e-05\n",
" * time: 0.12715888023376465\n",
" 41 1.144037e+00 4.752105e-06\n",
" * time: 0.12846589088439941\n",
" 42 1.144037e+00 2.917656e-06\n",
" * time: 0.12978696823120117\n",
" 43 1.144037e+00 1.570055e-06\n",
" * time: 0.1311349868774414\n",
" 44 1.144037e+00 8.902491e-07\n",
" * time: 0.13243699073791504\n",
" 45 1.144037e+00 5.143540e-07\n",
" * time: 0.13385891914367676\n",
" 46 1.144037e+00 4.679790e-07\n",
" * time: 0.13517189025878906\n",
" 47 1.144037e+00 3.188825e-07\n",
" * time: 0.13647985458374023\n",
" 48 1.144037e+00 2.872447e-07\n",
" * time: 0.13780593872070312\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown (in Ha):\n Kinetic 0.2682057 \n ExternalFromReal 0.4707475 \n LocalNonlinearity 0.4050836 \n\n total 1.144036852755 "
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"basis = PlaneWaveBasis(model, Ecut=500, 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.ρ)[:, 1, 1, 1] # converged density, first spin component\n",
"ψ_fourier = scfres.ψ[1][:, 1]; # first k-point, 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": [
"ψ = ifft(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.0508545227560829e-15"
},
"metadata": {},
"execution_count": 9
}
],
"cell_type": "code",
"source": [
"norm(scfres.ρ - 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+gvaeTAAAgAElEQVR4nOzdd3wVVfo/8Gfm3pve600vBAjphYSSBAgdUYS1rIpl13V1i6uufWXX+l137euq2F1XFBELKtKkQygBEkgngYSQenPTe7n3zvn9EX9ZpCZw5tbP+68kTJ6ZvJjJJ+ecOecIjDECAACwVaKpLwAAAMCUEIQAAGDTEIQAAGDTEIQAAGDTEIQAAGDTEIQAAGDTEIQAAGDTEIQAAGDTEIQAAGDTEIQAAGDTTBCEhYWFH3/88eiP1+l0sl0LjJZerzf1JQDp9XqsiWhyBoNBkiRTX4WtkyTJYDDwqmaCICwuLt66devojx8YGJDvYmCUBgYG8CvY5HQ6HX4Fm5xer+f4Kxguj8Fg4NhGQtcoAADYNAQhAADYNAQhAADYNAQhAADYNAQhAADYNAQhAADYNAQhAADYNAQhAADYNKWpLwAALoERFbaxik5W0SY627EodxbvRWEugqmvC8BKIAgBzFfbIH1QLn1ULjGiOE8h2JH0A7S1wZDXwqI9hLsmijdFikp06wBcGQQhgJlaWyU9cNCwMFj8eKZiqp9ARP39Q3Z2SoVCoZPohxrpjRLp1SLpwxmKZG+0Di1Aa2vrvffei+XZLtsjjzySlpYmR2UEIYDZ6dXTHbsN5R1s3VzlFL/zhJxKpGXh4rJw8ZMT0qLN+vtiFU8koWFo7hobG/ft2/fKK6+Y+kIs0jvvvFNYWIggBLAJrYO0eIs+1kPIW6a0u1S63T5enB8sXr1FX9PD3spQKNAyNG/u7u433HCDqa/CIm3ZskW+4vgrEsCMNPSxzPX62QHCBzMUl0zBYWpH2rlYWdnNfrnDoMfWFABjhyAEMBc9Orp6i+G2KPH5tLE17VxVtGGBckDP/rgf408AY4YgBDALEqNbdxlSfITLG+2zE2ntHOWxVvZiIVqFAGODIAQwCw/nGvr17J0MxWVXcFLSN3MVb5VK351GFgKMAYIQwPQ21bJvqtkXc5RXOCkwyFlYO1txT46hvpdxujQA64cgBDCxlgH6bY7hPzMVHnYcqk3xE/4Qo/j1HgOSEGCUEIQAJvabvYZbo4TsAG5TH55IFHv19K9idJCCjH7729/+97//Hf64q6srPT29r6/vQgcbDIaMjIyGhgZjXd3YIAgBTOmLKul0N3s29fKHBs+lFGnVLMXzxwzV3WgWgly6urr6+/uHP3755ZevvvpqJyenCx2sUChuu+22Z5991lhXNzYIQgCT6dHRw7nSv6ePdsrg6EW6Cg/EKf58EI1CuIQ1a9bU1dW9++67zz33HBENDg5+8sknTz311Jo1ayTpp/snPz//5Zdf/tvf/nbmF0fodLr33nvv9ttvJyKtVrtmzRoiysvLO3DgABF98803dXV1RHTTTTetXbu2s7PTmD/dKCEIAUzm2aOGuUHCDLUs68E8kiCWdbANtWgUwsWsWLHi6quvrq2t9fDwGBwczMjI2LVrV1hY2OrVq2+99dbhY1auXCmKYkRExAcffHDXXXedVSEvL8/JySk8PJyIqqur//rXvxLR+vXr165dS0TPP/98aWkpEXl4eEycOHH37t3G/OlGCUusAZhGWQf7uEIquk4lU307kV6fprh3v2FOoNKBZ88rcKOX6BfbDP1GfLHpD5PEZeFnt39+9atfPfDAA0S0cuXKwMDAjz76iIhuu+22yMjIioqKCRMmfPDBB8NH3nzzzd7e3itXrnRwcBj59qKiovHjx4/m7BMmTCgsLFyyZAmfH4YfBCGAaTx2SFqRpPB3lPEUC4KFOE9hZan0YDz6fsyRUqTHEsV+vfHOmOB1nu6HkZWs8/LySkpK5s2bN/xpb29veXn5hAkTXn/99ZUrVwqC4OLiotPpGhoaIiMjR769p6fH2dl5NGd3cXHp7u6+4h+CPwQhgAkcamZF7ezLObK31P6ZLs76QX9XtOgmV8sTrkiGv+kXSre3tx/+wNHR8dprr/3b3/428k/Ozs6FhYUvvvhiQUGBj4+PwWBwdnY+aycpHx+ftra24Y8dHBwGBgbO/Nf+/v6R5mNbW1tERISMP8nlwt+JACbwl8OGp5JFe/l7LCe6CwuCRUylgNFYuHDht99+S0Senp6enp6CIAiC0NTU5OHh4e3tTURr1qwZHBw867vS0tKKiooYY0QUFRXV3d2dn58//E/Hjx+vqamZNGnS8KfHjh1LT0833s8zamgRAhjbj/WsrpdujTLSn6HPpIqp6/S/ixb95OyGBStw9dVX5+XlxcTEpKend3V1VVZWFhYWZmVlOTk5TZkyxdPT08nJycvL66zvio6ODgoKOnz4cHp6upOT0wcffHDNNdeIomgwGNauXfvGG2/4+voS0enTp7u6uqZNm2aKn+wSEIQAxva3I4a/TxavcDW10QtzEW4aJ75cZHgxHe/MwNny8/NdXFxGPn3qqaceeOCB8vJyFxeXiRMnKhQKItq/f39hYaGjo2NMTExHR4ebmxsRvf/++3Z2Py2G9OCDD7733nvDrb0bbrhh2bJl999/f3d394cffqhS/dQp/9577913330jn5oVdI0CGNW2etarp+sijProPZ4oflQutZ/dpwVA7u7uw2l35lfS09NjYmJGvq5SqVJTU2NiYojIw8NDFEUicnNzGxn8u/322x0cHEZWllEqlb6+vt7e3iOxZzAYOjo67rvvPuP8UGOFFiGAUf2zwPBYomjkFySCnYUlYeLKMmnFZe3xBHBxoii++eabZ34lMzNzZN0ZIlIoFG+99ZbRr2u0EIQAxnOslVV00i8jTZBGjyWKs37QPxgnOuKhB/nNnTvX1JcwBvjzEMB4/n5Meihe5L6g2mhMdBem+Ikfn8Dro/AzsbGxR48elaOyTqebPn16U1PTRY65/vrrDx06JMfZxwRBCGAkJ7vY7kbpNxNN9tA9lii+UiRhfyY409NPPx0aGipH5Q8//DA+Pt7f3/8ix/z+97//y1/+IsfZxwRBCGAkb5ZKv40WXUz30tw0P8HPgX6oQaMQ/qenp2d4gvymTZtyc3NXr179wAMPfPXVV0S0a9euBx988N133x1ZaHvr1q1PPfXUQw899MUXXwxPHCQixtiqVav+/Oc/r1mzJicnZ/v27cNff+edd4ZXK+3t7X3xxReJqKioaN26dUS0bt26goICIsrOzq6oqDh+/Lixf+yfQxACGEOPjladkO6ONvET98cY8a1SBCH8z3PPPVdbW0tEX3/99fLlywsKCmJjY++999577rnno48+SkhIePvtt1944YXhgz/77LOwsLCkpKQ33njjscceG/7igw8+uHLlysmTJx88ePCWW2754YcfiKihoaGiomLq1KlE1NPT88QTTxBRXl7eqlWriOjTTz/Nzc0lIlEUs7KyNm3aZIof/X8wbg5gDJ+elGYFiGEuJl5P64ZI8eFcw/EOFu1h+pW9gCRD88rH2eDApY/kQ3Cde4NjYtaF/nnatGnDmVdbW/vNN98UFRUJguDh4fHSSy8Nd2B+/PHHRNTb25uSkjJz5swXX3yxs7Pz7bffrqysDAoKWr58eUlJyXCpoqKi8PDw0cwanDBhQnFxMZcf77IhCAGM4e0y6bWppp/PbifSXdHi22XS69NMfzFAosLjuj8yvc5oJ1T6BFzkX+Pi4oY/8Pf3j42NFQRh+OPW1lYikiTpkUce+eqrr7y8vERRbG1tHRgYqKqq8vT0DAoKGv7GlJSUoaEhIurt7XV0HNVSRk5OTj09PVfyQ105BCGA7HY3Mp1E2YFm0Qj7/SQx/mv9c5MVWIbbHKgCwk19Cf8zPFN+2Fmz7Ilo48aNO3bsqKiosLe312g0AQEBkiR5eHh0dXUZDIbh49vb24d3ovDz8xuOTyJydHSUJOnMRUr7+vpGtrNvaWm5+As1RoAxQgDZrSyT/hhj7En0FxLoJGQHiqtPYqQQxqatrc3Dw2N4q4r33ntv+IthYWHjxo17++23iaiqquqbb74Z/npycrJGo2lvbyciNze3qKio4eW8iai+vv7QoUPJycnDnxYWFpp8JW4EIYC8WgfpxzppubGW2B6Nu6PFjyoQhDA2S5Ys0Wq16enpU6dO1Wg0w18URXH16tUffvhhYGDgrbfeumjRouEeUWdn50WLFo28BfPxxx8/+eSTK1as2LlzZ1pa2pNPPhkbG0tEXV1dR44cWbx4sal+qGHoGgWQ16oT0tWhooedqa/jDHMDhbv7qaCNJZ5vm1awKSUlJcONvH//+98jXaN33XXXr371q+GP09LShue8e3h4FBQUlJeXe3p6BgYG/v3vfx/OvLi4uKNHj0qSJIri/Pnz58yZM/yNDz300IoVK2655RYimj59enl5+WuvvbZp06aNGzcqlT9Fz6effnrLLbd4enoa9Wc+B4IQQF4fV0ivmdmbKaJAd4wXPq4wi/d3wLRGXmkZGbQjInt7+5HdepVK5fB2E8MfD7fkiGgkvT7//PPi4uLg4OCcnJzTp0/feOONw1/PzMxctGiRRqNRq9Uj3+Li4jKSgkTU1tb25JNPyvWzjRqCEEBGR1pYl45mqs2u4fXrCWL6d/p/pimMsDkwWLeZM2d2d3e3tLQsWrTo/fffPzNQH3/88TOPjI+PP+sdnL/+9a9GusqLQhACyOg/FdKvJ4jm8p7MGcJdhXgvYX2NdL1xN4QC6xMYGHj33XeP5sjU1NTU1FS5r+cy4BkAkMuAgb6olO4Yb34xSEREv5mIV2YAiBCEAPJZXyMl+wihpl5N5kKWhYkHtUzTf+kjAawbghBALqtPsuXjzPcRc1TSNaHiF5VoFNq6goKC+++//4EHHtizZ4+pr8U0zPcpBbBo7YO0s1FaFm7Wj9jyKHE1gtC2HTlyZNGiRfHx8YmJiTfeeOPWrVtNfUUmgJdlAGTx1SlpfpDobk7TB881J1C4o4dVdLIJ7mbaf2vdDMxw/9YnhgzGW2v0xknXzg2feeZXXnzxxccee+yuu+4iooGBgZdeemnevHlGux4zgSAEkMVnldIDcWbdHCQihUA3RoqfV7KnUhCEJqAQFI9MuXfAMHjpQ3kQSAh1Czrri6WlpX/605+GP05NTf3HP/5hnIsxKwhCAP4a+lhxG1sUbO5BSETLo8TlOw1PpVjApVqlMPcQ015Af39/d3f38MddXV0jc+dtCu5+AP5WV7Jl4aJFzFVP9xUEgQ43M1NfCJjMRx99JEkSY+zDDz+0wX5RQosQQA5fVEovpFtCDBIR0S8jhS9PSWm+FnPBwJeXl1diYqJer/f29n7zzTdNfTkmgBYhAGfV3ex0D5thfsuqXcj1EeKXpxiahDbr7rvv3rt377Zt23Jycry9vU19OSaAIATg7MtTbFm4qLScZyvRS7AXKb8FUWi7PDw8RnaZt0GW87ACWIivTlneAp6/CBe+OoUJhbbo0UcfDQ4ONvVVmJiFPa4AZq62l1V1s+wAi+kXHXZ9hLi2Ci1CW3TPPfeMbJNksxCEADx9WcWWhllSv+iwFB9BEOhYK7IQbJGlPa8A5u3rauk6S+sXHXZduPAlekfBJlnkEwtgnhr76HgHmxNoYf2iw66LEL+tRosQbBGCEICb72ukq0JElWU+VWm+QpeOKjqRhWBzMKEegJvvTkt3TrDMGCQSiK4JFb47zR5JsMgWrfmzs7OrqKgYN26cqS/EIjU3N8+cOfPSx10WBCEAHz062t/EvphtqUFIRNeGic8dNTySYME/gjmbMGFCRUWFwWAw9YVYqtDQUJkqIwgB+NhYK2X4C64qU1/HFZgdKNyyk2n6Se1o6kuxUmFhYaa+BDgP/OkHwMd3p9m1YZb9QKlEmh8sbqjBu6NgWyz7uQUwEzqJttRJ14Ra/AN1bajw3Wm8LwO2xeKfWwBzsEfDxrsLAU6mvo4rtihE3KORevWmvg4AI0IQAnCwvsYamoNE5G5Hk32E7fXoHQUbYg2PLoDJbaxli0OtZNbB4lBxQy16R8GGIAgBrlR5J+vVUYKXlQThklBhfY2EJATbMYbpE1VVVT/++KOnp+eSJUscHc//enV+fv7BgwddXV2zs7OxtQfYiA017OpQwUpikGicm+CiEgpaWZK31fxMABcz2hZhTk5OampqSUnJhx9+mJWVNTg4eO4xDz/88JIlSwoLC7du3free+9xvU4A87WhVlocYlWZsThEQO8o2I7RtgifeeaZFStWPPzwwwaDYfLkyV999dXy5cvPPGDLli2ffvppcXGxj4+PDNcJYKa6dHS4mc0OtKpRhsUh4pN5hhVJVvVDAVzIqG70wcHB7du3L126lIgUCsU111yzadOms4758ssvb7/99tbW1u+//766upr7hQKYp611Uoa/4GLJC8qca0aAUNrBtP2mvg4AoxhVi1Cj0TDGAgMDhz8NDAzMyck565iqqqri4uI9e/ZER0ffeeedr7766u23337eau3t7SUlJc8///zIV26++eaLDCjqdDqdTjea6wT5DP8vCNYzEMbNhhqaH0jGuUWH/wskSfa5DQJRtpo21uiWR8p9Ksuj0+lEEW1lE9PpdAaDQaFQXPJIpVJ5yV9cowpCxhgRjdQSRfHcR3FoaGhoaOjIkSOiKG7cuPGWW25Zvnz5ea9Sr9cPDQ21t7ePfGVgYOAiz7YkSUZ48uHihv8XEIRnYURb6sVHYo10hxrzf2FhkLClTrg5HI/e2Yb/s/FLybSk/49LtVEFoVqtFgShqakpPDyciBobG0dahyMCAwO9vb2H/1DKyMjo7OzUaDRBQUHnVvP19U1OTn7ppZdGeYlDQ0P29vajPBhkotPp7O3tEYRnKWxjzirDJB8j3Z+SJNnZ2Y3mr+Artzic/e2oXmVnL+L//ByiKKpU1tUbbmlEUTQYDLyiYVQNfAcHh4yMjI0bNxIRY2zz5s1z5swhIr1eX1NTM5zJ8+fPLysrGz6+tLTU0dHRz8+PyyUCmK3NdWxhsHUGRbCz4Oco5Lfi3VGwfqN9a/Svf/3rzTffrNVqy8rKOjo6brrpJiI6derUhAkTmpubfXx8brnllldfffW2226Lj49/5513nn76afzFBFZvS530YLwx2mcmsTBY2FzLJvtYZ9IDjBjtkO+CBQt27typUqmys7MPHDjg7OxMRGq1+rPPPnN1dSUiJyen3NzcrKwsSZJWrVr16KOPynjVAGagV09HmtlMtdXmxIJgcQsWHQUbMIaVZRITExMTE8/8iqur6y233HLmp3fffTe3SwMwbzsbWJqvtU2cONMMtVDcxjqGyMPO1JcCICe8BAxwmbbUSQuCrfkJslfQNH9hRwMahWDlrPkxBpCVFb8pM2JBkLilDu/LgJVDEAJcjsou1qenOGvZceJCFoYImxGEYO0QhACX48d6Ni/I+qdVTnQXRIHKO5GFYM0QhACXY1s9mxdk9TlIRDQ7QNhajyAEa4YgBBgzA6NdjZKV7ThxIXODhG0IQrBqNvEkA/B1pJkFOQsBTqa+DqOYFyTuapR0eHUUrBeCEGDMttpMvygR+ThQhKtwuBmNQrBaCEKAMdvWIM0LsqFnZ14QhgnBmtnQwwzARa+e8ltYlvWurHauuUHiNkyrB+uFIAQYm92NLNVHcB7D6oQWb4ZaKGhlXdgeG6wUghBgbLbVS3NtqV+UiBwUNMVP2N2IRiFYJ1v6sxaAhx0N7N1MowbhkGFoe/We0taK8paTTnaOE72ikvzjpwVNNuY1zA4Udzawa0KNeU4AI0EQAoxB6yBV97BUI27Rt7f2wFv5H0W4h00OSJoRMNUgSic7qj8oWLWm9Js/Tf5tlGeEcS5jdoBwdw5ahGCdEIQAY7C9XpqhFpVGaRBKTHrh4L/L2yofnfKnFHUCEfX399vZ2U0NmnxLzC9+OPnjwzue+k3i8muiFhjhYib7CrW9rKmf/B2NcDYAo0IQAozBzkY2O9AYzUGJSf848K+Owa53F75irzh7P0BREJeMX5gemPLnbX8d0A/eEL1E7utRCJTpL+5ulG6MtK3xUbAFuKcBxmB7gzGCUGLsqb0v9Or6np+x4twUHKF29vvX3L9/W7Hxq+Pr5b4kIpodKOxowGxCsEIIQoDRqu9l7YMszlP2IFxT9k37QOezWX9RKVQXP9Lf2fe1uf/3WelXRc2lcl/V7EBhRyOCEKwQghBgtLY1sNmBoihzDpa2VKwt++5vGQ8qRcVojvdz8nl86n3P5rzcOdgl64XFewldQ+x0D7IQrA2CEGC0dsjfL9o11P10zouPTb3P39lv9N81JTB1dljWPw78S74LIyKBaFaAuAuNQrA6CEKA0drZwGYHyBuEHxz7NCM47TLmCP426faWvrYdp/fKcVUjZgcK27HoKFgdBCHAqJzsYoxovLuMQXiyvWpP7f5fx99yGd+rFBUPpP3u7fz/DOgHuF/YiNmBAlqEYH0QhACjsrOBZcvcHHzjyAd3Jd7mZu96ed8e5xud6Be3uvRrvld1pig3QSKq7EIWglVBEAKMym4NmylnEG6v3tOvH7hq3LwrKfK75Du+rdjU2NPE66rONSsAjUKwNghCgFHZ1ciyZXtTxsAMHxR8em/qb0Thik7h4+R93cSrPy76nNeFnQtBCNYHQQhwaSc6mUgU6SpXEG47tVvt4pfgF3vlpa6PXnKg/kh9d+OVlzqvWQGYVg/WBkEIcGmyNgclxj4v/eb2uF9yqeascrp2/KI1Zeu4VDtXlJugFOkkhgnBiiAIAS5tV6OMA4S7anIcVY7J/vG8Ct4wacmu0/uaerW8Cp5lphq9o2BVEIQAl7ZbI9cro4zYp8Vf/ir+Jo413excF0fNW1P2LceaZ5qJYUKwLghCgEuo6GQiUYQ8A4RHGo8JgjglMJVv2Rujr912anfPUC/fssOyA4UdDdibEKwHghDgEnbLOUC4rmLDLyYu5l7Wy9FzSmDqpqrt3CsTUaSroBKFik40CsFKIAgBLmG3hs2Sp1+0qbe5uPn4nLAsOYovm3jVuooNEpMlrmaqhT0aBCFYCQQhwCXs1bAZalmC8LsTmxZEznZQOshRPNYn2kXlnKc5JkfxGQHCHgwTgrVAEAJcTFU3MzAa58Y/CHUG3aaq7UuiFnCvPOLaCYvWVWyQozJeHAVrgiAEuJjdjWymPM3BnTX7ojwiQtyC5Cg+bE7YjOLm43LMoxjvLhgYVXcjC8EaIAgBLmaPbEuMbqrads14GZuDROSgtJ8TPmNz1U45is8IwDAhWAkEIcDF7JZnKn1Tb3Nle/W0wDHvOzhWCyNmbzm1gxH/xJqhFnajdxSsAoIQ4ILqelmvnk2QYQ/CH0/tnB2WpVKouFc+y0TvKDuFXUlzOffKM9EiBGuBIAS4oF2NbKZalKNj9MdTuxZEZstQ+DzmR8zaLMOEwkkeQpeO1fciC8HiIQgBLkimftGSluOM2CTvCdwrn9f8iOzdNfsHDUN8ywpEmf4iGoVgBRCEABck0wzCLVU7F0bO4V72QnwcvSZ6R+2vO8S98swADBOCNUAQApxfUz81D7BYT85BqJcMu2r2zY+YxbfsxS2IzP7x1C7uZbPUQk4TghAsHoIQ4Pz2aKRMtch9hDBPcyzULdjPyYdz3YvKDJ5aoC3uHurhWzbRS2joY9p+vlUBjA1BCHB+ezUsS4Z+0Z01+7LDMriXvThHpUOKOuFA/WG+ZUWBpvkJ+7XYiQIsG4IQ4PzkGCDUS4b9dYdmhkznW3Y0ZoVm7Dy9j3vZLLW4F+/LgIVDEAKcR+cQVXWxZG/OQTjcL+rj5M237GhkBE8pbC7p1fXxLZulxurbYPEQhADnkdPEpvgJKt7Ph0n6RYc5Kh2S/eP31eXyLZvmK1R0sm4d36oARoUgBDiPvRopS8356TBhv+gwOXpH7URK8RH2491RsGQIQoDz2NPI/02Z/KaCULcgk/SLDpselF6gLebeOzpDLezV4H0ZsGAIQoCz9eupqJ2l+3IOwpza3KyQaXxrjomTyjHBL+ZQQz7fsnhfBiwdghDgbAebWYKX4KTkWZMR219/eHpQGs+iY5cRPGV/PeclZqb5C/mtbMDAtyqA8SAIAc4mR79oRWuls8pR1m14R2N6UPqB+iN6iWdqOStpkodwuBmNQrBUCEKAs+U0SZn+nB+NffW5GcFT+Na8DN6OnsGugUXNpXzLZvoL6B0Fy4UgBPgZvUSHtGy6P/8BwozgdL41L09GcDr3SRRZaiGnCe/LgKVCEAL8zNFWFuYqeNnzrNnUq20baJ/kPZFn0cuVETxlb+1BvjWz1OL+JmZAmxAsE4IQ4Gdymlgm9+ZgXe70oHRRkGOL3zGL9AgTBfFUx2mONX0cKMBJKGpDEoJFQhAC/Iwca23vrz883Tz6RYdND07bx3sBbgwTguVCEAL8DyPa1yRlcg3CAf1AWUtFin8Cx5pXaEpgam5DHt+a2JsQLBeCEOB/KjqZo0IIceYZhPlNhdHe451UjhxrXqEkv7iT7VVdQ90ca2aphT2NeF8GLBKCEOB/5OgXPVifNyUwlW/NK2SnsEvwjc3XFHKsGeEqKEWhsguNQrA8CEKA/9nXxPj2ixLR4caj5haEJE/vaKYaw4RgkRCEAP+zV8P5ldGarjq9pA93D+FYk4tpQZNzG/IY8cytDH9hH4YJwQIhCAF+oumn9kE2yYNnEB5syJsaNJljQV4CXPwdlY6V7dUca+J9GbBQCEKAn+zVSJlqUeTaM3qoIT89IIVnRX64947GewqaPqbt51gSwBgQhAA/ydGwDK79ogP6wdKW8lR1IseaHE0NSj3INQhFgab5C/uw1hpYGgQhwE+4rylToC0e7xlpVhMnzpToF3eyvYrvPr0Z/iKGCcHiIAgBiIh6dFTRyVJ9eAbhkcZjkwOSOBbky15hN8l7QoG2mGNNDBOCJUIQAhAR7deyVB/BXsGz5mGNWQchEU0OSDrceIxjwSm+Qkk769NzLAkgOwQhABFRjkbi2y/a1t/e0tc60SuKY51c7nIAACAASURBVE3uJquTjnANQnsFJXgJudikFywKghCA6Kc3ZXg+Dkc0x1LUCaJg1o9YlGdk12C3tq+FY02svg0Wx6yfUgDj0El0pIXzZrxHNAWT1WbdL0pEoiCkqBPyNAUca2aqhRwNXhwFS4IgBKC8FjbOTXC341kzX1No5gOEw1LViXx7RzP9xVwt0yMKwXIgCAFoH++JE6c6axSCGOii5lhTJmkBSXmaYxzXWvO0pxAXoRCb9ILlQBAC0L4mzlPp8xoL0gKSORaUj7+zn4udS2X7KY41M/0xiQIsCYIQbN3wZrwZXDedyNMUmO2CMudKVSfy3ZIpQ43Vt8GSIAjB1lV0Miclz814JSYVNZcm+cfxKii3ZP/4o01FHAtm+GOTXrAkCEKwdTm8t14qbzvp5+Tj6eDBsaasUvwTCrQlBmbgVTDSVVCKQlU3GoVgGRCEYOu4DxDmawpT1AkcC8rNzd5V7exX3lrJseZ0fyEHswnBQiAIwdbl8N6V/mhTUbK/JQUhEaWoE442cR0mxCa9YDkQhGDTmvqpZYDF8NuMVyfpS1vKE/xieBU0jmT/BL7DhJloEYLlGEMQrlmz5q677lqxYoVGo7nIYWvXrn3jjTeu+MIAjCFHI2X4Cxw34y1rKQ9xC3K1c+FW0SiS/ONKWo7rDDpeBRO9hYY+1jrIqx6AjEYbhP/+979XrFgxc+bM1tbWzMzMwcHz3+C5ubn33nvvc889x+8KAWS0r4lN57rEaH5TYbJ/PMeCxuGscgp1Cy5rreBVUCFQuq+wH5v0giUY1a8Ag8HwyiuvvPXWW7fddts777zj6ur69ddfn3vY4ODgH/7wh2eeeYb3RQLIhftmvEc1RSmWNkA4LNk/Pp/rMOF0bNILFmJUQVhXV1dbW5udnT386axZs/bt23fuYc8999zSpUtjYixsdARsVq+eyjpYmi+3IBwyDJW3nYy3tAHCYSnqhKNNvDfpxTAhWALlaA5qbGx0cXGxt7cf/tTX1/fIkSNnHVNQUPD9998fPnz44MGDF69WV1e3d+/e6667buQrjz76aHz8BXuT+vv7FQqu+6XC2PX19QmCIAg8G08mt7tJTPAUpcG+Pk4FC1tKw91C2ZDUN8Sr5M/09/fr9XqZHodxzuHlrSc6ujvsFHxWH09woYI2VVt3n4N1Pb6Dg4OiKKpUKlNfiE3T6XQGg0GSLt337uDgIIqXaPKNKggdHByGhoZGPh0aGnJ0dDzzAL1ef9ddd73zzjsjYXkRXl5eYWFhv/zlL0e+EhUV5eDgcKHjdTrdRf4VjEOv1zs4OFhZEB5qY1lqxvHuKm2vSFHHy3e7Msbs7OxkCkIHcojwCD3VW5voF8upIMV4SEXd9llcZ6eYnCAICEKTUygUBoNhNM/aJVOQRhmEQUFBg4ODzc3Nvr6+RFRbWxsUFHTmAVVVVQUFBbfddhsRDQwMtLW1jRs3btu2bREREedWc3JyCg0NvfHGG0dzaiISRXE0PwnIavh/wcqCcJ9Wf3+sQuT3zmhBc8ny2Ovku13F/0+m+ol+cYXNJclqbi/7ZKrZ/mZhZqBVPb9y/y/AaIiiyBjj9b8wqiq+vr5ZWVmfffYZEXV0dGzcuHHZsmVE1NbWtm7dOiKKiIg4fvz41q1bt27d+sorr7i7u2/dujU4OJjLJQLIQS/RIS3PzXh1Bl1564k4n0m8Chpfkn/cMS3PYcIMf2zSCxZgVC1CIvrnP/+5dOnSnTt3lpaWzp8/f9q0aUR0/PjxX/ziF4wxlUoVGRk5fGRtba1CoRj5FMA8HWtjoS6C16X78kertLUiwiPUSeV46UPNVYJv7DM5L+kMOpWCT79fllq8c4/BwEhhVV0JYG1GG4TTp08/fvz4wYMHAwICkpN/2mgtJSWlvLz8rCOnTJly+PBhntcIIIMcDeeV1Y41FSX6WcyOE+flpHIMcQs63nYi3pfPi6++DuTvKJS0swQvJCGYrzF0sHp5eV111VUjKUhEDg4OEyZMOOswBweHsLAwPlcHIBvua20f0xZb0NZLF5LsF3+M6yQKrLUG5g/jvWCj9jVJHKfS6yT98VZuDSkTSuQ9TJipxm71YO4QhGCLTnYxURDCXbkF4fHWihC3IGeVE6+CppLoF1vWUqGXuO1NmOkv7EWLEMwbghBs0V4Nm8F1gLBAW5pk4QOEw5xVTkGuAeVtJ3gVHO8uGBg73YMsBPOFIARbtI/3EqMF2mJLf1NmRKJfbKG2lGPBDH8Rw4RgzhCEYIv2cn1lVGJSaUt5vK8FzyA8U4JfbKG2hGPBTH8ME4JZQxCCzWkeoKZ+FuvJLQhPtFf5Ofm42bvyKmhaCX4xhc2lEuM2ET5TjWFCMGsIQrA5ezVShr/AcYp3gbYkgdP6nObAw97d29GrquM0r4JJ3kJtDzbpBfOFIASbk6NhmWqed36httSagpCIEv1iC/j1jioEmuKHTXrBfCEIwebkNLEsfm/KMGLFzaW8dmwwEwl+MZyHCdV4XwbMF4IQbMvwZryT+W3Ge7qzzlHp6OPoxaugOUjyiyvQljDiFl1ZGCYEM4YgBNtyUMuSvAWOW8UWaksSLX9ltbP4OvnYK+zquhp4FZziKxS1s349r3oAPCEIwbbs1Uh8p9IXaksTLX9ltXMl+sdxHCZ0UlK8p3CoGY1CMEcIQrAtezUsi+ubMgXa4ng/KwzCBN9JHIOQ0DsKZgxBCDZEJ9HhZjbNj1uLUNOr1Un6YNdAXgXNR7xvTFEzz/VlstTiXmzSC2YJQQg2JL+FRbkJ7nbcChZqre190RGh7sED+gFtXwuvghn+wkEt0yMKwfwgCMGG7NGwLM4DhCUJ1tgvSkQCCXG+k4qay3gV9LSnMFfhWBt6R8HsIAjBhuzlHYRFzaVWsAfhhcT7xhRxXX0bw4RgnhCEYCsY0f4mKcOf2z3fOdil7WsZ5xnOq6C5SfDjPUyIvQnBLCEIwVaUtDNvByGA39a5Rc1lcT6TFAK/OYlmZoLXuMaepp6hXl4FZwaIORoJSQjmBkEItoJ/v6i2NN7PSrZeOi+FoJjoHVXcwm2YMMCJXFVCeQeiEMwLghBsxV4N5814C5utatOJ80rwjeU7TDgjQNiD3lEwMwhCsBU5GjYzgFsQDugHqzpqor3G8yponhL8Ygq4BmEmhgnB/CAIwSZUdjFGFOHKLQjLWivGeYQ7KO15FTRPMT4TT7ZX6Qw6XgVnqIVdjQhCMC8IQrAJe7g2B8mqZxCeyVHpEOoefLztBK+C490FiVF1N7IQzAiCEGwC96n0Rc1lthCERJTgG1OIYUKwaghCsAl7GhnHTSckJpW1VsT5WPMroyPifWOK+a0vQ5hNCOYHQQjWr76XdetYtAe3IDzZfsrX0dvN3pVXQXOW4BdT1FwmMW7RhRYhmBsEIVi/3Ro2I0Dk2DFa1FxqlVsvnZeng4e7vdvpzhpeBWM9hbZB1tCHLARzgSAE67dXw7NflIgKtaXxvjbRLzos3i+mkN9aawJRpr+Yg0YhmA0EIVi/3VwHCImouLnMitfaPle876QiLddhQjV6R8GMIAjByjUPkKafJXhxC8KGHg0JQoCLP6+C5i/Bl2eLkIhmBgi7MZsQzAaCEKzc7kYp01/kOEJYqC1NtKXmIBGFuAXpJB3HTXqTvIX6PtYywKsewBVBEIKV293IeSp9UXNpnI0FIRHF+URzXHRUIdA0P2GPBtvVg1lAEIKV28U7CAu1pQlWvenEecX7TuLdOyqidxTMBIIQrFnrINX1smRvbkHYOdjV2t8W6RHOq6CliPfjvFv9LAwTgtlAEII1290oZfgLCq4DhLG+0aJgcw/OBK9xml5t11A3r4Ip3kJ1D2sd5FUP4PLZ3PMMNmV3I5sZwPMmL24uS7C9AUL6/5v0lraU8yqoFGman7AXw4RgBhCEYM34DxDa0poyZ0nwjS3iuugohgnBTCAIwWq1DVJ1N0vhN0A4aBiq6jht9ZvxXki87yS+21BgmBDMBIIQrNYejZShFpT87vHSlnJb2Iz3QmJ9oyvaKocMQ7wKTvYRqrpZG4YJwdQQhGC1djWymWqed3ihttRG9iA8L0elQ5hbcHlbJa+CSpGmYpgQzACCEKzWzgaWHch5idE4W1pr+1zxfpOKuM4mnBUg7kLvKJgaghCsU+sgne7hOUAoMamk5bhNbTpxrnhfzrMJswOEHQ0IQjAxBCFYp50NUqY/zwHCyo5qXydvd3s3bhUtUIIv5016U32Emh4sOgomhiAE67SzkWUH8h0gLLGprZfOy8vR083eleMmvUqRMvyF3Y0YJgRTQhCCddrZwLL5LzFq60FIRAl+sXwXHc0OFHdimBBMCkEIVkjbT5p+lshvgJCG15Txi+VY0EIl+MYUaEs4FswOEHZimBBMCkEIVmhHgzRTLXJcYrS+u5EEQe3sx62ixUrw4xyESd5CUz/T9HMsCTA2CEKwQjsbOU+cKGwuTURzkIiIgl0DGZOaerW8CooCZarFXQ0YJgSTQRCCFdrVyGZxHiAssc21ts8r3jemgOskitmBAoYJwYQQhGBt6ntZ2yCL8+S7Kz0GCP8n3jeG77T62YHCdgwTgukgCMHabK1ncwJFkV8Odgx2dgx0hruHcqto4RL8Ygq5DhPGegq9OlbdjSwE00AQgrXZ3sDm8B0g1JbG+UaLAs+aFi3KM6K5r7VzsItXQYEoO1BEoxBMBUEI1mZnI5sTxHuAEP2iZxAFMdY3mu/ehHPQOwqmgyAEq1LawZQCRbryDMICbQleGT1Lom8s397RuUHC9gYJSQgmgSAEq7K9ns3j2hzs1fXVdTdM8IriWNMKJPrH8t2kN8xFcFUJJe2IQjABBCFYFe4DhMXNZdHe41WikmNNKxDtNb66s6ZPx3Ma/JxAYXs9ghBMAEEI1sPAaK9GmsV7re0EX/SLnk2lUE30iiptKedYE8OEYCoIQrAeR5pZkLOgduRZs0CLNWXOL8EvtrCZ5zDh7EBxj0bSYYUZMDoEIViPbQ1sLtd+0SHDUGXHqRifiRxrWo0Ev5iCJp5B6ONAka5CrhaNQjA2BCFYjx/rpPnBPG/p0pbyCPcwB6U9x5pWI853UnnbySHDEMea84KEbVh0FIwOQQhWoltH+a0sS8134gT6RS/IUekQ5h5S3lbJsea8IPHHOrQIwdgQhGAldjVKU3wFZ65vdxY2l8RjM94LS/SL5bslU6ZaKG5n7YMcSwJcGoIQrMTWejYviOf9rJcMZS0V2HTiIhL8YguaijkWdFDQdH9hVyN6R8GoEIRgJX6s4zyVvrztRJBrgIudM8eaVibRL7ak5bheMnCsOS9I3IrZhGBcCEKwBnW9rHWQJXnzDMJjTcVJfnEcC1ofVzuXABf/E+18hwmFHxGEYFwIQrAGm+vYvCCeWy8R0TFtcaI/gvASEv3ijnHtHY33Evr0rApbMoERIQjBGmzjvcSogRlKmo9jgPCSkvw5B6FANDcQ746CUSEIweIZGG2rlxYE8wzCirbKABd/N3tXjjWtUpJ/XHFLmcR4vt6yIBi9o2BUCEKweIeaWbCzEOjEeYAwEQOEo+Bm5+rr5HOivYpjzQXB4s4GaQivjoKxIAjB4m2ulRaFcN4+vkBbnIQBwtFJ4j1M6ONAE9yF/U1oFIKRIAjB4m2qYwu5rqwmMam4+Xg8BghHh/swIREtDBY216FJCEaCIATL1jJAJzrZdH+eLcIT7VU+jl6eDu4ca1qxRL+4ouZSvsOEC0PEzbVoEYKRIAjBsm2pk7IDRRXXGzlfU5isjudZ0ap5Orh7O3ryHSZM9xXq+1h9L7IQjAFBCJZtcx1bxPV9USI62lSU7J/At6Z1S1EnHG0q4lhQIdDcIHEL3h0Fo0AQggWTGG3lPXHCwAzFzWXYdGJMkv0Tjmp4BiENDxOidxSMAkEIFuxIC/NxEEJduC4x2noywEXtbu/GsabVS/KPK2ou5bvo6MJgcVsDNqwHY0AQggX7oUa6OpRzv2i+pjDFHwOEY+Nm5xrgoq5oO8mxpr8jRbkJOZhEAfJDEIIF+6GGLQ7hfA/nN+FNmcuR4h+frynkW3NxiLihBk1CkB2CECxVQx+r6WHT/Hi2CHWS/njriQRfDBCOWbI6Pr+JcxBeHSr8UIMWIchuDEG4b9+++fPnJyUlPfbYY4ODZ+8hrdVq//KXv8yYMSMtLe3ee+9tamriep0AZ/uhhi0IFpVc/5YrbSkPcQvCHoSXIdEvrqy1QmfQcayZ4iP06OlEJ7IQ5DXa3yItLS2LFy++6aabVq9efeDAgaeffvqsAyoqKgYGBp555pn333+/vr5+2bJlnK8U4Oc21LLFvAcIjzUVJWOA8LI4q5zC3EJKWys41hSIrgoRNuDdUZDZaINw1apVaWlpd955Z0xMzAsvvPD+++/rdD/70y8zM/O1117Lzs5OSkp65ZVXDhw40N3dLcMFAxARDRpod6O0gOvKakSUpylIUWMG4WVKUSfkaQr41lwcImyoxTAhyGu0v0eKiorS0tKGP548eXJbW1t9ff2FDs7Pzw8MDHR1xRY2IJcdDSzRS/C251lzQD9wor0KexBetlR1Yj7vIJwbJOZqWecQ36oAP6Mc5XFarTY6Onr4Y5VK5eLi0tTUFB4efu6RDQ0N99133+uvv36hUidPnvzuu+8iIiJGvvLuu+9Onz79Qsf39vYKAuceMBir3t5expj5/Ed8W6Wc6089PQMcax5uOhrlHqEf0PdQD8eyHPX399vZ2SkUClNfyPlFOIacbD+l7Wh2UjpyLDvVR7W+sm9piLm0CwcHB0VRVKlUpr4Qm6bT6QwGg16vv+SRTk5OoniJJt9og9Dd3b23t3f4Y0mS+vr6PDw8zj1Mq9XOnTv3vvvuu+GGGy5Uaty4cXPmzHnttddGvhIWFnaRZ5sx5uLiMsrrBPk4OzubSRAyok0N+q1XKVxcHDiWLS2vSA9KMeebTaFQmHMQElGMz8STvaemB6VzrPmLSGlLk/LWSebyU6tUKgShyQ0HoYMDn98Ao+0ajYiIqKj4aRi8srJSoVAEBwefdUxra+u8efNuvPHGxx9//CKlBEFwcXGJPIM5P9hghvJamIuKJrpzTuUjmmOp6iS+NW1NqjrxSCPn3tElYcKGWiwxAzIabRAuX75848aNlZWVRPT6668vXbrU2dmZiFatWrVt2zYiamtrmzdvXkZGxv3339/e3t7e3m4w8FxvCWDEd6ela8M4p2D7QEdzX2u0dxTfsrZmsjopT3OMb81AJ2GCu7BHg3dHQS6jDcJJkyY9+eSTqampwcHBBw8efPnll4e/vn79+gMHDhDRnj17qqur16xZM+7/q66ulumiwcZ9W82uDeP8vugRzbEk/3hRwBITV2S817j2gc7mvha+Za8NE787jSYhyGW0Y4RE9NBDD/3xj3/s7u729fUd+eLatWuHP1i6dOnSpUs5Xx3AOSq7WMsAS/fl3CLM0xSmYuLEFRMFIdk/Pl9TuCByNseyS8OEhZul16eRWYxRg9UZ29+/Dg4OZ6YggPF9e5pdGyaKvH8j5msKMUDIRYo6IY/3WmuTPARHBR1rRe8oyAIdQWBhvjstce8XremqJ2KhbkF8y9qmtIDkI41HGXEOrSVhAnpHQSYIQrAkTf1U3M7mBHFuDx5qyE8PSOFb02YFuqgdlY5VHaf5ll0WLn59Ci1CkAWCECzJN9XSVSGiHe/b9lBjfnoggpCb9MDkQw35fGtO9RM6h6isA1kI/CEIwZJ8dUq6Lpxzc3DIMFTcXJaqTuRb1palBaQcbjzKt6ZAtDRcQKMQ5IAgBIvRMkD5LWwh74W2C7Ql4zwjnFVOfMvashT/+LLWin49zwXwiOj6CPHragwTAn8IQrAY31RLC4JFxzFM+RmVQ41H0wOSORe1bQ5Kh0neE442FfEtm+kvaPvpZBcahcAZghAsxtenpOsj+E8kO4w3ZWSQFpB8uJHzMKEo0NJw4Sv0jgJvCEKwDO2DlNvMFoVwvmOb+1raBzrHe43jWxbSA1MONXAeJiSi68LFr0+hdxQ4QxCCZVh3WpobJDrz7hfNbcifHJAkmseuGtYk0iNsQD/Q0KPhW3ZmgFDTy051o1EIPCEIwTJ8XindHMk/rg42HJkalMq9LAgkTAlMPVB/mG9ZhUDXR4hrqhCEwBOCECyAtp/yWthVvPtFdQbd0aaiKQEIQllMC5p8oP4I97I3R4qfnkDvKPCEIAQLsKZKuiaU//uiR5uKItzD3OxdOdcFIiJKC0gubSnv0/XzLZuhFvoMVNyORiFwgyAEC/B5pXTzOP736oGGw9OCJnMvC8MclA4xPhO5b08oEN0YIXxeiUYhcIMgBHN3uodVdbO5gTIMENbnTQ9K414WRsjVOzpOXFPJ0CQEXhCEYO5WV7LrI0Ql71v1VGeNgRkiPMI414UzTA9K319/SOKdWUnegoOCDmkRhcAHghDM3eeV0s2R/G/Ug/VHpgelcy8LZwpw8XezdzvRVsm98k3jxM/QOwqcIAjBrOW3sG4dZajl6BfFxAljmB6Utp/3JAoiui1K+LxSGkIUAg8IQjBr/z0h/Wo89+3oqXOw62THqRT/BN6F4WzTg9Jy6nK5lw13FWI8hE21SELgAEEI5ksv0doq6dYo/s3BfXWH0gKS7RR23CvDWeJ8Y9oG2rkvMUNEd0wQ/3sCw4TAAYIQzNeGWmmCuzDOjX8Q7q09mBU8lXtZOJcoCDI1Cm+MEHc1Ss2c93oCW4QgBPP1yQl2xwT+t+iAfrBAWzwlEAOERpIVMjWn9iD3si4quipEXINXZuCKIQjBTLUO0o4G6foI/rdobkNerE+0i50z98pwXqn+iZUd1W397dwr/2q8+F8stwZXDEEIZmrVCemaUNFNxb/y3tqDWSHoFzUelUKVFpB8oIH/zPrZgULLAB1rxUghXBEEIZipD8ql30bzvz/1kiG3MS8jeAr3ynARWSHT9tYe4F5WFOjOieL75WgUwhVBEII5ytEwA6NMGaYPHmsqCnUL8nb05F4ZLmJqYGpRc1mvro975d9MENZUSn167oXBhiAIwRy9Xy7dHc1/+iAR7Ti9d1ZopgyF4WKcVU6JfrH76g5xrxzkLGT4i2ur0CiEy4cgBLPTOUTfn5ZujZKlXzSnLndGyDTuleGSZodl7TydI0fl30YL6B2FK4EgBLOz6qS0KET0deBf+YjmaJh7iL+zL//ScCkZwVMKtMVdQ93cK18VItb0YIdCuHwIQjAvjGhlqfS7SbLcmTtO52SjX9REHJUOqepEOXpHFQLdNVFcWYpGIVwmBCGYl231TCHQDBlek9EZdPvrDs0IRb+oycwOy9pxeq8cle+ZJH5eKbUPylEbrB+CEMzLGyXSA3Gy3Ja5jflRnhE+jl5yFIfRmBaUVtpS3jHYyb2y2pGuChE/xuR6uCwIQjAj1d3sgFa6eZxM/aJ7s8PQL2pKDkr79ICUvTIst0ZEf4oV3yqVJAwUwtghCMGMvFEq3TlBdFLyr9yr68ttyJsVmsG/NIzFvIiZP57aJUflqX6Clz1tqkMSwpghCMFc9OrpkxPS7+V5TWZ3zf4U/wR3ezc5isPoTQlMretukGNXJiK6N0b8d7FBjspg3RCEYC7ePy7NCRTDXeWYRk9bqnbMj8iWozKMiUJQzA7L/PHUTjmK3zxOPN5JR7H0KIwRghDMgl6ifxVLf5bnNZmmXu2pzpqp2HfJPMyPyN5ctYMR/7hSiXRvjPhKEV6ZgbFBEIJZ+KJKinKjKX6yNAc3V+2cE56lUsiwkwWM3USvKEelQ3HzcTmK/26SuKVOOt2DRiGMAYIQzMIrRdIjCQqZim+r3rUgYrZMxeEyzIuYJVPvqKuKfjNR/FcxGoUwBghCML0f65mB0fxgWZqDhdoSURCjvcfLURwuz/zwWbtO7xvQD8hR/L5Y8ZMTUism18OoIQjB9J7NNzyWKMteE0T0/Ykt10QtlKc2XCYfJ+94v0k75FmDO9BJuD5CfK0Ir4/CaCEIwcS21jPtAP0yUpZbsWuw+2DDkfmRs+QoDldiyfiF35/YLFPxFUni22VSiywNTrBCCEIwseeOGp5OERXytAc3V22fHpzuZucqS3W4AukBqe0DHRVtlXIUD3URrosQXy9BoxBGBUEIprS1nmn65WoOEtGGyq1L0C9qlkRBWBw174eTP8pU/69J4jtlUhtGCmEUEIRgSs/ky9gcPNZUxIjifKNlqQ5X7Opx83ec3tun65ejeKiLsCxcfAUjhTAKCEIwme9OS906ukm25uC6io3LJlwlU3G4cl6OninqhC2ndshU/6lk8d0yqa4XcwrhEhCEYBoGRiuOSC+mK2R6W1TTqz3aVLQwco4s1YGTG6Kv/er4eonJklVBzsKdE8XnjmJOIVwCghBM46NyydeBFsgzd5CIvjr+/eJx8xyVDjLVBy7ifSe52bscqD8sU/0nkhTfnpZKO9AohItBEIIJ9Onp2aPSS1PkWkqmT9e/5dTOZRMXy1QfOLo++tq1x7+TqbiHHT2SoFhxGI1CuBgEIZjAPwsMM9TCZB+5moM/nNySHpDi5+QjU33gaFbo9MYezfHWEzLVvzdGLGhjOxrQKIQLQhCCsVV1s7fLpBfT5br39JLh6/IfboheIlN94EshKJZNWPzV8fUy1XdQ0KtTxT/tN+jQLIQLQBCCsT1wQHokQRHkLFdz8MdTO4NcA7C4qAVZMn7h4cajtV31MtVfGiaGu9KbpUhCOD8EIRjVxlpW0ckekGffQSKSmLS69Os74m+SqT7IwVnltGziVZ+Vfi3fKV6bqvjHMYNGlimLYPEQhGA8vXq6d7/hjekKO9nuu62ndnk7eiX6xcp1ApDHDdHX7q87VNfdIFP9Ce7CXRPF+w9gfj2cB4IQjOeJw4ZZAcK8ILk6RSUmfVry5R1xv5SpPsjHWeW0dMKiz0u/ke8UqlQfbQAAFptJREFUT6UoitvYN9XoIIWzIQjBSA5q2dfV7GXZpkwQ0bbq3e727inqBPlOAfK5PnrJntoDDT0amerbK+iDGYr7DkjtWIAUfg5BCMYwYKA79xhenyp62ct1iiHD0AcFn92TfLtcJwCZudm53hh97XvHPpHvFNP8hGVhwoO56CCFn0EQgjE8kmtI9Baui5DxfvuqfH20d1S8b4x8pwC53ThpaWlLeXHzcflO8Y80xb4m9tUpdJDC/yAIQXabatn3NWzldBk7RbuHetaWfXd3EpqDls1eYffr+JtX5n/ESK757y4q+jxb8cf9htM9mGIPP0EQgrw0/fSbvfrPZik8ZesUJaL/FH6eHZYR7Boo4znAKBZEzh7QD+TUHpTvFKk+wp/jFHfsNhgQhUBECEKQlV6im3bofzdJkamW601RIjrRXrXj9N5fxd8s3ynAaERBvDf1rjfyPhzQD8h3lkcTRKVAT+VhsBCIEIQgq4dyDfYirUiS8TaTGHv10Nv3JN3ubu8m31nAmFLUCYl+sf8p/Fy+U4gCrZ2j/LwSg4VAhCAE+ayulDbVsS/mKGXagH7YtxUblaJy4TjsO2hV/pR614/Vu060V8l3Ci97+nqu4o/7DSXt6CG1dQhCkMW+Jvbng4Zv5io87GQ8i7av5b9Fax5O/4NAcoYtGJ2bvetvEpa/krvSwGTsvUzyFl5IUyzdamiWsRcWLACCEPg72cVu3G74ZKYyzlPGfJIY+8eBf90QvSTMPUS+s4CpLI6a52LnvKp4raxn+dUE8ZZxwlWb9b16Wc8DZg1BCJw19tH8TYbn00T5dp8f9kXZOp1Bd0vsdbKeBUxFIOHxafd/W7GppEXGaYVE9HSqIs5LuHmHQY/hQluFIASetP00Z6P+nmjxjvHy3lon26vWlK77a8ZDooB72Gr5OHo9lP77/9v3ap9Oxm0jBKL3MhWM2K27MKHCRuGXCHDTMkBzN+lvGic+lijvfdU11P3k3hfum/xbtbOfrCcCk8sKmTY5IOmfB1+Xb4o9EalE+nquslfP0C60TQhC4KO+l2Vv0F8bJjyZLO9NJTHp//a9OiNk2pzwGbKeCMzE/ZPv6RjolHuw0E6kL+coO4bYrbsMQ8hCG4MgBA4qu9isDYbrIoTnUmVcR23Y20c/lpiE1dRsh1JUPJ356PoTW/bVHZL1RA4KWj9fyYgWbtJ36WQ9FZgXBCFcqd2NLHO9/m/J4tMpsqfg2rJvD9YfeSrzEQwN2hQvR8/nZvzlpdw3iprLZD2RvYJWZysmegjZG/S1vRgwtBX4bQJX5K1S6aYd+k+zlbfL/HYMEW2u2vFV+fqXZz/tauci97nA3ER7j38y45G/7Xm+vO2krCdSCPR2huLmceLU7ww5GmShTUAQwmXq1tHtuwzvHZf2LVHOCZR9PvuPp3a9X7Dq1TnP+eMFGVuVok54KP0Pf9n1XGVHtdznejhe/M9MxfXb9a8USQhDq4cghMuRq2XJ6/SOSjqwRBnpKnsKfnV8/QcFn746+1nsL2HjskKm3T/5noe2P1moLZH7XPODhNxrleuqpYWb9I19cp8NTAlBCGPTq6eHcw1Lt+pfShffzVQ4KeU9ncTYu0f/+/3JzW/O/ydWkAEimhk6/cmMh5/c+889tQfkPleYi7BrsXK6v5i0TvdhOZqGVgtBCGOwrlqK/1qv7aei61TLwmW/eboGux/b9UxZa8Wb8/7p5+Qj9+nAUqSoE16a/cybeR++e/S/EpN3roNSpKdSxK2LlB+US9kb9AVtSEMrhCCEUclrE2dtMDydL32QpfhklsLHQfYzFjWX3r35wUj3sFfmPOtm7yr7+cCijPeMfH/RqxXtlQ/teLKpt1nu0yV4CfuuUd4UKS7cpL9nv1CPnlLrgiCES9jfxK7aor9tn+q2KCF/mXK2/O/FDOgH/n3k/adzXrpv8t2/T/m1QpB9VgZYInd7t5eyn56sTrp704Pfn9gs69IzRCQK9LtJYvmNKn9Hmvw9/X6fobobrUMrgSCE8xuS6LOT0tTv9bfvNiwLE48tHvzNRFHWnQWJSGJsc9WO29b/oVfX+/HiN6YHpcl7PrBwoiAuj73+9XnPb6ra/octjxQ1l8p9RjcVPZvMipaStz1N/lZ//XbDrkbEocWT+VUHsECHm9mqk9KaSinRW3giUbw6VBQF6umR96QSk3bX7P+05EtHpcMzWY/F+EyU93xgRcLdQ1YueHFb9e7n9r0a5Rm+PPb6WJ9oWc/o40D/N1nxeKLikxPSH/cZDIxuGy/eGiWEuWBfTIs0tiDU6XQqlerKjwFzMyTRvib2Q430TTVTiXRrlJh7rTJC/nkRRNQ+0Ln11M51FRu9Hb1+k3grWoFwGQQS5oXPmhkyfWPltv/b96qPk/e14xdmhUyzV8i4MbSLiv4QI/4hRjyoZZ+ckCZ/awh3Ea6LEBeFCAleAiLRggiMjapd39PTc8cdd2zdulWhUDz22GOPP/74uce8+OKL//jHPwwGw5w5cz755BNX1/O/4LB69eoNGzZ89tlno7zE7u7uC5WCKzFooPxWtlfD9mqkvRoW7SFcFSIuCxPivc7zCPf09Dg7Owv8nu7Owa799Yf31h4o0JZkBk9ZMn6h3H/FW4H+/n47OzuFAoOmFyMxaU/tgQ0ntx5vPZEVMjUrZFqqOsGOXyIODg6Konjun/sGRrsb2benpS11rFvH5gSKWWohUy1EuwsiUpE3nU5nMBgcHPi8tjfaIHziiSfy8vLWr1/f0NAwZcqUdevWTZ8+/cwDcnNzFy9enJubGxISsmzZsri4uBdeeOG8pRCEplLbyyo6qbiNFbez/FZ2vINN8hAy1cIMtTArQPSyv9j3cgnC9oHO460nCrTFx5qKa7rq0gKSM4KnZAZPcVI5XklZ24EgHJPmvpZdNftz6nJPtFXG+ExM9o+P95003muco/KKfnteKAjPdKqb7WpkezRsXxNr6mPJPkKStxDnKcR5ChPchYs/azAapgnCoKCgjz/+eN68eUT0wAMP9PX1vffee2ce8Pvf/16hULz55ptEtHPnzptuuqmpqem8pRCE8tFL1DpILQOsqZ8a+lhDH9X3slPdVN3DqrqYq4qiPYRYTyHeS0jyEhK9BYdR/0YdaxAOGYaa+1o1vdrGHk1dd+OpjppTnaf7dP3R3uPjfKNT/BMm+UxUiRiiHhsE4eXpGeot0JYcbSosbamo7KhWO/uGu4eGuYcEuwYEuKgDnP08HT1G/3LyaILwTB1DdKSZFbSx4nZW0s4qOplSoEg3IcJVCHGmEGch0JkCHAVfR/J1QEaOFt8gHNVvor6+voaGhpiYmOFPY2Ji1qxZc9YxJ0+eXLZs2cgBWq22q6vLzc3tvAWHhoba29tHPvXw8LjIb9juvt72fsubtqNn1Ku/9B8ZvXo6ayPQbh0N/3EiMerWMyIySNSjJyLq0hEx6tQxvUS9OurVU5+BeodYl466dNQxxHr15GVHng6Crz35Ogr+jhTpKGQFU7CzEOIsuPz8yR3soUEiIurT9RsuMCu5T99vYAYiGugfEFSCnhmIqGeojxHr1w/omaFb1zugH+jXD3brent0fV1D3Z2D3e2DnQOGIR8HDz9Hn0Bnv2DXgCVhs8JcgwPPXCZ0YACbvo0VGxiQ9EMCgnCMnIimecVM84ohIgMznO5uON1VV91df7DmkKZP29TX0jnU7Wbn4m7n6mbn6m7v6qJyclU5OyodHJT2TkpHe4XKTrRTiSoHpR0RiZKgUqiUSiUR2Sns7MXzd7raK+zsFCoiEojSXCjNhSj0p39qHaDaHlbXyxr6qbaFHa6h5gHWNkjNA6xPTx4qcrMXXJXkakeuKsFJSfYiuduTgshFJSgFclIREbnbERGJRK6qn355qhTk+PNbw0FB9qO4WVxVlrSZi0qlDPD05ltzVEE4HFouLj8t+e/m5tb2/9q7+5imrjYA4M9t62gRBNF2FLKOGj/AL0SYmwaHImKIEo3WuZG6CEi2bGZLwG0m24xvWMIfc1uMWzZdcNEsM1V0C3XLlqFTy9wUNoSkGbVYIp9llH4oFGjvvef94243fSlg99oP1vv8/jqePtz7KL0+7bnnnmO3+8fwAdwXOIfDMWkhNJvN9fX1DQ0NfM/Zs2dzc3OnOvt/6g91xTimelUIKACY9hmpOIA4gBQxAPemHwdmHKz3wQrQGsDxZSwlmuLwMhbE5K/TxxBKwhKuU0RASigxC7MZiCEgZWAuS81mqHgG4mhIYKh4GgAcAA4AM3+0/kD+tgiFmAwgHcD3jjRLUS7JyAPx8H1x34iEGhGTERG4xGRARI1RxCMCL0VoEYxTAABjImD//tzuocAzRQ3xiIg34DsJYgA5gJy/hL0AXqBH/7qEJkVAoDcexSwcKzieNDueYRiaph8aHxsbKxI9pNAHVAjnz59PUZTL5UpISAAAh8Px+OOPT4iRy+Uul4trO51OrmfSoy1atEij0QQ+NHr0+eM4NBpxQZ8sg/4PODQaOoGvY/tPh0ZRKAR3aDSgL8QxMTELFixoaWnh/nj79u309IkT/DIyMnwD0tLSYmNjg5IiQgghFDqBjgxXVFS899573d3d169fP3fu3P79+wFgaGhoy5Ytg4ODAFBeXn7hwoWffvqpp6enurq6oqIihFkjhBBCQRLotL2qqiqbzbZ+/frExMRPP/10xYoVAEAIGR8f5+adLlu27OTJk5WVlQ6HY9euXW+88UYIs0YIIYSCJNDHJ4Lonz4+8cEHH5SVlc2dOzekWaHpnThxoqioSKVSPTwUhYxOp1u6dCn3MRRFyvfffz979uz169dHOhFB++WXX2w2W3FxcVCO9i+YNHvmzJl79+5FOguhu3DhgtEY8j3B0fS+/fbbpqamSGchdJcvXzYYDJHOQuh+/vnnH3/8MVhH+xcUQoQQQih0sBAihBASNCyECCGEBC0Ck2VqampqamqmetzeX29vr0KhwMdXI8tqtSYkJMhkuDp2JA0ODspkMn4JJxQRdrtdLBZzq4ugSHG5XAzDJCUlPTSypKSkurp6+pgIFEKWZc1mc+CFbXx8PCYGV6KNMPwtzARer1csFj90vSgUUjRNUxSF6/tEFsMwhBBuxdfpKZXKh36Cj0AhRAghhGYO/GiJEEJI0LAQIoQQEjQshAghhAQNCyFCCCFBC3TR7fAwmUxffvklTdMlJSWTLqg4PDz8+eefd3d35+bm7ty5M/wZRj273a7X641GY3x8/M6dO5ctW+YfU1tbyzAM116yZEleXl54c4x+FovFd+fqbdu2paSkTIjxer21tbV37tzJzMzcu3cvTiUNura2tl9//dW3R6vVTthd7ty5c9z2qwCgVCqDtfQl6uzsbG5udjgce/bs8X1S5bffftPpdDKZrLS0NC0tzf8HBwcHa2trBwcHt27dmp+fH+DpZtDF09HR8fTTTxNC4uLicnNzW1snbq5OCCksLLx69erChQsPHTr0/vvvRyTP6Hbw4MH6+nqFQuFyudasWTPpan6vvvpqW1ubxWKxWCzcJlwouJqbm2tqaix/Gx0d9Y/RarVnz55dtGjR8ePHX3/99fAnGfWcTif/K/jmm2/effdd/4e+jhw5YjAYuJj+/v6I5Bl9BgYGsrOzT5w48dJLL/355598/40bN/Lz8+fPn+92u3Nycnp7eyf8oNvtXrt2rclkevLJJ0tKSnQ6XaCnJDPGa6+9tn//fq795ptvvvjiixMCGhoaUlNTPR4PIaSxsVGhUHCbQKEgGh0d5dsHDx7UaDT+MTExMX19fWFMSnB0Ol1BQcE0ASaTSSaTORwOQsi9e/ekUunAwEC4shMijUZTVVXl35+RkdHY2Bj+fKIb94wgTdMAcOfOHb5/+/bt1dXVXHvPnj3vvPPOhB88derUU089xbIsIeSrr75auXJlgGecQd8Ir127VlhYyLU3b9587do1/4ANGzZwH8rWrl3rdrv/+OOPcGcZ7aRSKd8eGxubahGTL7744tixYzdv3gxXXoLT399/9OjR2tpaq9Xq/6rBYMjJyUlMTAQAlUqlVqsnDOKhIBoaGtLr9WVlZZO++vXXX3/44Ye+Q9noEU01zn/9+vXNmzdz7UlrBBdAURQX0NbWZrfbAzrjI2QbZP39/fy6awqFwmq1kv992N9qtfIBIpFILpf39fWFO0vBaG1tPX36dGVlpf9LeXl5Dx48MJvNRUVFhw8fDn9uUS8+Pj4zM9PpdF66dCkjI6O5uXlCgO+1AAAKhQKvhdA5c+ZMVlbW0qVL/V/KysoCgL6+vn379mm12rCnJiBjY2MOh8O3RviPRfsWkXnz5kkkkgDHq2fQZJlZs2Zx34UBgKZpiUTCFXaeRCLh52gAgNfrfeyxx8KaomB0dXXt2LHjo48+mnTK0g8//MA1ysvLc3JyDhw4oFAowptglCsqKioqKuLaVVVVhw8f/u6773wD8FoIp9OnTx84cGDSl/gNxisrKxcvXnzr1q01a9aEMTUB4RYX9K0R/u95iUTCBzAMwzBMgNfFDPpGmJqayn+q7e3tTU1N9Q/g7456PB6bzeY/lQ49uq6urg0bNrz11lvl5eXTR2ZlZUmlUtw2OaTWrVtnsVgmdPpeCwDQ29uL10KI3Lp1q6Oj47nnnps+LCUlRa1Wd3Z2hicrAZo1a5ZcLuff9pO+532LCNdQKpWBHHwGFcLi4uLz589z7fPnz/MTkQ0Gw9DQEBdw5coVbsxXr9c/8cQT6enpkco2WvX09OTn57/yyisvv/yyb39LS0tXVxcAjI2N8Z2XL19mGGbhwoXhzjLa8f/IhJBLly4tX76c+2NTUxP3H0FhYaHRaLx79y7X6XK5nn322UhlG91OnTq1e/fuOXPm8D3t7e0mkwkAPB4Py7J8p9lsnvRxIxQsxcXFdXV1AEAIqaur42oEy7JXrlwZGRnhAvR6PXf51NXV5efnB7pVy6PO7wkem822ePHiLVu2bN++XaVS9fT0cP1JSUl6vZ5rl5aWpqen79u3Ty6XX7x4MXLJRq3du3dLpdLsv2m1Wq4/Ly+vpqaGEKLT6ZYsWfLCCy9s27YtLi7us88+i2i+0WnHjh0bN27UarWrVq1Sq9Vms5nrz8rK+uSTT7j2kSNHVCpVWVlZcnLyxx9/HLlko5nb7U5MTDQYDL6dpaWlFRUVhJCmpiaVSqXRaHbt2jVnzpy33347QmlGoU2bNq1evRoAli9fnp2dPTw8TAjp6OhITk7WaDQbN27MzMy8f/8+IWR4eBgA2traCCE0TRcWFmZnZ+/du3fevHk3btwI8HQza/cJt9vd0NDAMExBQUF8fDzX2dLSolaruQlyANDY2Njd3f3MM8+o1erIZRq17t69yz8gDACxsbEZGRkA0N7enpCQoFQqaZq+ffu22WyOi4vLyckJcOQB/SNOp/PmzZt2u12pVK5bt46/z2E0GuVyOX9H9vfffzeZTCtXrsQvIiEyMjLS3t6+evVq3/kKnZ2dFEWlpaWxLGs0Gtvb2yUSSWZm5oIFCyKYapRpbW3l7/YBwKpVq7h9r5xOZ0NDQ2xs7KZNm7iN4ViWbW5uXrFiBbfXEsMwV69etdlseXl5ycnJAZ5uZhVChBBCKMxm0D1ChBBCKPywECKEEBI0LIQIIYQEDQshQgghQcNCiBBCSNCwECKEEBI0LIQIIYQEDQshQgghQcNCiBBCSNCwECKEEBI0LIQIIYQE7b+jVSXpXP//tgAAAABJRU5ErkJggg==",
"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 $k$-points. Here, we just\n",
"have one $k$-point:"
],
"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 k-point, 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": "4.33442190517645e-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.00022340105025458898"
},
"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.8.2"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.2",
"language": "julia"
}
},
"nbformat": 4
}