{
"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.597445e+02 1.635877e+02\n",
" * time: 0.0006971359252929688\n",
" 1 1.540207e+02 1.529545e+02\n",
" * time: 0.002201080322265625\n",
" 2 1.168600e+02 1.706955e+02\n",
" * time: 0.003958940505981445\n",
" 3 4.719435e+01 1.074110e+02\n",
" * time: 0.005813121795654297\n",
" 4 1.990498e+01 6.233995e+01\n",
" * time: 0.0075130462646484375\n",
" 5 1.446872e+01 5.470987e+01\n",
" * time: 0.008764028549194336\n",
" 6 1.028559e+01 3.463054e+01\n",
" * time: 0.009814977645874023\n",
" 7 9.239520e+00 1.018030e+01\n",
" * time: 0.010827064514160156\n",
" 8 6.585533e+00 1.100819e+01\n",
" * time: 0.011834144592285156\n",
" 9 3.948380e+00 1.142955e+01\n",
" * time: 0.012946128845214844\n",
" 10 2.122278e+00 2.391355e+00\n",
" * time: 0.01401209831237793\n",
" 11 1.416845e+00 1.553289e+00\n",
" * time: 0.015066146850585938\n",
" 12 1.369622e+00 1.573589e+00\n",
" * time: 0.01585698127746582\n",
" 13 1.250320e+00 7.940685e-01\n",
" * time: 0.01666712760925293\n",
" 14 1.184392e+00 8.198353e-01\n",
" * time: 0.017460107803344727\n",
" 15 1.158106e+00 4.092584e-01\n",
" * time: 0.018252134323120117\n",
" 16 1.148319e+00 1.811436e-01\n",
" * time: 0.019036054611206055\n",
" 17 1.145951e+00 1.127054e-01\n",
" * time: 0.019865989685058594\n",
" 18 1.145023e+00 7.041936e-02\n",
" * time: 0.02069711685180664\n",
" 19 1.144439e+00 4.518698e-02\n",
" * time: 0.021502017974853516\n",
" 20 1.144231e+00 6.553943e-02\n",
" * time: 0.022083044052124023\n",
" 21 1.144195e+00 6.369244e-02\n",
" * time: 0.022686004638671875\n",
" 22 1.144130e+00 3.549207e-02\n",
" * time: 0.02433609962463379\n",
" 23 1.144091e+00 1.907598e-02\n",
" * time: 0.02518606185913086\n",
" 24 1.144058e+00 1.124082e-02\n",
" * time: 0.0260009765625\n",
" 25 1.144040e+00 7.059753e-03\n",
" * time: 0.026803970336914062\n",
" 26 1.144040e+00 4.578568e-03\n",
" * time: 0.027383089065551758\n",
" 27 1.144038e+00 4.604402e-03\n",
" * time: 0.028222084045410156\n",
" 28 1.144037e+00 2.416074e-03\n",
" * time: 0.029050111770629883\n",
" 29 1.144037e+00 1.055905e-03\n",
" * time: 0.02986311912536621\n",
" 30 1.144037e+00 6.336158e-04\n",
" * time: 0.03067493438720703\n",
" 31 1.144037e+00 7.362609e-04\n",
" * time: 0.03147411346435547\n",
" 32 1.144037e+00 2.860822e-04\n",
" * time: 0.03205394744873047\n",
" 33 1.144037e+00 1.964097e-04\n",
" * time: 0.03290200233459473\n",
" 34 1.144037e+00 1.304958e-04\n",
" * time: 0.03381705284118652\n",
" 35 1.144037e+00 1.078917e-04\n",
" * time: 0.0346379280090332\n",
" 36 1.144037e+00 4.127986e-05\n",
" * time: 0.03550004959106445\n",
" 37 1.144037e+00 6.999419e-05\n",
" * time: 0.036366939544677734\n",
" 38 1.144037e+00 3.934010e-05\n",
" * time: 0.03724193572998047\n",
" 39 1.144037e+00 2.635852e-05\n",
" * time: 0.03808093070983887\n",
" 40 1.144037e+00 1.961053e-05\n",
" * time: 0.03891801834106445\n",
" 41 1.144037e+00 1.166020e-05\n",
" * time: 0.03971505165100098\n",
" 42 1.144037e+00 7.174118e-06\n",
" * time: 0.04052400588989258\n",
" 43 1.144037e+00 6.365032e-06\n",
" * time: 0.0413510799407959\n",
" 44 1.144037e+00 4.842403e-06\n",
" * time: 0.04215693473815918\n",
" 45 1.144037e+00 3.860060e-06\n",
" * time: 0.04297900199890137\n",
" 46 1.144037e+00 3.154550e-06\n",
" * time: 0.044181108474731445\n",
" 47 1.144037e+00 2.319926e-06\n",
" * time: 0.0450289249420166\n",
" 48 1.144037e+00 1.045393e-06\n",
" * time: 0.04588508605957031\n",
" 49 1.144037e+00 6.560508e-07\n",
" * time: 0.04670095443725586\n",
" 50 1.144037e+00 2.975209e-07\n",
" * time: 0.04751706123352051\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": [
"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 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": "8.338746851893258e-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+gvaeTAAAgAElEQVR4nOzdd3xUVfo/8OfemUnvddJIoYV0EhJKEnoVRViVVbHsuq5ucdW1rLq41u+6a19Xxe66ooiIoiJNOoQSSiCdBBJC6mTSe5mZe35/jL8sUidw7tTP+68kXJ5788q9+eScc885AmOMAAAAHJVo6QsAAACwJAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAsEYWtr6+nTp00/3mAwyHcxYCJJkix9CUCSJGFNRItjjOGnYHGMMY6/lCwQhJs3b37yySdNP763t1e+iwET9fb24uG3uIGBAfxFYnGDg4N6vd7SV+Ho9Hr94OAgr2roGgUAAIeGIAQAAIeGIAQAAIeGIAQAAIeGIAQAAIeGIAQAAIeGIAQAAIeGIAQAAIemtPQFAMBlMKKCVlbewcpbRXcnNsqbJfpRpIdg6esCsBMIQgDr1TpAH5ZJH5dJjCjBVwh3JX0/ba03HG1msT7C3WPFm2NEJbp1AK4OghDASq2plB48aJgfLn4yTTEpSCCivr5BJyelQqHQSfRDtfRmsfRaofTRVMV4f7QObUBLS8t9992HxZOv2KOPPpqeni5HZQQhgNXp0dOduw1l7WzdbOXEoAuEnEqkJVHikijx05PSgs36++MVf01Bw9DaNTQ07Nu379VXX7X0hdikd999t6CgAEEI4BBaBmjhFn28j3B0idLpcul2x2hxbrh47RZ9dTd7O1OhQMvQunl7e990002WvgqbtGXLFvmK469IACtS38uy1utnhggfTlVcNgWN1K60c6Gyoov9codBj60pAIYPQQhgLbp1dO0Ww+2jxBfSh9e081TRhnnKfj37436MPwEMG4IQwCpIjG7bZUgNEK5stM9JpDWzlMdb2EsFaBUCDA+CEMAqPJJr6NOzdzMVV1zBTUnfzFa8XSJ9dwZZCDAMCEIAy9tUw76pYl/OUl7lpMAwd2HNTMW9OYa6Hsbp0gDsH4IQwMKa++m3OYb/TFP4OHGoNjFI+EOc4td7DEhCABMhCAEs7Dd7DbeNEmaEcJv68NdksUdP/ypCBynI6Le//e1///tf48ednZ0ZGRm9vb0XO9hgMGRmZtbX15vr6oYHQQhgSV9WSme62HNpVz40eD6lSCunK144bqjqQrMQ5NLZ2dnX12f8+JVXXrn22mvd3NwudrBCobj99tufe+45c13d8CAIASymW0eP5Er/nmLqlEHTxXgKDyYo/nwQjUK4jNWrV9fW1r733nvPP/88EQ0MDHz66adPP/306tWrJemn+ycvL++VV17529/+dvYXh+h0uvfff/+OO+4gIq1Wu3r1aiI6evTogQMHiOibb76pra0loptvvnnNmjUdHR3m/O5MhCAEsJjnjhlmhwlT1bKsB/NokljazjbUoFEIl7J8+fJrr722pqbGx8dnYGAgMzNz165dkZGRq1atuu2224zHrFixQhTF6OjoDz/88O677z6nwtGjR93c3KKiooioqqrqySefJKL169evWbOGiF544YWSkhIi8vHxGTt27O7du8353ZkIS6wBWEZpO/ukXCq8QSVTfSeR3pisuG+/YVao0oVnzyvw9HaJ9K0Zp7tsmHeBdft+9atfPfjgg0S0YsWK0NDQjz/+mIhuv/32mJiY8vLyMWPGfPjhh8Yjb7nlFn9//xUrVri4uAz998LCwtGjR5ty9jFjxhQUFCxatIjL98IRghDAMh47JC1PUQS7yniKeeFCgq+wokR6KBF9P1Zqdpgw1tt8f6coL9T7MLSS9dGjR4uLi+fMmWP8tKenp6ysbMyYMW+88caKFSsEQfDw8NDpdPX19TExMUP/vbu7293d3ZSze3h4dHV1Xe33IAMEIYAFHGpihW3sq1my/wb8Z4Y4/Qf93bGil1wtT7gqY72Fsd4WvgZnZ2fjB66urtdff/3f/va3oX9yd3cvKCh46aWX8vPzAwICDAaDu7v7OTtJBQQEtLa2Gj92cXHp7+8/+1/7+vqGmo+tra3R0dEyfidXCn8nAljAE4cNT48XneVvCYz1FuaFi5hKAaaYP3/+t99+S0S+vr6+vr6CIAiC0NjY6OPj4+/vT0SrV68eGBg453+lp6cXFhYyxoho1KhRXV1deXl5xn86ceJEdXX1uHHjjJ8eP348IyPDfN+PydAiBDC3H+tYbQ/dNspMf4Y+myamrdP/LlYMkrMbFuzAtddee/To0bi4uIyMjM7OzoqKioKCguzsbDc3t4kTJ/r6+rq5ufn5+Z3zv2JjY8PCwg4fPpyRkeHm5vbhhx9ed911oigaDIY1a9a8+eabgYGBRHTmzJnOzs7Jkydb4ju7DAQhgLn97Yjh7xPEq1xNzXSRHsLNI8VXCg0vZeCdGThXXl6eh4fH0KdPP/30gw8+WFZW5uHhMXbsWIVCQUT79+8vKChwdXWNi4trb2/38vIiog8++MDJ6afFkB566KH333/f2Nq76aablixZ8sADD3R1dX300Ucq1U+d8u+///79998/9KlVQdcogFltq2M9eroh2qyP3uPJ4sdlUtu5fVoA5O3tbUy7s7+SkZERFxc39HWVSpWWlhYXF0dEPj4+oigSkZeX19Dg3x133OHi4jK0soxSqQwMDPT39x+KPYPB0N7efv/995vnmxoutAgBzOqf+YbHkkUz7yQf7i4sihRXlErLr2iPJ4BLE0XxrbfeOvsrWVlZQ+vOEJFCoXj77bfNfl2mQhACmM/xFlbeQb+MsUAaPZYsTv9B/1CC6IqHHuQ3e/ZsS1/CMODPQwDz+ftx6eFEkfuCaqYY6y1MDBI/OYnXR+Fn4uPjjx07JkdlnU43ZcqUxsbGSxxz4403Hjp0SI6zDwuCEMBMTnWy3Q3Sb8Za7KF7LFl8tVDC/kxwtmeeeWbEiBFyVP7oo48SExODg4Mvcczvf//7J554Qo6zDwuCEMBM3iqRfhsreljupbnJQUKQC/1QjUYh/E93d7dxgvymTZtyc3NXrVr14IMPrl27loh27dr10EMPvffee0MLbW/duvXpp59++OGHv/zyS+PEQSJijK1cufLPf/7z6tWrc3Jytm/fbvz6u+++a1yttKen56WXXiKiwsLCdevWEdG6devy8/OJaMaMGeXl5SdOnDD3t/1zCEIAc+jW0cqT0j2xFn7i/hgnvl2CIIT/ef7552tqaojo66+/XrZsWX5+fnx8/H333Xfvvfd+/PHHSUlJ77zzzosvvmg8+PPPP4+MjExJSXnzzTcfe+wx4xcfeuihFStWTJgw4eDBg7feeusPP/xARPX19eXl5ZMmTSKi7u7uv/71r0R09OjRlStXEtFnn32Wm5tLRKIoZmdnb9q0yRLf+v9g3BzAHD47JU0PESM9zPy66LluihEfyTWcaGexPha+EjDq2rG275j5NmQIfPB1QXHRX/uTJ082Zl5NTc0333xTWFgoCIKPj8/LL79s7MD85JNPiKinpyc1NXXatGkvvfRSR0fHO++8U1FRERYWtmzZsuLiYmOpwsLCqKgoU2YNjhkzpqioiMt3d8UQhADm8E6p9Poky89ndxLp7ljxnVLpjcmWvxggIrfUac6jk812OkG81M89ISHB+EFwcHB8fLwgCMaPW1paiEiSpEcffXTt2rV+fn6iKLa0tPT391dWVvr6+oaFhRn/Y2pq6uDgIBH19PS4upq0lJGbm1t3d/fVfFNXD0EIILvdDUwn0YxQq2iE/X6cmPi1/vkJCizDbQ0UPoEKn0BLX8VPjDPljc6ZZU9EGzdu3LFjR3l5ubOzs0ajCQkJkSTJx8ens7PTYDAYj29razPuRBEUFGSMTyJydXWVJOnsRUp7e3uHtrNvbm6+9As1ZoAxQgDZrSiV/hhn7kn0FxPqJswIFVedwkghDE9ra6uPj49xq4r333/f+MXIyMiRI0e+8847RFRZWfnNN98Yvz5+/HiNRtPW1kZEXl5eo0aNMi7nTUR1dXWHDh0aP3688dOCggKLr8SNIASQV8sA/VgrLTPXEtumuCdW/LgcQQjDs2jRIq1Wm5GRMWnSJI1GY/yiKIqrVq366KOPQkNDb7vttgULFhh7RN3d3RcsWDD0Fswnn3zy1FNPLV++fOfOnenp6U899VR8fDwRdXZ2HjlyZOHChZb6pozQNQogr5UnpWtHiD5Olr6Os8wOFe7po/xWluxnJc1UsJji4mJjI+/f//73UNfo3Xff/atf/cr4cXp6unHOu4+PT35+fllZma+vb2ho6N///ndj5iUkJBw7dkySJFEU586dO2vWLON/fPjhh5cvX37rrbcS0ZQpU8rKyl5//fVNmzZt3LhRqfwpej777LNbb73V19fXrN/zeRCEAPL6pFx63creTBEFunO08Em5Vby/A5Y19ErL0KAdETk7Ow/t1qtUKo3bTRg/NrbkiGgovb744ouioqLw8PCcnJwzZ84sXbrU+PWsrKwFCxZoNBq1Wj30Xzw8PIZSkIhaW1ufeuopub43kyEIAWR0pJl16mia2uoaXr8eI2Z8p/9nusIMmwODfZs2bVpXV1dzc/OCBQs++OCDswP18ccfP/vIxMTEc97BefLJJ810lZeEIASQ0X/KpV+PEa3lPZmzRHkKiX7C+mrpRvNuCAX2JzQ09J577jHlyLS0tLS0NLmv5wrgGQCQS7+BvqyQ7hxtfTFIRES/GYtXZgCIEIQA8llfLY0PEEZYejWZi1kSKR7UMk3f5Y8EsG8IQgC5rDrFlo203kfMVUnXjRC/rECj0NHl5+c/8MADDz744J49eyx9LZZhvU8pgE1rG6CdDdKSKKt+xJaNElchCB3bkSNHFixYkJiYmJycvHTp0q1bt1r6iiwAL8sAyGLtaWlumOhtTdMHzzcrVLizm5V3sDHeVtp/a/dWl67bUbXXbKd7e95LKvFnv/Zfeumlxx577O677yai/v7+l19+ec6cOWa7HiuBIASQxecV0oMJVt0cJCKFQEtjxC8q2NOpCELLmBmZPT440WynUwjnTpcpKSn505/+ZPw4LS3tH//4h9kuxnogCAH4q+9lRa1sQbi1ByERLRslLttpeDrVBi7VLgW5BQS5BVjwAvr6+rq6uowfd3Z2Ds2ddyi4+wH4W1XBlkSJNjFXPSNQEAQ63MQsfSFgMR9//LEkSYyxjz76yAH7RQktQgA5fFkhvZhhCzFIRES/jBG+Oi2lB9rMBQNffn5+ycnJer3e39//rbfesvTlWABahACcVXWxM91sqvUtq3YxN0aLX51maBI6rHvuuWfv3r3btm3Lycnx9/e39OVYAIIQgLOvTrMlUaLSdp6tZD/BWaS8ZkSh4/Lx8RnaZd4B2c7DCmAj1p62vQU8fxElrD2NCYWO6C9/+Ut4eLilr8LCbOxxBbByNT2ssovNCLGZflGjG6PFNZVoETqie++9d2ibJIeFIATg6atKtjjSlvpFjVIDBEGg4y3IQnBEtva8Ali3r6ukG2ytX9TohijhK/SOgkOyyScWwDo19NKJdjYr1Mb6RY1uiBa/rUKLEBwRghCAm++rpWsiRJVtPlXpgUKnjso7kIXgcDChHoCb785Id42xzRgkEoiuGyF8d4Y9mmSTLVrr5+3t3dvbO3LkSEtfiE1qamqaNm2aTMURhAB8dOtofyP7cqatBiERXR8pPn/M8GiSDX8L1iwiImL//v09PT2WvhBbNWLECJkqIwgB+NhYI2UGC54qS1/HVZgZKty6k2n6SO1q6UuxU8HBwZa+BLgA/OkHwMd3Z9j1kbb9QKlEmhsubqjGu6PgWGz7uQWwEjqJttRK142w+Qfq+hHCd2fwvgw4Fpt/bgGswR4NG+0thLhZ+jqu2oIIcY9G6tFb+joAzAhBCMDB+mp7aA4SkbcTTQgQttehdxQciD08ugAWt7GGLRxhJ7MOFo4QN9SgdxQcCIIQ4GqVdbAeHSX52UkQLhohrK+WkITgOIYxfaKysvLHH3/09fVdtGiRq+uFX6/Oy8s7ePCgp6fnjBkzsLUHOIgN1ezaEYKdxCDRSC/BQyXkt7AUf7v5ngAuxdQWYU5OTlpaWnFx8UcffZSdnT0wMHD+MY888siiRYsKCgq2bt36/vvvc71OAOu1oUZaGGFXmbEwQkDvKDgOU1uEzz777PLlyx955BGDwTBhwoS1a9cuW7bs7AO2bNny2WefFRUVBQQEyHCdAFaqU0eHm9jMULsaZVgYIT511LA8xa6+KYCLMelGHxgY2L59++LFi4lIoVBcd911mzZtOueYr7766o477mhpafn++++rqqq4XyiAddpaK2UGCx62vKDM+aaGCCXtTNtn6esAMAuTWoQajYYxFhoaavw0NDQ0JyfnnGMqKyuLior27NkTGxt71113vfbaa3fcccfFCp46deqFF144+yuRkZFLly694ME6nU6n05lynSAf409BsJ+BMG42VNPcUDLPLWr8EUiS7HMbBKIZatpYrVsWI/epbI9OpxNFtJUtTKfTGQwGhUJx2SOVSuVlf3GZFISMMSIaqiWK4vmP4uDg4ODg4JEjR0RR3Lhx46233rps2bKLXeXg4GBbW9vZX/Hz87vY4y1JkhmefLg0408BQXgORrSlTnw03kx3qDl/CvPDhC21wi1RePTOZfxh45eSZUn/H5dqJgWhWq0WBKGxsTEqKoqIGhoahlqHQ0JDQ/39/Y1/KGVmZnZ0dGg0mrCwsAsWjIuLe/nll028xMHBQWdnZxMPBpnodDpnZ2cE4TkKWpm7yjAuwEz3pyRJTk5OpvwVfPUWRrG/HdOrnJxF/MzPI4qiSmVfveG2RhRFg8HAKxpMauC7uLhkZmZu3LiRiBhjmzdvnjVrFhHp9frq6mpjJs+dO7e0tNR4fElJiaura1BQEJdLBLBam2vZ/HD7DIpwdyHIVchrwbujYP9MfWv0ySefvOWWW7RabWlpaXt7+80330xEp0+fHjNmTFNTU0BAwK233vraa6/dfvvtiYmJ77777jPPPIO/mMDubamVHko0R/vMIuaHC5tr2IQA+0x6gCGmDvnOmzdv586dKpVqxowZBw4ccHd3JyK1Wv355597enoSkZubW25ubnZ2tiRJK1eu/Mtf/iLjVQNYgR49HWli09R2mxPzwsUtWHQUHMAwVpZJTk5OTk4++yuenp633nrr2Z/ec8893C4NwLrtrGfpgfY2ceJsU9VCUStrHyQfJ0tfCoCc8BIwwBXaUivNC7fnJ8hZQZODhR31aBSCnbPnxxhAVnb8psyQeWHillq8LwN2DkEIcCUqOlmvnhLsZceJi5kfIWxGEIK9QxACXIkf69icMPufVjnWWxAFKutAFoI9QxACXIltdWxOmN3nIBHRzBBhax2CEOwZghBg2AyMdjVIdrbjxMXMDhO2IQjBrjnEkwzA15EmFuYuhLhZ+jrMYk6YuKtB0uHVUbBfCEKAYdvqMP2iRBTgQtGewuEmNArBbiEIAYZtW700J8yBnp05YRgmBHvmQA8zABc9esprZtn2u7La+WaHidswrR7sF4IQYHh2N7C0AMF9GKsT2rypaiG/hXVie2ywUwhCgOHZVifNdqR+USJyUdDEIGF3AxqFYJ8c6c9aAB521LP3sswahIOGwe1Ve0paysuaT7k5uY71G5USnDg5bII5r2FmqLiznl03wpznBDATBCHAMLQMUFU3SzPjFn17aw68nfdxtHfkhJCUqSGTDKJ0qr3qw/yVq0u++dOE347yjTbPZcwMEe7JQYsQ7BOCEGAYttdJU9Wi0iwNQolJLx78d1lrxV8m/ilVnUREfX19Tk5Ok8Im3Br3ix9O/fjIjqd/k7zsulHzzHAxEwKFmh7W2EfBrmY4G4BZIQgBhmFnA5sZao7moMSkfxz4V/tA53vzX3VWnLsfoCiIi0bPzwhN/fO2J/v1AzfFLpL7ehQCZQWLuxukpTGONT4KjgD3NMAwbK83RxBKjD2998UeXe8LU5efn4JD1O5B/5r992/LN649sV7uSyKimaHCjnrMJgQ7hCAEMFVdD2sbYAm+sgfh6tJv2vo7nst+QqVQXfrIYPfA12f/3+clawubSuS+qpmhwo4GBCHYIQQhgKm21bOZoaIocw6WNJevKf3ub5kPKUWFKccHuQU8Pun+53Je6RjolPXCEv2EzkF2phtZCPYGQQhgqh3y94t2DnY9k/PSY5PuD3YPMv1/TQxNmxmZ/Y8D/5LvwohIIJoeIu5CoxDsDoIQwFQ769nMEHmD8MPjn2WGp1/BHMHfptzR3Nu648xeOa5qyMxQYTsWHQW7gyAEMMmpTsaIRnvLGISn2ir31Oz/deKtV/B/laLiwfTfvZP3n359P/cLGzIzVECLEOwPghDAJDvr2QyZm4NvHvnw7uTbvZw9r+y/JwTGJgclrCr5mu9VnW2UlyARVXQiC8GuIAgBTLJbw6bJGYTbq/b06fuvGTnnaor8bvyd35Zvauhu5HVV55segkYh2BsEIYBJdjWwGbK9KWNghg/zP7sv7TeicFWnCHDzv2HstZ8UfsHrws6HIAT7gyAEuLyTHUwkivGUKwi3nd6t9ghKCoq/+lI3xi46UHekrqvh6ktd0PQQTKsHe4MgBLg8WZuDEmNflHxzR8IvuVRzV7ldP3rB6tJ1XKqdb5SXoBTpFIYJwY4gCAEub1eDjAOEu6pzXFWu44MTeRW8adyiXWf2NfZoeRU8xzQ1ekfBriAIAS5vt0auV0YZsc+KvvpV4s0ca3o5eS4cNWd16bcca55tGoYJwb4gCAEuo7yDiUTR8gwQHmk4LgjixNA0vmWXxl6/7fTu7sEevmWNZoQKO+qxNyHYDwQhwGXslnOAcF35hl+MXci9rJ+r78TQtE2V27lXJqIYT0ElCuUdaBSCnUAQAlzGbg2bLk+/aGNPU1HTiVmR2XIUXzL2mnXlGyQmS1xNUwt7NAhCsBMIQoDL2KthU9WyBOF3JzfNi5nponSRo3h8QKyHyv2o5rgcxaeGCHswTAj2AkEIcCmVXczAaKQX/yDUGXSbKrcvGjWPe+Uh149ZsK58gxyV8eIo2BMEIcCl7G5g0+RpDu6s3jfKJzrCK0yO4kazIqcWNZ2QYx7FaG/BwKiqC1kI9gBBCHApe2RbYnRT5bbrRsvYHCQiF6XzrKipmyt3ylF8agiGCcFOIAgBLmW3PFPpG3uaKtqqJocOe9/B4ZofPXPL6R2M+CfWVLWwG72jYBcQhAAXVdvDevRsjAx7EP54eufMyGyVQsW98jnG+o9yUjgVN5VxrzwNLUKwFwhCgIva1cCmqUU5OkZ/PL1rXswMGQpfwNzo6ZtlmFA4zkfo1LG6HmQh2DwEIcBFydQvWtx8ghEb5z+Ge+ULmhs9Y3f1/gHDIN+yAlFWsIhGIdgBBCHARck0g3BL5c75MbO4l72YAFe/sf6j9tce4l55WgiGCcEeIAgBLqyxj5r6Wbwv5yDUS4Zd1fvmRk/nW/bS5sXM+PH0Lu5ls9VCTiOCEGweghDgwvZopCy1yH2E8Kjm+Aiv8CC3AM51LykrfFK+tqhrsJtv2WQ/ob6Xafv4VgUwNwQhwIXt1bBsGfpFd1bvmxGZyb3spbkqXVLVSQfqDvMtKwo0OUjYr8VOFGDbEIQAFybHAKFeMuyvPTQtYgrfsqaYPiJz55l93Mtmq8W9eF8GbByCEOACOgapspON9+cchMZ+0QA3f75lTZEZPrGgqbhH18u3bLYaq2+DzUMQAlxATiObGCSoeD8fFukXNXJVuowPTtxXm8u3bHqgUN7BunR8qwKYFYIQ4AL2aqRsNeenw4L9okZy9I46iZQaIOzHu6NgyxCEABewp4H/mzJ5jfkjvMIs0i9qNCUsI19bxL13dKpa2KvB+zJgwxCEAOfq01NhG8sI5ByEOTW52RGT+dYcFjeVa1JQ3KH6PL5l8b4M2DoEIcC5DjaxJD/BTcmzJiO2v+7wlLB0nkWHLzN84v46zkvMTA4W8lpYv4FvVQDzQRACnEuOftHylgp3laus2/CaYkpYxoG6I3qJZ2q5K2mcj3C4CY1CsFUIQoBz5TRKWcGcH419dbmZ4RP51rwC/q6+4Z6hhU0lfMtmBQvoHQXbhSAE+Bm9RIe0bEow/wHCzPAMvjWvTGZ4BvdJFNlqIacR78uArUIQAvzMsRYW6Sn4OfOs2dijbe1vG+c/lmfRK5UZPnFvzUG+NbPV4v5GZkCbEGwTghDgZ3IaWRb35mBt7pSwDFGQY4vfYYvxiRQF8XT7GY41A1woxE0obEUSgk1CEAL8jBxrbe+vOzzFOvpFjaaEp+/jvQA3hgnBdiEIAf6HEe1rlLK4BmG/vr+0uTw1OIljzas0MTQtt/4o35rYmxBsF4IQ4H/KO5irQohw5xmEeY0Fsf6j3VSuHGtepZSghFNtlZ2DXRxrZquFPQ14XwZsEoIQ4H/k6Bc9WHd0Ymga35pXyUnhlBQYn6cp4Fgz2lNQikJFJxqFYHsQhAD/s6+R8e0XJaLDDcesLQhJnt7RLDWGCcEmIQgB/mevhvMro9WdtXpJH+UdwbEmF5PDJuTWH2XEM7cyg4V9GCYEG4QgBPiJpo/aBtg4H55BeLD+6KSwCRwL8hLiEeyqdK1oq+JYE+/LgI1CEAL8ZK9GylKLItee0UP1eRkhqTwr8sO9dzTRV9D0Mm0fx5IA5oAgBPhJjoZlcu0X7dcPlDSXpamTOdbkaFJY2kGuQSgKNDlY2Ie11sDWIAgBfsJ9TZl8bdFo3xirmjhxtuSghFNtlXz36c0MFjFMCDYHQQhARNSto/IOlhbAMwiPNByfEJLCsSBfzgqncf5j8rVFHGtimBBsEYIQgIhov5alBQjOCp41D2usOgiJaEJIyuGG4xwLTgwUittYr55jSQDZIQgBiIhyNBLfftHWvrbm3paxfqM41uRugjrlCNcgdFZQkp+Qi016waYgCAGIfqyKdWAAACAASURBVHpThufjcERzPFWdJApW/YiN8o3pHOjS9jZzrInVt8HmWPVTCmAeOomONHPejPeIJn+C2qr7RYlIFIRUddJRTT7HmllqIUeDF0fBliAIAehoMxvpJXg78ayZpymw8gFCozR1Mt/e0axgMVfL9IhCsB0IQgDax3vixOmOaoUghnqoOdaUSXpIylHNcY5rrfk6U4SHUIBNesF2IAgBaF8j56n0Rxvy00PGcywon2D3IA8nj4q20xxrZgVjEgXYEgQhODrjZryZXDedOKrJt9oFZc6Xpk7muyVTphqrb4MtQRCCoyvvYG5KnpvxSkwqbCpJCU7gVVBu44MTjzUWciyYGYxNesGWIAjB0eXw3nqprPVUkFuAr4sPx5qySg1OytcWG5iBV8EYT0EpCpVdaBSCbUAQgqPjPkCYpylIVSdxLCg3L2dPtXtQWUsFx5pTgoUczCYEG4EgBEeXw3tX+mONheODbSkIiShVnXSskeswITbpBduBIASH1thHzf0sjt9mvDpJX9JclhQUx6ugeYwPTuI7TJiFFiHYjmEE4erVq+++++7ly5drNJpLHLZmzZo333zzqi8MwBxyNFJmsMBxM97S5rIIrzBPJw9uFc0iJTihuPmEzqDjVTDZX6jvZS0DvOoByMjUIPz3v/+9fPnyadOmtbS0ZGVlDQxc+AbPzc297777nn/+eX5XCCCjfY1sCtclRvMaC8YHJ3IsaB7uKrcRXuGlLeW8CioEyggU9mOTXrAFJv0KMBgMr7766ttvv3377be/++67np6eX3/99fmHDQwM/OEPf3j22Wd5XySAXLhvxntMU5hqawOERuODE/O4DhNOwSa9YCNMCsLa2tqampoZM2YYP50+ffq+ffvOP+z5559fvHhxXJyNjY6Aw+rRU2k7Sw/kFoSDhsGy1lOJtjZAaJSqTjrWyHuTXgwTgi1QmnJQQ0ODh4eHs7Oz8dPAwMAjR46cc0x+fv73339/+PDhgwcPXrbgkSNHbrjhhrO/EhcX98QTT1zw4L6+PoWC636pMHy9vb2CIAgCz8aTxe1uFJN8RWmgt5dTwYLmkiivEWxQ6h3kVfJn+vr69Hq9TI/DSPeospaT7V3tTgo+q48neVB+q6q1q9fFvh7fgYEBURRVKpWlL8Sh6XQ6g8EgSZfve3dxcRHFyzT5TApCFxeXwcHBoU8HBwddXV3PPkCv1999993vvvvuUFheWlhY2C9/+cuzvxIYGOji4nLBg3U63cX+CcxGr9e7uLjYWRAeamXZasbx7ippK09VJ8p3uzLGnJycZApCF3KJ9hlxuqcmOSieU0GK85EKu5yzuc5OsThBEBCEFqdQKAwGgynP2mVTkEwMwrCwsIGBgaampsDAQCKqqakJCws7+4DKysr8/Pzbb7+diPr7+1tbW0eOHLlt27bo6OgLFgwJCVm6dKkppyYiURRN+U5AVsafgp0F4T6t/oF4hcjvndH8puJl8TfId7uK/59M9ZODEgqaiserub3sk6Vm+5uEaaF29fzK/VMAU4iiyBjj9VMwqUpgYGB2dvbnn39ORO3t7Rs3blyyZAkRtba2rlu3joiio6NPnDixdevWrVu3vvrqq97e3lu3bg0PD+dyiQBy0Et0SMtzM16dQVfWcjIhYByvguaXEpxwXMtzmDAzGJv0gg0wqUVIRP/85z8XL168c+fOkpKSuXPnTp48mYhOnDjxi1/8gjGmUqliYmKMR9bU1CgUiqFPAazT8VY2wkPwM6kv3yQlLeXRPiPcVK6XP9RaJQXGP5vzss6gUyn49Ptlq8W79hgMjBR21ZUA9sbUIJwyZcqJEycOHjwYEhIyfvxPG62lpqaWlZWdc+TEiRMPHz7M8xoBZJCj4byy2vHGwuQgm9lx4oLcVK4RXmEnWk8mBvJ58TXQhYJdheI2luSHJATrNYwOVj8/v2uuuWYoBYnIxcVlzJgx5xzm4uISGRnJ5+oAZMN9re3j2iIb2nrpYsYHJR7nOokCa62B9cN4LziofY0Sx6n0Okl/ooVbQ8qCknkPE2apsVs9WDsEITiiU51MFIQoT25BeKKlPMIrzF3lxqugpSQHxZc2l+slbnsTZgULe9EiBOuGIARHtFfDpnIdIMzXlqTY+AChkbvKLcwzpKz1JK+Co70FA2NnupGFYL0QhOCI9vFeYjRfW2Trb8oMSQ6KL9CWcCyYGSximBCsGYIQHNFerq+MSkwqaS5LDLThGYRnSwqKL9AWcyyYFYxhQrBqCEJwOE391NjH4n25BeHJtsogtwAvZ09eBS0rKSiuoKlEYtwmwmepMUwIVg1BCA5nr0bKDBY4TvHO1xYncVqf0xr4OHv7u/pVtp/hVTDFX6jpxia9YL0QhOBwcjQsS83zzi/QlthTEBJRclB8Pr/eUYVAE4OwSS9YLwQhOJycRpbN700ZRqyoqYTXjg1WIikojvMwoRrvy4D1QhCCYzFuxjuB32a8ZzpqXZWuAa5+vApag5SghHxtMSNu0ZWNYUKwYghCcCwHtSzFX+C4VWyBtjjZ9ldWO0egW4Czwqm2s55XwYmBQmEb69PzqgfAE4IQHMtejcR3Kn2BtiTZ9ldWO19ycALHYUI3JSX6Coea0CgEa4QgBMeyV8Oyub4pk68tSgyywyBMChzHMQgJvaNgxRCE4EB0Eh1uYpODuLUINT1anaQP9wzlVdB6JAbGFTbxXF8mWy3uxSa9YJUQhOBA8prZKC/B24lbwQKtvb0vOmSEd3i/vl/b28yrYGawcFDL9IhCsD4IQnAgezQsm/MAYXGSPfaLEpFAQkLguMKmUl4FfZ0p0lM43oreUbA6CEJwIHt5B2FhU4kd7EF4MYmBcYVcV9/GMCFYJwQhOApGtL9Rygzmds93DHRqe5tH+kbxKmhtkoJ4DxNib0KwSghCcBTFbczfRQjht3VuYVNpQsA4hcBvTqKVGeM3sqG7sXuwh1fBaSFijkZCEoK1QRCCo+DfL6otSQyyk62XLkghKMb6jypq5jZMGOJGniqhrB1RCNYFQQiOYq+G82a8BU12tenEBSUFxvMdJpwaIuxB7yhYGQQhOIocDZsWwi0I+/UDle3VsX6jeRW0TklBcflcgzALw4RgfRCE4BAqOhkjivbkFoSlLeUjfaJclM68ClqnuICxp9oqdQYdr4JT1cKuBgQhWBcEITiEPVybg2TXMwjP5qp0GeEdfqL1JK+Co70FiVFVF7IQrAiCEBwC96n0hU2ljhCERJQUGFeAYUKwawhCcAh7GhjHTSckJpW2lCcE2PMro0MSA+OK+K0vQ5hNCNYHQQj2r66HdelYrA+3IDzVdjrQ1d/L2ZNXQWuWFBRX2FQqMW7RhRYhWBsEIdi/3Ro2NUTk2DFa2FRil1svXZCvi4+3s9eZjmpeBeN9hdYBVt+LLARrgSAE+7dXw7NflIgKtCWJgQ7RL2qUGBRXwG+tNYEoK1jMQaMQrAaCEOzfbq4DhERU1FRqx2ttny8xcFyhluswoRq9o2BFEIRg55r6SdPHkvy4BWF9t4YEIcQjmFdB65cUyLNFSETTQoTdmE0IVgNBCHZud4OUFSxyHCEs0JYkO1JzkIgivMJ0ko7jJr0p/kJdL2vu51UP4KogCMHO7W7gPJW+sKkkwcGCkIgSAmI5LjqqEGhykLBHg+3qwSogCMHO7eIdhAXakiS73nTighIDx/HuHRXROwpWAkEI9qxlgGp72Hh/bkHYMdDZ0tca4xPFq6CtSAzivFv9dAwTgtVAEII9290gZQYLCq4DhPGBsaLgcA/OGL+Rmh5t52AXr4Kp/kJVN2sZ4FUP4Mo53PMMDmV3A5sWwvMmL2oqTXK8AUL6/5v0ljSX8SqoFGlykLAXw4RgBRCEYM/4DxA60poy50gKjC/kuugohgnBSiAIwW61DlBVF0vlN0A4YBisbD9j95vxXkxi4Di+21BgmBCsBIIQ7NYejZSpFpT87vGS5jJH2Iz3YuIDY8tbKwYNg7wKTggQKrtYK4YJwdIQhGC3djWwaWqed3iBtsRB9iC8IFelS6RXeFlrBa+CSpEmYZgQrACCEOzWzno2I5TzEqMJjrTW9vkSg8YVcp1NOD1E3IXeUbA0BCHYp5YBOtPNc4BQYlJx8wmH2nTifImBnGcTzggRdtQjCMHCEIRgn3bWS1nBPAcIK9qrAt38vZ29uFW0QUmBnDfpTQsQqrux6ChYGIIQ7NPOBjYjlO8AYbFDbb10QX6uvl7Onhw36VWKlBks7G7AMCFYEoIQ7NPOejaD/xKjjh6ERJQUFM930dEZoeJODBOCRSEIwQ5p+0jTx5L5DRCScU2ZoHiOBW1UUmBcvraYY8EZIcJODBOCRSEIwQ7tqJemqUWOS4zWdTWQIKjdg7hVtFlJQZyDMMVfaOxjmj6OJQGGB0EIdmhnA+eJEwVNJcloDhIRUbhnKGNSY4+WV0FRoCy1uKsew4RgMQhCsEO7Gth0zgOExY651vYFJQbG5XOdRDEzVMAwIVgQghDsTV0Pax1gCb58d6XHAOH/JAbG8Z1WPzNU2I5hQrAcBCHYm611bFaoKPLLwfaBjvb+jijvEdwq2rikoLgCrsOE8b5Cj45VdSELwTIQhGBvttezWXwHCLUlCYGxosCzpk0b5Rvd1NvSMdDJq6BANCNURKMQLAVBCPZmZwObFcZ7gBD9omcRBTE+MJbv3oSz0DsKloMgBLtS0s6UAsV48gzCfG0xXhk9R3JgPN/e0dlhwvZ6CUkIFoEgBLuyvY7N4doc7NH11nbVj/EbxbGmHUgOjue7SW+kh+CpEorbEIVgAQhCsCvcBwiLmkpj/UerRCXHmnYg1m90VUd1r47nNPhZocL2OgQhWACCEOyHgdFejTSd91rbSYHoFz2XSqEa6zeqpLmMY00ME4KlIAjBfhxpYmHugtqVZ818LdaUubCkoPiCJp7DhDNDxT0aSYcVZsDsEIRgP7bVs9lc+0UHDYMV7afjAsZyrGk3koLi8ht5BmGAC8V4CrlaNArB3BCEYD9+rJXmhvO8pUuay6K9I12Uzhxr2o2EwHFlracGDYMca84JE7Zh0VEwOwQh2IkuHeW1sGw134kT6Be9KFelS6R3RFlrBceac8LEH2vRIgRzQxCCndjVIE0MFNy5vt1Z0FSciM14Ly45KJ7vlkxZaqGojbUNcCwJcHkIQrATW+vYnDCe97NeMpQ2l2PTiUtICorPbyziWNBFQVOChV0N6B0Fs0IQgp34sZbzVPqy1pNhniEeTu4ca9qZ5KD44uYTesnAseacMHErZhOCeSEIwR7U9rCWAZbizzMIjzcWpQQlcCxofzydPEI8gk+28R0mFH5EEIJ5IQjBHmyuZXPCeG69RETHtUXJwQjCy0gOSjjOtXc00U/o1bNKbMkEZoQgBHuwjfcSowZmKG46gQHCy0oJ5hyEAtHsULw7CmaFIASbZ2C0rU6aF84zCMtbK0I8gr2cPTnWtEspwQlFzaUS4/l6y7xw9I6CWSEIweYdamLh7kKoG+cBwmQMEJrAy8kz0C3gZFslx5rzwsWd9dIgXh0Fc0EQgs3bXCMtiOC8fXy+tigFA4SmSeE9TBjgQmO8hf2NaBSCmSAIweZtqmXzua6sJjGpqOlEIgYITcN9mJCI5ocLm2vRJAQzQRCCbWvup5MdbEowzxbhybbKAFc/XxdvjjXtWHJQQmFTCd9hwvkR4uYatAjBTBCEYNu21EozQkUV1xs5T1MwXp3Is6Jd83Xx9nf15TtMmBEo1PWyuh5kIZgDghBs2+ZatoDr+6JEdKyxcHxwEt+a9i1VnXSssZBjQYVAs8PELXh3FMwCQQg2TGK0lffECQMzFDWVYtOJYRkfnHRMwzMIyThMiN5RMAsEIdiwI80swEUY4cF1idGWUyEeam9nL4417V5KcEJhUwnfRUfnh4vb6rFhPZgDghBs2A/V0rUjOPeL5mkKUoMxQDg8Xk6eIR7q8tZTHGsGu9IoLyEHkyhAfghCsGE/VLOFEZzv4bxGvClzJVKDE/M0BXxrLowQN1SjSQiyQxCCrarvZdXdbHIQzxahTtKfaDmZFIgBwmEbr07Ma+QchNeOEH6oRosQZDeMINy3b9/cuXNTUlIee+yxgYFz95DWarVPPPHE1KlT09PT77vvvsbGRq7XCXCuH6rZvHBRyfVvuZLmsgivMOxBeAWSgxJKW8p1Bh3HmqkBQreeTnYgC0Fepv4WaW5uXrhw4c0337xq1aoDBw4888wz5xxQXl7e39//7LPPfvDBB3V1dUuWLOF8pQA/t6GGLeQ9QHi8sXA8BgiviLvKLdIroqSlnGNNgeiaCGED3h0FmZkahCtXrkxPT7/rrrvi4uJefPHFDz74QKf72Z9+WVlZr7/++owZM1JSUl599dUDBw50dXXJcMEAREQDBtrdIM3jurIaER3V5KeqMYPwCqWqk45q8vnWXBghbKjBMCHIy9TfI4WFhenp6caPJ0yY0NraWldXd7GD8/LyQkNDPT2xhQ3IZUc9S/YT/J151uzX959sq8QehFcsTZ2cxzsIZ4eJuVrWMci3KsDPKE08TqvVxsbGGj9WqVQeHh6NjY1RUVHnH1lfX3///fe/8cYbl6i2ffv26Ojos78yadKkDz744IIH9/T0CALnHjAYrp6eHsaY9fwgvq1Uzg6m7u5+jjUPNx4b5R2t79d3UzfHshz19fU5OTkpFApLX8iFRbtGnGo7rW1vclO6ciw7KUC1vqJ3cYS1tAsHBgZEUVSpVJa+EIem0+kMBoNer7/skW5ubqJ4mSafqUHo7e3d09Nj/FiSpN7eXh8fn/MP02q1s2fPvv/++2+66aZLVJs8efKrr7569ldcXV09PDwueDBj7GL/BObk7u5uJUHIiDbV67deo/DwcOFYtqSsPCMs1ZpvNoVCYc1BSERxAWNP9ZyeEpbBseYvYqQtjcrbxlnLd61SqRCEFmcMQhcXPr8BTA3C6Ojo8vKfhsErKioUCkV4ePg5x7S0tMyZM2fp0qWPP/74pau5ubnFxMQM91oBjI42Mw8VjfXmnMpHNMcfnfgnvjUdTZo6+UhDPt8gXBQp/PWIQScp+C6tDjDE1Dtr2bJlGzdurKioIKI33nhj8eLF7u7uRLRy5cpt27YRUWtr65w5czIzMx944IG2tra2tjaDged6SwBDvjsjXR/JOQXb+tubelti/UfxLetoJqhTjmqO860Z6iaM8Rb2aPDuKMjF1CAcN27cU089lZaWFh4efvDgwVdeecX49fXr1x84cICI9uzZU1VVtXr16pH/X1VVlUwXDQ7u2yp2fSTn1sERzfGU4ERRQKPjqoz2G9nW39HU28y37PWR4ndnrGWMEOyPqV2jRPTwww//8Y9/7OrqCgwMHPrimjVrjB8sXrx48eLFnK8O4DwVnay5n2UEcm4RHtUUpGHixFUTBWF8cGKepmBezEyOZRdHCvM3S29MJqsYowa7M7y/f11cXM5OQQDz+/YMuz5SFHn/RszTFKSpUzgXdUip6qSjvNdaG+cjuCroeAt6R0EW6AgCG/PdGYl7v2h1Zx0RG+EVxresY0oPGX+k4RgjzqG1KFJA7yjIBEEItqSxj4ra2Kwwzu3BQ/V5GSGpfGs6rFAPtavStbL9DN+yS6LEr0+jRQiyQBCCLfmmSromQnTifdseasjLCEUQcpMROv5QfR7fmpOChI5BKm1HFgJ/CEKwJWtPSzdEcW4ODhoGi5pK09TJfMs6svSQ1MMNx/jWFIgWRwloFIIcEIRgM5r7Ka+Zzee90Ha+tnikb7S7yo1vWUeWGpxY2lLep+e5AB4R3Rgtfl2FYULgD0EINuObKmleuOg6jCk/JjnUcCwjZDznoo7NRekyzn/MscZCvmWzggVtH53qRKMQOEMQgs34+rR0YzT/iWSH8aaMDNJDxh9u4DxMKAq0OEpYi95R4A1BCLahbYBym9iCCM53bFNvc1t/x2i/kXzLQkZo6qF6zsOERHRDlPj1afSOAmcIQrAN685Is8NEd979orn1eRNCUkTr2FXDnsT4RPbr++u7NXzLTgsRqnvY6S40CoEnBCHYhi8qpFti+MfVwfojk8LSuJcFgYSJoWkH6g7zLasQ6MZocXUlghB4QhCCDdD20dFmdg3vflGdQXessXBiCIJQFpPDJhyoO8K97C0x4mcn0TsKPCEIwQasrpSuG8H/fdFjjYXR3pFezp6c6wIREaWHjC9pLuvV9fEtm6kWeg1U1IZGIXCDIAQb8EWFdMtI/vfqgfrDk8MmcC8LRi5Kl7iAsdy3JxSIlkYLX1SgUQjcIAjB2p3pZpVdbHaoDAOEdUenhKVzLwtD5OodHSmurmBoEgIvCEKwdqsq2I3RopL3rXq6o9rADNE+kZzrwlmmhGXsrzsk8c6sFH/BRUGHtIhC4ANBCNbuiwrplhj+N+rBuiNTwjK4l4WzhXgEezl7nWyt4F755pHi5+gdBU4QhGDV8ppZl44y1XL0i2LihDlMCUvfz3sSBRHdPkr4okIaRBQCDwhCsGr/PSn9ajT37eipY6DzVPvp1OAk3oXhXFPC0nNqc7mXjfIU4nyETTVIQuAAQQjWSy/RmkrptlH8m4P7ag+lh4x3UjhxrwznSAiMa+1v477EDBHdOUb870kMEwIHCEKwXhtqpDHewkgv/kG4t+Zgdvgk7mXhfKIgyNQoXBot7mqQmjjv9QSOCEEI1uvTk+zOMfxv0X79QL62aGIoBgjNJDtiUk7NQe5lPVR0TYS4Gq/MwFVDEIKVahmgHfXSjdH8b9Hc+qPxAbEeTu7cK8MFpQUnV7RXtfa1ca/8q9Hif7HcGlw1BCFYqZUnpetGiF4q/pX31hzMjkC/qPmoFKr0kPEH6vnPrJ8ZKjT30/EWjBTCVUEQgpX6sEz6bSz/+1MvGXIbjmaGT+ReGS4hO2Ly3poD3MuKAt01VvygDI1CuCoIQrBGORpmYJQlw/TB442FI7zC/F19uVeGS5gUmlbYVNqj6+Ve+TdjhNUVUq+ee2FwIAhCsEYflEn3xPKfPkhEO87snT4iS4bCcCnuKrfkoPh9tYe4Vw5zFzKDxTWVaBTClUMQgtXpGKTvz0i3jZKlXzSnNndqxGTuleGyZkZm7zyTI0fl38YK6B2Fq4EgBKuz8pS0IEIMdOFf+YjmWKR3RLB7IP/ScDmZ4RPztUWdg13cK18TIVZ3Y4dCuHIIQrAujGhFifS7cbLcmTvO5MxAv6iFuCpd0tTJcvSOKgS6e6y4ogSNQrhCCEKwLtvqmEKgqTK8JqMz6PbXHpo6Av2iFjMzMnvHmb1yVL53nPhFhdQ2IEdtsH8IQrAubxZLDybIclvmNuSN8o0OcPWToziYYnJYeklzWftAB/fKale6JkL8BJPr4YogCMGKVHWxA1rplpEy9YvunRGJflFLclE6Z4Sk7pVhuTUi+lO8+HaJJGGgEIYPQQhW5M0S6a4xopuSf+UeXW9u/dHpIzL5l4bhmBM97cfTu+SoPClI8HOmTbVIQhg2BCFYix49fXpS+r08r8nsrt6fGpzk7ewlR3Ew3cTQtNquejl2ZSKi++LEfxcZ5KgM9g1BCNbigxPSrFAxylOOafS0pXLH3OgZclSGYVEIipmRWT+e3ilH8VtGiic66BiWHoVhQhCCVdBL9K8i6c/yvCbT2KM93VE9CfsuWYe50TM2V+5gxD+uVCLdFye+WohXZmB4EIRgFb6slEZ50cQgWZqDmyt3zorKVilk2MkChm+s3yhXpUtR0wk5iv9unLilVjrTjUYhDAOCEKzCq4XSo0kKmYpvq9o1L3qmTMXhCsyJni5T76inin4zVvxXERqFMAwIQrC8H+uYgdHccFmagwXaYlEQY/1Hy1EcrszcqOm7zuzr1/fLUfz+ePHTk1ILJteDyRCEYHnP5RkeS5Zlrwki+v7klutGzZenNlyhADf/xKBxO+RZgzvUTbgxWny9EK+PgqkQhGBhW+uYtp9+GSPLrdg50HWw/sjcmOlyFIersWj0/O9Pbpap+PIU8Z1SqVmWBifYIQQhWNjzxwzPpIoKedqDmyu3TwnP8HLylKU6XIWMkLS2/vby1go5io/wEG6IFt8oRqMQTIIgBEvaWsc0fXI1B4loQ8XWRegXtUqiICwcNeeHUz/KVP/JFPHdUqkVI4VgAgQhWNKzeTI2B483FjKihMBYWarDVbt25NwdZ/b26vrkKD7CQ1gSJb6KkUIwAYIQLOa7M1KXjm6WrTm4rnzjkjHXyFQcrp6fq2+qOmnL6R0y1X96vPheqVTbgzmFcBkIQrAMA6PlR6SXMhQyvS2q6dEeayycHzNLlurAyU2x1689sV5ismRVmLtw11jx+WOYUwiXgSAEy/i4TAp0oXnyzB0korUnvl84co6r0kWm+sBFYuA4L2ePA3WHZar/1xTFt2ekknY0CuFSEIRgAb16eu6Y9PJEuZaS6dX1bTm9c8nYhTLVB45ujL1+zYnvZCru40SPJimWH0ajEC4FQQgW8M98w1S1MCFArubgD6e2ZISkBrkFyFQfOJo+YkpDt+ZEy0mZ6t8XJ+a3sh31aBTCRSEIwdwqu9g7pdJLGXLde3rJ8HXZDzfFLpKpPvClEBRLxixce2K9TPVdFPTaJPFP+w06NAvhIhCEYG4PHpAeTVKEucvVHPzx9M4wzxAsLmpDFo2ef7jhWE1nnUz1F0eKUZ70VgmSEC4MQQhmtbGGlXewB+XZd5CIJCatKvn6zsSbZaoPcnBXuS0Ze83nJV/Ld4rXJyn+cdygkWXKItg8BCGYT4+e7ttveHOKwkm2+27r6V3+rn7JQfFynQDkcVPs9ftrD9V21ctUf4y3cPdY8YEDmF8PF4AgBPP562HD9BBhTphcnaISkz4r/urOhF/KVB/k465yWzxmwRcl38h3iqdTFUWt7JsqdJDCuRCEYCYHtezrKvaKbFMmiGhbUDBtVgAAFXdJREFU1W5vZ+9UdZJ8pwD53Bi7aE/NgfpujUz1nRX04VTF/QekNixACj+HIARz6DfQXXsMb0wS/ZzlOsWgYfDD/M/vHX+HXCcAmXk5eS6Nvf7945/Kd4rJQcKSSOGhXHSQws8gCMEcHs01JPsLN0TLeL+tLVsf6z8qMTBOvlOA3JaOW1zSXFbUdEK+U/wjXbGvka09jQ5S+B8EIchuUw37vpqtmCJjp2jXYPea0u/uSUFz0LY5K5x+nXjLiryPGck1/91DRV/MUPxxv+FMN6bYw08QhCAvTR/9Zq/+8+kKX9k6RYnoPwVfzIjMDPcMlfEcYBbzYmb26/tzag7Kd4q0AOHPCYo7dxsMiEIgIgQhyEov0c079L8bp8hSy/WmKBGdbKvccWbvrxJvke8UYDaiIN6XdvebRz/q1/fLd5a/JIlKgZ4+isFCIEIQgqwezjU4i7Q8RcbbTGLstUPv3Jtyh7ezl3xnAXNKVSclB8X/p+AL+U4hCrRmlvKLCgwWAhGCEOSzqkLaVMu+nKWUaQN6o2/LNypF5fyR2HfQrvwp7e4fq3adbKuU7xR+zvT1bMUf9xuK29BD6ugQhCCLfY3szwcN38xW+DjJeBZtb/N/C1c/kvEHgeQMWzA7L2fP3yQtezV3hYHJ2HuZ4i+8mK5YvNXQJGMvLNgABCHwd6qTLd1u+HSaMsFXxnySGPvHgX/dFLso0jtCvrOApSwcNcfDyX1l0RpZz/KrMeKtI4VrNut79LKeB6waghA4a+iluZsML6SL8u0+b/Rl6TqdQXdr/A2yngUsRSDh8ckPfFu+qbhZxmmFRPRMmiLBT7hlh0GP4UJHhSAEnrR9NGuj/t5Y8c7R8t5ap9oqV5esezLzYVHAPWy3Alz9Hs74/f/te61XJ+O2EQLR+1kKRuy2XZhQ4aDwSwS4ae6n2Zv0N48UH0uW977qHOx6au+L90/4rdo9SNYTgcVlR0yeEJLyz4NvyDfFnohUIn09W9mjZ2gXOiYEIfBR18NmbNBfHyk8NV7em0pi0v/te21qxORZUVNlPRFYiQcm3Nve3yH3YKGTSF/NUrYPstt2GQaRhQ4GQQgcVHSy6RsMN0QLz6fJuI6a0TvHPpGYhNXUHIdSVDyT9Zf1J7fsqz0k64lcFLR+rpIRzd+k79TJeiqwLghCuFq7G1jWev3fxovPpMqegmtKvz1Yd+TprEcxNOhQ/Fx9n5/6xMu5bxY2lcp6ImcFrZqhGOsjzNigr+nBgKGjwG8TuCpvl0g379B/NkN5h8xvxxDR5soda8vWvzLzGU8nD7nPBdYm1n/0U5mP/m3PC2Wtp2Q9kUKgdzIVt4wUJ31nyNEgCx0CghCuUJeO7thleP+EtG+Rclao7PPZfzy964P8la/Nej4YL8g4qlR10sMZf3hi1/MV7VVyn+uRRPE/0xQ3bte/WighDO0eghCuRK6WjV+nd1XSgUXKGE/ZU3DtifUf5n/22sznsL+Eg8uOmPzAhHsf3v5UgbZY7nPNDRNyr1euq5Lmb9I39Mp9NrAkBCEMT4+eHsk1LN6qfzlDfC9L4aaU93QSY+8d++/3pza/NfefWEEGiGjaiClPZT7y1N5/7qk5IPe5Ij2EXQuVU4LFlHW6j8rQNLRbCEIYhnVVUuLXem0fFd6gWhIl+83TOdD12K5nS1vK35rzzyC3ALlPB7YiVZ308sxn3zr60XvH/isxeec6KEV6OlXcukD5YZk0Y4M+vxVpaIcQhGCSo63i9A2GZ/KkD7MVn05XBLjIfsbCppJ7Nj8U4x356qznvJw9ZT8f2JTRvjEfLHitvK3i4R1PNfY0yX26JD9h33XKm2PE+Zv09+4X6tBTal8QhHAZ+xvZNVv0t+9T3T5KyFuinCn/ezH9+v5/H/ngmZyX759wz+9Tf60QZJ+VAbbI29nr5RnPTFCn3LPpoe9PbpZ16RkiEgX63TixbKkq2JUmfE+/32eo6kLr0E4gCOHCBiX6/JQ06Xv9HbsNSyLF4wsHfjNWlHVnQSKSGNtcueP29X/o0fV8svDNKWHp8p4PbJwoiMvib3xjzgubKrf/YcujhU0lcp/RS0XPjWeFi8nfmSZ8q79xu2FXA+LQ5sn8qgPYoMNNbOUpaXWFlOwv/DVZvHaEKArU3S3vSSUm7a7e/1nxV65Kl2ezH4sLGCvv+cCORHlHrJj30raq3c/ve22Ub9Sy+BvjA2JlPWOAC/3fBMXjyYpPT0p/3GcwMLp9tHjbKCHSA/ti2qThBaFOp1OpVFd/DFibQYn2NbIfqqVvqphKpNtGibnXK6PlnxdBRG39HVtP71xXvtHf1e83ybehFQhXQCBhTtT0aRFTNlZs+799rwW4+V8/en52xGRnhYwbQ3uo6A9x4h/ixINa9ulJacK3higP4YZocUGEkOQnIBJtiMCYSe367u7uO++8c+vWrQqF4rHHHnv88cfPP+all176xz/+YTAYZs2a9emnn3p6XvgFh1WrVm3YsOHzzz838RK7urouVgquxoCB8lrYXg3bq5H2alisj3BNhLgkUkj0u8Aj3N3d7e7uLvB7ujsGOvfXHd5bcyBfW5wVPnHR6Ply/xVvB/r6+pycnBQKDJpeisSkPTUHNpzaeqLlZHbEpOyIyWnqJCd+iTgwMCCK4vl/7hsY7W5g356RttSyLh2bFSpmq4UstRDrLYhIRd50Op3BYHBx4fPanqlB+Ne//vXo0aPr16+vr6+fOHHiunXrpkyZcvYBubm5CxcuzM3NjYiIWLJkSUJCwosvvnjBUghCS6npYeUdVNTKitpYXgs70c7G+QhZamGqWpgeIvo5X+r/cgnCtv6OEy0n87VFxxuLqjtr00PGZ4ZPzAqf6KZyvZqyjgNBOCxNvc27qvfn1OaebK2ICxg7PjgxMXDcaL+Rrsqr+u15sSA82+kutquB7dGwfY2ssZeNDxBS/IUEXyHBVxjjLVz6WQNTWCYIw8LCPvnkkzlz5hDRgw8+2Nvb+/777599wO9//3uFQvHWW28R0c6dO2+++ebGxsYLlkIQykcvUcsANfezxj6q72X1vVTXw053UVU3q+xkniqK9RHifYVEPyHFT0j2F1xM/o063CAcNAw29bZoerQN3ZrarobT7dWnO8706vpi/UcnBMamBieNCxirEjFEPTwIwivTPdiTry0+1lhQ0lxe0V6ldg+M8h4R6R0R7hkS4qEOcQ/ydfUx/eVkU4LwbO2DdKSJ5beyojZW3MbKO5hSoBgvIdpTiHCnCHch1J1CXIVAVwp0QUaaim8QmvSbqLe3t76+Pi4uzvhpXFzc6tWrzznm1KlTS5YsGTpAq9V2dnZ6eXldsKBOp2tra/vZdSiVF0u73oGBtj7bm7ajZ9Sjv/wfGT16Omcj0C4dGf84kRh16RkRGSTq1hMRdeqIGHXomF6iHh316KnXQD2DrFNHnTpqH2Q9evJzIl8XIdCZAl2FYFeKcRWywyncXYhwFzx+/uQOdNMAERH16voMF5mV3KvvMzADEfX39QsqQc8MRNQ92MuI9en79czQpevp1/f36Qe6dD3dut7Owa6Oga62gY5+w2CAi0+Qa0Coe1C4Z8iiyOmRnuGhZy8T2t+PTd+Gi/X3S/pBAUE4TG5Ek/3iJvvFEZGBGc501Z/prK3qqjtYfUjTq23sbe4Y7PJy8vB28vRy8vR29vRQuXmq3F2VLi5KZzelq7NC5SQ6qUSVi9KJiERJUClUSqWSiJwUTs7ihTtdnRVOTgoVEQlE6R6U7kE04qd/aumnmm5W28Pq+6immR2upqZ+1jpATf2sV08+KvJyFjyV5OlEnirBTUnOInk7k4LIQyUoBXJTERF5OxERiUSeqp/+PFUpyPXnt4aLgpxNuFk8Vba0mYu7i7O/lw/fmiYFoTG0PDx+WvLfy8urtbX1/GOGDjBGWltb28WCcPv27TExMWd/JTMz8/xwNfpq33cbOr435TrtlUD0/9q785gmtjUA4KctAgXKXnnVaqQoUq0YEUnQEOOLC0YgUcnlBYnPPVGfolYTDIobxuC+BiOb0YRAgUgMwUZMEEOUAEbArSmWYCtYlLYMtZbSZd4fc224FLXe23Zq5/v9NT3zdeZjmXztOTPnoB8+IxWAUABCUxgIEf/0BmQ2IOUwUiLUacfxmRYa/TuHZ1oQA//z9D44zcuCE410HPniNIYF+ZuRD458zSjEQvM301hmFGBCQWYay4QQ0iCkQajberSP9vy0ADgZE6EYhMaOSFtoNMxLp2V8GWb067xoOgauoyOMgQ/QaSM0fJSOjDTcREcGGkIIjdCR5VvPyCgNjX6nhozScaPdIwkMhNgIsa2XsBEhIzLp/7yEJoQjig488kb/VfCfC8Q3QpPJ9NN4Pz8/Ov0nhd6uQhgeHk6j0TAMCwoKQghpNJqIiIhxMWw2G8MwYntoaIho+d4Bk5OT7e8a/e+///gfa6udwcBJHH6zDPgboGvUeeyfx/ZXu0aBMzi2a9SuL8Q+Pj48Hu/FixfEy46OjpiY8Tf48fn8sQEzZszw8/NzSIoAAACA89jbM7x9+/b8/HyFQvHkyRORSLRt2zaEkEqlWrVq1efPnxFCW7durampaWxs/PDhw6lTp7Zv3+7ErAEAAAAHsfe2PaFQODg4mJSUFBwcXFhYOG/ePIQQjuMGg4G473Tu3Lm3bt06cOCARqNZv379oUOHnJg1AAAA4CD2Pj7hQPD4xO8IxgjdAYwRugMYI3QHJIwRkmvfvn0ymYzsLKju+PHjbW1tZGdBdZcuXWpoaCA7C6orKyurqqoiOwuqu3fvXlFRkaOO9hsUwvb2dq1WS3YWVNfV1aVSqcjOguokEkl/fz/ZWVCdTCbr7e0lOwuqk8vl3d3dP4+zz29QCAEAAADngUIIAACA0ki4Waa1tXXTpk0Gg8HOeI1Gw2KxiAmNAFkwDGMymd7eTlzUBvyUVqudNGmSo24QAH/Ply9fGAwGkwkzxZNJr9ebzWbrdGY/UFdXx+fzfxxDQiFECCkUCqPR6PrzAgAAoBQul/vTT/DkFEIAAADATcAYIQAAAEqDQggAAIDSoBACAACgNCiEAAAAKM2NCiGO4zk5OaGhocHBwXv37jWbzbYxdXV1PB7P399/2bJlCoXC9Ul6vIqKivj4eF9fXzabvWvXrpGREdsYPp8f9U1OTo7rk/R49fX1UWO0tLTYxjx79kwgEDCZzPj4+FevXrk+SY9XXFwc9VfESqtjrVy50ro3MzOTlDw9j06n27dvX1JSUlRU1NhJfEZHR7dt2xYUFBQeHp6fnz/he0tLSzkcDovFWrt2re3f67twtyESiXg8Xn9/v0qlmjdv3q1bt8YFqNVqFov14MEDo9G4Z8+elJQUUvL0bGVlZY2NjXq9/v3797GxsXl5ebYxPj4+z549k8lkMpns06dPrk/S41VWVi5ZskT2zdevX8cFGI1GLpdbXFxsNpsLCgpiY2NJydOzDQ0NWf8EBw8eXLJkiW0Mn88XiURETF9fn+uT9EgYhh07dqy2thYhJJVKre3nzp1LSEjAMKy3t3fq1KlisXjcGyUSCYvFev78uV6vX7du3a5du+w8oxsVwtWrV589e5bYLioqSkxMHBdQWFho/V9UKpVeXl4fP350aYoUk5+fn5qaatvu4+PT39/v+nyoo7Kycvny5T8IqK+vnz59usViwXHcYDAEBga2tbW5KjsqmjNnTmlpqW07n89vbm52fT5UYDKZxhXCOXPmVFRUENtHjhzJyMgY95bDhw9nZmYS221tbSwWa3R01J5zuVHXqFQqFQgExLZAILCdULW7u5tYBxEhFBERERoaCqtSOI/FYhGLxYmJiRPuXbx4cWRk5IYNG6CD2klaW1unTJkyf/78goIC22GC7u5ugUBArIrl7e0dHR0tlUrJSJMSWlpa5HJ5enr6hHszMzO5XG5aWhp0UDsVjuMymeynNWJsgFarVSqV9hzcjQqhRqOxzpfDYrHUajX+14f9NRqNv7+/9WVgYKBarXZpilRy8uTJ4eHh/fv32+4SiURPnz5taGiwWCwpKSnEBzfgQAsWLGhoaOjs7Lx8+fLNmzfPnTs3LgCuBVcqKSnJyMiYcFXUixcvNjc3P336dObMmStWrPiFQSnwi3Q6ncFgGFsjbNfDGVtEfH19J02aZOd14UaFMDw8HMMwYntoaIjNZo9bBjYsLGx4eNj6kohxaYqUcenSpfLycrFYPOG0lmlpaRwOZ+bMmaWlpW/fvoXvIg43a9ashIQENpu9bNmyvLw829Xv4FpwGZ1OJxKJtmzZMuHe5OTkadOmTZ8+/cKFC97e3s3NzS5OjzoCAgKYTObYGjF58uRxMWFhYdYAnU5nNBrtvC7cqBDGxMR0dXUR211dXbNnz/5BQF9fH4ZhUVFRLk2RGoqKiq5cufLw4UMOh2NPPA6z9DnTuI+DhJiYmJcvXxK/+ZGREalUanu9AIeoqqricDjfGyMYi0aDGSudKzo62v4a0dXVFRISYlssJ/YPxzMd6P79+1wuVyKR9Pb2RkdH37lzh2jPysoibgTAMCw4OLi8vBzDsM2bN6enp5Oar2cqKSkJCAiorq5ub29vb29/8+YN0Z6bm1tTU4PjeGdnp1gsViqVPT09GzZsEAgERqOR1JQ9UE1NTUdHh0qlampq4vF4p0+fJtr3799fX1+P47jZbObxeOfPn9dqtbm5uYsWLSI1X0+WlJRkvYmPcPXq1Rs3buA4LpfLq6ur+/r6FArFoUOHIiIiiAEd8M91dHS0trYihO7du9fe3m4ymXAcv3btWmxsbF9fX2dn5+TJkx8/fozj+MjISFpaWm9vL47jPT09AQEBjx49UqlUycnJBw4csPN0brS2UWpqqkQiSU5ONpvNW7ZsycrKIto/ffpErNkUGBhYW1srFAqzs7OTkpJu3rxJar6e6fXr17Nnzz5z5gzxks/n3717FyGkVqt1Oh1CyGAwnDhxoqenx8/Pb/HixXV1dbBClsO9e/fu6NGjSqWSy+Xu2LFDKBQS7Wq1Wq/XI4TodHptbe3u3bvPnDkTFxdXUVFBar4ea2BgQK/Xb9y4cWzj8PAwg8FACFksluvXr+/Zs8fLyysuLq6hoSEkJISkTD2NUCgcGhpauHAh8bxgU1OTv7//zp075XJ5QkKCr6/vsWPHli5dihDCcXxgYIC4UyEyMvL27dvZ2dmDg4Nr1qw5deqUnaeD7/IAAAAozY3GCAEAAADXg0IIAACA0qAQAgAAoDQohAAAACgNCiEAAABKg0IIAACA0qAQAgAAoDQohAAAACgNCiEAAABKg0IIAACA0qAQAgAAoDQohAAAACjt/5WjI6/qqWu0AAAAAElFTkSuQmCC",
"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": "3.326664185386434e-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.00022348749912162638"
},
"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.2"
},
"kernelspec": {
"name": "julia-1.6",
"display_name": "Julia 1.6.2",
"language": "julia"
}
},
"nbformat": 4
}