{
"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.659236e+02 1.273413e+02\n",
" * time: 0.0005431175231933594\n",
" 1 1.629640e+02 1.076857e+02\n",
" * time: 0.0017230510711669922\n",
" 2 1.189965e+02 1.233141e+02\n",
" * time: 0.0031621456146240234\n",
" 3 6.382899e+01 8.841720e+01\n",
" * time: 0.0046901702880859375\n",
" 4 5.285382e+01 8.996111e+01\n",
" * time: 0.006206035614013672\n",
" 5 3.672449e+01 6.707361e+01\n",
" * time: 0.007630109786987305\n",
" 6 9.657292e+00 2.334245e+01\n",
" * time: 0.009122133255004883\n",
" 7 5.714108e+00 5.499301e+00\n",
" * time: 0.010296106338500977\n",
" 8 5.460833e+00 6.627074e+00\n",
" * time: 0.011226177215576172\n",
" 9 4.599566e+00 4.976298e+00\n",
" * time: 0.012139081954956055\n",
" 10 3.679981e+00 3.796006e+00\n",
" * time: 0.013181209564208984\n",
" 11 3.566766e+00 9.160392e+00\n",
" * time: 0.013978004455566406\n",
" 12 2.671424e+00 4.426654e+00\n",
" * time: 0.01498723030090332\n",
" 13 1.940267e+00 1.816700e+00\n",
" * time: 0.015896081924438477\n",
" 14 1.709309e+00 2.622563e+00\n",
" * time: 0.016662120819091797\n",
" 15 1.482670e+00 2.091593e+00\n",
" * time: 0.01723504066467285\n",
" 16 1.291693e+00 1.296827e+00\n",
" * time: 0.01830005645751953\n",
" 17 1.200778e+00 7.728255e-01\n",
" * time: 0.019262075424194336\n",
" 18 1.156592e+00 5.778878e-01\n",
" * time: 0.02006220817565918\n",
" 19 1.152121e+00 4.035080e-01\n",
" * time: 0.02081918716430664\n",
" 20 1.149436e+00 2.725012e-01\n",
" * time: 0.021610021591186523\n",
" 21 1.146264e+00 1.498671e-01\n",
" * time: 0.022431135177612305\n",
" 22 1.144894e+00 6.952439e-02\n",
" * time: 0.02320718765258789\n",
" 23 1.144380e+00 4.165815e-02\n",
" * time: 0.02394700050354004\n",
" 24 1.144285e+00 5.653299e-02\n",
" * time: 0.024468183517456055\n",
" 25 1.144139e+00 2.574190e-02\n",
" * time: 0.025199174880981445\n",
" 26 1.144084e+00 1.957682e-02\n",
" * time: 0.026006221771240234\n",
" 27 1.144044e+00 1.027799e-02\n",
" * time: 0.026789188385009766\n",
" 28 1.144040e+00 4.777109e-03\n",
" * time: 0.02761220932006836\n",
" 29 1.144038e+00 2.389800e-03\n",
" * time: 0.028562068939208984\n",
" 30 1.144037e+00 1.594167e-03\n",
" * time: 0.02924323081970215\n",
" 31 1.144037e+00 1.154946e-03\n",
" * time: 0.03010106086730957\n",
" 32 1.144037e+00 9.267397e-04\n",
" * time: 0.030907154083251953\n",
" 33 1.144037e+00 5.423600e-04\n",
" * time: 0.03166913986206055\n",
" 34 1.144037e+00 4.033488e-04\n",
" * time: 0.032476186752319336\n",
" 35 1.144037e+00 3.783241e-04\n",
" * time: 0.033306121826171875\n",
" 36 1.144037e+00 2.676693e-04\n",
" * time: 0.03416800498962402\n",
" 37 1.144037e+00 1.821987e-04\n",
" * time: 0.03498411178588867\n",
" 38 1.144037e+00 1.222264e-04\n",
" * time: 0.035822153091430664\n",
" 39 1.144037e+00 7.712673e-05\n",
" * time: 0.036646127700805664\n",
" 40 1.144037e+00 5.254643e-05\n",
" * time: 0.037484169006347656\n",
" 41 1.144037e+00 2.344311e-05\n",
" * time: 0.038346052169799805\n",
" 42 1.144037e+00 2.212903e-05\n",
" * time: 0.03913617134094238\n",
" 43 1.144037e+00 1.542135e-05\n",
" * time: 0.03997921943664551\n",
" 44 1.144037e+00 9.354198e-06\n",
" * time: 0.04082822799682617\n",
" 45 1.144037e+00 6.800788e-06\n",
" * time: 0.041671037673950195\n",
" 46 1.144037e+00 4.920860e-06\n",
" * time: 0.04252910614013672\n",
" 47 1.144037e+00 3.296588e-06\n",
" * time: 0.04340004920959473\n",
" 48 1.144037e+00 1.754460e-06\n",
" * time: 0.04426121711730957\n",
" 49 1.144037e+00 1.175943e-06\n",
" * time: 0.04513120651245117\n",
" 50 1.144037e+00 8.524324e-07\n",
" * time: 0.04598808288574219\n",
" 51 1.144037e+00 5.671743e-07\n",
" * time: 0.04689311981201172\n",
" 52 1.144037e+00 3.244964e-07\n",
" * time: 0.047750234603881836\n",
" 53 1.144037e+00 1.671850e-07\n",
" * time: 0.048693180084228516\n",
" 54 1.144037e+00 9.836016e-08\n",
" * time: 0.04957318305969238\n",
" 55 1.144037e+00 1.091886e-07\n",
" * time: 0.05047917366027832\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": "9.290622208760838e-16"
},
"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+gvaeTAAAgAElEQVR4nOzdeXxU1dkH8OfemUkmO2SdbGRhy76xk4R9FUGo4opo3arWV63aqtVa0bbWBZdqq6321YILVSkigiA7JEAgBEI2kpAQsk5msu/JzNzz/jG+FFkTOHfW3/evJFyee/PJvfnlnHPPOQJjjAAAAJyVaO0LAAAAsCYEIQAAODUEIQAAODUEIQAAODUEIQAAODUEIQAAODUEIQAAODUEIQAAODUEIQAAODUEIQAAODUrBGFbW9uZM2cGf7zJZJLvYmCQJEmy9iUASZKENRGtjjGGn4LVMcY4/lKyQhBu2bLl2WefHfzxPT098l0MDFJPTw8efqvr7+/HXyRWNzAwYDQarX0Vzs5oNA4MDPCqhq5RAABwaghCAABwaghCAABwaghCAABwaghCAABwaghCAABwaghCAABwaghCAABwakprXwAAXAEjOtHCyttZWYvo7sJG+bBEX4rwFKx9XQAOAkEIYLta+umjUul/SyVGlDBcCHMjYx9trzcdbWIxw4R7x4q3RYtKdOsAXBsEIYCN+nel9KtDpgVh4ifTFZMDBSLq7R1wcVEqFAqDRN9VS+8WSW8VSP+cpkj1Q+vQDlRXVz/yyCPWvgo79uijj86ZM0eOyghCAJvTbaSVe0xl7WzDHOWkwIuEnEqkZZHiskhxTbm0cKvxf+IVz6WgYWjrOjo68vPz3333XWtfiF169913h7Rbw5AgCAFsS3M/LdpmjB8mHF2mdLlSuq0cLc4LE6/fZqzpYn9NVyjQMrRtXl5eS5YssfZV2KWNGzfKVxx/RQLYkPoelrnJOCtY+Gia4oopaKZxo92LlBWd7JZdJiO2pgAYOgQhgK3oMtD120wrRol/mjC0pp2XijbPV/YZ2S8PYPNOgCFDEALYBInRij2mNH/ht1c12uci0pezlceb2esn0CoEGBoEIYBNeCrH1GtkH6QrrrqCu5L+M0fxXrG08QyyEGAIEIQA1vd9DftPFfv3bOU1TgoM9RC+nKX4RZapvodxujQAx4cgBLCypj66P8v08XTFMBcO1SYFCg/HKe7ea0ISAgwSghDAyu7db1oxSpgZzG3qw2+TxW4jvVOIDlKAQUEQAljTvyulM53spXFXPzR4IaVIa2co/njcVNWJZiHAlSEIAaymy0BP5Uh/mTrYKYODF+0lPJ6g+NUhNApBLi+++OLZSe4DAwN33nlnX1/fpQ6WJOmuu+7q6Oiw1NUNDYIQwGpeOmaaEypM08iyHsyvk8SSNra5Bo1CkEVJSUlDQ4P54w8++MDPz0+tVl/qYFEUY2JiXnvtNUtd3dAgCAGso6SNfVIm/XkCz07Rc7mI9M4UxeMHTX2YZA+Xlp2d3dHRsW/fvnXr1pm/cuTIkc8//zw/P//sMW1tbVu2bFm7du2JEycurMAYe/fdd++9914i6ujoyM7OJqLKysqysjIiOnz4cHNzMxHdddddH3744WVajVaEIASwjqcPS8+lKILcZDzF/DAhYbjwt2J0kMIlrVy5csmSJe+88052drYkSbfccsvjjz9+9OjRlStXvvDCC+Zjbrzxxo0bNx4/fnz58uUvv/zyeRUKCgq6uroSExOJ6OTJk3fddRcRrVmz5v333yeihx9++OjRo0QUEhISHByclZVl0W9vcLDoNoAVHNazglb21Wy5moNn/XmiOOM7430xordK7lPB1Vh/WtrTYLnu67cmKy6cq5qZmWmOt88//7yuri4rK0sQhM7OztGjR99///3h4eE7d+40H/nMM89ERUU988wzKtV/76e8vLz4+PjBnD0hISE3N1emrZSuBYIQwAqePWL6faroKnsO0lgfYX6Y+Hah9EIqun9sUYCbMMbHcqcTLjYePXfuXPMH+/fv7+npefDBB8/+U2FhYXh4+JYtW957770zZ8709fX19fXV1dVFRkaePaa1tdXb23swZ/fx8WlpabmW65cJghDA0n6oY7XdtGKUhZJp1Thx3AbjgzFioJzdsHB1pmnkeltq8Nzd3c0fmEym+Pj45cuXmz9dvnx5UlJSRUXF3XffvWnTprS0NKVS6eXlZTAYzv3vPj4+nZ2d5o9VKtV5/2owGM42Hzs6OsLCwuT9Zq4K/kgEsLTf5Zr+OF68xtXUBi/CU7h1pPhGAd6ZgSvIyMgoKiqaMWPGnDlz5syZM3v27ICAgPLy8hEjRkyaNEmlUu3du7e7u/u8/5WSklJUVGT+ODo6WqfT1dbWmj/V6/UVFRWjRo0yf1pUVJSammqxb2fw0CIEsKgddazbSDdGWfRv0GeSxZT/GJ9NVgx3teRpwc6sWLFiy5YtEyZMmDNnTk9Pz/bt2w8cODBp0iStVrty5Up/f/8jR45c2AuakpLi4uJSUlISGxvr4+OzatWq9PR0X1/fgYGBb7/99le/+lV4eDgR6XS6mpqaadOmWeM7uwK0CAEs6s/5pqeTRQv3hYV5CEsixL+V4PVRON+aNWvGjBlj/lgUxXXr1v3zn/9MSEiYNWvWnj17/P39hw8fnp+fP3/+/PT09G3btm3atCk0NJSIXnzxxRtuuMH8vx555JGPP/7YXOQ3v/nNvn37xowZExUVtWPHjrNvma5du/buu+8+2w1rU9AiBLCc482srJ1uibbCH6BPJ4szvjM+kSC64aGHc6Snp5/3lbS0tLS0tHO/4ufnd8cdd5g/Ptuki42NPXvAI488csstt/T29rq5uRFRREREbGxse3t7VFSU+QBJkg4cOPDhhx/K9F1cIzwTAJbzx+PSk4ki9wXVBmOsjzApUPykXHooFv1AwJlarT673JpZdHT0uaOJoiiuX7/e4tc1WAhCAAs51cH2NkgfT7PahL6nk8U795geiBEVVn5LERzfypUrrX0JQ4C/DQEs5L1i6f4Y0dN6E9unBAqBavquGiOFAD+BIASwhC4DrS2XHoix8hP3yzjxr1hxDeCnEIQAlvDpKWlGsBjhaeVOyeXRYkELO9mGLSngR0888URVVZVMxf/0pz/pdLrLHPDHP/7x8gdYBoIQwBLeL5F+GWf9x81FpPtixPcxjwL+X0dHh9FolKPy3r17d+7cGRgYeJlj/P39V61aJcfZh8T6TyaAw9vbwAwSzQyxiXdUHooVPz0ldRiufCQ4g5deeikiIoKI9Hp9W1tbU1PTvn37zBsn9fb2ZmVl1dXVnT3YaDQeP348Ozu7tbX13CJNTU379+9vaWlpbW01/18iWr169f33309EJpOpsrKSiDo7OxsbG4mosbHRvEnvihUrvvrqK6svQIogBJDd30qkX8ZZehL9pYS4CzNDxM9PoVEIREQZGRnmXQafffbZFStWLFy48JVXXhkzZsyWLVtmzpz52muvJSYmbtu2zXxwQkLC888/v3r16ri4uG+//db8xQ0bNsTFxb3zzjvz58+/5ZZb/vCHPxBRV1fXDz/8sGDBAiJqamoyz9lfv379Qw89REQPP/yweftDDw+PcePGbd682Rrf+n9h+gSAvJr76Yda6e8ZNrQN0gMx4vO5pgcxodAG9JcfH6g5ZbHTec34GYmX/Lm3tLQcPHhQqVT+4he/uPfee4uKinx9ff/+97+//fbb8+fPJ6Ljx4+bt6E/cODAz3/+8yVLlgwMDDz88MNffvnljBkz+vr6kpKSzFsyHTt2LDg4eNiwYVe8pKSkpMOHD9955538vsshQxACyGttuXT9CHGYi7Wv4xxzQoQHeim/hSX72kgz1XlJfb1SV5sFT3i596QWLlyoVCqJKD4+XqfT+fr6mj/+y1/+Yj4gJydnzZo19fX1AwMDp06dMhgMFRUVAwMDM2bMICK1Wr1w4ULzkXq93vzfr8jPz0++t3UGCUEIIK9PyqS3psi/8eBQiALdNVr4pEx6a7JtXZgTckuc4pY4xdpX8SNza4+IFAqFebE088cmk4mI8vPzb7/99n/84x+xsbGdnZ0pKSkGg0GSJPGcJqZCoWCMEZGnp+fZlWXMX5Sk//bGm0wmheLHe6+7u9vLy0v+b+5y0DcCIKPcJtZhoOnW3nDuQj8fI352SurH1kwwaLm5uZMnT160aFF0dPTZjZZGjhxpMpmOHj1KRCaTaceOHeavJyYmVlVV9ff3E5Gfn5+Hh0dJSYn5n0wmU1FR0dllSEtLSxMTEy39zfwUWoQAMvq4TPr5GNFW3pM5R6SXkOgrbKqWbrLshlBgv6ZMmfLEE0+88soroih+99135i+q1eo33nhj8eLFS5cuLSoq8vT0NDf1goOD4+LiDh06NH36dEEQfv/73y9dujQ2NraiomLhwoUjRoyYPn06EUmSlJWV9dJLL1nzG0OLEEA+fSb6d4V012jbi0EiIrp3rPi/ZXh31Nm9+eab5sbZgw8+uGTJEvMX58+f/9hjj5k/HjNmzNtvv01EcXFxe/bsMRqNXl5eGzdu/PLLL11dXYnonnvu2bt37+zZsz/++OOxY8ee3YP+4YcfXrNmjfnjJ5988uuvv/bx8RFF8Yknnti+fbs5L3fs2BEXF3d2HyhrQYsQQC6bqqVUf2GEtVeTuZRlEeIjB0zaXtK4WftSwHqWLl1q/mD8+PFnvzhq1Kiz28r7+fmZZ0EQUWpq6tkt5pcvX27+4MiRI6IoZmZm7tu3b8OGDS+88IL563fffffHH3/c0NAQHBxMRMnJybNnz+7u7j5bjYjef//9V155RcZvb3AQhABy+fwUu2Ok7Xa6uClp8Qjx3xXSYwm2e5Fg+3p6el599VW9Xj9ixIhNmzadHfxTKpXZ2dnnHunn5zdy5Mhzv7JhwwbLXeilIQgBZNHaT7sbpE+m29D0wQvdMUr8Xa4JQQjXYvr06eYBvytavHjx4sWL5b6eq4AHAEAWX5+W5oWKPrY0ffBCs0OE6i5W1o41uMGpIQgBZPFZhXT7KBsdHTxLIdDN0eIXFQhCcGoIQgD+6ntYYQtbGGYHz9cdo8RPse4oODc7eFAB7M7nFWxZpOhqD8u2TAwQBIGO6NEoBOeFIATg798V0m02/L7oeW6JFr46jUahU2toaDBvkOSc7OZZBbAXVZ3sTBebZnvLql3KTVHiV6cZmoTOqb+//6abbsrMzMzMzJwxY8Z5Gw06CUyfAODsq9NsWaSotJ8/MpN9BVeR8prYOH+7CW+HUd5aWdtRb7HTTR+RLgo/+Sm/8cYbRqOxrKyMiMy7Ca5evdpi12MjEIQAnH19WvrDeHsYHjzHzyKFr09L4/zt7LIdwKnW04fqci12uszwKaLwk5/yd9999/zzz5t3kHjggQeeeOIJi12M7UAQAvBU080qO9nMYDtrWt0UJS7faXplgrWvw/ksjJ69MHq2FS9Ar9d7e3ubP/b29m5vb7fixViL/XTfANiDryrZ0gh76hc1S/MXBIGON2Og0Okwxo4fP27++NixY2PHjrXu9VgFWoQAPK2vkn6XapcdjDdGCl+dllL87PLi4Vp88MEHXl5eoiiuWrXq448/tvblWIG9/eEKYMMaeuhkG5sdYmf9omY3RonfVKFF6IzefPPN8vLygwcPfvbZZ+duDeE80CIE4Obbaum6cFFln39eTggQOgxU1s7G+NhlkMNV8/Pz++Mf/2jtq7Am+3xkAWzSxjPSDRH2miIC0eIRwsYzaBQ6Fx8fH6XS2VtECEIAProMdKCRzbeH9UUv5YYIceMZLDHjXPLy8lJSUqx9FVZmxw8tgE3ZUiOlBwleNr3/4BXMChFK2pi219rXAWBZCEIAPjaeYTdE2PcDpRJpXpi4uRqNQnAu9v3cAtgIg0TbaqXFI+z+gboBw4TgfOz+uQWwBfu0bLSPEOxu7eu4ZgvDxX1aqdto7esAsCAEIQAHm6odoTlIRD4uNN5f2FmH3lFwIo7w6AJY3ZYatmiEvU6cOM+iEeLmGvSOghNx9ukjANeutJ11GyjJ10GCcMkI4fUTpg9I4SDfjy3p7e3Nycmx9lXYJb1eL1/xIQRhZWXlDz/8MHz48CVLlri5uV30mLy8vEOHDnl5ec2cOTMsLIzTRQLYtM3V7PoRgsPExkhvwVMl5DezFD+H+Z5sQmRk5LRp0x599FFrX4i9CgwMlKnyYIMwKytr8eLFK1asKC0tff3117Ozs11dXc875qmnnlq3bt3111/f09NTXl7+0ksv8b5aAFu0uUZ6LN6hRhkWhQubaxCEnHl6ejrnkta2b7BBuGrVqueee+6pp54ymUzjx4//+uuv77jjjnMP2LZt26efflpYWOjv7y/DdQLYqA4DHdGzWSEOFoTiC0dNz6U41DcFcCmDutH7+/t37ty5dOlSIlIoFIsXL/7+++/PO+arr75auXJlc3Pzt99+W1VVxf1CAWzT9lopPUjwtOcFZS40LVgobmM6LDEDzmFQLUKtVssYCwkJMX8aEhKSlZV13jGVlZWFhYX79u2LiYm555573nzzzZUrV16q4KlTp85b7DwyMvLmm2++6MEGg8FgMAzmOkE+5p+C4DgDYdxsrqZ5IWSZW9T8I5Ak2ec2CEQzNbSl2nBHtNynsj8Gg0EU0Va2MoPBYDKZFIorb5+pVCqv+ItrUEHIGCOis7VEUbzwURwYGBgYGMjNzRVFccuWLbfffvsdd9xxqas0GAxtbW3nfqWzs/NSj7ckSRZ48uHyzD8FBOF5GNG2OvHX8Ra6Qy35U1gQKmyrFW6LxKN3PvMPG7+UrEv6f1yqDSoINRqNIAiNjY2RkZFE1NDQcLZ1eFZISIifn5/5D6X09PT29natVhsaGnrRgrGxsa+//vogL3FgYODCF3PAwgwGg6urK4LwPCdamIfKFOtvoftTkiQXF5fB/BV87RZFst8dM6pcXEX8zC8giqJK5Vi94fZGFEWTycQrGgbVwFer1enp6Vu2bCEixtjWrVtnz55NREajsbq62pzJ8+bNKykpMR9fXFzs5uYm36uuADZiay1bEOaYQRHmIQS6CXnNmFkPjm+wb40+//zzt912m06nKykpaWtru/XWW4no9OnTY8aM0ev1/v7+t99++5tvvnnnnXcmJiZ+8MEHL774Iv5iAoe3rVZ6ItES7TOrWBAmbK1h4/0dM+kBzhrskO/8+fN3796tUqlmzpx58OBBDw8PItJoNJ999pmXlxcRubu75+TkZGZmSpK0du3a3/zmNzJeNYAN6DZSrp5N1zhsTswPE7dh0VFwAkNYWSY5OTk5Ofncr3h5ed1+++3nfvrAAw9wuzQA27a7nk0IcLSJE+eaphEKW1jbAA1zsfalAMgJLwEDXKVttdL8MEd+glwVNCVI2FWPRiE4OEd+jAFk5cBvypw1P1TcVov3ZcDBIQgBrkZFB+sxUoKj7DhxKQvCha0IQnB0CEKAq/FDHZsb6vjTKsf6CKJApe3IQnBkCEKAq7Gjjs0NdfgcJCKaFSxsr0MQgiNDEAIMmYnRngbJwXacuJQ5ocIOBCE4NKd4kgH4ytWzUA8h2N3a12ERc0PFPQ2SAa+OguNCEAIM2Xan6RclIn81RXkJR/RoFILDQhACDNmOemluqBM9O3NDMUwIjsyJHmYALrqNlNfEMh13ZbULzQkVd2BaPTguBCHA0OxtYOP8BY8hrE5o96ZphPxm1oHtscFBIQgBhmZHnTTHmfpFiUitoEmBwt4GNArBMTnTn7UAPOyqZ3/PsGgQDpgGdlTtK2kuO9lU7uHiPtZ3VEpQ4pTQ8Za8hlkh4u56tniEJc8JYCEIQoAhaO6nqi42zoJb9O2vOfjXvP+N8okYH5wyLXiySZROtVV9lL92XcmGR8ffP3JYpGUuY1aw8EAWWoTgmBCEAEOws06aphGVFmkQSkx69dBfSlsqfjPpf9I0SUTU29vr4uIyOXT87XE/++7UD0/ufOG+5DuuHzXfAhczPkCo6WaNvRTkZoGzAVgUghBgCHY3sFkhlmgOSkx65eDbbf0df1+w2lVx/n6AoiAuGb1gYkjar3Y832vsXx6zRO7rUQiUESTubZBujnau8VFwBrinAYZgZ70lglBi7Pf7X+029Pxp2nMXpuBZGo/At+f88ZuyLetLN8l9SUQ0K0TYVY/ZhOCAEIQAg1XXzVr7WcJw2YNwXcl/WvvaX8p8VqVQXf7III+At+b84dOirwv0xXJf1awQYVcDghAcEIIQYLB21LNZIaIocw4WN5V9WbLxd+lPKEXFYI4PdPd/ZvKjL2W90d7fIeuFJfoKHQPsTBeyEBwNghBgsHbJ3y/aMdD5YtarT09+NMgjcPD/a1LIuFkRma8cfFu+CyMigWhGsLgHjUJwOAhCgMHaXc9mBcsbhB8d/zQ9bOJVzBG8P2VlU0/LrjP75biqs2aFCDux6Cg4HAQhwKCc6mCMaLSPjEF4qrVyX82BnyfefhX/VykqHp/w4Pt5H/cZ+7hf2FmzQgS0CMHxIAgBBmV3PZspc3Pw3dyP7ku+09vV6+r+e0JATHJgwufF6/le1blGeQsSUUUHshAcCoIQYFD2atl0OYNwR9XeXmPfdSPnXkuRB1Pv+qbs+4auRl5XdaEZwWgUgqNBEAIMyp4GNlO2N2VMzPTP/M8eGXevKFzTKfzd/W4ce/0nBV/wurALIQjB8SAIAa6svJ2JRNFecgXh9tN7NZ6BSYHx117qppglB+ty6zobrr3URc0IxrR6cDQIQoArk7U5KDHps6KvVibcwqWah8r9htELvyj+D5dqFxrlLShFOoVhQnAgCEKAK9vTIOMA4e7qbE8Xz9SgRF4Fl8cu2Vt9oLFbx6vgeaZr0DsKDgVBCHBle7VyvTLKiH1a+OXdibdyrOnt4rVo1Nx1Jd9wrHmu6RgmBMeCIAS4grJ2JhJFyTNAmNtwXBQUk0LG8S17c8wN20/v6Rro5lvWbGaIsKseexOC40AQAlzBXjkHCDeUbf7Z2EXcy/q6DZ8cMv77yp3cKxNRtJegEoWydjQKwUEgCAGuYK+WzZCnX7SxW1+oPzk7IlOO4svGXrehbLPEZImr6RphnxZBCA4CQQhwBfu1bJpGliDcWP79/OhZaqVajuLx/jGeKo+j2uNyFJ8WLOzDMCE4CgQhwOVUdjITo5He/IPQYDJ8X7lzyaj53CufdcOYhRvKNstRGS+OgiNBEAJczt4GNl2e5uCemuxRw6LCvUPlKG42O2Jaof5kY7eee+XRPoKJUVUnshAcAYIQ4HL2ybbE6JaKHYtHy9gcJCK10nV25LStlbvkKD4tGMOE4CAQhACXs1eeqfT6nqaK1qopIUPed3CoFkTN2nZ6FyP+iTVNI+xF7yg4BAQhwCXVdrNuIxsjwx6EWyt3zYrIVClU3CufZ6zfKBeFS5G+lHvl6WgRgqNAEAJc0p4GNl0jytEx+sPpPfOjZ8pQ+CLmRc3Ydpp/72jsMKHDwOq6kYVg9xCEAJckU79ocVMpIxbrN4Z75YuaFzVzz5nsftMA37ICUUaQiEYhOAAEIcAlyTSDcGvlrgXRs7mXvRR/N9+xfqMO1B7mXnl6MIYJwREgCAEurrGX9H0sfjjnIDRKpj3V2fOiZvAte3nzo2f+cHoP97KZGiGrEUEIdg9BCHBx+7RShkbkPkJ4VHt8hHdYoLs/57qXlRE2OV9X2DnQxbdssq9Q38N0vXyrAlgaghDg4vZrWaYM/aK7q7NnRqRzL3t5bkp1mibpYN0RvmVFgaYECgd02IkC7BuCEODi5BggNEqmA7WHp4dP5Vt2MGaMSN99Jpt72UyNuB/vy4CdQxACXET7AFV2sFQ/zkFo7hf1d/fjW3Yw0sMmndAXdRt6+JbN1GD1bbB7CEKAi8hqZJMCBRXv58Mq/aJmbkp1alBidm0O37ITAoSydtZp4FsVwKIQhAAXsV8rZWo4Px1W7Bc1k6N31EWkNH/hAN4dBXuGIAS4iH0N/N+UyWvMH+EdapV+UbOpoRPzdYXce0enaYT9WrwvA3YMQQhwvl4jFbSyiQGcgzCrJiczfArfmkPirnJLCow7XJ/HtyzelwF7hyAEON8hPUvyFdyVPGsyYgfqjkwNncCz6NClh006UMd5iZkpQUJeM+sz8a0KYDkIQoDzydEvWtZc4aFyk3Ub3sGYGjrxYF2uUeKZWh5Kih0mHNGjUQj2CkEIcL6sRikjiPOjkV2Xkx42iW/Nq+DnNjzMK6RAX8y3bEaQgN5RsF8IQoCfMEp0WMemBvEfIEwPm8i35tVJD5vIfRJFpkbIasT7MmCvEIQAP3GsmUV4Cb6uPGs2duta+lpj/cbyLHq10sMm7a85xLdmpkY80MhMaBOCfUIQAvxEViPL4N4crM2ZGjpRFOTY4nfIoodFiIJ4uu0Mx5r+agp2FwpakIRglxCEAD8hx1rbB+qOTLWNflGzqWETsnkvwI1hQrBfCEKA/2JE2Y1SBtcg7DP2lTSVpQUlcax5jSaFjMupP8q3JvYmBPuFIAT4r7J25qYQwj14BmFe44kYv9HuKjeONa9RSmDCqdbKjoFOjjUzNcK+BrwvA3YJQQjwX3L0ix6qOzopZBzfmtfIReGSFBCfpz3BsWaUl6AUhYoONArB/iAIAf4ru5Hx7RcloiMNx2wtCEme3tEMDYYJwS4hCAH+a7+W8yuj1R21RskY6RPOsSYXU0LH59QfZcQzt9KDhGwME4IdQhAC/EjbS639LHYYzyA8VH90cuh4jgV5CfYMclO6VbRWcayJ92XATiEIAX60XytlaESRa8/o4fq8icFpPCvyw713NHG4oO1hul6OJQEsAUEI8KMsLUvn2i/aZ+wvbiodp0nmWJOjyaHjDnENQlGgKUFCNtZaA3uDIAT4Efc1ZfJ1haOHR9vUxIlzJQcmnGqt5LtPb3qQiGFCsDsIQgAioi4DlbWzcf48gzC34fj44BSOBflyVbjE+o3J1xVyrIlhQrBHCEIAIqIDOjbOX3BV8Kx5RGvTQUhE44NTjjQc51hwUoBQ1Mp6jBxLAsgOQQhARJSllfj2i7b0tjb1NI/1HcWxJnfjNSm5XIPQVY5VP+oAACAASURBVEFJvkIONukFu4IgBCD68U0Zno9DrvZ4miZJFGz6ERs1PLqjv1PX08SxJlbfBrtj008pgGUYJMpt4rwZb642f7zGpvtFiUgUhDRN0lFtPseaGRohS4sXR8GeIAgB6GgTG+kt+LjwrJmnPWHjA4Rm4zTJfHtHM4LEHB0zIgrBfiAIASib98SJ0+3VCkEM8dRwrCmTCcEpR7XHOa61NtyVwj2FE9ikF+wHghCAshs5T6U/2pA/ITiVY0H5BHkEerp4VrSe5lgzIwiTKMCeIAjB2Zk3403nuunEUW2+zS4oc6FxmmS+WzKla7D6NtgTBCE4u7J25q7kuRmvxKQCfXFKUAKvgnJLDUo81ljAsWB6EDbpBXuCIARnl8V766XSllOB7v7D1cM41pRVWlBSvq7IxEy8CkZ7CUpRqOxEoxDsA4IQnB33AcI87Yk0TRLHgnLzdvXSeASWNldwrDk1SMjCbEKwEwhCcHZZvHelP9ZYkBpkT0FIRGmapGONXIcJsUkv2A8EITi1xl5q6mNx/DbjNUjG4qbSpMA4XgUtIzUoie8wYQZahGA/hhCE69atu++++5577jmtVnuZw7788st33333mi8MwBKytFJ6kMBxM96SptJw71AvF09uFS0iJSihqOmkwWTgVTDZT6jvYc39vOoByGiwQfiXv/zlueeemz59enNzc0ZGRn//xW/wnJycRx555OWXX+Z3hQAyym5kU7kuMZrXeCI1KJFjQcvwULmP8A4raS7jVVAh0MQA4QA26QV7MKhfASaTafXq1X/961/vvPPODz74wMvLa/369Rce1t/f//DDD69atYr3RQLIhftmvMe0BWn2NkBolhqUmMd1mHAqNukFOzGoIKytra2pqZk5c6b50xkzZmRnZ1942Msvv7x06dK4ODsbHQGn1W2kkjY2IYBbEA6YBkpbTiXa2wChWZom6Vgj7016MUwI9kA5mIMaGho8PT1dXV3NnwYEBOTm5p53TH5+/rfffnvkyJFDhw5dsWBubu7Pfvazc78SHx//7LPPXvTg3t5ehYLrfqkwdD09PYIgCALPxpPV7W0Uk4aLUn9PD6eCJ5qKI71HsAGpZ4BXyZ/o7e01Go0yPQ4jPSJLm8vbOttcFHxWH0/ypPwWVUtnj9qxHt/+/n5RFFUqlbUvxKkZDAaTySRJV+57V6vVoniFJt+gglCtVg8MDJz9dGBgwM3N7dwDjEbjfffd98EHH5wNy8sLDQ299dZbz/1KYGCgWq2+6MEGg+FS/wQWYzQa1Wq1gwXh4RaWqWEc767i1rI0TaJ8tytjzMXFRaYgVJM6atiI0901yYHxnApS3DCpoNM1k+vsFKsTBAFBaHUKhcJkMg3mWbtiCtIggzA0NLS/v1+v1wcEBBBRTU1NaGjouQdUVlbm5+ffeeedRNTX19fS0jJy5MgdO3ZERUVdtGBwcPDNN988mFMTkSiKg/lOQFbmn4KDBWG2zvhYvELk985ovr7ojvgb5btdxf8nU/3kwIQT+qJUDbeXfTI07IBemB7iUM+v3D8FGAxRFBljvH4Kg6oSEBCQmZn52WefEVFbW9uWLVuWLVtGRC0tLRs2bCCiqKiokydPbt++ffv27atXr/bx8dm+fXtYWBiXSwSQg1Giwzqem/EaTIbS5vIE/1heBS0vJSjhuI7nMGF6EDbpBTswqBYhEf35z39eunTp7t27i4uL582bN2XKFCI6efLkz372M8aYSqWKjo42H1lTU6NQKM5+CmCbjrewEZ6C76D68geluLksatgId5XblQ+1VUkB8auyXjeYDCoFn36/TI14zz6TiZHCoboSwNEMNginTp168uTJQ4cOBQcHp6b+uNFaWlpaaWnpeUdOmjTpyJEjPK8RQAZZWs4rqx1vLEgOtJsdJy7KXeUW7h16sqU8MYDPi68BagpyE4paWZIvkhBs1xA6WH19fa+77rqzKUhEarV6zJgx5x2mVqsjIiL4XB2AbLivtX1cV2hHWy9dSmpg4nGukyiw1hrYPoz3gpPKbpQ4TqU3SMaTzdwaUlaUzHuYMEOD3erB1iEIwRmd6mCiIER6cQvCk81l4d6hHip3XgWtJTkwvqSpzChx25swI0jYjxYh2DYEITij/Vo2jesAYb6uOMXOBwjNPFTuoV7BpS3lvAqO9hFMjJ3pQhaC7UIQgjPK5r3EaL6u0N7flDkrOTD+hK6YY8H0IBHDhGDLEITgjPZzfWVUYlJxU2ligB3PIDxXUmD8CV0Rx4IZQRgmBJuGIASno++jxl4WP5xbEJa3Vga6+3u7evEqaF1JgXEn9MUS4zYRPkODYUKwaQhCcDr7tVJ6kMBxine+riiJ0/qctmCYq4+fm29l2xleBVP8hJoubNILtgtBCE4nS8syNDzv/BO6YkcKQiJKDozP59c7qhBoUiA26QXbhSAEp5PVyDL5vSnDiBXqi3nt2GAjkgLjOA8TavC+DNguBCE4F/NmvOP5bcZ7pr3WTenm7+bLq6AtSAlMyNcVMeIWXZkYJgQbhiAE53JIx1L8BI5bxZ7QFSXb/8pq5wlw93dVuNR21PMqOClAKGhlvUZe9QB4QhCCc9mvlfhOpT+hK062/5XVLpQclMBxmNBdSYnDhcN6NArBFiEIwbns17JMrm/K5OsKEwMdMAiTAmI5BiGhdxRsGIIQnIhBoiN6NiWQW4tQ260zSMYwrxBeBW1HYkBcgZ7n+jKZGnE/NukFm4QgBCeS18RGeQs+LtwKntA52vuiZ43wCesz9ul6mngVTA8SDumYEVEItgdBCE5kn5Zlch4gLEpyxH5RIhJISAiILdCX8Co43JUivITjLegdBZuDIAQnsp93EBboix1gD8JLSQyIK+C6+jaGCcE2IQjBWTCiA41SehC3e769v0PX0zRyeCSvgrYmKZD3MCH2JgSbhCAEZ1HUyvzUQjC/rXML9CUJ/rEKgd+cRBszxndkQ1dj10A3r4LTg8UsrYQkBFuDIARnwb9fVFecGOggWy9dlEJQjPUbVdjEbZgw2J28VEJpG6IQbAuCEJzFfi3nzXhP6B1q04mLSgqI5ztMOC1Y2IfeUbAxCEJwFllaNj2YWxD2Gfsr26pjfEfzKmibkgLj8rkGYQaGCcH2IAjBKVR0MEYU5cUtCEuay0YOi1QrXXkVtE1x/mNPtVYaTAZeBadphD0NCEKwLQhCcAr7uDYHyaFnEJ7LTake4RN2sqWcV8HRPoLEqKoTWQg2BEEIToH7VPoCfYkzBCERJQXEncAwITg0BCE4hX0NjOOmExKTSprLEvwd+ZXRsxID4gr5rS9DmE0ItgdBCI6vrpt1GljMMG5BeKr1dICbn7erF6+CtiwpMK5AXyIxbtGFFiHYGgQhOL69WjYtWOTYMVqgL3bIrZcuarh6mI+r95n2al4F44cLLf2svgdZCLYCQQiOb7+WZ78oEZ3QFScGOEW/qFliYNwJfmutCUQZQWIWGoVgMxCE4Pj2ch0gJKJCfYkDr7V9ocSA2AId12FCDXpHwYYgCMHB6ftI28uSfLkFYX2XlgQh2DOIV0HblxTAs0VIRNODhb2YTQg2A0EIDm5vg5QRJHIcITyhK052puYgEYV7hxokA8dNelP8hLoe1tTHqx7ANUEQgoPb28B5Kn2BvjjByYKQiBL8YzguOqoQaEqgsE+L7erBJiAIwcHt4R2EJ3TFSQ696cRFJQbE8u4dFdE7CjYCQQiOrLmfartZqh+3IGzv72jubYkeFsmroL1IDOS8W/0MDBOCzUAQgiPb2yClBwkKrgOE8QExouB0D84Y35Habl3HQCevgml+QlUXa+7nVQ/g6jnd8wxOZW8Dmx7M8yYv1JckOd8AIf3/Jr3FTaW8CipFmhIo7McwIdgABCE4Mv4DhM60psx5kgLiC7guOophQrARCEJwWC39VNXJ0vgNEPabBirbzjj8ZryXkhgQy3cbCgwTgo1AEILD2qeV0jWCkt89XtxU6gyb8V5KfEBMWUvFgGmAV8Hx/kJlJ2vBMCFYG4IQHNaeBjZdw/MOP6ErdpI9CC/KTamO8A4rbangVVAp0mQME4INQBCCw9pdz2aGcF5iNMGZ1tq+UGJgbAHX2YQzgsU96B0Fa0MQgmNq7qczXTwHCCUmFTWddKpNJy6UGMB5NuHMYGFXPYIQrAxBCI5pd72UEcRzgLCirSrA3c/H1ZtbRTuUFMB5k95x/kJ1FxYdBStDEIJj2t3AZobwHSAscqqtly7K1224t6sXx016lSKlBwl7GzBMCNaEIATHtLuezeS/xKizByERJQXG8110dGaIuBvDhGBVCEJwQLpe0vayZH4DhGReUyYwnmNBO5UUEJevK+JYcGawsBvDhGBVCEJwQLvqpekakeMSo3WdDSQIGo9AbhXtVlIg5yBM8RMae5m2l2NJgKFBEIID2t3AeeLECX1xMpqDREQU5hXCmNTYreNVUBQoQyPuqccwIVgNghAc0J4GNoPzAGGRc661fVGJAXH5XCdRzAoRMEwIVoQgBEdT181a+lnCcL670mOA8L8SA+L4TqufFSLsxDAhWA+CEBzN9jo2O0QU+eVgW397W197pM8IbhXtXFJg3Amuw4Txw4VuA6vqRBaCdSAIwdHsrGez+Q4Q6ooTAmJEgWdNuzZqeJS+p7m9v4NXQYFoZoiIRiFYC4IQHM3uBjY7lPcAIfpFzyEKYnxADN+9CWejdxSsB0EIDqW4jSkFivbiGYT5uiK8Mnqe5IB4vr2jc0KFnfUSkhCsAkEIDmVnHZvLtTnYbeip7awf4zuKY00HkBwUz3eT3ghPwUslFLUiCsEKEITgULgPEBbqS2L8RqtEJceaDiDGd3RVe3WPgec0+Nkhws46BCFYAYIQHIeJ0X6tNIP3WttJAegXPZ9KoRrrO6q4qZRjTQwTgrUgCMFx5OpZqIegceNZM1+HNWUuLikw/oSe5zDhrBBxn1YyYIUZsDgEITiOHfVsDtd+0QHTQEXb6Tj/sRxrOoykwLj8Rp5B6K+maC8hR4dGIVgaghAcxw+10rwwnrd0cVNplE+EWunKsabDSAiILW05NWAa4FhzbqiwA4uOgsUhCMFBdBoor5llavhOnEC/6CW5KdURPuGlLRUca84NFX+oRYsQLA1BCA5iT4M0KUDw4Pp25wl9USI247205MB4vlsyZWiEwlbW2s+xJMCVIQjBQWyvY3NDed7PRslU0lSGTScuIykwPr+xkGNBtYKmBgl7GtA7ChaFIAQH8UMt56n0pS3loV7Bni4eHGs6mOTA+KKmk0bJxLHm3FBxO2YTgmUhCMER1Haz5n6W4sczCI83FqYEJnAs6Hi8XDyDPYPKW/kOEwo/IAjBshCE4Ai21rK5oTy3XiKi47rC5CAE4RUkByYc59o7mugr9BhZJbZkAgtCEIIj2MF7iVETMxXpT2KA8IpSgjgHoUA0JwTvjoJFIQjB7pkY7aiT5ofxDMKylopgzyBvVy+ONR1SSlBCYVOJxHi+3jI/DL2jYFEIQrB7h/UszEMIcec8QJiMAcJB8HbxCnD3L2+t5Fhzfpi4u14awKujYCkIQrB7W2ukheGct4/P1xWmYIBwcFJ4DxP6q2mMj3CgEY1CsBAEIdi972vZAq4rq0lMKtSfTMQA4eBwHyYkogVhwtZaNAnBQhCEYN+a+qi8nU0N4tkiLG+t9HfzHa724VjTgSUHJhToi/kOEy4IF7fWoEUIFoIgBPu2rVaaGSKquN7IedoTqZpEnhUd2nC1j5/bcL7DhBMDhLoeVteNLARLQBCCfdtayxZyfV+UiI41FqQGJfGt6djSNEnHGgs4FlQINCdU3IZ3R8EiEIRgxyRG23lPnDAxU6G+BJtODElqUNIxLc8gJPMwIXpHwSIQhGDHcpuYv1oY4cl1idHmU8GeGh9Xb441HV5KUEKBvpjvoqMLwsQd9diwHiwBQQh27Ltq6foRnPtF87Qn0oIwQDg03i5ewZ6aspZTHGsGudEobyELkyhAfghCsGPfVbNF4Zzv4bxGvClzNdKCEvO0J/jWXBQubq5GkxBkhyAEe1Xfw6q72JRAni1Cg2Q82VyeFIABwiFL1STmNXIOwutHCN9Vo0UIshtCEGZnZ8+bNy8lJeXpp5/u7z9/D2mdTvfss89OmzZtwoQJjzzySGNjI9frBDjfd9Vsfpio5Pq3XHFTabh3KPYgvArJgQklzWUGk4FjzTR/octI5e3IQpDXYH+LNDU1LVq06NZbb/38888PHjz44osvnndAWVlZX1/fqlWrPvzww7q6umXLlnG+UoCf2lzDFvEeIDzeWJCKAcKr4qFyj/AOL24u41hTILouXNiMd0dBZoMNwrVr106YMOGee+6Ji4t79dVXP/zwQ4PhJ3/6ZWRkvPXWWzNnzkxJSVm9evXBgwc7OztluGAAIqJ+E+1tkOZzXVmNiI5q89M0mEF4ldI0SUe1+XxrLgoXNtdgmBDkNdjfIwUFBRMmTDB/PH78+JaWlrq6uksdnJeXFxIS4uWFLWxALrvqWbKv4OfKs2afsa+8tRJ7EF61cZrkPN5BOCdUzNGx9gG+VQF+QjnI43Q6XUxMjPljlUrl6enZ2NgYGRl54ZH19fWPPvroO++8c5lqO3fujIiIOPcrU6ZM+eijjy56cHd3tyBw7gGDoeru7maM2c4P4ptK5Zwg6urq41jzSOOxUT5Rxj5jF3VxLMtRb2+vi4uLQqGw9oVcXJRb+KnW07o2vbvSjWPZyf6qTRU9S8NtpV3Y398viqJKpbL2hTg1g8FgMpmMRuMVj3R3dxfFKzT5BhuEPj4+3d3d5o8lSerp6Rk2bNiFh+l0ujlz5jz66KPLly+/TLUpU6a89dZb537Fzc3N09Pzogczxi71T2BJHh4eNhKEjOj7euP26xSenmqOZYtLyyaGptnyzaZQKGw5CIkozn/sqe7TU0Mncqz5s2hpW6NyRaytfNcqlQpBaHXmIFSr+fwGGGwQRkVFlZX9OAxeUVGhUCjCwsLOO6a5uXnu3Lk333zzM888c/lq7u7uUVFRQ71WALOjTcxTRWN9OKdyrvb4ryf9D9+azmacJjm3IZ9vEC6JEH6bazJICr5LqwOcNdg764477tiyZUtFRQURvfPOO0uXLvXw8CCitWvX7tixg4haWlrmzp2bnp7+2GOPtba2tra2mkw811sCOGvjGemGCM4p2NrXpu9pjvEbxbessxmvSTmqPc63Zoi7MMZH2KfFu6Mgl8EGYWxs7AsvvDBu3LiwsLBDhw698cYb5q9v2rTp4MGDRLRv376qqqp169aN/H9VVVUyXTQ4uW+q2A0RnFsHudrjKUGJooBGxzUZ7Tuyta9d39PEt+wNEeLGM7YyRgiOZ7Bdo0T05JNP/vKXv+zs7AwICDj7xS+//NL8wdKlS5cuXcr56gAuUNHBmvrYxADOLcKj2hPjMHHimomCkBqUmKc9MT96FseySyOEBVuld6aQTYxRg8MZ2t+/arX63BQEsLxvzrAbIkSR92/EPO2JcZoUzkWdUpom6SjvtdZihwluCjrejN5RkAU6gsDObDwjce8Xre6oI2IjvEP5lnVOE4JTcxuOMeIcWksiBPSOgkwQhGBPGnupsJXNCuHcHjxcnzcxOI1vTacV4qlxU7pVtp3hW3ZphLj+NFqEIAsEIdiT/1RJ14WLrrxnlB1uyJsYgiDkZmJI6uH6PL41pwQJ7QNU0oYsBP4QhGBPvj4t3RjJuTk4YBoo1JeM0yTzLevMJgSnHWk4xremQLQ0UkCjEOSAIAS70dRHeU1sAe+FtvN1RSOHR3mo3PmWdWZpQYklzWW9Rp4L4BHRTVHi+ioMEwJ/CEKwG/+pkuaHiW5DmPIzKIcbjk0MTuVc1LmplepYvzHHGgv4ls0IEnS9dKoDjULgDEEIdmP9aemmKP4TyY7gTRkZTAhOPdLAeZhQFGhppPA1ekeBNwQh2IfWfsrRs4XhnO9YfU9Ta1/7aN+RfMvCxJC0w/WchwmJ6MZIcf1p9I4CZwhCsA8bzkhzQkUP3v2iOfV544NTRNvYVcORRA+L6DP21Xdp+ZadHixUd7PTnWgUAk8IQrAPX1RIt0Xzj6tD9bmTQ8dxLwsCCZNCxh2sO8K3rEKgm6LEdZUIQuAJQQh2QNdLR5vYdbz7RQ0mw7HGgknBCEJZTAkdf7Aul3vZ26LFT8vROwo8IQjBDqyrlBaP4P++6LHGgiifCG9XL851gYiIJgSnFjeV9hh6+ZZN1wg9JipsRaMQuEEQgh34okK6bST/e/Vg/ZEpoeO5lwUztVId5z+W+/aEAtHNUcIXFWgUAjcIQrB1Z7pYZSebw3t9USI6VHd0augE7mXhLLl6R0eK6yoYmoTAC4IQbN3nFeymKFHJ+1Y93V5tYqaoYRGc68I5poZOPFB3WOKdWSl+glpBh3WIQuADQQi27osK6bZo/jfqobrcqaETuZeFcwV7Bnm7epe3VHCvfOtI8TP0jgInCEKwaXlNrNNA6Ro5+kUxccISpoZOOMB7EgUR3TlK+KJCGkAUAg8IQrBp/yqX7h7NfTt6au/vONV2Oi0oiXdhON/U0AlZtTncy0Z6CXHDhO9rkITAAYIQbJdRoi8rpRWj+DcHs2sPTwhOdVG4cK8M50kIiGvpa+W+xAwR3TVG/Fc5hgmBAwQh2K7NNdIYH2GkN/8g3F9zKDNsMveycCFREGRqFN4cJe5pkPSc93oCZ4QgBNu1ppzdNYb/Ldpn7M/XFU4KwQChhWSGT86qOcS9rKeKrgsX1+GVGbhmCEKwUc39tKteuimK/y2aU3803j/G08WDe2W4qHFByRVtVS29rdwr3z1a/BeWW4NrhiAEG7W2XFo8QvRW8a+8v+ZQRvgk/nXhElQK1YTg1IP1/GfWzwoRmvroeDNGCuGaIAjBRn1UKt0fw//+NEqmnIaj6WEIQovKDJ+yv+Yg97KiQPeMFT8sRaMQrgmCEGxRlpaZGGXIMH3weGNBmFeIv5sv98pwGZNDxp3QFXcberhXvneMsK5C6jFyLwxOBEEItujDUumBGP7TB4lo15n9MyMyZCgMl+Ohck8JSsiuPcy9cqiHkB4kflmJRiFcPQQh2Jz2Afr2jLRilCz9olm1OdPDp3KvDFc0KyJz95ksOSrfHyOgdxSuBYIQbM7aU9LCcDFAzb9yrvbYCO+wII8A/qXhStLDJuXrCjsGOrlXvi5crO7CDoVw9RCEYFsY0d+KpQdjZbkzd53JmhWRKUdluCI3pXqcJlmO3lGFQPeNFf9WjEYhXCUEIdiWHXVMIdA0GV6TMZgMB2oPTxsxhXtlGKRZEZm7zuyXo/IvYsUvKqTWfjlqg+NDEIJtebdIejxBltsypyFv1PAovC9qRVNCJxQ3lbb1t3OvrHGj68LFTzC5Hq4KghBsSFUnO6iTbhspU78o3he1MrXSdWJw2n4Zllsjov+JF/9aLEkYKIShQxCCDXm3WLpnjOiu5F+529CTU390xoh0/qVhKOZGTf/h9B45Kk8OFHxd6ftaJCEMGYIQbEW3kdaUSw/J85rM3uoDaUFJPq7echSHwZsUMq62s16OXZmI6JE48S+FJjkqg2NDEIKt+PCkNDtEjPSSYxo9bavcNS9qphyVYUgUgmJWRMYPp3fLUfy2keLJdjqGpUdhiBCEYBOMEr1dKP1KntdkGrt1p9urJ2PfJdswL2rm1spdjPjHlUqkR+LE1QV4ZQaGBkEINuHfldIob5oUKEtzcGvl7tmRmSqFDDtZwNCN9R3lplQX6k/KUfzBWHFbrXSmC41CGAIEIdiE1QXSr5MUMhXfUbVnftQsmYrDVZgbNUOm3lEvFd07Vny7EI1CGAIEIVjfD3XMxGhemCzNwRO6IlEQY/xGy1Ecrs68yBl7zmT3GfvkKP5ovLimXGrG5HoYNAQhWN9Leaank2XZa4KIvi3ftnjUAnlqw1Xyd/dLDIzdJc8a3CHuwk1R4lsFeH0UBgtBCFa2vY7p+uiWaFluxY7+zkP1ufOiZ8hRHK7FktELvi3fKlPx51LE90ukJlkanOCAEIRgZS8fM72YJirkaQ9urdw5NWyit4uXLNXhGkwMHtfa11bWUiFH8RGewo1R4jtFaBTCoCAIwZq21zFtr1zNQSLaXLF9CfpFbZIoCItGzf3u1A8y1X8+RfygRGrBSCEMAoIQrGlVnozNweONBYwoISBGlupwza4fOW/Xmf09hl45io/wFJZFiqsxUgiDgCAEq9l4Ruo00K2yNQc3lG1ZNuY6mYrDtfN1G56mSdp2epdM9X+fKv69RKrtxpxCuAIEIViHidFzudJrExUyvS2q7dYdayxYED1blurAyfKYG74+uUlismRVqIdwz1jx5WOYUwhXgCAE6/jfUilATfPlmTtIRF+f/HbRyLluSrVM9YGLxIBYb1fPg3VHZKr/2xTFN2ek4jY0CuFyEIRgBT1GeumY9PokuZaS6TH0bju9e9nYRTLVB45uirnhy5MbZSo+zIV+naR47ggahXA5CEKwgj/nm6ZphPH+cjUHvzu1bWJwWqC7v0z1gaMZI6Y2dGlPNpfLVP+RODG/he2qR6MQLglBCJZW2cneL5FemyjXvWeUTOtLv1ses0Sm+sCXQlAsG7Po65ObZKqvVtCbk8X/OWAyoFkIl4AgBEt7/KD06yRFqIdczcEfTu8O9QrG4qJ2ZMnoBUcajtV01MlUf2mEGOlF7xUjCeHiEIRgUVtqWFk7e1yefQeJSGLS58Xr70q8Vab6IAcPlfuysdd9VrxevlO8NVnxynGTVpYpi2D3EIRgOd1GeuSA6d2pChfZ7rvtp/f4ufkmB8bLdQKQx/KYGw7UHq7trJep/hgf4f4Y8bGDmF8PF4EgBMv57RHTjGBhbqhcnaISkz4t+uquhFtkqg/y8VC5Lx2z8Ivi/8h3ihdSFYUt7D9V6CCF8yEIwUIO6dj6vBPm8QAAFYpJREFUKvaGbFMmiGhH1V4fV580TZJ8pwD53BSzZF/NwfourUz1XRX0z2mKRw9KrViAFH4KQQiW0Geie/aZ3pks+rrKdYoB08BH+Z/9InWlXCcAmXm7eN0cc8M/jq+R7xSTA4VlEcITOegghZ9AEIIl/DrHlOwn3Bgl4/32demmGL9RiQFx8p0C5HZz7NLiptJC/Un5TvHKBEV2I/v6NDpI4b8QhCC772vYt9Xsb1Nl7BTtHOj6smTjAyloDto3V4XLzxNv+1ve/zKSa/67p4q+mKn45QHTmS5MsYcfIQhBXtpeune/8bMZiuGydYoS0ccnvpgZkR7mFSLjOcAi5kfP6jP2ZdUcku8U4/yFXyUo7tprMiEKgYgQhCAro0S37jI+GKvI0Mj1pigRlbdW7jqz/+7E2+Q7BViMKIiPjLvv3aP/7DP2yXeW3ySJSoF+fxSDhUCEIARZPZljchXpuRQZbzOJsTcPv/+LlJU+rt7ynQUsKU2TlBwY//GJL+Q7hSjQl7OVX1RgsBCIEIQgn88rpO9r2b9nK2XagN7sm7ItSlG5YCT2HXQo/zPuvh+q9pS3Vsp3Cl9XWj9H8csDpqJW9JA6OwQhyCK7kf3qkOk/cxTDXGQ8i66n6V8F656a+LBAcoYtWJy3q9e9SXeszvmbicnYe5niJ7w2UbF0u0kvYy8s2AEEIfB3qoPdvNO0ZroyYbiM+SQx9srBt5fHLInwCZfvLGAti0bN9XTxWFv4paxnuWu0ePtI4bqtxm6jrOcBm4YgBM4aemje96Y/TRDl233e7N8lGwwmw+3xN8p6FrAWgYRnpjz2Tdn3RU0yTiskohfHKRJ8hdt2mYwYLnRWCELgSddLs7cYfxEj3jVa3lvrVGvluuINz6c/KQq4hx2Wv5vvkxMf+kP2mz0GGbeNEIj+kaFgxFbswYQKJ4VfIsBNUx/N+d5460jx6WR576uOgc4X9r/66Pj7NR6Bsp4IrC4zfMr44JQ/H3pHvin2RKQSaf0cZbeRoV3onBCEwEddN5u52bg0QnghVd6bSmLSH7LfnBY+ZXbkNFlPBDbisfG/aOtrl3uw0EWkr2Yr2wbYij2mAWShk0EQAgcVHWzGZtONUcJL42RcR83s/WOfMMawmprzUIqKFzN+s6l8W3btYVlPpFbQpnlKRrTge2OHQdZTgW1BEMK12tvAMjYZf5cqvpgmewp+WfJNTv3RFzKewtCgU/F1G/7ytGdfz3m3QF8i64lcFfT5TMUYH2HmZmNNNwYMnQV+m8A1+WuxdOsu46czlStlfjuGiLZW7vq6dNPrM3/v5eIp97nA1sT4jX4h/de/2/en0pZTsp5IIdAHGYrbRoqTN5qytMhCp4AghKvUaaCVe0z/OCllL1HODpF9PvsPp/d8mL/2zdkvB+EFGWeVpkl6cuLDz+55uaKtSu5zPZUofjxdcdNO4+oCCWHo8BCEcDVydCx1g9FNSQeXKKO9ZE/Br09u+ij/0zdnvYT9JZxcZviUx8b/4smdL5zQFcl9rnmhwqElyg1V0oLvjQ09cp8NrAlBCEPTbaSnckxLtxtfnyj+PUPhrpT3dBJjfz/2r29PbX1v3p+xggwQ0fQRU19If+qF/X/eV3NQ7nNFegl7FimnBokpGwz/LEXT0GEhCGEINlRJieuNul4quFG1LFL2m6ejv/PpPatKmsvem/vnQHd/uU8H9iJNk/T6rFXvHf3n34/9S2LyznVQivT7NHH7QuVHpdLMzcb8FqShA0IQwqAcbRFnbDa9mCd9lKlYM0Phr5b9jAX64ge2PjFyWOTq2S95u3rJfj6wK6OHR3+48M2y1oond73Q2K2X+3RJvkL2YuWt0eKC742/OCDUoafUsSAI4QoONLLrthlXZKnuHCXkLVPOkv+9mD5j319yP3wx6/VHx9//YOrdCkH2WRlgj3xcvV+f+eK4oOQHvn/i2/Ktsi49Q0SiQA/GiqU3q4LcaNy39FC2qaoTrUMHgSCEixuQ6LNT0uRvjSv3mpZFiPnX9987VpR1Z0EikhjbWrnrzk0Pdxm6P1n07tTQifKeD+ycKIgrEpa/M/dP31fufHjbrwv0xXKf0VtFL6WywqXk60rjvzHetNO0pwFxaPdkftUB7NARPVt7SlpXIaX4Cb9NFq8fIYoCdXXJe1KJSXurD3xa9JWbUr0q8+k4/7Hyng8cSKRP+N/mv7ajau/L2W+OGh55R/xN8f4xsp7RX01/HK94Nlmx9pT0y2yTidGdo8UVo4QIT+yLaZeGFoQGg0GlUl37MWBrBiTKbmTfVUv/qWIqke4YKebcoIySf14EEbX2tW8/vXtD2RY/N997k1dMDZ1ggZOCgxFImBs5Y3r41M0V2/+Q/aa/u98Noxdkhk9xVci4MbSnih6KFR+KFQ/p2Jpyafw3pkhP4cYocWG4kOQrIBLtiMDYoNr1XV1dd9111/bt2xUKxdNPP/3MM89ceMxrr732yiuvmEym2bNnr1mzxsvr4i84fP755999993nn38+yEvs7Oy8VCm4Fv0mymtm+7Vsv1bar2Wxw4SF4eKyCCHR9yKPcFdXl4eHh8Dv6W7v7zhQd2R/zcF8XVFG2KQloxfI/Ve8A+jt7XVxcVEoMGh6ORKT9tUc3Hxq+8nm8szwyZnhU8Zpklz4JWJ/f78oihf+uW9itLeBfXNG2lbLOg1sdoiYqREyNEKMjyAiFXkzGAwmk0mt5vPa3mCD8Le//e3Ro0c3bdpUX18/adKkDRs2TJ069dwDcnJyFi1alJOTEx4evmzZsoSEhFdfffWipRCE1lLTzcraqbCFFbayvGZ2so3FDhMyNMI0jTAjWPR1vdz/5RKErX3tJ5vL83WFxxsLqztqJwSnpodNygib5K5yu5ayzgNBOCT6nqY91QeyanPKWyri/MemBiUmBsSO9h3pprym356XCsJzne5kexrYPi3LbmSNPSzVX0jxExKGCwnDhTE+wuWfNRgM6wRhaGjoJ598MnfuXCJ6/PHHe3p6/vGPf5x7wEMPPaRQKN577z0i2r1796233trY2HjRUghC+Rglau6npj7W2Ev1Pay+h+q62elOqupilR3MS0Uxw4T44UKir5DiKyT7CepB/0YdahAOmAb0Pc3abl1Dl7a2s+F0W/Xp9jM9ht4Yv9EJATFpQUmx/mNVIoaohwZBeHW6BrrzdUXHGk8UN5VVtFVpPAIifUZE+ISHeQUHe2qCPQKHuw0b/MvJgwnCc7UNUK6e5bewwlZW1MrK2plSoGhvIcpLCPegcA8hxIOC3YQANwpQIyMHi28QDuo3UU9PT319fVxcnPnTuLi4devWnXfMqVOnli1bdvYAnU7X0dHh7e190YIGg6GlpeXcr6hUqkulXU9/f2tv92Cu06YYGXUbr/xHRreRTD+dENxpIPMKFhKjTiMjIpNEXUYiog4DEaMOAzNK1G2gbhN1G6l7gHUaqWOA2gZYt5GGu5CvWvB3pUB3IUhNI92EaWEU6iGEewieP31y+7uon4iIegx9Jma66OX1GHtNTCKivt5eQSUYmYmIugZ6GLFeY5+RmToN3X3G/l5jX6ehp8vQ0zHQ2d7f2drf3mca8FcPC3TzD/EIDPPSLImYHuEVFnLuMqF9vdj0bahYX59kVAkIwiFyJ5riGzvFN5aIjJLpTFdddUddVWf9oeocbY++saepfaDT28XTx8XL28XLx9XLU+XupfJwU7qqlWp3pdpV4eIiqlSiSq10ISJRElQKF6VSSUQuCpWrePFOV1eFi4tCRUQC0QRPmuBJNOLHf2ruo5puVtvF6nupponl1pC+l7X0k76P9RhpmIq8XAUvFXmryFMluCtJLZKPK4kCeSoFpUDuKiIiHxciIpHIS/Xjn6cqBbkrzrsGch3EzeKlsqfNXDzUaj/vYXxrDioIW1tbicjT88cl/729vc+LMfMxZw8wR1pra+ulgnDXrl0jR4489ysZGRlffPHFRQ/+Knvj5vZvB3OdjkogosvOkfIg8iAKVhCZb/p+MvVTQwc1DK7+/7V35zFNbGsAwE8rsrVlURar1UhBpFowopIoIYYXvWAEEpU8X3B5EdFEfYhQTDQoRsUQ3BUNBLQaSQwWUKIGqxi3EFEBBdwIWKwgWNCuvNqWLvP+mGfDLXit17ZTO9/vr+mZrzNfKJOvPefMHC8Thfqdw3uZ0Djs/6f3wChuJgxvpGLIE6OMMyGaEXlgyNOI/E0UmpHCMCK6AfkaKQwDQkiOkByhLvPRrMwHALvyRigCoZEj0iYKRemmHhr3X9W4frUbRT0OU1ORchw2QKVoKdgwFekpmIGKdBSEENJSkelbz8gwBQ1/p4YMUzG91SMJVIQCEAowX8J6hPRI/+0SGhOGSDrwyB6eVPSvY/gvQoPB8MN4b29vKvUHhd6qQhgQEEChUJRKpa+vL0JILpcHBwdbxAQGBiqVSnxboVDgLd87YEJCgvVdo//+xz//w9hoZTCwE5tPlgF/A3SN2o/1z7H92a5RYA+27Rq16gexh4cHm81+8eIF/rK1tTUiwnKCH4fDGRkwffp0b29vm6QIAAAA2I+1PcObNm0qKCjo7e199OiRQCDIyMhACEml0oSEhM+fPyOENm7cWFNTc//+/Y8fPx48eHDTpk12zBoAAACwEWun7fF4vC9fvsTFxfn5+ZWUlERGRiKEMAzT6XT4vNPZs2eXlZXl5OTI5fJVq1bt3LnTjlkDAAAANmLt7RM2BLdP/I5gjNAZwBihM4AxQmdAwBghsXJzc9+/f090FmRXUFDw/PlzorMgu+Li4nv37hGdBdlVVFRcvXqV6CzI7saNG3w+31ZH+w0KYWNjo3k+KiBKS0vL4OAg0VmQ3cuXL3t7e4nOguw6OjpEIhHRWZBdd3f327dvbXW036AQAgAAAPYDhRAAAACpETBZ5unTp+vWrdNoNFbGq1QqGo0GEwSINTQ05OnpCRMEiKVWq8ePH+/ubselhcAPaTQaCoViq2ka4O/RarUmk8mau9Vv375tfj7o9xBQCBFCEolEr9c7/rwAAABIZdKkST/8Bk9MIQQAAACcBIwRAgAAIDUohAAAAEgNCiEAAABSg0IIAACA1JyoEGIYtmvXrgkTJvj5+W3fvt1oHGPN9Js3b7LZbBqNFh8fD4/YsIfKysr58+d7enoGBgZu3bpVq9WOjuFwOKHf7Nq1y/FJury6urrQEZ48eTI6prGxkcvlenl5zZ8//9WrV45P0uWdO3cu9M/wlVZH+uOPP8x709LSCMnT9ajV6h07dsTFxYWGhorFYnP78PBwRkaGr69vQEBAQUHBmO/l8/lMJpPBYKxYsWL05/VdmNMQCARsNru/v18qlUZGRpaVlVkEyGQyBoNx69YtvV6fmZmZlJRESJ6u7cKFC/fv39doNB8+fIiKisrPzx8d4+Hh0djYKBKJRCLR4OCg45N0eVeuXImNjRV98/XrV4sAvV7PYrHOnTtnNBqLioqioqIIydO1KRQK80eQm5sbGxs7OobD4QgEAjymr6/P8Um6JKVSuW/fvtraWoRQZ2enuf3IkSMxMTFKpVIsFk+ZMkUoFFq8saOjg8FgtLS0aDSalStXbt261cozOlEhXLZs2eHDh/Ht8vLyhQsXWgSUlJSY/xclEombm9unT58cmiLJFBQUJCcnj2738PDo7+93fD7kceXKlSVLlvxFQF1d3bRp00wmE4ZhOp3Ox8enqanJUdmR0axZs/h8/uh2DofT0NDg+HzIwGAwWBTCWbNmVVZW4tt79uxZvXq1xVt2796dlpaGbzc1NTEYjOHhYWvO5URdo52dnVwuF9/mcrldXV0WAV1dXfg6iAih4ODgCRMmwKNv7cdkMgmFwoULF465d9GiRSEhIWvWrIEOajt59uzZ5MmT58yZU1RUNHqYoKuri8vl4qtiubu7h4eHd3Z2EpEmKTx58qSnpyc1NXXMvWlpaSwWKyUlBTqo7QrDMJFI9MMaMTJgaGhIIpFYc3AnKoRyuZxOp+PbDAZDJpNhf77ZXy6X02g080sfHx+ZTObQFMnkwIEDKpUqOzt79C6BQPD48eP6+nqTyZSUlIR/cQM2NHfu3Pr6+ra2tpMnT5aWlh45csQiAK4FRzp//vzq1avHXBX1+PHjDQ0Njx8/DgsLW7p06U8MSoGfpFardTrdyBohlUotYkYWEfyRkFZeF05UCAMCAszLLSkUisDAQItlYCdOnKhSqcwv8RiHpkgaJ06cuHz5slAoHPOBiikpKUwmMywsjM/nv337Fn6L2NyMGTNiYmICAwPj4+Pz8/OrqqosAuBacBi1Wi0QCNLT08fcm5iYOHXq1GnTph07dszd3b2hocHB6ZEHnU738vIaWSOCgoIsYiZOnGgOUKvVer3eyuvCiQphREREe3s7vt3e3j5z5sy/COjr61MqlaGhoQ5NkRzKy8tPnTp1584dJpNpTTwGT+mzJ4uvg7iIiIiXL1/if3mtVtvZ2Tn6egE2UVVVxWQyvzdGMBKFAk+stK/w8HDra0R7e7u/v//oYjm2XxzPtKHr16+zWKyOjg6xWBweHn7p0iW8fe3atfhEAKVS6efnd/nyZaVSuWHDhtTUVELzdU3nz5+n0+nV1dXNzc3Nzc1v3rzB2/Py8mpqajAMa2trEwqFEomku7t7zZo1XC5Xr9cTmrILqqmpaW1tlUqlDx8+ZLPZhw4dwtuzs7Pr6uowDDMajWw2++jRo0NDQ3l5eQsWLCA0X1cWFxdnnsSHO3369NmzZzEM6+npqa6u7uvr6+3t3blzZ3BwMD6gA35da2vrs2fPEELXrl1rbm42GAwYhhUXF0dFRfX19bW1tQUFBT148ADDMK1Wm5KSIhaLMQzr7u6m0+l3796VSqWJiYk5OTlWns7NZsX6lyUnJ3d0dCQmJhqNxvT09LVr1+Ltg4ODOp0OIeTj41NbW8vj8bKysuLi4kpLSwnN1zW9fv165syZhYWF+EsOh1NRUYEQkslkarUaIaTT6fbv39/d3e3t7b1o0aKbN2+6uTnRf5FrePfu3d69eyUSCYvF2rx5M4/Hw9tlMhm+fhmVSq2trd22bVthYWF0dHRlZSWh+bqsgYEBjUazfv36kY0qlQpfFc5kMp05cyYzM9PNzS06Orq+vt7f35+gTF0Nj8dTKBTz5s3D7xd8+PAhjUbbsmVLT09PTEyMp6fnvn37Fi9ejBDCMGxgYACfqRASEnLx4sWsrKwvX74sX7784MGDVp4OfssDAAAgNScaIwQAAAAcDwohAAAAUoNCCAAAgNSgEAIAACA1KIQAAABIDQohAAAAUoNCCAAAgNSgEAIAACA1KIQAAABIDQohAAAAUoNCCAAAgNSgEAIAACC1/wECiFt0bsnxEgAAAABJRU5ErkJggg==",
"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.8411086231969123e-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.00022342725151536436"
},
"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.6.0"
},
"kernelspec": {
"name": "julia-1.6",
"display_name": "Julia 1.6.0",
"language": "julia"
}
},
"nbformat": 4
}