{
"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 pick a harmonic\n",
"potential. We use the function `ExternalFromReal` which uses\n",
"Cartesian coordinates ( see Lattices and lattice vectors)."
],
"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": [
"n Energy log10(ΔE) log10(Δρ) Δtime\n",
"--- --------------- --------- --------- ------\n",
" 1 +159.8469026645 -1.57 320ms\n",
" 2 +122.5658803513 1.57 -0.90 1.58ms\n",
" 3 +68.42271237955 1.73 -0.65 1.56ms\n",
" 4 +12.04805716751 1.75 -0.45 1.53ms\n",
" 5 +9.191318452175 0.46 -0.31 1.09ms\n",
" 6 +7.832911040811 0.13 -0.60 835μs\n",
" 7 +4.887964593339 0.47 -0.33 850μs\n",
" 8 +2.748013272908 0.33 -0.14 822μs\n",
" 9 +2.671390247657 -1.12 -0.18 701μs\n",
" 10 +2.370654782410 -0.52 -0.30 664μs\n",
" 11 +1.451754614743 -0.04 -0.15 828μs\n",
" 12 +1.297595520654 -0.81 -0.66 657μs\n",
" 13 +1.165546209566 -0.88 -0.96 660μs\n",
" 14 +1.153364585952 -1.91 -1.49 722μs\n",
" 15 +1.148183746230 -2.29 -1.36 910μs\n",
" 16 +1.145903970534 -2.64 -1.75 657μs\n",
" 17 +1.144903477262 -3.00 -2.32 660μs\n",
" 18 +1.144441760142 -3.34 -2.53 647μs\n",
" 19 +1.144122953493 -3.50 -2.55 475μs\n",
" 20 +1.144077686419 -4.34 -2.93 671μs\n",
" 21 +1.144049380286 -4.55 -2.79 697μs\n",
" 22 +1.144045255958 -5.38 -2.48 672μs\n",
" 23 +1.144041266372 -5.40 -2.87 839μs\n",
" 24 +1.144038883350 -5.62 -2.77 655μs\n",
" 25 +1.144038062934 -6.09 -3.16 656μs\n",
" 26 +1.144037435939 -6.20 -4.01 646μs\n",
" 27 +1.144037101964 -6.48 -3.90 673μs\n",
" 28 +1.144036920624 -6.74 -4.11 654μs\n",
" 29 +1.144036880220 -7.39 -3.85 644μs\n",
" 30 +1.144036861308 -7.72 -4.06 646μs\n",
" 31 +1.144036856415 -8.31 -4.68 814μs\n",
" 32 +1.144036854693 -8.76 -4.65 649μs\n",
" 33 +1.144036854153 -9.27 -4.67 674μs\n",
" 34 +1.144036853308 -9.07 -4.93 650μs\n",
" 35 +1.144036853094 -9.67 -5.32 644μs\n",
" 36 +1.144036852928 -9.78 -4.93 651μs\n",
" 37 +1.144036852849 -10.11 -5.31 650μs\n",
" 38 +1.144036852793 -10.25 -5.97 652μs\n",
" 39 +1.144036852766 -10.57 -5.99 868μs\n",
" 40 +1.144036852761 -11.26 -6.51 663μs\n",
" 41 +1.144036852758 -11.50 -6.35 649μs\n",
" 42 +1.144036852756 -11.92 -6.52 656μs\n",
" 43 +1.144036852756 -12.25 -6.36 654μs\n",
" 44 +1.144036852756 -12.90 -6.84 650μs\n",
" 45 +1.144036852756 -12.68 -7.13 680μs\n",
" 46 +1.144036852755 -13.14 -7.20 835μs\n",
" 47 +1.144036852755 -13.57 -6.87 650μs\n",
" 48 +1.144036852755 -13.46 -7.36 652μs\n",
" 49 +1.144036852755 -13.69 -7.42 648μs\n",
" 50 +1.144036852755 -14.45 -7.79 678μs\n",
" 51 +1.144036852755 -14.51 -8.11 658μs\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": "8.639717975914011e-16"
},
"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+gvaeTAAAgAElEQVR4nOzdd3wUx/k/8Gf37qRTR/3UJRAghDoSTaKIbmNjCLZjG7cYO4l74thxIW5xfrZjx04cl9jBSYhxb9hgekcSHYEKEpJQQ+106r3c3c7vj/NXwVQJZq9+3n9JYnl29dKuPpqZnRmBMUYAAACOSrT0BQAAAFgSghAAABwaghAAABwaghAAABwaghAAABwaghAAABwaghAAABwaghAAABwaghAAABwaghAAAByaBYIwPz9/zZo1wz9er9fLdi0wXAaDwdKXAGQwGLAmosUZjUZJkix9FY5OkiSj0cirmgWCsLCwcPv27cM/vr+/X76LgWHq7+/Hr2CL0+v1+BVscQaDgeOvYLgyRqORYxsJXaMAAODQEIQAAODQEIQAAODQEIQAAODQEIQAAODQEIQAAODQEIQAAODQEIQAAODQlJa+AAC4DEaU38pKO1hpq+jmxKK9WLwPRbgLlr4uADuBIASwXq0D9GGJ9O8SiRHFeQuhLmTop+31xmPNLGaUcO948ZbRohLdOgBXB0EIYKW+rJB+c9C4KFRcM0sxNUAgor6+QScnpUKh0Ev0wxnp7ZPSmwXSv2Yqkn3ROrRJer3+nnvuGRgYsPSF2IYnnngiLS1NjsoIQgCr02Ogu/YaS9rZunnKKQEXCDmVSMsixWWR4kdl0jVbDI9MVDyThIah7enq6vruu+/+/e9/W/pCbMD777+fn5+PIARwCC0DtHirYeIo4dgypdPl0u3OseKCUPG6rYYz3ezddIUCLUNb4+TkdNNNN1n6KmzA1q1b5SuOvyIBrEh9L8vYYJgTJHw4U3HZFDTRuNDuxcryLvbzXUYDtqYAGDkEIYC16NbTdVuNd0SLL6eNrGnnoaKNC5X9BvbgfmwPBDBiCEIAqyAxun2PMcVPuLLRPieRvpyrPNHCXstHqxBgZBCEAFbh8UPGPgN7P11xxRVclfTtPMW7RdL31chCgBFAEAJY3uYa9m0V+2Ku8ionBYa4CV/OUfwq21jXwzhdGoD9QxACWFhzP92XbfzPLMUoJw7VpgQID8QqfrHPiCQEGCYEIYCFrcwy3h4tZAZxm/rwTKLYY6C/FaKDFMzq2Weffe2110wf6/X66dOnNzY2XuL4xYsXFxQUmOXSLgNBCGBJX1RI1V3sj5OufGjwfEqR1s5WvHzCWNWFZiGYT3d3d09Pj+njDz/8MCEhITAw8BLH/+pXv1q1apVZLu0yEIQAFtOtp8cPSX+fPtwpg8M32kP4TZzitwfRKIQRW79+fUVFxUcfffTUU08ZjUaj0fjll18+//zza9asGRwcNB1TXFz81ltvrVq16uwvnu2999678847iai7u9u0dE5xcfHOnTuJaOvWraWlpUR07bXXHj58uKqqymzf2sUgCAEs5o/HjfNChJkaWdaDeSJBLG5nG2vQKISRef3116+77roTJ074+/sbjcbFixd/8cUX4eHhO3fuXLRoEWOMiD788MOBgYHo6Oj169ffcMMN51SoqamprKycPHkyEbW1tT300ENElJWVtXr1aiJ699139+/fT0RKpXL69Olbtmwx93d4HiyxBmAZxe1sTalUsFwlU30nkd6apnhov3FusFLNs+cV5HXHHqO2z3x/vvx8tHjv+HNbRNddd51ptO/bb79tbW3dvHmzIAgrV65MSUnZu3fv7Nmz33jjDdORd911V1hYWFVVVWRk5NB/z8vLGz16tFJ5+XwZP358fn4+t2/mSiEIASzjycPSqiRFoIuMp1gYKsR5C+8VSY/Fo+/HZjwaJ7abcTuK8aMu8MWhta1zc3Pr6+sXLFhg+rS+vr6kpGT27Nlr16595ZVXJElyc3Nra2urra09Owi7u7vd3NyGc3Y3N7fa2tqr/BauHoIQwAION7GCNvbVXNlbaq9OFmf/YLg3RvSUq+UJnKX6WX7pdGdnZ9MHLi4uM2fOfPfdd4f+ydXVtbGx8aGHHsrPz4+IiCCiyMhIo/Ena/v5+fm1tbWZPlar1Xq9/uwD+vr6XFx+/AOwtbXV399f1u9lOPB3IoAFPH3E+Hyy6Cx/j+V4L2FhqIipFHBlFixYsH379u7ubm9vb29vbycnJ6PR2Nzc7OzsHBQURES7d++urq4+53+lpaVVVlZ2d3cTkb+/f2Bg4NDeEVqt9tixYwkJCaZP8/LyTEOJloUWIYC5batjtT10e7SZ/gx9cZI4aZ3h1zFigJzdsGCX0tLSXnjhhcmTJ6ekpAwODpaUlOzatWvChAkpKSmpqalhYWH9/f1jx4495395eXnNnj17x44dS5cuJaI1a9asXLlSFMXu7u6UlJRnnnlmwoQJRNTZ2ZmXl7dw4UILfGM/hSAEMLdnjxr/X6p4laupDV+Eu3DLGPEvBcbXJuOdGbi8TZs2qdXqoU8ffPDBu+6669SpU2q1euzYsaZe002bNpnmwsfHx3d3d7u6uhLRSy+9JAg/9us+9thjb731likI582bV1FR8corr+Tk5Hz77bdD/aJr165dsWKFt7e3mb/B8yEIAcxqRx3rMdDyKLOOSjyVKCZ9a3g6UeHtbM7Tgk3y8PA45yvu7u6pqalnf0UUxcTERNPHnp6eQ4cNHbBw4cKdO3fqdLqAgAAiUigUAQEBXl5eQylIRJWVlc8++6wc38JIIQgBzOrVPOOTiaKZX4cIdROWRIjvFUurrmiPJ4ArMLTcmklSUtLZSUlEf/nLX8x7RReFpwLAfE60sNIO+vloCzx3TyaK75w09hnMf2YAIqLJkyffdtttlr6KC0MQApjP/zsh/S5e5L6g2nCM9xKmBIhryvD6KFzGvHnzNm/eLEdlxtjChQuLi4svccxDDz20fv16Oc5+CQhCADM53cn2Nkgrz1vFw2yeTBTfKJCwPxNc2m9/+9u4uDg5Km/YsMHJycn0yujFPPDAA08//bRpITezwRghgJm8UyTdFyO6W25i+7QAIUBNP5yRbojAX8BwUb29vQaDgYiysrL0en1PT8/OnTsTExPvvvvuY8eOffbZZ2FhYffff7/p9dGcnJydO3e2tbUlJibefvvtQ8uqffvttzk5ObGxscnJyTU1Nab1SN9///27777bdMBrr7324IMP6nS67OzsO+64Y/v27S4uLhkZGbGxsa6urnv27MnMzDTbt4znAcAcuvW0tkz6ZYyFn7gHY8V3i9A7Cpfy1ltvnTx5koi2bdu2cuXKzZs3JyQk/OlPf7rvvvteeeWVhISE77777vHHHzcd/Nlnn2k0mrS0tK+++moo5F5//fVnn302MTGxurp6+fLln376KRH19/fv3r171qxZpmOeeeaZ7u7usrIy07I169ev37Fjh+mfZs2aJVPf7MWgRQhgDh+flmYHiRHuFl4966bR4uOHjKfaWcwoy6/jBRfU8q8/GtubzHY6t6mL3NIXX+xfx4wZ89577xHR4ODgiy++WF1dberbvOWWW95++20ieuedd4ior68vMzMzKirqww8/dHZ2fvnll3fu3JmSkkJEVVVVAwMDRFRWVubs7HzpHQpNxo4d+9133/H6BocDQQhgDv8olv461fLz2Z1EujdG/Eex9NY0y18MXJDX9fdIA31mO53S+1JLfQ4NFgYEBIwdO9bJycn0cUtLi+nrL7/88urVq728vJRKpdForK+vd3Nz6+zsHJplmJKScuDAASLq6ek5exLhJbi5uZmWZzMbBCGA7PY2ML1EmcFW0Qi7f4IY/43hpVQFluG2TsqAUEtfwv8oFIoLfmxy7NixDz74oKCgwNPT02AwuLq6SpLk5eVFRB0dHT4+PkQ0tPp2QEBAW1ubJEmiKBKRi4tLX9//8r63tzcsLMz0cXNz83AajhxhjBBAdu8VSw/GmnsS/cUEuwqZweKnpzFSCFerpaXFzc3NNE1+zZo1er2eiNRq9dy5c02T5RsbG9euXWs6OCoqysvLq6SkxPRpSkrKunXrTB+3t7cPdaUSUV5e3pQpU8z5jaBFCCCvlgHaVit9kGFF7a9fxoh/OGr89QT8HQxXZdasWT4+PomJiT4+PiEhIUOrhq5evXrFihVBQUGBgYELFy40tfwEQbjxxhs3btxomj7x3nvvrVixorm5ub29PTY2duXKlfPmzSMiSZJ27tz53HPPmfMbQRACyGttmXRduDjKydLXcZZ5wcIv+yivlSX6WEkzFazIjh07VCoVET3zzDNDX7z++uvnz59v+jg8PNy09ZKzs/O+ffvKysrUanVERER7e7tp3dGIiIjs7GxTL+jdd989btw403989NFHly9f/tvf/lahUEycOPHEiRNffPHFyy+/fOjQoaFlvtevXz9lypSoqChzfssIQgB5rSmV/mplb6aIAt01VlhTahXv74C1Gcqks99tUalUpnQkIlEUR40aNfTx+PHjTR8PfXHr1q3bt2+Pjo4+fvz4jh07hhYdHTdu3AMPPFBWVhYTE2P6ire3t4uLy9mbXdTU1Lz66qtyfW8Xgb4RABkdbWadepqlsbqG1y/GiZ+clgaMlz8SYKRSU1MnTpzY2to6derUgoIC0wYUJvfff/9QChJRVFTUHXfccfb/ffjhh8eMGWO+ayUitAgBZPWfUukX40RreU/mLJEeQryPsOGMdKN5N4QCR+Dr6/uLX/xiOEeOHTv2/H19zQ/PAIBc+o30Rbl011jri0EiIlo5Xvx3Kd4dBUAQAshmwxkp2U8It/RqMhezLEI8qGNa803dBrBSCEIAuXx6mq0YY72PmIuSrg8XvyhHoxDOVV5e/thjjz388MObNm2y9LWYg/U+pQA2rW2AdjdIyyKt+hFbES1+iiCEnyovL585c2ZkZOT06dMfeOCBjz/+2NJXJDu8LAMgi68rpQUhopc1TR8839xg4a5uVtrBxnlZaf+tA3p6z0stfW1mO92CqMwbY64/+ytvv/327bff/sgjjxCRq6vrqlWrbr/9drNdj0UgCAFk8Um59Js4q24OEpFCoJtHi5+Vs+dTEITW4uHU+7oGzbfkdKhH8DlfKSoqWrlypenjSZMmlZaWDi0Qaq8QhAD81feywlZ2TagN/O5YES2u2G18PsUGLtVBBLtrLHsBfX19XV1dpo87Ozvd3d3tOwUJY4QAcvi0nC2LFJ1tYdmWyf6CINCRJmbpCwEr8tFHHw0ODhLR6tWrFyxYYOnLkR1ahAD8fVEu/XmyLcQgERH9fLTwVaWU5m8zFwxyGzNmTHJyskqlMhqNGzdutPTlyA4tQgDOqrpYdTebaX3Lql3MjVHiV5UMTUIYsmzZsoMHD65fv76goCA8PNzSlyM7BCEAZ19VsmWRotJ2nq1EH8FZpNxmRCH8j4eHhyNEoIntPKwANuLrSttbwPNnkcLXlZhQCEREv/71r2NjYy19FWZlY48rgJWr6WEVXSwzyGb6RU1ujBK/rECLEIiIVqxYER0dbemrMCsEIQBPX1WwpRG21C9qkuInCAKdaEEWgiOytecVwLp9UyUtt7V+UZPlkcJX6B0Fh2STTyyAdWropVPtbG6wjfWLmiyPEr+rQosQHBGCEICb9Weka8NElW0+VWn+QqeeSjuQheBwMKEegJvvq6V7xtlmDBIJRNeHC99XsycSbLJFa4uUSmVnZ+eYMWMsfSE2oKmpadasWTIVRxAC8NGtp/2N7Is5thqERHRDhPjSceMTCTb8LdgWT0/PqqqqgYEBS1+IbZBvXiOCEICPTTVSeqDgobL0dVyFOcHCbbuZto80Lpa+FIcREhJi6UsAjBECcPJ9NbshwrYfKJVIC0LFjWfw7ig4Ftt+bgGshF6irbXS9eE2/0DdEC58X433ZcCx2PxzC2AN9mnZWC8hyNXS13HVrgkT92mlHoOlrwPAjBCEABxsOGMPzUEi8nKiVD9hZx16R8GB2MOjC2Bxm2rY4nA7mXWwOFzcWIPeUXAgCEKAq1XSwXr0lOBjJ0G4JFzYcEZCEoLjGMH0iYqKim3btnl7ey9ZssTF5cKvV+fm5h48eNDDwyMzMzM0NJTTRQJYtY1n2HXhgp3EINEYT8FdJeS1sCRfu/meAC5luC3C7OzsSZMmnTx58l//+teMGTMuOAP08ccfX7JkSX5+/vbt2//5z39yvU4A67WxRlocZleZsThMQO8oOI7htghffPHFVatWPf7440ajMTU19euvv16xYsXZB2zduvXjjz8uLCz08/OT4ToBrFSnno40sTnBdjXKsDhMfO6YcVWSXX1TABczrBt9YGBg586dS5cuJSKFQnH99ddv3rz5nGO++uqrO++8s6WlZf369VVVVdwvFMA6ba+V0gMFd1teUOZ8M4OEonam67P0dQCYxbBahFqtljEWHBxs+jQ4ODg7O/ucYyoqKgoLC/ft2xcTE3PPPfe8+eabd9555wWrtbW1nTx58uWXXx76yq233nqJAUW9Xq/X64dznSAf009BsJ+BMG42nqEFwWSeW9T0I5Ak2ec2CESZGtp0Rr9itNynsj16vV4U0Va2ML1ebzQaFQrFZY9UKpWX/cU1rCBkjBHRUC1RFM9/FAcHBwcHB48ePSqK4qZNm2677bYVK1Zc8CoNBsPg4GBbW9vQV/r7+y/xbEuSZIYnHy7N9FNAEJ6DEW2tE5+YaKY71Jw/hUUhwtZa4dZIPHrnMv2w8UvJsqT/w6XasIJQo9EIgtDY2BgZGUlEDQ0NQ63DIcHBwb6+vqY/lNLT0zs6OrRa7QXXk/X3909OTn799deHeYmDg4POzs7DPBhkotfrnZ2dEYTnyG9lbirjBD8z3Z+SJDk5OQ3nr+CrtziSPXvcoHJyFvEzP48oiiqVffWG2xpRFI1GI69oGFYDX61Wp6enb9q0iYgYY1u2bJk7dy4RGQyGM2fOmDJ5wYIFxcXFpuOLiopcXFwCAgK4XCKA1dpSyxaF2mdQhLoJAS5CbgveHQX7N9y3Rv/whz/ceuutOp2uuLi4vb39lltuIaLKyspx48Y1NTX5+fnddtttb7755h133BEfH//++++/8MIL+IsJ7N7WWumxeHO0zyxiUaiwpYal+tln0gMMGe6Q78KFC3fv3q1SqTIzMw8cOODm5kZEGo3mk08+8fDwICJXV9dDhw7NmDFDkqS1a9f+/ve/l/GqAaxAj4GONrFZGrvNiYWh4lYsOgoOYAQryyQmJiYmJp79FQ8Pj9tuu+3sT3/5y19yuzQA67a7nqX529vEibPN1AiFrax9kEY5WfpSAOSEl4ABrtDWWmlhqD0/Qc4KmhYo7KpHoxDsnD0/xgCysuM3ZYYsDBG31uJ9GbBzCEKAK1HeyXoNFGcvO05czKIwYQuCEOwdghDgSmyrY/ND7H9a5XgvQRSopANZCPYMQQhwJXbUsfkhdp+DRERzgoTtdQhCsGcIQoARMzLa0yDZ2Y4TFzMvRNiBIAS75hBPMgBfR5tYiJsQ5Grp6zCL+SHingZJj1dHwX4hCAFGbLvD9IsSkZ+aojyEI01oFILdQhACjNiOeml+iAM9O/NDMEwI9syBHmYALnoMlNvMZtjvymrnmxci7sC0erBfCEKAkdnbwCb5CW4jWJ3Q5s3UCHktrBPbY4OdQhACjMyOOmmeI/WLEpFaQVMChL0NaBSCfXKkP2sBeNhVzz7IMGsQDhoHd1btK2opLWk+7erkMt4nOikwflpIqjmvYU6wuLueXR9uznMCmAmCEGAEWgaoqptNMuMWfVk1B97N/XeUV0RqUNLMoKlGUTrdXvVh3trPi759OPW+aO8o81zGnCDhl9loEYJ9QhACjMDOOmmmRlSapUEoMenPB/9e0lr++ykPp2gSiKivr8/JyWlqSOptsT/74fS2x3c9vzJxxfXRC81wMan+Qk0Pa+yjQBcznA3ArBCEACOwu4HNCTZHc1Bi0isH/tY+0PnBojecFefuBygK4pKxiyYHp/x2xx/6DQM3xSyR+3oUAmUEinsbpJtHO9b4KDgC3NMAI7Cz3hxBKDH2fNafe/S9L89cdX4KDtG4Bfxt3v/7rnTT16c2yH1JRDQnWNhVj9mEYIcQhADDVdfD2gZYnLfsQfh58bdt/R1/nPG0SqG69JGBbv5/nfenT4q+Lmgqkvuq5gQLuxoQhGCHEIQAw7Wjns0JFkWZc7CoufTL4u+fTX9MKSqGc3yAq99TUx/5Y/ZfOgY6Zb2weB+hc5BVdyMLwd4gCAGGa5f8/aKdg10vZL/25NRHAt0Chv+/pgRPmhMx45UDf5PvwohIIJodJO5BoxDsDoIQYLh217M5QfIG4YcnPk4PTbuCOYL3Jd3Z3Nu6qzpLjqsaMidY2IlFR8HuIAgBhuV0J2NEY71kDMLTbRX7avb/Iv62K/i/SlHxm7Rf/yP3P/2Gfu4XNmROsIAWIdgfBCHAsOyuZ5kyNwffPvrhvYl3eDp7XNl/j/OPSQyI+7ToG75XdbZoT0EiKu9EFoJdQRACDMteLZslZxDurNrXZ+i/dsz8qyny6+S7vivd3NDdyOuqzjc7CI1CsDcIQoBh2dPAMmV7U8bIjB/mffzQpJWicFWn8HP1XT7+ujUFn/G6sPMhCMH+IAgBLq+sg4lEoz3kCsIdlXs17gEJAROvvtSNMUsO1B2t62q4+lIXNDsI0+rB3iAIAS5P1uagxNhnRd/eGfdzLtXcVK43jL3m8+J1XKqdL9pTUIp0GsOEYEcQhACXt6dBxgHCPWeyXVQuyYHxvAreNGHJnuqcxh4dr4LnmKVB7yjYFQQhwOXt1cr1yigj9nHhV3fH38KxpqeTx+Lo+Z8Xf8ex5tlmYZgQ7AuCEOAySjuYSBQlzwDh0YYTgiBOCZ7Et+zNMTfsqNzbPdjDt6xJZrCwqx57E4L9QBACXMZeOQcI15Vu/Nn4xdzL+rh4TwmetLliJ/fKRDTaQ1CJQmkHGoVgJxCEAJexV8tmy9Mv2tjTVNh0am7EDDmKLxt/7brSjRKTJa5maYR9WgQh2AkEIcBlZGnZTI0sQfh92eaFo+eolWo5ik/0i3FXuR3TnpCj+MwgYR+GCcFeIAgBLqWiixkZjfHkH4R6o35zxc4l0Qu5Vx5yw7hr1pVulKMyXhwFe4IgBLiUvQ1sljzNwd1ncqJHRYV5hshR3GRuxMzCplNyzKMY6yUYGVV1IQvBHiAIAS5ln2xLjG6u2HH9WBmbg0SkVjrPjZy5pWK3HMVnBmGYEOwEghDgUvbKM5W+saepvK1qWvCI9x0cqUVRc7ZW7mLEP7FmaoS96B0Fu4AgBLio2h7WY2DjZNiDcFvl7jkRM1QKFffK5xjvG+2kcDrZVMK98iy0CMFeIAgBLmpPA5ulEeXoGN1WuWfh6EwZCl/AgqjZW2SYUDhhlNCpZ3U9yEKweQhCgIuSqV/0ZPMpRmyC7zjulS9oQVTm3jP7B4yDfMsKRBmBIhqFYAcQhAAXJdMMwq0VuxeNnsu97MX4ufiM943eX3uYe+VZQRgmBHuAIAS4sMY+aupnE705B6FBMu45k7Mgajbfspe2cHTmtso93MvO0AjZjQhCsHkIQoAL26eVMjQi9xHCY9oT4Z6hAa5+nOteUkbo1DxdYddgN9+yiT5CfS/T9fGtCmBuCEKAC8vSshky9IvuPpOTGZHOveyluSjVKZqEA3VH+JYVBZoWIOzXYScKsG0IQoALk2OA0CAZ99cenhU2nW/Z4Zgdnr67Ood72RkaMQvvy4CNQxACXEDHIFV0smRfzkFo6hf1c/XlW3Y40kOn5Ded7NH38i07Q4PVt8HmIQgBLiC7kU0JEFS8nw+L9IuauCjVyYHxObWH+JZN8xdKO1iXnm9VALNCEAJcQJZWmqHh/HRYsF/URI7eUSeRUvyE/Xh3FGwZghDgAvY18H9TJrcxL9wzxCL9oibTQybn6Qq5947O1AhZWrwvAzYMQQhwrj4DFbSxyf6cgzC75tCMsGl8a46Iq8olISD2cH0u37J4XwZsHYIQ4FwHm1iCj+Cq5FmTEdtfd2R6SBrPoiOXHjplfx3nJWamBQq5LazfyLcqgPkgCAHOJUe/aGlLuZvKRdZteIdjesjkA3VHDRLP1HJT0oRRwpEmNArBViEIAc6V3ShlBHJ+NHLqDqWHTuFb8wr4uniHegQXNBXxLZsRKKB3FGwXghDgJwwSHdax6YH8BwjTQyfzrXll0kMnc59EMUMjZDfifRmwVQhCgJ843sIiPAQfZ541G3t0rf1tE3zH8yx6pdJDp2TVHORbc4ZG3N/IjGgTgm1CEAL8RHYjy+DeHKw9ND1ksijIscXviI0eFSEKYmV7NceafmoKchUKWpGEYJMQhAA/Icda2/vrjky3jn5Rk+mhaTm8F+DGMCHYLgQhwP8wopxGKYNrEPYb+oubS1MCEzjWvEpTgicdqj/Gtyb2JgTbhSAE+J/SDuaiEMLceAZhbmN+jO9YV5ULx5pXKSkg7nRbRedgF8eaMzTCvga8LwM2CUEI8D9y9IserDs2JXgS35pXyUnhlOA/MVebz7FmlIegFIXyTjQKwfYgCAH+J6eR8e0XJaIjDcetLQhJnt7RDA2GCcEmIQgB/idLy/mV0TOdtQbJEOkVxrEmF9NCUg/VH2PEM7fSA4UcDBOCDUIQAvxI20dtA2zCKJ5BeLD+2NSQVI4FeQlyD3RRupS3VXGsifdlwEYhCAF+lKWVMjSiyLVn9HB97uSgFJ4V+eHeOxrvLWh7ma6PY0kAc0AQAvwoW8vSufaL9hsGippLJmkSOdbkaGrIpINcg1AUaFqgkIO11sDWIAgBfsR9TZk8XeFY79FWNXHibIkBcafbKvju05seKGKYEGwOghCAiKhbT6UdbJIfzyA82nAiNSiJY0G+nBVOE3zH5ekKOdbEMCHYIgQhABHRfh2b5Cc4K3jWPKK16iAkotSgpCMNJzgWnOIvnGxjvQaOJQFkhyAEICLK1kp8+0Vb+9qae1vG+0RzrMldqibpKNcgdFZQgo9wCJv0gk1BEAIQ/fimDM/H4SKukK4AACAASURBVKj2RIomQRSs+hGL9h7dOdCl623mWBOrb4PNseqnFMA89BIdbea8Ge9RbV6qxqr7RYlIFIQUTcIxbR7HmhkaIVuLF0fBliAIAehYMxvjKXg58ayZq8238gFCk0maRL69oxmB4iEdMyAKwXYgCAEoh/fEicqOMwpBDHbXcKwpk7SgpGPaExzXWvN2pjB3IR+b9ILtQBACUE4j56n0xxry0oKSORaUT6BbgLuTe3lbJceaGYGYRAG2BEEIjs60GW86100njmnzrHZBmfNN0iTy3ZIpXYPVt8GWIAjB0ZV2MFclz814JSYVNBUlBcbxKii35MD4440FHAumB2KTXrAlCEJwdNm8t14qaT0d4OrnrR7FsaasUgIT8nQnjczIq+BoD0EpChVdaBSCbUAQgqPjPkCYq81P0SRwLCg3T2cPjVtASUs5x5rTA4VszCYEG4EgBEeXzXtX+uONBcmBthSERJSiSTjeyHWYEJv0gu1AEIJDa+yj5n4Wy28zXr1kKGouSQiI5VXQPJIDE/gOE2agRQi2YwRB+Pnnn997772rVq3SarWXOOzLL798++23r/rCAMwhWyulBwocN+Mtbi4J8wzxcHLnVtEskgLjTjaf0hv1vAom+gr1vaxlgFc9ABkNNwj//ve/r1q1atasWS0tLRkZGQMDF77BDx069NBDD7300kv8rhBARjmNbDrXJUZzG/OTA+M5FjQPN5VruGdocUspr4IKgSb7C/uxSS/YgmH9CjAajW+88ca77757xx13vP/++x4eHt988835hw0MDDzwwAMvvvgi74sEkAv3zXiPawtSbG2A0CQ5MD6X6zDhdGzSCzZiWEFYW1tbU1OTmZlp+nT27Nk5OTnnH/bSSy8tXbo0NtbGRkfAYfUYqLidpflzC8JB42BJ6+l4WxsgNEnRJBxv5L1JL4YJwRYoh3NQQ0ODu7u7s7Oz6VN/f/+jR4+ec0xeXt769euPHDly8ODBS1erra3Nyspavnz50Fd+//vfx8dftDepr69PoeC6XyqMXG9vryAIgsCz8WRxexvFBG9RGujt5VQwv7ko0jOcDUq9g7xK/kRfX5/BYJDpcRjjFlnSUtbe1e6k4LP6eII75bWqWrt61fb1+A4MDIiiqFKpLH0hDk2v1xuNRkm6fN+7Wq0Wxcs0+YYVhGq1enBwcOjTwcFBFxeXsw8wGAz33nvv+++/PxSWl+Dj4xMREfHzn/986CvR0dFqtfpix+v1+kv8K5iHwWBQq9V2FoSHW9kMDeN4dxW1laZo4uW7XRljTk5OMgWhmtRRo8Ire2oSAyZyKkixo6SCLucZXGenWJwgCAhCi1MoFEajcTjP2mVTkIYZhCEhIQMDA01NTf7+/kRUU1MTEhJy9gEVFRV5eXl33HEHEfX397e2to4ZM2bHjh1RUVHnV3N1dQ0PD7/55puHc2oiEkVxON8JyMr0U7CzIMzRGR6dqBD5vTOa13RyxcTl8t2u4v+RqX5iQFx+08lkDbeXfTI0bH+TMCvYrp5fuX8KMByiKDLGeP0UhlXF399/xowZn3zyCRG1t7dv2rRp2bJlRNTa2rpu3ToiioqKOnXq1Pbt27dv3/7GG294eXlt3749NDSUyyUCyMEg0WEdz8149UZ9SUtZnN8EXgXNLykw7oSO5zBheiA26QUbMKwWIRG9+uqrS5cu3b17d1FR0YIFC6ZNm0ZEp06d+tnPfsYYU6lUo0ePNh1ZU1OjUCiGPgWwTidaWbi74HP5vvzhKmopjRoV7qpyufyh1irBf+KL2a/rjXqVgk+/3wyNeM8+o5GRwq66EsDeDDcIp0+ffurUqYMHDwYFBSUn/7jRWkpKSklJyTlHTpky5ciRIzyvEUAG2VrOK6udaCxIDLCZHScuyFXlEuYZcqq1LN6fz4uv/moKdBFOtrEEHyQhWK8RdLD6+Phce+21QylIRGq1ety4ceccplarIyIi+FwdgGy4r7V9QldoQ1svXUxyQPwJrpMosNYaWD+M94KDymmUOE6l10uGUy3cGlIWlMh7mDBDg93qwdohCMERne5koiBEenALwlMtpWGeIW4qV14FLSUxYGJxc6lB4rY3YUagkIUWIVg3BCE4oiwtm8l1gDBPV5Rk4wOEJm4q1xCPoJLWMl4Fx3oJRsaqu5GFYL0QhOCIcngvMZqnK7T1N2WGJAZMzNcVcSyYHihimBCsGYIQHFEW11dGJSYVNZfE+9vwDMKzJQRMzNed5FgwIxDDhGDVEITgcJr6qbGPTfTmFoRlbRUBrn6ezh68ClpWQkBsflORxLhNhM/QYJgQrBqCEBxOllZKDxQ4TvHO051M4LQ+pzUY5ezl6+JT0V7Nq2CSr1DTjU16wXohCMHhZGtZhobnnZ+vK7KnICSixICJefx6RxUCTQnAJr1gvRCE4HCyG9kMfm/KMGKFTUW8dmywEgkBsZyHCTV4XwasF4IQHItpM95UfpvxVnfUuihd/Fx8eBW0BkkBcXm6k4y4RdcMDBOCFUMQgmM5qGNJvgLHrWLzdScTbX9ltXP4u/o5K5xqO+t5FZziLxS0sT4Dr3oAPCEIwbFkaSW+U+nzdUWJtr+y2vkSA+M4DhO6KineWzjchEYhWCMEITiWLC2bwfVNmTxdYXyAHQZhgv8EjkFI6B0FK4YgBAeil+hIE5sWwK1FqO3R6SVDqEcwr4LWI94/tqCJ5/oyMzRiFjbpBauEIAQHktvMoj0FLyduBfN19va+6JBwr9B+Q7+ut5lXwfRA4aCOGRCFYH0QhOBA9mnZDM4DhCcT7LFflIgEEuL8JxQ0FfMq6O1MER7CiVb0joLVQRCCA8niHYQFTUV2sAfhxcT7xxZwXX0bw4RgnRCE4CgY0f5GKT2Q2z3fMdCp620e4x3Jq6C1SQjgPUyIvQnBKiEIwVGcbGO+aiGI39a5BU3FcX4TFAK/OYlWZpzPmIbuxu7BHl4FZwWJ2VoJSQjWBkEIjoJ/v6iuKD7ATrZeuiCFoBjvG13YzG2YMMiVPFRCSTuiEKwLghAcRZaW82a8+U12tenEBSX4T+Q7TDgzSNiH3lGwMghCcBTZWjYriFsQ9hsGKtrPxPiM5VXQOiUExOZxDcIMDBOC9UEQgkMo72SMKMqDWxAWt5SOGRWpVjrzKmidYv3Gn26r0Bv1vArO1Ah7GhCEYF0QhOAQ9nFtDpJdzyA8m4tSHe4Veqq1jFfBsV6CxKiqC1kIVgRBCA6B+1T6gqZiRwhCIkrwj83HMCHYNQQhOIR9DYzjphMSk4pbSuP87PmV0SHx/rGF/NaXIcwmBOuDIAT7V9fDuvQsZhS3IDzdVunv4uvp7MGroDVLCIgtaCqWGLfoQosQrA2CEOzfXi2bGSRy7BgtaCqyy62XLshbPcrL2bO64wyvghO9hdYBVt+LLARrgSAE+5el5dkvSkT5uqJ4f4foFzWJD4jN57fWmkCUEShmo1EIVgNBCPZvL9cBQiIqbCq247W2zxfvP6FAx3WYUIPeUbAiCEKwc039pO1jCT7cgrC+W0uCEOQeyKug9Uvw59kiJKJZQcJezCYEq4EgBDu3t0HKCBQ5jhDm64oSHak5SERhniF6Sc9xk94kX6GulzX386oHcFUQhGDn9jZwnkpf0FQU52BBSERxfjEcFx1VCDQtQNinxXb1YBUQhGDn9vAOwnxdUYJdbzpxQfH+E3j3joroHQUrgSAEe9YyQLU9LNmXWxB2DHS29LWOHhXJq6CtiA/gvFv9bAwTgtVAEII929sgpQcKCq4DhBP9Y0TB4R6ccT5jtD26zsEuXgVTfIWqbtYywKsewJVzuOcZHMreBjYriOdNXthUnOB4A4T0f5v0FjWX8CqoFGlagJCFYUKwAghCsGf8BwgdaU2ZcyT4TyzguugohgnBSiAIwW61DlBVF0vhN0A4YBysaK+2+814LybefwLfbSgwTAhWAkEIdmufVkrXCEp+93hRc4kjbMZ7MRP9Y0pbyweNg7wKpvoJFV2sFcOEYGkIQrBbexrYLA3POzxfV+QgexBekItSHeEZWtJazqugUqSpGCYEK4AgBLu1u55lBnNeYjTOkdbaPl98wIQCrrMJZweJe9A7CpaGIAT71DJA1d08BwglJp1sPuVQm06cL96f82zCzCBhVz2CECwMQQj2aXe9lBHIc4CwvL3K39XXy9mTW0UblODPeZPeSX7CmW4sOgoWhiAE+7S7gWUG8x0gPOlQWy9dkI+Lt6ezB8dNepUipQcKexswTAiWhCAE+7S7nmXyX2LU0YOQiBICJvJddDQzWNyNYUKwKAQh2CFdH2n7WCK/AUIyrSkTMJFjQRuV4B+bpzvJsWBmkLAbw4RgUQhCsEO76qVZGpHjEqN1XQ0kCBq3AG4VbVZCAOcgTPIVGvuYto9jSYCRQRCCHdrdwHniRH5TUSKag0REFOoRzJjU2KPjVVAUKEMj7qnHMCFYDIIQ7NCeBjab8wDhScdca/uC4v1j87hOopgTLGCYECwIQQj2pq6HtQ6wOG++u9JjgPB/4v1j+U6rnxMs7MQwIVgOghDszfY6NjdYFPnlYPtAR3t/R6RXOLeKNi4hIDaf6zDhRG+hR8+qupCFYBkIQrA3O+vZXL4DhLqiOP8YUeBZ06ZFe0c19bZ0DHTyKigQZQaLaBSCpSAIwd7sbmBzQ3gPEKJf9CyiIE70j+G7N+Fc9I6C5SAIwa4UtTOlQKM9eAZhnu4kXhk9R6L/RL69o/NChJ31EpIQLAJBCHZlZx2bz7U52KPvre2qH+cTzbGmHUgMnMh3k94Id8FDJZxsQxSCBSAIwa5wHyAsbCqO8R2rEpUca9qBGJ+xVR1nevU8p8HPDRZ21iEIwQIQhGA/jIyytNJs3mttJ/ijX/RcKoVqvE90UXMJx5oYJgRLQRCC/TjaxELcBI0Lz5p5Oqwpc2EJARPzm3gOE84JFvdpJT1WmAGzQxCC/dhRz+Zx7RcdNA6Wt1fG+o3nWNNuJATE5jXyDEI/NY32EA7p0CgEc0MQgv3YVistCOV5Sxc1l0R5RaiVzhxr2o04/wklracHjYMca84PEXZg0VEwOwQh2IkuPeW2sBkavhMn0C96US5KdYRXWElrOcea80PEbbVoEYK5IQjBTuxpkKb4C25c3+7MbzoZj814Ly4xYCLfLZkyNEJhG2sb4FgS4PIQhGAnttex+SE872eDZCxuLsWmE5eQEDAxr7GQY0G1gqYHCnsa0DsKZoUgBDuxrZbzVPqS1rIQjyB3JzeONe1MYsDEk82nDJKRY835IeJ2zCYE80IQgj2o7WEtAyzJl2cQnmgsTAqI41jQ/ng4uQe5B5a18R0mFLYhCMG8EIRgD7bUsvkhPLdeIqITusLEQAThZSQGxJ3g2jsa7yP0GlgFtmQCM0IQgj3YwXuJUSMznmw6hQHCy0oK5ByEAtG8YLw7CmaFIASbZ2S0o05aGMozCEtby4PcAz2dPTjWtEtJgXGFzcUS4/l6y8JQ9I6CWSEIweYdbmKhbkKwK+cBwkQMEA6Dp5OHv6tfWVsFx5oLQ8Xd9dIgXh0Fc0EQgs3bUiNdE8Z5+/g8XWESBgiHJ4n3MKGfmsZ5Cfsb0SgEM0EQgs3bXMsWcV1ZTWJSYdOpeAwQDg/3YUIiWhQqbKlFkxDMBEEItq25n8o62PRAni3CsrYKPxcfb7UXx5p2LDEgrqCpiO8w4aIwcUsNWoRgJghCsG1ba6XMYFHF9UbO1eYna+J5VrRr3movXxdvvsOEk/2Ful5W14MsBHNAEIJt21LLruH6vigRHW8sSA5M4FvTvqVoEo43FnAsqBBoXoi4Fe+OglkgCMGGSYy28544YWTGwqZibDoxIsmBCce1PIOQTMOE6B0Fs0AQgg072sz81EK4O9clRltOB7lrvJw9Oda0e0mBcQVNRXwXHV0UKu6ox4b1YA4IQrBhP5yRrgvn3C+aq81PCcQA4ch4OnkEuWtKW09zrBnoQtGeQjYmUYD8EIRgw344wxaHcb6HcxvxpsyVSAmMz9Xm8625OEzceAZNQpAdghBsVX0vO9PNpgXwbBHqJcOplrIEfwwQjliyJj63kXMQXhcu/HAGLUKQ3QiCMCcnZ8GCBUlJSU8++eTAwLl7SOt0uqeffnrmzJlpaWkPPfRQY2Mj1+sEONcPZ9jCUFHJ9W+5ouaSMM8Q7EF4BRID4opbSvVGPceaKX5Ct4HKOpCFIK/h/hZpbm5evHjxLbfc8umnnx44cOCFF14454DS0tL+/v4XX3xx9erVdXV1y5Yt43ylAD+1sYYt5j1AeKKxIBkDhFfETeUa4RlW1FLKsaZAdG2YsBHvjoLMhhuEa9euTUtLu+eee2JjY//85z+vXr1ar//Jn34ZGRl//etfMzMzk5KS3njjjQMHDnR1dclwwQBERANG2tsgLeS6shoRHdPmpWgwg/AKpWgSjmnz+NZcHCZsrMEwIchruL9HCgoK0tLSTB+npqa2trbW1dVd7ODc3Nzg4GAPD2xhA3LZVc8SfQRfZ541+w39ZW0V2IPwik3SJObyDsJ5IeIhHesY5FsV4CeUwzxOp9PFxMSYPlapVO7u7o2NjZGRkecfWV9f/8gjj7z11lsXK3X69Onvv/8+Kipq6CsffPDB9OnTL3Z8T0+PIHDuAYOR6unpYYxZzw/iuwrlvEDq7u7nWPNI4/ForyhDv6GbujmW5aivr8/JyUmhUFj6Qi4syiXsdFulrr3JVenCsexUP9WG8t6lYdbSLhwYGBBFUaVSWfpCHJperzcajQaD4bJHurq6iuJlmnzDDUIvL6+enh7Tx5Ik9fb2jho16vzDdDrdvHnzHnnkkZtuuulipcaMGTN37ty//vWvQ1+JiIi4xLPNGHN3dx/mdYJ83NzcrCQIGdHmesP2axXu7mqOZYtKSieHpFjzzaZQKKw5CIko1m/86Z7K6SGTOdb82Whpa6Py9gnW8l2rVCoEocWZglCt5vMbYLhdo1FRUaWlPw6Dl5eXKxSK0NDQc45paWmZP3/+zTff/NRTT12ilCAI7u7uo89izQ82WKFjzcxdReO9OKfyUe2JSZokvjUdzSRN4tEGzr2jSyKEjTVYYgZkNNwgXLFixaZNm8rLy4norbfeWrp0qZubGxGtXbt2x44dRNTa2jp//vz09PRHH320ra2tra3NaOS53hLAkO+rpRsiOKdgW397U29LjG8037KOJlWTdEx7gm/NYFdhnJewT4t3R0Euww3CCRMmPPfcc5MmTQoNDT148OBf/vIX09c3bNhw4MABItq3b19VVdXnn38+5v9UVVXJdNHg4L6rYjdEcH5f9Kj2RFJgvChgiYmrMtZnTFt/R1NvM9+yN0SI31ejSQhyGe4YIRH97ne/e/DBB7u6uvz9/Ye++OWXX5o+WLp06dKlSzlfHcB5yjtZcz+b7M+5RXhMmz8JEyeumigIyYHxudr8haPncCy7NEJYtEV6axpZxRg12J2R/f2rVqvPTkEA8/uumt0QIYq8fyPmavMxQMhFiibhGO+11iaMElwUdKIFvaMgC3QEgY35vlri3i96prOOiIV7hvAt65jSgpKPNhxnxDm0lkQI6B0FmSAIwZY09lFhG5sTzLk9eLg+d3JQCt+aDivYXeOidKlor+ZbdmmE+E0lWoQgCwQh2JJvq6Rrw0Rn3tNtDjfkTg5GEHIzOTj5cH0u35rTAoWOQSpuRxYCfwhCsCVfV0rLIzk3BweNg4VNxZM0iXzLOrK0oJQjDcf51hSIlkYKaBSCHBCEYDOa+ym3mS3ivdB2nu7kGO8oN5Ur37KOLCUwvriltM/AcwE8IroxSvymCsOEwB+CEGzGt1XSwlDRZQRTfoblcMPxyUHJnIs6NrVSPcF33PHGAr5lMwIFXR+d7kSjEDhDEILN+KZSujGK/0SyI3hTRgZpQclHGjgPE4oCLY0UvkbvKPCGIATb0DZAh5rYNWGc79im3ua2/o6xPmP4loXJwSmH6zkPExLR8kjxm0r0jgJnCEKwDeuqpXkhohvvftFD9bmpQUmideyqYU9Gj4roN/TXd2v5lp0VJJzpYZVdaBQCTwhCsA2flUu3juYfVwfrj04NmcS9LAgkTAmedKDuCN+yCoFujBI/r0AQAk8IQrABuj461syu5d0vqjfqjzcWTAlCEMpiWkjqgbqj3MveOlr8uAy9o8ATghBswOcV0vXh/N8XPd5YEOUV4enswbkuEBFRWlByUXNJr76Pb9l0jdBrpMI2NAqBGwQh2IDPyqVbx/C/Vw/UH5kWksq9LJiolepYv/HctycUiG6OEj4rR6MQuEEQgrWr7mYVXWwe7/VFiehg3bHpIWncy8IQuXpHx4iflzM0CYEXBCFYu0/L2Y1RopL3rVrZccbIjFGjIjjXhbNMD5m8v+6wxDuzknwFtYIO6xCFwAeCEKzdZ+XSraP536gH645OD5nMvSycLcg90NPZs6y1nHvlW8aIn6B3FDhBEIJVy21mXXpK18jRL4qJE+YwPSRtP+9JFER0R7TwWbk0iCgEHhCEYNX+WybdPZb7dvTUMdB5ur0yJTCBd2E41/SQtOzaQ9zLRnoIsaOEzTVIQuAAQQjWyyDRlxXS7dH8m4M5tYfTgpKdFE7cK8M54vxjW/vbuC8xQ0R3jRP/W4ZhQuAAQQjWa2ONNM5LGOPJPwizag7OCJ3KvSycTxQEmRqFN0eJexqkJs57PYEjQhCC9fqojN01jv8t2m8YyNMVTgnGAKGZzAibml1zkHtZdxVdGyZ+jldm4KohCMFKtQzQrnrpxij+t+ih+mMT/WLcndy4V4YLmhSYWN5e1drXxr3y3WPF/2K5NbhqCEKwUmvLpOvDRU8V/8pZNQdnhKFf1HxUClVaUPKBev4z6+cEC839dKIFI4VwVRCEYKU+LJHui+F/fxok46GGY+mhU7hXhkuYETYtq+YA97KiQPeMF1eXoFEIVwVBCNYoW8uMjDJkmD54orEgzCPE18Wbe2W4hKnBk/J1RT36Xu6VV44TPi+Xeg3cC4MDQRCCNVpdIv0yhv/0QSLaVZ2VGZEhQ2G4FDeVa1JgXE7tYe6VQ9yE9EDxywo0CuHKIQjB6nQM0vpq6fZoWfpFs2sPzQybxr0yXNaciBm7q7PlqHxfjIDeUbgaCEKwOmtPS9eEif5q/pWPao9HeIUFuvnzLw2Xkx46JU9X2DnYxb3ytWHimW7sUAhXDkEI1oURvVck/XqCLHfmrurszHD0i1qGi1I9SZMoR++oQqB7x4vvFaFRCFcIQQjWZUcdUwg0U4bXZPRG/f7awzPD0S9qMXMiZuyqzpKj8q8miJ+VS20DctQG+4cgBOvy9knpN3Gy3JaHGnKjvaP8XHzkKA7DMS0krai5pH2gg3tljQtdGyauweR6uCIIQrAiVV3sgE66dYxM/aJ4X9TC1ErnyUEpWTIst0ZED08U3y2SJAwUwsghCMGKvF0k3TNOdFXyr9yj7z1Uf2x2eDr/0jAS86NmbavcI0flqQGCjzNtrkUSwoghCMFa9BjoozLpfnlek9l7Zn9KYIKXs6ccxWH4pgRPqu2ql2NXJiJ6KFb8e6FRjspg3xCEYC1Wn5LmBouRHnJMo6etFbsWRGXKURlGRCEo5kRkbKvcLUfxW8eIpzroOJYehRFCEIJVMEj0t0Lpt/K8JtPYo6vsODMV+y5ZhwVRmVsqdjHiH1cqkR6KFd8owCszMDIIQrAKX1RI0Z40JUCW5uCWit1zI2eoFDLsZAEjN94n2kWpLmw6JUfxX08Qt9ZK1d1oFMIIIAjBKrxRID2RoJCp+I6qPQuj5shUHK7A/KjZMvWOeqho5Xjxb4VoFMIIIAjB8rbVMSOjBaGyNAfzdSdFQYzxHStHcbgyCyJn76nO6Tf0y1H8kYniR2VSCybXw7AhCMHy/phrfDJRlr0miGh92dbroxfJUxuukJ+rb3zAhF3yrMEd7CrcGCX+tQCvj8JwIQjBwrbXMV0//Xy0LLdi50DXwfqjC0bPlqM4XI0lYxetL9siU/FVSeI/iqVmWRqcYIcQhGBhLx03vpAiKuRpD26p2Dk9dLKnk4cs1eEqTA6a1NbfXtpaLkfxcHdheZT41kk0CmFYEIRgSdvrmLZPruYgEW0s374E/aJWSRSExdHzfzi9Tab6f0gS3y+WWjFSCMOAIARLejFXxubgicYCRhTnHyNLdbhq141ZsKs6q1ffJ0fxcHdhWaT4BkYKYRgQhGAx31dLXXq6Rbbm4LrSTcvGXStTcbh6Pi7eKZqErZW7ZKr/fLL4QbFU24M5hXAZCEKwDCOjVUel1yYrZHpbVNujO95YsGj0XFmqAyc3xdzw9akNEpMlq0LchHvGiy8dx5xCuAwEIVjGv0skfzUtlGfuIBF9fWr94jHzXZRqmeoDF/H+Ezyd3Q/UHZGp/jNJiu+qpaJ2NArhUhCEYAG9Bvrjcen1KXItJdOr79tauXvZ+MUy1QeOboy54ctT38tUfJQTPZGgWHUEjUK4FAQhWMCrecaZGiHVT67m4A+nt04OSglw9ZOpPnA0O3x6Q7f2VEuZTPUfihXzWtmuejQK4aIQhGBuFV3sH8XSa5PluvcMkvGbkh9uilkiU33gSyEolo1b/PWpDTLVVyvozaniw/uNejQL4SIQhGBuvzkgPZGgCHGTqzm4rXJ3iEcQFhe1IUvGLjrScLyms06m+ksjxEgPeqcISQgXhiAEs9pUw0o72G/k2XeQiCQmfVr0zV3xt8hUH+TgpnJdNv7aT4q+ke8Uf52qeOWEUSvLlEWweQhCMJ8eAz203/j2dIWTbPfd9so9vi4+iQET5ToByOOmmBv21x6u7aqXqf44L+He8eKjBzC/Hi4AQQjm88wR4+wgYX6IXJ2iEpM+PvnVXXE/l6k+yMdN5bp03DWfFX0r3ymeT1EUtrJvq9BBCudCEIKZHNSxb6rYX2SbMkFEO6r2J+0+3gAAFoZJREFUejl7pWgS5DsFyOfGmCX7ag7Ud2tlqu+soA9nKh45ILVhAVL4KQQhmEO/ke7ZZ3xrqujjLNcpBo2DH+Z98qvkO+U6AcjM08nj5pgb/nniI/lOMS1AWBYhPHYIHaTwEwhCMIcnDhkTfYXlUTLeb1+XbIjxjY73j5XvFCC3mycsLWouKWw6Jd8pXklT5DSyryvRQQr/gyAE2W2uYevPsPemy9gp2jXY/WXx979MQnPQtjkrnH4Rf+t7uf9mJNf8d3cVfZapeHC/sbobU+zhRwhCkJe2j1ZmGT6ZrfCWrVOUiP6T/1lmRHqoR7CM5wCzWDh6Tr+hP7vmoHynmOQn/DZOcddeoxFRCESEIARZGSS6ZZfh1xMUGRq53hQlorK2il3VWXfH3yrfKcBsREF8aNK9bx/7V7+hX76z/D5BVAr0/DEMFgIRghBk9btDRmeRViXJeJtJjL15+B+/SrrTy9lTvrOAOaVoEhIDJv4n/zP5TiEK9OVc5WflGCwEIgQhyOfTcmlzLftirlKmDehNvivdpBSVi8Zg30G78vCke7dV7Slrq5DvFD7O9M08xYP7jSfb0EPq6BCEIIucRvbbg8Zv5ylGOcl4Fl1v838LPn988gMCyRm2YHaezh4rE1a8ceg9I5Ox9zLJV/hzmmLpdmOTjL2wYAMQhMDf6U52807jR7OUcd4y5pPE2CsH/nZTzJIIrzD5zgKWsjh6vruT29rCL2U9y93jxNvGCNduMfQYZD0PWDUEIXDW0EsLNhtfThPl233e5IvidXqj/raJy2U9C1iKQMJT0x79rnTzyWYZpxUS0QuTFHE+wq27jAYMFzoqBCHwpOujuZsMv4oR7xor7611uq3i86J1f0j/nSjgHrZbfi4+v5t8/59y3uzVy7hthED0zwwFI3b7HkyocFD4JQLcNPfTvM2GW8aITybKe191DnY9l/XnR1Lv07gFyHoisLgZYdNSg5JePfiWfFPsiUgl0jfzlD0GhnahY0IQAh91PSxzo+GGCOG5ZHlvKolJf8p5c2bYtLmRM2U9EViJR1N/1d7fIfdgoZNIX81Vtg+y2/cYB5GFDgZBCByUd7LZG43Lo4SXJsm4jprJP46vkZiE1dQch1JUvJDx+w1lW3NqD8t6IrWCNixQMqJFmw2dellPBdYFQQhXa28Dy9hgeDZZfCFF9hT8svi7g3VHn894AkODDsXHxfulmU+/fujtgqZiWU/krKBPMxXjRwmZGw01PRgwdBT4bQJX5d0i6ZZdho8zlXfK/HYMEW2p2PV1yYa/zHnBw8ld7nOBtYnxHftc+hPP7nu5pPW0rCdSCPSPdMWtY8Sp3xuztchCh4AghCvUpac79xj/eUrKWaKcGyz7fPZtlXtW5619c+5LgXhBxlGlaBJ+N/mBp/e8VN5eJfe5Ho8X/zNLceNOwxsFEsLQ7iEI4Uoc0rHkdQYXJR1YohztIXsKfn1qw4d5H78554/YX8LBzQib9mjqr36387l83Um5z7UgRDi4RLmuSlq02dDQK/fZwJIQhDAyPQZ6/JBx6XbD65PFDzIUrkp5Tycx9sHx/64/veWdBa9iBRkgolnh059Lf/y5rFf31RyQ+1yRHsKexcrpgWLSOv2/StA0tFsIQhiBdVVS/DcGXR8VLFcti5T95ukc6Hpyz4vFLaXvzH81wNVP7tOBrUjRJLw+58V3jv3rg+P/lZi8cx2UIj2fIm6/RvlhiZS50ZDXijS0QwhCGJZjreLsjcYXcqUPZyg+mq3wU8t+xoKmol9ueWy0V8Qbc//o6ewh+/nApoz1Hr36mjdL28p/t+u5xp4muU+X4CPkXK+8ZbS4aLPhV/uFOvSU2hcEIVzG/kZ27VbDHTmqO6KF3GXKOfK/F9Nv6P/70dUvZL/+SOp996f8QiHIPisDbJGXs+frmS9M0iT9cvNj68u2yLr0DBGJAv16glhysyrQhVLX0/05xqoutA7tBIIQLmxQok9OS1PXG+7ca1wWIZ5YPLByvCjrzoJEJDG2pWLXHRse6NH3rFn89vSQyfKeD2ycKIi3T7zxrfkvb67Y+cDWJwqaiuQ+o6eK/pjMCpaSrzOlfme4cadxTwPi0ObJ/KoD2KAjTWztaenzcinRV3gmUbwuXBQF6u6W96QSk/ae2f/xya9clOoXZzwZ6zde3vOBHYn0Cntv4Ws7qva+lPNmtHfkiok3TvSLkfWMfmr6U6riqUTFR2XSgzlGI6M7xoq3RwsR7tgX0yaNLAj1er1Kpbr6Y8DaDEqU08h+OCN9W8VUIt0eLR66QRkl/7wIImrr79heuXtd6SZfF5+VibdPD0kzw0nBzggkzI+cPSts+qbyHX/KedPP1feGsYtmhE1zVsi4MbS7ih6IFR+IFQ/q2EdlUup3xkh3YXmUeE2YkOAjIBJtiMDYsNr13d3dd9111/bt2xUKxZNPPvnUU0+df8xrr732yiuvGI3GuXPnfvTRRx4eF37B4dNPP924ceMnn3wyzEvs6uq6WCm4GgNGym1hWVqWpZWytCxmlHBtmLgsQoj3ucAj3N3d7ebmJvB7ujsGOvfXHcmqOZCnO5kROmXJ2EVy/xVvB/r6+pycnBQKDJpeisSkfTUHNp7efqqlbEbY1Blh0yZpEpz4JeLAwIAoiuf/uW9ktLeBfVctba1lXXo2N1icoREyNEKMlyAiFXnT6/VGo1Gt5vPa3nCD8Jlnnjl27NiGDRvq6+unTJmybt266dOnn33AoUOHFi9efOjQobCwsGXLlsXFxf35z3++YCkEoaXU9LDSDipsZYVtLLeFnWpnE0YJGRphpkaYHST6OF/q/3IJwrb+jlMtZXm6whONhWc6a9OCktNDp2SETnFVuVxNWceBIByRpt7mPWf2Z9ceKmstj/UbnxwYH+8/YazPGBflVf32vFgQnq2yi+1pYPu0LKeRNfayZD8hyVeI8xbivIVxXsKlnzUYDssEYUhIyJo1a+bPn09Ev/nNb3p7e//5z3+efcD999+vUCjeeecdItq9e/ctt9zS2Nh4wVIIQvkYJGoZoOZ+1thH9b2svpfqelhlF1V1s4pO5qGimFHCRG8h3kdI8hESfQX1sH+jjjQIB42DTb0t2h5dQ7e2tquhsv1MZUd1r74vxndsnH9MSmDCBL/xKhFD1CODILwy3YM9ebqTxxvzi5pLy9urNG7+kV7hEV5hoR5BQe6aILcAb5dRw385eThBeLb2QTraxPJaWWEbO9nGSjuYUqDRnkKUhxDmRmFuQrAbBbkI/i7kr0ZGDhffIBzWb6Le3t76+vrY2FjTp7GxsZ9//vk5x5w+fXrZsmVDB+h0us7OTk9PzwsWHBwcbGtrG/p01KhRl/gN29Xb09Zne9N2DIx6DJf/I6PHQOdsBNqlJ9MfJxKjLgMjIqNE3QYiok49EaMOPTNI1KOnHgP1GqlnkHXqqVNP7YOsx0A+TuStFvydyd9FCHSh0S7CjFAKdRPC3AT3nz65A900QEREvfo+40VmJfca+ozMSET9ff2CSjAwIxF1D/YyYn2GfgMzdul7+g39fYaBLn1Pt763c7CrY6CrbaCj3zjopx4V4OIX7BYQ6hG0JGJ2hEdo8NnLhPb3Y9O3kWL9/ZJhUEAQjpAr0TSf2Gk+sURkZMbqrvrqztqqrrqDZw5re3WNvc0dg12eTu5eTh6eTh5ezh7uKlcPlZuLUq1WOrsqXZwVKifRSSWq1EonIhIlQaVQKZVKInJSODmLF+50dVY4OSlURCQQpblTmjtR+I//1NJPNd2stofV91FNMztyhpr6WesANfWzXgONUpGns+ChJA8n8lAJrkpyFsnLmRRE7ipBKZCriojIy4mISCTyUP34y1OlIJef3hpqBTkP42bxUNnSZi4qlTLI25dvzWEFoSm03N1/XPLf09OztbX1/GOGDjA14Nra2i4YhGVlZevXr9+xY8fQVz777LOMjIyLnf3F9U+dcW672L86AoGILjlHyp3+f3t3H9PU1QYA/LltHS2CINqOQlap8QP8QoS5aXQoIoZMotE635G6CEi2bGZLwH0k24wJS/hDt8U448cCi2bZUsVtWd2yZejUMjeFDSFpRi2WyGcZhbYKBdp773n/uNtNX4pYX9te1vv8/jqePtz7KL0+7bnnngNxAClSAO5NPw7MONjvgR2gJYjjK1hK8oDDK1iQkr9PH0MoGUu4TgkBOaGkLMxkIIaAnIHZLDWToeIZiKMhgaHiaQBwAjgBrPzR+oL52yIUZgqAdAD/O9IsRbllI/elw/ekvSMyakRKRiTglpJ+CTVGEa8EfBShJTBOAQCMSYD953O7lwLvA2qIV0J8Qd9JkAIoAZT8JewD8AE9+vclNCkCIr3xKGXhaP6xpJnxDMPQNP3Q+NjYWInkIYU+qEI4d+5ciqLcbndCQgIAOJ3OJ598ckKMUql0u91c2+VycT2THm3hwoU6nS74odEj/zmGQ6OCC/lkGfR/wKHR8Al+HdtHHRpF4RDaodGgvhDHxMTMnz+/ubmZ++OtW7fS0ydO8MvIyPAPSEtLi42NDUmKCCGEUPgEOzJcXl7+wQcfdHV1Xbt27dy5c/v27QOAwcHBLVu2DAwMAEBZWdmFCxd+/vnn7u7uqqqq8vLyMGaNEEIIhUiw0/YqKysdDsf69esTExNPnDixfPlyACCEjI+Pc/NOly5devr06YqKCqfTuXPnzjfffDOMWSOEEEIhEuzjEyH0qI9PfPjhh6WlpbNnzw5rVmhqp06dKiws1Gg0Dw9FYWMwGJYsWcJ9DEVC+eGHH2bOnLl+/XqhExG1X3/91eFwFBUVheRo/4JJs2fPnr17967QWYjdhQsXzOaw7wmOpvbdd981NjYKnYXYXbp0yWQyCZ2F2P3yyy8//fRTqI72LyiECCGEUPhgIUQIISRqWAgRQgiJmgCTZaqrq6urqx/0uH2gnp4elUqFj68Ky263JyQkKBS4OraQBgYGFAoFv4QTEsTQ0JBUKuVWF0FCcbvdDMMkJSU9NLK4uLiqqmrqGAEKIcuyVqs1+MI2Pj4eE4Mr0QoMfwvTgc/nk0qlD10vCoUVTdMUReH6PsJiGIYQwq34OjW1Wv3QT/ACFEKEEEJo+sCPlgghhEQNCyFCCCFRw0KIEEJI1LAQIoQQErVgF92ODIvF8vnnn9M0XVxcPOmCisPDw59++mlXV9e6det27NgR+Qyj3tDQkNFoNJvN8fHxO3bsWLp0aWBMTU0NwzBce/Hixbm5uZHNMfrZbDb/nau3bt2akpIyIcbn89XU1Ny+fTszM3PPnj04lTTkWltbf/vtN/8evV4/YXe5c+fOcduvAoBarQ7V0peoo6OjqanJ6XTu3r3b/0mV33//3WAwKBSKkpKStLS0wB8cGBioqakZGBh4/vnn8/LygjzdNLp42tvbn3nmGUJIXFzcunXrWlombq5OCCkoKLhy5cqCBQveeeedw4cPC5JndDtw4MC3336rUqncbvfq1asnXc3vtddea21ttdlsNpuN24QLhVZTU1N1dbXtH6Ojo4Exer3+yy+/XLhw4bFjx954443IJxn1XC4X/yv45ptv3n///cCHvg4dOmQymbiYvr4+QfKMPv39/dnZ2adOnXr55Zf/+usvvv/69et5eXlz5871eDw5OTk9PT0TftDj8axZs8ZiscybN6+4uNhgMAR7SjJtvP766/v27ePab7311ksvvTQhoL6+PjU11ev1EkIaGhpUKhW3CRQKodHRUb594MABnU4XGBMTE9Pb2xvBpETHYDDk5+dPEWCxWBQKhdPpJITcvXtXLpf39/dHKjsx0ul0lZWVgf0ZGRkNDQ2Rzye6cc8I0jQNALdv3+b7t23bVlVVxbV379793nvvTfjB2trap59+mmVZQsgXX3yxYsWKIM84jb4RXr16taCggGtv3rz56tWrgQEbNmzgPpStWbPG4/H8+eefkc4y2snlcr49Njb2oEVMPvvss6NHj964cSNSeYlOX1/fkSNHampq7HZ74KsmkyknJycxMREANBqNVqudMIiHQmhwcNBoNJaWlk766tdff/3RRx/5D2Wjx/Sgcf5r165t3ryZa09aI7gAiqK4gNbW1qGhoaDO+BjZhlhfXx+/7ppKpbLb7eR/H/a32+18gEQiUSqVvb29kc5SNFpaWs6cOVNRURH4Um5u7v37961Wa2Fh4cGDByOfW9SLj4/PzMx0uVwXL17MyMhoamqaEOB/LQCASqXCayF8zp49m5WVtWTJksCXsrKyAKC3t3fv3r16vT7iqYnI2NiY0+n0rxGBY9H+RWTOnDkymSzI8eppNFlmxowZ3HdhAKBpWiaTcYWdJ5PJ+DkaAODz+Z544omIpiganZ2d27dv//jjjyedsvTjjz9yjbKyspycnP3796tUqsgmGOUKCwsLCwu5dmVl5cGDB7///nv/ALwWIunMmTP79++f9CV+g/GKiopFixbdvHlz9erVEUxNRLjFBf1rROB7XiaT8QEMwzAME+R1MY2+EaampvKfant6elJTUwMD+LujXq/X4XAETqVDj6+zs3PDhg1vv/12WVnZ1JFZWVlyuRy3TQ6rtWvX2my2CZ3+1wIA9PT04LUQJjdv3mxvb3/hhRemDktJSdFqtR0dHZHJSoRmzJihVCr5t/2k73n/IsI11Gp1MAefRoWwqKjo/PnzXPv8+fP8RGSTyTQ4OMgFXL58mRvzNRqNTz31VHp6ulDZRqvu7u68vLxXX331lVde8e9vbm7u7OwEgLGxMb7z0qVLDMMsWLAg0llGO/4fmRBy8eLFZcuWcX9sbGzk/iMoKCgwm8137tzhOt1u93PPPSdUttGttrZ2165ds2bN4nva2tosFgsAeL1elmX5TqvVOunjRihUioqK6urqAIAQUldXx9UIlmUvX748MjLCBRiNRu7yqaury8vLC3arlsed3xM6Dodj0aJFW7Zs2bZtm0aj6e7u5vqTkpKMRiPXLikpSU9P37t3r1Kp/Oqrr4RLNmrt2rVLLpdn/0Ov13P9ubm51dXVhBCDwbB48eIXX3xx69atcXFxJ0+eFDTf6LR9+/aNGzfq9fqVK1dqtVqr1cr1Z2VlHT9+nGsfOnRIo9GUlpYmJyd/8sknwiUbzTweT2Jioslk8u8sKSkpLy8nhDQ2Nmo0Gp1Ot3PnzlmzZr377rsCpRmFNm3atGrVKgBYtmxZdnb28PAwIaS9vT05OVmn023cuDEzM/PevXuEkOHhYQBobW0lhNA0XVBQkJ2dvWfPnjlz5ly/fj3I002v3Sc8Hk99fT3DMPn5+fHx8Vxnc3OzVqvlJsgBQENDQ1dX17PPPqvVaoXLNGrduXOHf0AYAGJjYzMyMgCgra0tISFBrVbTNH3r1i2r1RoXF5eTkxPkyAN6JC6X68aNG0NDQ2q1eu3atfx9DrPZrFQq+Tuyf/zxh8ViWbFiBX4RCZORkZG2trZVq1b5z1fo6OigKCotLY1lWbPZ3NbWJpPJMjMz58+fL2CqUaalpYW/2wcAK1eu5Pa9crlc9fX1sbGxmzZt4jaGY1m2qalp+fLl3F5LDMNcuXLF4XDk5uYmJycHebrpVQgRQgihCJtG9wgRQgihyMNCiBBCSNSwECKEEBI1LIQIIYREDQshQgghUcNCiBBCSNSwECKEEBI1LIQIIYREDQshQgghUcNCiBBCSNSwECKEEBK1/wJSIhsieC/viwAAAABJRU5ErkJggg==",
"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": "1.9793316300831036e-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.00022344388025419073"
},
"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.10.2"
},
"kernelspec": {
"name": "julia-1.10",
"display_name": "Julia 1.10.2",
"language": "julia"
}
},
"nbformat": 4
}