{
"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.796554e+02 1.621610e+02\n",
" * time: 0.00046706199645996094\n",
" 1 1.765557e+02 1.381380e+02\n",
" * time: 0.0013620853424072266\n",
" 2 1.322890e+02 1.575009e+02\n",
" * time: 0.0027070045471191406\n",
" 3 1.291856e+02 1.438058e+02\n",
" * time: 0.004106998443603516\n",
" 4 8.519556e+01 1.400220e+02\n",
" * time: 0.005408048629760742\n",
" 5 3.577844e+01 7.357446e+01\n",
" * time: 0.006927967071533203\n",
" 6 1.253144e+01 1.527640e+01\n",
" * time: 0.008049964904785156\n",
" 7 8.806660e+00 3.319643e+01\n",
" * time: 0.00855112075805664\n",
" 8 4.613286e+00 6.924283e+00\n",
" * time: 0.009604930877685547\n",
" 9 3.816805e+00 2.065872e+01\n",
" * time: 0.01052093505859375\n",
" 10 2.969985e+00 1.412308e+01\n",
" * time: 0.01155400276184082\n",
" 11 2.823231e+00 4.246204e+00\n",
" * time: 0.01222991943359375\n",
" 12 2.804328e+00 1.095390e+01\n",
" * time: 0.012892007827758789\n",
" 13 2.801366e+00 9.364745e+00\n",
" * time: 0.01359415054321289\n",
" 14 2.503603e+00 6.227278e+00\n",
" * time: 0.01429605484008789\n",
" 15 2.007359e+00 3.262992e+00\n",
" * time: 0.01511693000793457\n",
" 16 1.699088e+00 2.141388e+00\n",
" * time: 0.015716075897216797\n",
" 17 1.615603e+00 2.562667e+00\n",
" * time: 0.016309022903442383\n",
" 18 1.324555e+00 1.425976e+00\n",
" * time: 0.016967058181762695\n",
" 19 1.211161e+00 8.615679e-01\n",
" * time: 0.01774311065673828\n",
" 20 1.156973e+00 7.266140e-01\n",
" * time: 0.01848912239074707\n",
" 21 1.149926e+00 2.129798e-01\n",
" * time: 0.019238948822021484\n",
" 22 1.147934e+00 4.172996e-01\n",
" * time: 0.019976139068603516\n",
" 23 1.145721e+00 1.601387e-01\n",
" * time: 0.02072906494140625\n",
" 24 1.144693e+00 4.780865e-02\n",
" * time: 0.021558046340942383\n",
" 25 1.144553e+00 6.495186e-02\n",
" * time: 0.02232813835144043\n",
" 26 1.144236e+00 5.347165e-02\n",
" * time: 0.023051977157592773\n",
" 27 1.144113e+00 5.457617e-02\n",
" * time: 0.023749113082885742\n",
" 28 1.144065e+00 2.518325e-02\n",
" * time: 0.02445697784423828\n",
" 29 1.144053e+00 1.149139e-02\n",
" * time: 0.025213003158569336\n",
" 30 1.144042e+00 5.884988e-03\n",
" * time: 0.025954008102416992\n",
" 31 1.144039e+00 3.964863e-03\n",
" * time: 0.026658058166503906\n",
" 32 1.144038e+00 2.576994e-03\n",
" * time: 0.027341127395629883\n",
" 33 1.144037e+00 2.183956e-03\n",
" * time: 0.028028011322021484\n",
" 34 1.144037e+00 1.397237e-03\n",
" * time: 0.06602597236633301\n",
" 35 1.144037e+00 9.794050e-04\n",
" * time: 0.06675601005554199\n",
" 36 1.144037e+00 9.784813e-04\n",
" * time: 0.06742000579833984\n",
" 37 1.144037e+00 1.000341e-03\n",
" * time: 0.0681149959564209\n",
" 38 1.144037e+00 6.879353e-04\n",
" * time: 0.06882596015930176\n",
" 39 1.144037e+00 3.680227e-04\n",
" * time: 0.06939697265625\n",
" 40 1.144037e+00 3.060667e-04\n",
" * time: 0.0699150562286377\n",
" 41 1.144037e+00 1.928656e-04\n",
" * time: 0.07056808471679688\n",
" 42 1.144037e+00 9.736329e-05\n",
" * time: 0.07121396064758301\n",
" 43 1.144037e+00 3.826276e-05\n",
" * time: 0.0718541145324707\n",
" 44 1.144037e+00 3.344021e-05\n",
" * time: 0.0724949836730957\n",
" 45 1.144037e+00 3.469809e-05\n",
" * time: 0.07315897941589355\n",
" 46 1.144037e+00 3.149593e-05\n",
" * time: 0.07380795478820801\n",
" 47 1.144037e+00 2.086549e-05\n",
" * time: 0.07444000244140625\n",
" 48 1.144037e+00 1.217044e-05\n",
" * time: 0.07503604888916016\n",
" 49 1.144037e+00 9.643829e-06\n",
" * time: 0.0756380558013916\n",
" 50 1.144037e+00 3.833179e-06\n",
" * time: 0.07624197006225586\n",
" 51 1.144037e+00 2.246018e-06\n",
" * time: 0.07684206962585449\n",
" 52 1.144037e+00 1.864299e-06\n",
" * time: 0.07745003700256348\n",
" 53 1.144037e+00 1.484146e-06\n",
" * time: 0.07799005508422852\n",
" 54 1.144037e+00 9.633378e-07\n",
" * time: 0.07853198051452637\n",
" 55 1.144037e+00 3.599088e-07\n",
" * time: 0.07897615432739258\n",
" 56 1.144037e+00 3.506398e-07\n",
" * time: 0.0795290470123291\n",
" 57 1.144037e+00 2.020302e-07\n",
" * time: 0.08013701438903809\n",
" 58 1.144037e+00 1.125714e-07\n",
" * time: 0.08073902130126953\n",
" 59 1.144037e+00 1.207900e-07\n",
" * time: 0.0814211368560791\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.024068132162509e-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": "2.0013135757351903e-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.00022342068330310092"
},
"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
}