{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Gross-Pitaevskii equation in one dimension\n",
"In this example we will use DFTK to solve\n",
"the Gross-Pitaevskii equation, and use this opportunity to explore a few internals."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## The model\n",
"The [Gross-Pitaevskii equation](https://en.wikipedia.org/wiki/Gross%E2%80%93Pitaevskii_equation) (GPE)\n",
"is a simple non-linear equation used to model bosonic systems\n",
"in a mean-field approach. Denoting by ``ψ`` the effective one-particle bosonic\n",
"wave function, the time-independent GPE reads in atomic units:\n",
"$$\n",
" H ψ = \\left(-\\frac12 Δ + V + 2 C |ψ|^2\\right) ψ = μ ψ \\qquad \\|ψ\\|_{L^2} = 1\n",
"$$\n",
"where ``C`` provides the strength of the boson-boson coupling.\n",
"It's in particular a favorite model of applied mathematicians because it\n",
"has a structure simpler than but similar to that of DFT, and displays\n",
"interesting behavior (especially in higher dimensions with magnetic fields, see\n",
"Gross-Pitaevskii equation with external magnetic field)."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We wish to model this equation in 1D using DFTK.\n",
"First we set up the lattice. For a 1D case we supply two zero lattice vectors,"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"a = 10\n",
"lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"which is special cased in DFTK to support 1D models.\n",
"\n",
"For the potential term `V` we just pick a harmonic\n",
"potential. The real-space grid is in ``[0,1)``\n",
"in fractional coordinates( see\n",
"Lattices and lattice vectors),\n",
"therefore:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"pot(x) = (x - a/2)^2;"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"We setup each energy term in sequence: kinetic, potential and nonlinear term.\n",
"For the non-linearity we use the `LocalNonlinearity(f)` term of DFTK, with f(ρ) = C ρ^α.\n",
"This object introduces an energy term ``C ∫ ρ(r)^α dr``\n",
"to the total energy functional, thus a potential term ``α C ρ^{α-1}``.\n",
"In our case we thus need the parameters"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"C = 1.0\n",
"α = 2;"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"... and with this build the model"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using DFTK\n",
"using LinearAlgebra\n",
"\n",
"n_electrons = 1 # Increase this for fun\n",
"terms = [Kinetic(),\n",
" ExternalFromReal(r -> pot(r[1])),\n",
" LocalNonlinearity(ρ -> C * ρ^α),\n",
"]\n",
"model = Model(lattice; n_electrons, terms, spin_polarization=:spinless); # spinless electrons"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"We discretize using a moderate Ecut (For 1D values up to `5000` are completely fine)\n",
"and run a direct minimization algorithm:"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter Function value Gradient norm \n",
" 0 1.765089e+02 1.182079e+02\n",
" * time: 0.0006511211395263672\n",
" 1 1.708323e+02 9.753269e+01\n",
" * time: 0.0019850730895996094\n",
" 2 1.256698e+02 9.437160e+01\n",
" * time: 0.004372119903564453\n",
" 3 1.076115e+02 8.272244e+01\n",
" * time: 0.007039070129394531\n",
" 4 4.933072e+01 6.689338e+01\n",
" * time: 0.009698152542114258\n",
" 5 3.165150e+01 4.701766e+01\n",
" * time: 0.012061119079589844\n",
" 6 1.017234e+01 9.863353e+00\n",
" * time: 0.01418614387512207\n",
" 7 1.013702e+01 1.736950e+01\n",
" * time: 0.014991998672485352\n",
" 8 6.807519e+00 7.913527e+00\n",
" * time: 0.016422033309936523\n",
" 9 3.305763e+00 3.922085e+00\n",
" * time: 0.018199920654296875\n",
" 10 2.527361e+00 9.239986e+00\n",
" * time: 0.019628047943115234\n",
" 11 1.819008e+00 4.401821e+00\n",
" * time: 0.021079063415527344\n",
" 12 1.653479e+00 2.715344e+00\n",
" * time: 0.022576093673706055\n",
" 13 1.491605e+00 1.672313e+00\n",
" * time: 0.02370309829711914\n",
" 14 1.436884e+00 2.405348e+00\n",
" * time: 0.024821996688842773\n",
" 15 1.311629e+00 1.844782e+00\n",
" * time: 0.026051998138427734\n",
" 16 1.232842e+00 1.205995e+00\n",
" * time: 0.027163982391357422\n",
" 17 1.226469e+00 7.379742e-01\n",
" * time: 0.02796006202697754\n",
" 18 1.185173e+00 5.690920e-01\n",
" * time: 0.02908015251159668\n",
" 19 1.151745e+00 2.702299e-01\n",
" * time: 0.030213117599487305\n",
" 20 1.146948e+00 2.121817e-01\n",
" * time: 0.03132009506225586\n",
" 21 1.144795e+00 1.357340e-01\n",
" * time: 0.03242802619934082\n",
" 22 1.144300e+00 3.863231e-02\n",
" * time: 0.03355813026428223\n",
" 23 1.144113e+00 1.533640e-02\n",
" * time: 0.03465700149536133\n",
" 24 1.144099e+00 1.431367e-02\n",
" * time: 0.03576207160949707\n",
" 25 1.144072e+00 1.131954e-02\n",
" * time: 0.03686714172363281\n",
" 26 1.144058e+00 1.158915e-02\n",
" * time: 0.03799796104431152\n",
" 27 1.144042e+00 5.220822e-03\n",
" * time: 0.03909802436828613\n",
" 28 1.144039e+00 3.325529e-03\n",
" * time: 0.040216922760009766\n",
" 29 1.144038e+00 3.065392e-03\n",
" * time: 0.041355133056640625\n",
" 30 1.144037e+00 1.240680e-03\n",
" * time: 0.04250812530517578\n",
" 31 1.144037e+00 1.036939e-03\n",
" * time: 0.04366898536682129\n",
" 32 1.144037e+00 9.482615e-04\n",
" * time: 0.04481196403503418\n",
" 33 1.144037e+00 7.319115e-04\n",
" * time: 0.0459599494934082\n",
" 34 1.144037e+00 3.303911e-04\n",
" * time: 0.047070980072021484\n",
" 35 1.144037e+00 2.160884e-04\n",
" * time: 0.04818296432495117\n",
" 36 1.144037e+00 1.803334e-04\n",
" * time: 0.049317121505737305\n",
" 37 1.144037e+00 1.028342e-04\n",
" * time: 0.05043292045593262\n",
" 38 1.144037e+00 7.265890e-05\n",
" * time: 0.051553964614868164\n",
" 39 1.144037e+00 4.889936e-05\n",
" * time: 0.052670955657958984\n",
" 40 1.144037e+00 2.569856e-05\n",
" * time: 0.05381011962890625\n",
" 41 1.144037e+00 1.637467e-05\n",
" * time: 0.05491495132446289\n",
" 42 1.144037e+00 1.474808e-05\n",
" * time: 0.05602693557739258\n",
" 43 1.144037e+00 9.318347e-06\n",
" * time: 0.057138919830322266\n",
" 44 1.144037e+00 5.161005e-06\n",
" * time: 0.05825495719909668\n",
" 45 1.144037e+00 3.368800e-06\n",
" * time: 0.05934596061706543\n",
" 46 1.144037e+00 2.132506e-06\n",
" * time: 0.06044912338256836\n",
" 47 1.144037e+00 1.623603e-06\n",
" * time: 0.061231136322021484\n",
" 48 1.144037e+00 1.731472e-06\n",
" * time: 0.06238102912902832\n",
" 49 1.144037e+00 8.606646e-07\n",
" * time: 0.06348395347595215\n",
" 50 1.144037e+00 4.209374e-07\n",
" * time: 0.0645899772644043\n",
" 51 1.144037e+00 3.540226e-07\n",
" * time: 0.06574797630310059\n",
" 52 1.144037e+00 1.692595e-07\n",
" * time: 0.06686806678771973\n",
" 53 1.144037e+00 2.020204e-07\n",
" * time: 0.06765007972717285\n",
" 54 1.144037e+00 9.579822e-08\n",
" * time: 0.06843113899230957\n",
" 55 1.144037e+00 1.073623e-07\n",
" * time: 0.06955504417419434\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": [
"ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n",
"ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"Check whether ``ψ`` is normalised:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"x = a * vec(first.(DFTK.r_vectors(basis)))\n",
"N = length(x)\n",
"dx = a / N # real-space grid spacing\n",
"@assert sum(abs2.(ψ)) * dx ≈ 1.0"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"The density is simply built from ψ:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "1.161130389502071e-15"
},
"metadata": {},
"execution_count": 9
}
],
"cell_type": "code",
"source": [
"norm(scfres.ρ - abs2.(ψ))"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"We summarize the ground state in a nice plot:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=3}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ3wc1bkG8Hdmd9V7XfXiJqtbsuQiyb1iMHYoAUzJTQikECCUQGJCvSGhhwCmE4LBGGMwYNxwt+UiF9nqlmzJsupq1XvZ3Tn3w3IV47qyz2x9/p8kefzO6KcZPTrnzDlHYIwRAACAoxItfQEAAACWhCAEAACHhiAEAACHhiAEAACHhiAEAACHhiAEAACHhiAEAACHhiAEAACHhiAEAACHhiAEAACHZoEgbGtrO336tOnHGwwG+S4GTCRJkqUvAUiSJKyJaHGMMfwULI4xxvGXkgWCcPPmzU888YTpx/f19cl3MWCivr4+PPwWNzg4iL9ILG5oaEiv11v6KhydXq8fGhriVQ1dowAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NAQhAAA4NCUlr4AALgMRlTYxio6WUWb6O7ERnuzJD+K8hAsfV0AdgJBCGC92gbpg3Lpo3KJESX6CuGupB+grQ2Goy0szke4e5x4S6yoRLcOwNVBEAJYqTVV0oMHDQvCxY+nKyYHCUTU3z/k5KRUKBQ6ib6vkd4okV4tkj6cppjgj9ahDWhtbb3vvvuwePIVe/TRRzMyMuSojCAEsDq9erprt6G8g62bo5wUdIGQU4m0NFpcGi1+clJauFl/f4LiL6loGFq7xsbGffv2vfLKK5a+EJv0zjvvFBYWIggBHELrIC3aok/wEY4uVTpdLt3uHCPOCxev3aKv6WFvZSkUaBlaN29v75tuusnSV2GTtmzZIl9x/BUJYEUa+lj2ev2sEOGDaYrLpqCR2pV2LlJWdrOf7zDosTUFwMghCAGsRY+Ort1iuGO0+HzGyJp2niraMF85oGe/34/xJ4ARQxACWAWJ0e27DGkBwpWN9jmJtGa28ngre7EQrUKAkUEQAliFR/IM/Xr2Tpbiiiu4KenrOYq3SqVvzyALAUYAQQhgeZtq2dfV7IvZyqucFBjmLqyZpbg311DfyzhdGoD9QxACWFjLAP061/Dv6QofJw7VJgUJv4tX/M8eA5IQwEQIQgAL+9Vew+2jhZkh3KY+/CVF7NXTP4vRQQoy+vWvf/2f//zH+HFXV1dmZmZfX9/FDjYYDFlZWQ0NDea6upFBEAJY0hdV0plu9mz6lQ8Nnk8p0soZiuePG6q70SwEuXR1dfX39xs/fvnll6+99lo3N7eLHaxQKO64445nn33WXFc3MghCAIvp0dEjedK/ppo6ZdB0sZ7Cg4mKPx5EoxAuY/Xq1XV1de++++5zzz1HRIODg5988slTTz21evVqSfrx/snPz3/55Zf/+te/nv3FYTqd7r333rvzzjuJSKvVrl69moiOHj164MABIvr666/r6uqI6JZbblmzZk1nZ6c5vzsTIQgBLObZY4Y5YcI0tSzrwTyaLJZ1sA21aBTCpSxfvvzaa6+tra318fEZHBzMysratWtXVFTUqlWrbr/9duMxK1asEEUxJibmgw8+uPvuu8+pcPToUTc3t+joaCKqrq5+4okniGj9+vVr1qwhoueff760tJSIfHx8xo0bt3v3bnN+dybCEmsAllHWwT6ukIpuUMlU30mk16co7ttvmB2qdOHZ8wo8vVUqfWPG6S4b5l9g3b5f/OIXDz74IBGtWLEiNDT0o48+IqI77rgjNja2oqJi7NixH3zwgfHIW2+91d/ff8WKFS4uLsP/vaioaMyYMaacfezYsYWFhYsXL+byvXCEIASwjMcOSctTFcGuMp5ifriQ6CusKJUeSkLfj5WaEyaM8zbf3ynKC/U+DK9kffTo0ZKSkrlz5xo/7e3tLS8vHzt27Ouvv75ixQpBEDw8PHQ6XUNDQ2xs7PB/7+npcXd3N+XsHh4e3d3dV/s9yABBCGABh5pZUTv7crbsvwH/kSnO+F5/d5zoJVfLE67KOG9hnLeFr8HZ2dn4gaur6/XXX//Xv/51+J/c3d0LCwtffPHFgoKCgIAAg8Hg7u5+zk5SAQEBbW1txo9dXFwGBgbO/tf+/v7h5mNbW1tMTIyM38mVwt+JABbw58OGpyaIzvK3BMZ5C/PDRUylAFMsWLDgm2++ISJfX19fX19BEARBaGpq8vHx8ff3J6LVq1cPDg6e878yMjKKiooYY0Q0evTo7u7u/Px84z+dOHGipqZm/Pjxxk+PHz+emZlpvu/HZGgRApjbD/WsrpduH22mP0OfSRfT1+l/EycGydkNC3bg2muvPXr0aHx8fGZmZldXV2VlZWFhYU5Ojpub26RJk3x9fd3c3Pz8/M75X3FxcWFhYYcPH87MzHRzc/vggw+uu+46URQNBsOaNWveeOONwMBAIjpz5kxXV9eUKVMs8Z1dBoIQwNz+esTwt4niVa6mZrooD+GWUeLLRYYXM/HODJwrPz/fw8Nj+NOnnnrqwQcfLC8v9/DwGDdunEKhIKL9+/cXFha6urrGx8d3dHR4eXkR0fvvv+/k9ONiSA899NB7771nbO3ddNNNS5cufeCBB7q7uz/88EOV6sdO+ffee+/+++8f/tSqoGsUwKy21bNePd0QY9ZH7/EU8aNyqf3cPi0A8vb2Nqbd2V/JzMyMj48f/rpKpUpPT4+PjyciHx8fURSJyMvLa3jw784773RxcRleWUapVAYGBvr7+w/HnsFg6OjouP/++83zTY0UWoQAZvWPAsNjKaKZd5IPdxcWR4kryqTlV7THE8CliaL45ptvnv2V7Ozs4XVniEihULz11ltmvy5TIQgBzOd4K6vopJ/HWiCNHksRZ3yvfyhRdMVDD/KbM2eOpS9hBPDnIYD5/O249HCSyH1BNVOM8xYmBYkfn8Tro/ATCQkJx44dk6OyTqebOnVqU1PTJY658cYbDx06JMfZRwRBCGAmp7rY7kbpV+Ms9tA9liK+UiRhfyY429NPPx0ZGSlH5Q8//DApKSk4OPgSx/z2t7/985//LMfZRwRBCGAmb5ZKv44TPSz30tyUICHIhb6vQaMQ/qunp8c4QX7Tpk15eXmrVq168MEH165dS0S7du166KGH3n333eGFtrdu3frUU089/PDDX3zxhXHiIBExxlauXPnHP/5x9erVubm527dvN379nXfeMa5W2tvb++KLLxJRUVHRunXriGjdunUFBQVENHPmzIqKihMnTpj72/4pBCGAOfToaOVJ6Z44Cz9xv48X3ypFEMJ/Pffcc7W1tUT01VdfLVu2rKCgICEh4b777rv33ns/+uij5OTkt99++4UXXjAe/Nlnn0VFRaWmpr7xxhuPPfaY8YsPPfTQihUrJk6cePDgwdtuu+37778nooaGhoqKismTJxNRT0/PX/7yFyI6evToypUriejTTz/Ny8sjIlEUc3JyNm3aZIlv/b8wbg5gDp+ekmaEiFEeZn5d9Fw3xYqP5BlOdLA4HwtfCRh171jbf8x8GzIEPviaoLjor/0pU6YYM6+2tvbrr78uKioSBMHHx+ell14ydmB+/PHHRNTb25uWljZ9+vQXX3yxs7Pz7bffrqysDAsLW7ZsWUlJibFUUVFRdHS0KbMGx44dW1xczOW7u2IIQgBzeLtMem2y5eezO4l0d5z4dpn0+hTLXwwQkVvadOcxKWY7nSBe6ueemJho/CA4ODghIUEQBOPHra2tRCRJ0qOPPrp27Vo/Pz9RFFtbWwcGBqqqqnx9fcPCwoz/MS0tbWhoiIh6e3tdXU1aysjNza2np+dqvqmrhyAEkN3uRqaTaGaoVTTCfjteTPpK/9xEBZbhtgYKn0CFT6Clr+JHxpnyRufMsieijRs37tixo6KiwtnZWaPRhISESJLk4+PT1dVlMBiMx7e3txt3oggKCjLGJxG5urpKknT2IqV9fX3D29m3tLRc+oUaM8AYIYDsVpRJv4839yT6iwl1E2aGiqtOYaQQRqatrc3Hx8e4VcV7771n/GJUVNSoUaPefvttIqqqqvr666+NX58wYYJGo2lvbyciLy+v0aNHG5fzJqL6+vpDhw5NmDDB+GlhYaHFV+JGEALIq3WQfqiTlplriW1T3BMnflSBIISRWbx4sVarzczMnDx5skajMX5RFMVVq1Z9+OGHoaGht99++8KFC409ou7u7gsXLhx+C+bjjz9+8sknly9fvnPnzoyMjCeffDIhIYGIurq6jhw5smjRIkt9U0boGgWQ18qT0rWRoo+Tpa/jLHNChXv6qaCNpfhZSTMVLKakpMTYyPvXv/413DV69913/+IXvzB+nJGRYZzz7uPjU1BQUF5e7uvrGxoa+re//c2YeYmJiceOHZMkSRTFefPmzZ492/gfH3744eXLl992221ENHXq1PLy8tdee23Tpk0bN25UKn+Mnk8//fS2227z9fU16/d8HgQhgLw+rpBes7I3U0SB7hojfFxhFe/vgGUNv9IyPGhHRM7OzsO79SqVSuN2E8aPjS05IhpOr88//7y4uDg8PDw3N/fMmTM333yz8evZ2dkLFy7UaDRqtXr4v3h4eAynIBG1tbU9+eSTcn1vJkMQAsjoSAvr0tF0tdU1vP5nrJj5rf4fGQozbA4M9m369Ond3d0tLS0LFy58//33zw7Uxx9//Owjk5KSznkH54knnjDTVV4SghBARv+ukP5nrGgt78mcJdpTSPIT1tdIN5p3QyiwP6Ghoffcc48pR6anp6enp8t9PVcAzwCAXAYM9EWldNcY64tBIiL61Ti8MgNAhCAEkM/6GmlCgBBp6dVkLmZplHhQyzT9lz8SwL4hCAHksuoUWzbKeh8xVyVdFyl+UYlGoaMrKCh44IEHHnzwwT179lj6WizDep9SAJvWPkg7G6Wl0Vb9iC0bLa5CEDq2I0eOLFy4MCkpKSUl5eabb966daulr8gC8LIMgCzWnpbmhYne1jR98HyzQ4W7elhFJxvrbaX9t3Zvddm6HdV7zXa6t+a/qBJ/8mv/xRdffOyxx+6++24iGhgYeOmll+bOnWu267ESCEIAWXxWKT2YaNXNQSJSCHRzrPh5JXsqDUFoGbOiciYEJ5ntdArh3OkypaWlf/jDH4wfp6en//3vfzfbxVgPBCEAfw19rLiNLQy39iAkomWjxWU7DU+l2cCl2qUgt4AgtwALXkB/f393d7fx466uruG58w4Fdz8Af6sq2dJo0SbmqmcGCoJAh5uZpS8ELOajjz6SJIkx9uGHHzpgvyihRQgghy8qpRcybSEGiYjo57HCl6eljECbuWDgy8/PLyUlRa/X+/v7v/nmm5a+HAtAixCAs+pudqaHTbO+ZdUu5sYY8cvTDE1Ch3XPPffs3bt327Ztubm5/v7+lr4cC0AQAnD25Wm2NFpU2s6zleInOIuU34IodFw+Pj7Du8w7INt5WAFsxNrTtreA58+ihbWnMaHQEf3pT38KDw+39FVYmI09rgBWrraXVXWzmSE20y9qdGOMuKYKLUJHdO+99w5vk+SwEIQAPH1ZxZZE2VK/qFFagCAIdLwVWQiOyNaeVwDr9lW1dIOt9Ysa3RAtfIneUXBINvnEAlinxj460cFmh9pYv6jRDTHiN9VoEYIjQhACcPNdjXRNhKiyzacqI1Do0lFFJ7IQHA4m1ANw8+0Z6ZdjbTMGiQSi6yKFb8+wR5NtskVr/by9vfv6+kaNGmXpC7FJzc3N06dPl6k4ghCAjx4d7W9iX8yy1SAkouujxOeOGR5NtuFvwZpFRETs37+/t7fX0hdiqyIjI2WqjCAE4GNjrZQVLHiqLH0dV2FWqHDbTqbpJ7WrpS/FTgUHB1v6EuAC8KcfAB/fnmHXR9n2A6USaV64uKEG746CY7Ht5xbASugk2lInXRdp8w/U9ZHCt2fwvgw4Fpt/bgGswR4NG+MthLhZ+jqu2sIIcY9G6tVb+joAzAhBCMDB+hp7aA4SkbcTTQwQttejdxQciD08ugAWt7GWLYq0k1kHiyLFDbXoHQUHgiAEuFrlnaxXR8l+dhKEiyOF9TUSkhAcxwimT1RVVf3www++vr6LFy92db3w69X5+fkHDx709PScOXMmtvYAB7Ghhl0bKdhJDBKN8hI8VEJBK0v1t5vvCeBSTG0R5ubmpqenl5SUfPjhhzk5OYODg+cf88gjjyxevLiwsHDr1q3vvfce1+sEsF4baqVFEXaVGYsiBPSOguMwtUX4zDPPLF++/JFHHjEYDBMnTly7du2yZcvOPmDLli2ffvppcXFxQECADNcJYKW6dHS4mc0KtatRhkUR4pNHDctT7eqbArgYk270wcHB7du3L1myhIgUCsV11123adOmc4758ssv77zzztbW1u+++666upr7hQJYp611Ulaw4GHLC8qcb1qIUNrBtP2Wvg4AszCpRajRaBhjoaGhxk9DQ0Nzc3PPOaaqqqq4uHjPnj1xcXG//OUvX3311TvvvPNiBU+dOvX8888bP/b29r7nnnsucXadTqfT6Uy5TpCP8acg2M9AGDcbamheKJnnFjX+CCRJ9rkNAtFMNW2s0S2LlftUtken04ki2soWptPpDAaDQqG47JFKpfKyv7hMCkLGGBEN1xJF8fxHcWhoaGho6MiRI6Iobty48bbbblu2bNnFrnJoaKi9vX24+KUfbEmSzPDkw6UZfwoIwnMwoi314qMJZrpDzflTWBAmbKkTbo3Go3cu4w8bv5QsS/p/XKqZFIRqtVoQhKampujoaCJqbGwcbh0OCw0N9ff3N/6hlJWV1dnZqdFowsLCLlgwPj7+pZdeMvESh4aGnJ2dTTwYZKLT6ZydnRGE5yhsY+4qw/gAM92fkiQ5OTmZ8lfw1VsUzf56TK9ychbxMz+PKIoqlX31htsaURQNBgOvaDCpge/i4pKVlbVx40YiYoxt3rx59uzZRKTX62tqaoyZPG/evLKyMuPxpaWlrq6uQUFBXC4RwGptrmMLwu0zKMLdhSBXIb8V746C/TP1rdEnnnji1ltv1Wq1ZWVlHR0dt9xyCxGdPn167Nixzc3NAQEBt91226uvvnrHHXckJSW98847Tz/9NP5iAru3pU56KMkc7TOLWBAubK5lEwPsM+kBhpk65Dt//vydO3eqVKqZM2ceOHDA3d2diNRq9Weffebp6UlEbm5ueXl5OTk5kiStXLnyT3/6k4xXDWAFevV0pJlNV9ttTswPF7dg0VFwACNYWSYlJSUlJeXsr3h6et52221nf3rp9z8B7MnOBpYRaG8TJ842TS0Ut7GOIfJxsvSlAMgJLwEDXKEtddL8cHt+gpwVNCVY2NGARiHYOXt+jAFkZcdvygybHyZuqcP7MmDnEIQAV6Kyi/XpKdFedpy4mAURwmYEIdg7BCHAlfihns0Ns/9pleO8BVGg8k5kIdgzBCHAldhWz+aG2X0OEhHNChG21iMIwZ4hCAFGzMBoV6NkZztOXMycMGEbghDsmkM8yQB8HWlmYe5CiJulr8Ms5oaJuxolHV4dBfuFIAQYsa0O0y9KRAEuFOMpHG5GoxDsFoIQYMS2NUhzwxzo2ZkbhmFCsGcO9DADcNGrp/wWlmO/K6udb06YuA3T6sF+IQgBRmZ3I0sPENxHsDqhzZumFgpaWRe2xwY7hSAEGJlt9dIcR+oXJSIXBU0KEnY3olEI9smR/qwF4GFHA3s326xBOGQY2l69p7S1orzllJuT6zi/0anBSVPCJprzGmaFijsb2HWR5jwngJkgCAFGoHWQqntYuhm36Ntbe+Ct/I9ivKMmhqROC5lsEKVTHdUfFKxcXfr1Hyb+erRvjHkuY1aIcE8uWoRgnxCEACOwvV6aphaVZmkQSkx64eC/ytsq/zTpD2nqZCLq7+93cnKaHDbxtviffX/qh0d2PPWrlGXXjZ5vhouZGCjU9rKmfgp2NcPZAMwKQQgwAjsb2axQczQHJSb9/cA/Owa73l3wirPi3P0ARUFcPGZBZmjaH7c9MaAfvClusdzXoxAoO1jc3SjdHOtY46PgCHBPA4zA9gZzBKHE2FN7X+jV9T0/bfn5KThM7R70zzl/+6Zi49oT6+W+JCKaFSrsaMBsQrBDCEIAU9X3svZBlugrexCuLvu6faDz2Zw/qxSqSx8Z7B742pz//ax0bVFzqdxXNStU2NGIIAQ7hCAEMNW2BjYrVBRlzsHSloo1Zd/+Neshpagw5fggt4DHJ9//bO7LnYNdsl5Ykp/QNcTO9CALwd4gCAFMtUP+ftGuoe6nc198bPL9we5Bpv+vSaHps6Jy/n7gn/JdGBEJRDNCxF1oFILdQRACmGpnA5sVIm8QfnD806zwjCuYI/jr1Dtb+tp2nNkrx1UNmxUqbMeio2B3EIQAJjnVxRjRGG8Zg/BUe9We2v3/k3TbFfxfpah4MOM3b+f/e0A/wP3Chs0KFdAiBPuDIAQwyc4GNlPm5uAbRz64O+UOL2fPK/vviYFxKUGJq0q/4ntVZxvtJUhElV3IQrArCEIAk+zWsOlyBuH26j39+oFrRs29miK/mXDXNxWbGnuaeF3V+WaEoFEI9gZBCGCSXY1spmxvyhiY4YOCT+9L/5UoXNUpAtz8bxh37cdFn/O6sPMhCMH+IAgBLu9kJxOJYj3lCsJtp3erPYKSgxKuvtSNcYsP1B+p7268+lIXNCME0+rB3iAIAS5P1uagxNjnpV/fmfhzLtXcVW7Xj1m4umwdl2rnG+0lKEU6hWFCsCMIQoDL29Uo4wDhrppcV5XrhOAkXgVvGr9415l9Tb1aXgXPMV2N3lGwKwhCgMvbrZHrlVFG7NPiL3+RdAvHml5OnotGz11d9g3HmmebjmFCsC8IQoDLqOhkIlGMPAOERxqPC4I4KTSdb9mb467fdnp3z1Av37JGM0OFHQ3YmxDsB4IQ4DJ2yzlAuK5iw8/GLeJe1s/Vd1Jo+qaq7dwrE1Gsp6AShYpONArBTiAIAS5jt4bNkKdftKm3ubj5xOyoHDmKLx13zbqKDRKTJa6mq4U9GgQh2AkEIcBl7NWwaWpZgvDbk5vmx85yUbrIUTwhIM5D5X5Uc1yO4tNChD0YJgR7gSAEuJSqbmZgNMqLfxDqDLpNVdsXj57PvfKw68cuXFexQY7KeHEU7AmCEOBSdjey6fI0B3fW7BvtExPhFSZHcaPZUdOKm0/IMY9ijLdgYFTdjSwEe4AgBLiUPbItMbqpatt1Y2RsDhKRi9J5dvS0zVU75Sg+LQTDhGAnEIQAl7Jbnqn0Tb3Nle3VU0JHvO/gSC2ImbXl9A5G/BNrmlrYjd5RsAsIQoCLqutlvXo2VoY9CH84vXNWVI5KoeJe+Rzj/Ec7KZxKmsu5V56OFiHYCwQhwEXtamTT1aIcHaM/nN41P3amDIUvYF7MjM0yTCgc7yN06Vh9L7IQbB6CEOCiZOoXLWk5wYiN9x/LvfIFzYuZubtm/6BhiG9ZgSg7WESjEOwAghDgomSaQbilaueC2Nncy15MgKvfOP/R++sOca88PQTDhGAPEIQAF9bUT80DLMGXcxDqJcOumn3zYmbwLXtp82Nn/nB6F/eyOWohtwlBCDYPQQhwYXs0UrZa5D5CeFRzPNIrPMgtgHPdS8oOn1ygLe4e6uFbNsVPaOhj2n6+VQHMDUEIcGF7NSxHhn7RnTX7ZkZlcS97aa5KlzR18oH6w3zLigJNCRL2a7ETBdg2BCHAhckxQKiXDPvrDk2PmMq3rClmRGbtPLOPe9kctbgX78uAjUMQAlxA5xBVdbEJ/pyD0NgvGuDmz7esKbLCJxU2l/Tq+viWzVFj9W2weQhCgAvIbWKTggQV7+fDIv2iRq5KlwnBSfvq8viWzQgUKjpZt45vVQCzQhACXMBejZSj5vx0WLBf1EiO3lEnkdIChP14dxRsGYIQ4AL2NPJ/Uya/qSDSK8wi/aJGU8MyC7TF3HtHp6mFvRq8LwM2DEEIcK5+PRW1s8xAzkGYW5uXEzGFb80RcVO5JgfFH2rI51sW78uArUMQApzrYDNL9hPclDxrMmL76w9PDcvgWXTkssIn7a/nvMTMlGAhv5UNGPhWBTAfBCHAueToF61orXRXucq6Da8ppoZlHqg/opd4ppa7ksb7CIeb0SgEW4UgBDhXbpOUHcz50dhXn5cVPolvzSvg7+ob7hla1FzKt2x2sIDeUbBdCEKAn9BLdEjLpgbzHyDMCs/kW/PKZIVncp9EkaMWcpvwvgzYKgQhwE8ca2VRnoKfM8+aTb3atoH28f7jeBa9Ulnhk/bWHuRbM0ct7m9iBrQJwTYhCAF+IreJZXNvDtblTQ3LFAU5tvgdsVifKFEQT3ec4VgzwIVC3ISiNiQh2CQEIcBPyLHW9v76w1Oto1/UaGp4xj7eC3BjmBBsF4IQ4L8Y0b4mKZtrEA7oB8paKtKCkznWvEqTQtPzGo7yrYm9CcF2IQgB/quik7kqhAh3nkGY31QY5z/GTeXKseZVSg1KPNVe1TXUzbFmjlrY04j3ZcAmIQgB/kuOftGD9UcnhabzrXmVnBROyYEJ+ZpCjjVjPAWlKFR2oVEItgdBCPBf+5oY335RIjrceMzagpDk6R3NVmOYEGwSghDgv/ZqOL8yWtNVp5f00d4RHGtyMSVsYl7DUUY8cysrWNiHYUKwQQhCgB9p+ql9kI334RmEBxuOTg6byLEgLyEewa5K18r2ao418b4M2CgEIcCP9mqkbLUocu0ZPdSQnxmSxrMiP9x7R5N8BU0f0/ZzLAlgDghCgB/lalgW137RAf1gaUt5ujqFY02OJoelH+QahKJAU4KFfVhrDWwNghDgR9zXlCnQFo/xjbWqiRNnSwlKPNVexXef3qxgEcOEYHMQhABERD06quhk6QE8g/BI4/GJIakcC/LlrHAa7z+2QFvMsSaGCcEWIQgBiIj2a1l6gOCs4FnzsMaqg5CIJoakHm48zrHgpEChpJ316TmWBJAdghCAiChXI/HtF23rb2/pax3nN5pjTe4mqlOPcA1CZwUl+wl52KQXbAqCEIDoxzdleD4OR3IBpxsAACAASURBVDTH09TJomDVj9ho39iuwW5tXwvHmlh9G2yOVT+lAOahk+hIC+fNeI9oCiaqrbpflIhEQUhTJx/VFHCsma0WcjV4cRRsCYIQgI62sFFegrcTz5r5mkIrHyA0Slen8O0dzQ4W87RMjygE24EgBKB9vCdOnO6sUQhiqIeaY02ZZISkHtUc57jWmq8zRXgIhdikF2wHghCA9jVxnkp/tLEgI2QCx4LyCXYP8nDyqGw/zbFmdjAmUYAtQRCCozNuxpvFddOJo5oCq11Q5nzp6hS+WzJlqbH6NtgSBCE4uopO5qbkuRmvxKSi5tLU4EReBeU2ITjpWFMRx4JZwdikF2wJghAcXS7vrZfK204FuQX4uvhwrCmrtODkAm2JgRl4FYz1FJSiUNWNRiHYBgQhODruA4T5msI0dTLHgnLzcvZUuweVt1ZyrDk1WMjFbEKwEQhCcHS5vHelP9ZUNCHYloKQiNLUyceauA4TYpNesB0IQnBoTf3UMsDi+W3Gq5P0pS3lyUHxvAqax4TgZL7DhNloEYLtGEEQrl69+u67716+fLlGo7nEYWvWrHnjjTeu+sIAzCFXI2UFCxw34y1rKY/wCvN08uBW0SxSgxNLWk7oDDpeBVP8hYY+1jrIqx6AjEwNwn/961/Lly+fPn16a2trdnb24OCFb/C8vLz77rvvueee43eFADLa18Smcl1iNL+pcEJwEseC5uGucov0Ci9rreBVUCFQZqCwH5v0gi0w6VeAwWB45ZVX3nrrrTvuuOOdd97x9PT86quvzj9scHDwd7/73TPPPMP7IgHkwn0z3mOaojRbGyA0mhCclM91mHAqNukFG2FSENbV1dXW1s6cOdP46YwZM/bt23f+Yc8999ySJUvi421sdAQcVq+eyjpYRiC3IBwyDJW3nUqytQFCozR18rEm3pv0YpgQbIHSlIMaGxs9PDycnZ2NnwYGBh45cuScYwoKCr777rvDhw8fPHjwsgWPHDlyww03GD8OCAh47bXXLnFwf3+/QsF1v1QYub6+PkEQBIFn48nidjeJyb6iNNjXx6lgYUtptFckG5L6hniV/In+/n69Xi/T4zDKPbq89WRHd4eTgs/q48keVNCmauvuc7Gvx3dwcFAURZVKZekLcWg6nc5gMEjS5fveXVxcRPEyTT6TgtDFxWVoaGj406GhIVdX17MP0Ov1d9999zvvvDMclpcWFhb285//3Pixu7u7i4vLJQ7W6XSXPgDMQK/Xu7i42FkQHmpjOWrG8e4qba9IUyfJd7syxpycnGQKQhdyifGJPN1bmxKUwKkgxftIRd3OOVxnp1icIAgIQotTKBQGg8GUZ+2yKUgmBmFYWNjg4GBzc3NgYCAR1dbWhoWFnX1AVVVVQUHBHXfcQUQDAwNtbW2jRo3atm1bTEzMBQuGhITcfPPNppyaiERRNOU7AVkZfwp2FoT7tPoHEhQiv3dGC5pLliXcIN/tKv4/meqnBCUWNpdMUHN72SdbzfY3C9ND7er5lfunAKYQRZExxuunYFKVwMDAnJyczz77jIg6Ojo2bty4dOlSImpra1u3bh0RxcTEnDhxYuvWrVu3bn3llVe8vb23bt0aHh7O5RIB5KCX6JCW52a8OoOuvPVkYsB4XgXNLzU48biW5zBhVjA26QUbYFKLkIj+8Y9/LFmyZOfOnaWlpfPmzZsyZQoRnThx4mc/+xljTKVSxcbGGo+sra1VKBTDnwJYp+NtLNJD8DOpL98kpa0VMT6RbirXyx9qrZIDE57JfUln0KkUfPr9ctTiL/cYDIwUdtWVAPbG1CCcOnXqiRMnDh48GBISMmHCjxutpaWllZeXn3PkpEmTDh8+zPMaAWSQq+G8strxpqKUIJvZceKC3FSuEV5hJ9pOJgXyefE10IWCXYWSdpbshyQE6zWCDlY/P79rrrlmOAWJyMXFZezYsecc5uLiEhUVxefqAGTDfa3t49piG9p66WImBCUd5zqJAmutgfXDeC84qH1NEsep9DpJf6KVW0PKglJ4DxNmq7FbPVg7BCE4olNdTBSEaE9uQXiitSLCK8xd5caroKWkBCWUtVToJW57E2YHC3vRIgTrhiAER7RXw6ZxHSAs0Jam2vgAoZG7yi3MM6S87SSvgmO8BQNjZ3qQhWC9EITgiPbxXmK0QFts62/KDEsJSijUlnIsmBUsYpgQrBmCEBzRXq6vjEpMKm0pTwq04RmEZ0sOSijUlnAsmB2MYUKwaghCcDjNA9TUzxJ8uQXhyfaqILcAL2dPXgUtKzkovrC5VGLcJsJnqzFMCFYNQQgOZ69GygoWOE7xLtCWJHNan9Ma+Dh7+7v6VXWc4VUw1V+o7cEmvWC9EITgcHI1LFvN884v1JbaUxASUUpQQgG/3lGFQJOCsEkvWC8EITic3CaWw+9NGUasuLmU144NViI5KJ7zMKEa78uA9UIQgmMxbsY7kd9mvGc661yVrgGufrwKWoPUoMQCbQkjbtGVg2FCsGIIQnAsB7Us1V/guFVsobYkxfZXVjtHoFuAs8KprquBV8FJgUJRO+vX86oHwBOCEBzLXo3Edyp9obY0xfZXVjtfSnAix2FCNyUl+QqHmtEoBGuEIATHslfDcri+KVOgLU4KssMgTA4czzEICb2jYMUQhOBAdBIdbmZTgri1CDW9Wp2kD/cM5VXQeiQFxhc181xfJkct7sUmvWCVEITgQPJb2GgvwduJW8FCrb29Lzos0jt8QD+g7WvhVTArWDioZXpEIVgfBCE4kD0alsN5gLAk2R77RYlIICExcHxRcxmvgr7OFOUpHG9D7yhYHQQhOJC9vIOwqLnUDvYgvJikwPgirqtvY5gQrBOCEBwFI9rfJGUFc7vnOwe7tH0to3yjeRW0NslBvIcJsTchWCUEITiKknbm7yKE8Ns6t6i5LDFgvELgNyfRyoz1G9XY09Qz1Mur4PQQMVcjIQnB2iAIwVHw7xfVliYF2cnWSxekEBTj/EcXt3AbJgxxI0+VUN6BKATrgiAER7FXw3kz3sJmu9p04oKSAxP4DhNOCxH2oHcUrAyCEBxFroZND+EWhAP6waqOmji/MbwKWqfkoPgCrkGYjWFCsD4IQnAIlV2MEcV4cgvCstaKUT7RLkpnXgWtU3zAuFPtVTqDjlfBaWphVyOCEKwLghAcwh6uzUGy6xmEZ3NVukR6h59oO8mr4BhvQWJU3Y0sBCuCIASHwH0qfVFzmSMEIRElB8YXYpgQ7BqCEBzCnkbGcdMJiUllrRWJAfb8yuiwpMD4Yn7ryxBmE4L1QRCC/avvZd06FufDLQhPtZ8OdPX3cvbkVdCaJQfFFzWXSYxbdKFFCNYGQQj2b7eGTQsROXaMFjWX2uXWSxfk6+Lj7ex1prOGV8EEX6FtkDX0IQvBWiAIwf7t1fDsFyWiQm1pUqBD9IsaJQXFF/Jba00gyg4Wc9EoBKuBIAT7t5vrACERFTeX2fFa2+dLChxfpOU6TKhG7yhYEQQh2LnmAdL0s2Q/bkHY0KMhQQjxCOZV0PolB/JsERLR9BBhN2YTgtVAEIKd290oZQeLHEcIC7WlKY7UHCSiCK8wnaTjuElvqr9Q38daBnjVA7gqCEKwc7sbOU+lL2ouTXSwICSixIA4jouOKgSaEiTs0WC7erAKCEKwc7t4B2GhtjTZrjeduKCkwPG8e0dF9I6ClUAQgj1rHaS6XjbBn1sQdg52tfa3xfpE8ypoK5KCOO9WPwPDhGA1EIRgz3Y3SlnBgoLrAGFCYJwoONyDM9ZvlKZX2zXUzatgmr9Q3cNaB3nVA7hyDvc8g0PZ3cimh/C8yYuby5Idb4CQ/n+T3tKWcl4FlSJNCRL2YpgQrACCEOwZ/wFCR1pT5hzJgQlFXBcdxTAhWAkEIdittkGq7mZp/AYIBw1DVR1n7H4z3otJChzPdxsKDBOClUAQgt3ao5Gy1IKS3z1e2lLuCJvxXkxCYFxFW+WQYYhXwYkBQlU3a8MwIVgaghDs1q5GNl3N8w4v1JY6yB6EF+SqdInyCi9vq+RVUCnSZAwTghVAEILd2tnAZoZyXmI00ZHW2j5fUtD4Iq6zCWeEiLvQOwqWhiAE+9Q6SGd6eA4QSkwqaTnhUJtOnC8pkPNswpkhwo4GBCFYGIIQ7NPOBik7mOcAYWVHdaCbv7ezF7eKNig5kPMmvekBQk0PFh0FC0MQgn3a2chmhvIdICxxqK2XLsjP1dfL2ZPjJr1KkbKChd2NGCYES0IQgn3a2cBm8l9i1NGDkIiSgxL4Ljo6M1TciWFCsCgEIdghbT9p+lkKvwFCMq4pE5TAsaCNSg6ML9CWcCw4M0TYiWFCsCgEIdihHQ3SdLXIcYnR+u5GEgS1exC3ijYrOYhzEKb6C039TNPPsSTAyCAIwQ7tbOQ8caKwuTQFzUEiIgr3DGVMaurV8iooCpStFnc1YJgQLAZBCHZoVyObwXmAsMQx19q+oKTA+AKukyhmhQoYJgQLQhCCvanvZW2DLNGX7670GCD8r6TAeL7T6meFCtsxTAiWgyAEe7O1ns0OFUV+Odgx2Nkx0BntHcmtoo1LDoov5DpMmOAr9OpYdTeyECwDQQj2ZnsDm813gFBbmhgYJwo8a9q00b4xzX2tnYNdvAoKRDNDRTQKwVIQhGBvdjay2WG8BwjRL3oWURATAuP47k04G72jYDkIQrArpR1MKVCsJ88gLNCW4JXRc6QEJvDtHZ0TJmxvkJCEYBEIQrAr2+vZXK7NwV5dX113w1i/0Rxr2oGU4AS+m/RGeQieKqGkHVEIFoAgBLvCfYCwuLkszn+MSlRyrGkH4vzGVHfW9Ol4ToOfHSpsr0cQggUgCMF+GBjt1UgzeK+1nRyIftFzqRSqcX6jS1vKOdbEMCFYCoIQ7MeRZhbmLqhdedYs0GJNmQtLDkoobOY5TDgrVNyjkXRYYQbMDkEI9mNbA5vDtV90yDBU2XE6PmAcx5p2IzkovqCJZxAGuFCsp5CnRaMQzA1BCPbjhzppXjjPW7q0pTzGO8pF6cyxpt1IDBxf3nZqyDDEsebcMGEbFh0Fs0MQgp3o1lF+K8tR8504gX7Ri3JVukR5R5S3VXKsOTdM/KEOLUIwNwQh2IldjdKkQMGd69udhc0lSdiM9+JSghL4bsmUrRaK21n7IMeSAJeHIAQ7sbWezQ3jeT/rJUNZSwU2nbiE5KCEgqZijgVdFDQ1WNjViN5RMCsEIdiJH+o4T6UvbzsZ5hni4eTOsaadSQlKKGk5oZcMHGvODRO3YjYhmBeCEOxBXS9rHWSp/jyD8HhTcWpQIseC9sfTySPEI/hkO99hQuEHBCGYF4IQ7MHmOjY3jOfWS0R0XFucEowgvIyUoMTjXHtHk/yEPj2rwpZMYEYIQrAH23gvMWpghpLmExggvKzUYM5BKBDNCcW7o2BWCEKweQZG2+ql+eE8g7CirTLEI9jL2ZNjTbuUGpxY3FImMZ6vt8wPR+8omBWCEGzeoWYW7i6EunEeIEzBAKEJvJw8A90CTrZXcaw5P1zc2SAN4dVRMBcEIdi8zbXSwgjO28cXaItTMUBomlTew4QBLjTWW9jfhEYhmAmCEGzepjq2gOvKahKTiptPJGGA0DTchwmJaEG4sLkOTUIwEwQh2LaWATrZyaYG82wRnmyvCnD183Xx5ljTjqUEJRY1l/IdJlwQIW6uRYsQzARBCLZtS500M1RUcb2R8zWFE9RJPCvaNV8Xb39XX77DhJmBQn0fq+9FFoI5IAjBtm2uYwu5vi9KRMeaiiYEJ/Otad/S1MnHmoo4FlQINCdM3IJ3R8EsEIRgwyRGW3lPnDAwQ3FzGTadGJEJwcnHNDyDkIzDhOgdBbNAEIINO9LCAlyESA+uS4y2ngrxUHs7e3GsafdSgxOLmkv5Ljq6IFzc1oAN68EcEIRgw76vka6N5Nwvmq8pTAvGAOHIeDl5hnioK9pOcawZ7EqjvYRcTKIA+SEIwYZ9X8MWRXC+h/Ob8KbMlUgLTsrXFPKtuShC3FCDJiHIDkEItqqhj9X0sClBPFuEOkl/ovVkciAGCEdsgjopv4lzEF4bKXxfgxYhyG4EQbhv37558+alpqY+9thjg4Pn7iGt1Wr//Oc/T5s2LSMj47777mtqauJ6nQDn+r6GzQ8XlVz/littKY/wCsMehFcgJSixrLVCZ9BxrJkWIPTo6WQnshDkZepvkZaWlkWLFt1yyy2rVq06cODA008/fc4BFRUVAwMDzzzzzPvvv19fX7906VLOVwrwUxtq2SLeA4THm4omYIDwirir3KK8IkpbKzjWFIiuiRA24N1RkJmpQbhy5cqMjIxf/vKX8fHxL7zwwvvvv6/T/eRPv+zs7Ndee23mzJmpqamvvPLKgQMHuru7ZbhgACKiQQPtbpTmc11ZjYiOagrS1JhBeIXS1MlHNQV8ay6KEDbUYpgQ5GXq75GioqKMjAzjxxMnTmxra6uvr7/Ywfn5+aGhoZ6e2MIG5LKjgaX4Cf7OPGsO6AdOtldhD8Irlq5OyecdhHPCxDwt6xziWxXgJ5QmHqfVauPi4owfq1QqDw+Ppqam6Ojo849saGi4//77X3/99UtU2759e0xMjPHjsLCwzZs3X+Lg3t5eQeDcAwYj1dvbyxiznh/EN1XKOcHU0zPAsebhpmOjvWP0A/oe6uFYlqP+/n4nJyeFQmHpC7mwGNeIU+2ntR3NbkpXjmUnB6jWV/YtibCWduHg4KAoiiqVytIX4tB0Op3BYNDr9Zc90s3NTRQv0+QzNQi9vb17e3uNH0uS1NfX5+Pjc/5hWq12zpw5999//0033XSJalOmTHnllVd+vAKl0sPD4xIHM8YufQCYh7u7u5UEISPa1KDfeo3Cw8OFY9nS8orMsDRrvtkUCoU1ByERxQeMO9V7empYJseaP4uVtjQpbx9vLd+1SqVCEFqcMQhdXPj8BjC1azQmJqai4sdh8MrKSoVCER4efs4xra2tc+fOvfnmmx9//PFLV3Nzc4v9f5GRkSO9aHBwR1uYh4rGeXNO5SOa4+nqVL41HU26OuVII+fe0cVRwoZaLDEDMjI1CJctW7Zx48bKykoiev3115csWeLu7k5EK1eu3LZtGxG1tbXNnTs3KyvrgQceaG9vb29vNxh4rrcEMOzbM9L1UZxTsH2go7mvNc5/NN+yjmaiOvWo5jjfmqFuwlhvYY8G746CXEwNwvHjxz/55JPp6enh4eEHDx58+eWXjV9fv379gQMHiGjPnj3V1dWrV68e9f+qq6tlumhwcN9Us+ujOL8vekRzPDU4SRSwxMRVGeM3qn2gs7mvhW/Z66PEb8+gSQhyMXWMkIgefvjh3//+993d3YGBgcNfXLNmjfGDJUuWLFmyhPPVAZynsou1DLDMQM4twqOawnRMnLhqoiBMCE7K1xTOj53FseySKGHBZun1KWQVY9Rgd0b296+Li8vZKQhgft+cYddHiSLv34j5mkIMEHKRpk4+ynuttfE+gquCjreidxRkgY4gsDHfnpG494vWdNUTsUivML5lHVNGyIQjjccYcQ6txVECekdBJghCsCVN/VTczmaHcW4PHmrIzwxJ41vTYYV6qF2VrlUdZ/iWXRotfnUaLUKQBYIQbMnX1dI1EaIT79v2UGN+ZiiCkJvM0AmHGvL51pwcJHQOUVkHshD4QxCCLVl7WrohmnNzcMgwVNxclq5O4VvWkWWEpB1uPMa3pkC0JFpAoxDkgCAEm9EyQPktbAHvhbYLtCWjfGPcVW58yzqytOCkstaKfj3PBfCI6MYY8atqDBMCfwhCsBlfV0vzw0XXEUz5McmhxmOZIRM4F3VsLkqX8f5jjzUV8S2bHSxo++lUFxqFwBmCEGzGV6elG2P4TyQ7jDdlZJARMuFwI+dhQlGgJdHCWvSOAm8IQrAN7YOU18wWRnC+Y5v7WtoHOsf4jeJbFjJD0w41cB4mJKIbosWvTqN3FDhDEIJtWHdGmhMmuvPuF81ryJ8Ykipax64a9iTWJ2pAP9DQo+FbdnqIUNPLTnejUQg8IQjBNnxeKd0ayz+uDjYcmRyWzr0sCCRMCk0/UH+Yb1mFQDfGiKurEITAE4IQbIC2n462sGt494vqDLpjTUWTQhCEspgSNvFA/RHuZW+NFT89id5R4AlBCDZgdZV0XST/90WPNRXFeEd5OXtyrgtERJQRMqG0pbxP18+3bJZa6DNQcTsahcANghBswOeV0q2j+N+rBxoOTwmbyL0sGLkoXeIDxnHfnlAgujlG+LwSjULgBkEI1u5MD6vqZnNCZRggrD86NSyDe1kYJlfv6ChxdSVDkxB4QRCCtVtVyW6MEZW8b9XTnTUGZojxieJcF84yNSxzf/0hiXdmpfoLLgo6pEUUAh8IQrB2n1dKt8byv1EP1h+ZGpbJvSycLcQj2MvZ62RbJffKt4wSP0PvKHCCIASrlt/CunWUpZajXxQTJ8xhaljGft6TKIjojtHC55XSEKIQeEAQglX7z0npF2O4b0dPnYNdpzpOpwUn8y4M55oalpFbl8e9bLSnEO8jbKpFEgIHCEKwXnqJ1lRJt4/m3xzcV3coI2SCk8KJe2U4R2JgfNtAO/clZojorrHif05imBA4QBCC9dpQK431FkZ58Q/CvbUHc8Incy8L5xMFQaZG4c0x4q5GqZnzXk/giBCEYL0+OcnuGsv/Fh3QDxZoiyeFYoDQTHIiJufWHuRe1kNF10SIq/HKDFw1BCFYqdZB2tEg3RjD/xbNaziaEBDn4eTOvTJcUHpwSmVHdVt/O/fKvxgj/gfLrcFVQxCClVp5UrouUvRS8a+8t/ZgTgT6Rc1HpVBlhEw40MB/Zv2sUKFlgI63YqQQrgqCEKzUB+XSr+P43596yZDXeDQrfBL3ynAJORFT9tYe4F5WFOiX48T3y9EohKuCIARrlKthBkbZMkwfPN5UFOkV5u/qy70yXMLk0PSi5rJeXR/3yr8aK6yulPr03AuDA0EQgjV6v1y6J47/9EEi2nFm74zIbBkKw6W4q9xSghL21R3iXjnMXcgKFtdUoVEIVw5BCFanc4i+OyPdPlqWftHcurxpEVO4V4bLmhWVs/NMrhyVfx0noHcUrgaCEKzOylPSwggx0IV/5SOaY1HeEcHugfxLw+VkhU8q0BZ3DXVzr3xNhFjTgx0K4cohCMG6MKIVpdJvxstyZ+44kzsT/aIW4qp0SVenyNE7qhDo7nHiilI0CuEKIQjBumyrZwqBpsnwmozOoNtfd2haJPpFLWZWVM6OM3vlqHzvePHzSql9UI7aYP8QhGBd3iiRHkyU5bbMa8wf7RsT4OonR3EwxZSwjNKW8o7BTu6V1a50TYT4MSbXwxVBEIIVqe5mB7TSraNk6hfdOzMK/aKW5KJ0zgxJ2yvDcmtE9IcE8a1SScJAIYwcghCsyBul0i/Him5K/pV7dX15DUdnRGbxLw0jMTdm+g+nd8lReXKQ4OdMm+qQhDBiCEKwFr16+uSk9Ft5XpPZXbM/LTjZ29lLjuJgukmh6XXdDXLsykRE98WL/yo2yFEZ7BuCEKzF+yek2aFitKcc0+hpS9WOeTEz5agMI6IQFLOisn84vVOO4reOEk900jEsPQojhCAEq6CX6J/F0h/leU2mqVd7urNmMvZdsg7zYmZurtrBiH9cqUS6L158pQivzMDIIAjBKnxRJY32oklBsjQHN1ftnB2do1LIsJMFjNw4v9GuSpfi5hNyFP/NeHFLnXSmB41CGAEEIViFV4qkR5MVMhXfVr1rfswsmYrDFZgbM0Om3lFPFf1qnPjPYjQKYQQQhGB5P9QzA6N54bI0Bwu1JaIgxvmPkaM4XJl50TN2ndk3oB+Qo/j9CeInJ6VWTK4HkyEIwfKezTc8liLLXhNE9N3JLdeNXiBPbbhCAW7+SUHjd8izBneom3BjjPhaEV4fBVMhCMHCttYz7QD9PFaWW7FrsPtgw5F5sTPkKA5XY/GYBd+d3CxT8eWp4ttlUossDU6wQwhCsLDnjhmeThMV8rQHN1dtnxqe6eXkKUt1uAqZIentAx0VbZVyFI/0EG6IEV8vQaMQTIIgBEvaWs80/XI1B4loQ+XWxegXtUqiICwaPff7Uz/IVP+JVPGdMqkNI4VgAgQhWNIz+TI2B483FTGixMA4WarDVbt21LwdZ/b26frlKB7pISyNFl/BSCGYAEEIFvPtGalbR7fI1hxcV7Fx6dhrZCoOV8/P1TdNnbzl9A6Z6j81QXy3TKrrxZxCuAwEIViGgdHyI9KLmQqZ3hbV9GqPNRUtiJ0tS3Xg5Ka469eeWC8xWbIqzF345TjxuWOYUwiXgSAEy/ioXAp0ofnyzB0korUnvls0aq6r0kWm+sBFUuB4L2ePA/WHZar/l1TFN2ek0g40CuFSEIRgAX16evaY9NIkuZaS6dP1bzm9c+m4RTLVB45ujLt+zYlvZSru40SPJiuWH0ajEC4FQQgW8I8CwzS1MDFArubg96e2ZIakBbkFyFQfOJoRObWxR3Oi9aRM9e+LFwva2I4GNArhohCEYG5V3eztMunFTLnuPb1k+Kr8+5viFstUH/hSCIqlYxetPbFepvouCnp1sviH/QYdmoVwEQhCMLcHD0iPJivC3OVqDv5wemeYZwgWF7Uhi8csONx4rLarXqb6S6LEaE96sxRJCBeGIASz2ljLKjrZg/LsO0hEEpNWlX51V9ItMtUHObir3JaOu+az0q/kO8VrkxV/P27QyDJlEWweghDMp1dP9+03vDFV4STbfbf19C5/V7+UoAS5TgDyuCnu+v11h+q6G2SqP9ZbuHuc+MABzK+HC0AQgvn85bBhRogwN0yuTlGJSZ+WfHlX4s9lqg/ycVe5LRm78PPSr+U7xVNpiuI29nU1OkjhXAhCMJODWvZVNXtZtikTRLStere3s3eaOlm+U4B8boxbvKf2N1vnKgAAFktJREFUQEOPRqb6zgr6YJri/gNSOxYghZ9CEII5DBjol3sMr08W/ZzlOsWQYeiDgs/unXCnXCcAmXk5ed4cd/17xz+R7xRTgoSlUcJDeegghZ9AEII5PJpnSPEXboiR8X5bW74+zn90UmC8fKcAud08fklpS3lx8wn5TvH3DMW+Jrb2NDpI4b8QhCC7TbXsuxq2YqqMnaLdQz1ryr69JxXNQdvmrHD6n6RbV+R/xEiu+e8eKvp8puL3+w1nejDFHn6EIAR5afrpV3v1n81Q+MrWKUpE/y78fGZUVrhnqIznALOYHztrQD+QW3tQvlOkBwh/TFTctdtgQBQCESEIQVZ6iW7Zof/NeEW2Wq43RYnoZHvVjjN7f5F0q3ynALMRBfG+9LvfOPrhgH5AvrP8KVlUCvTUUQwWAhGCEGT1cJ7BWaTlqTLeZhJjrx56+97UO72dveQ7C5hTmjo5JSjh34Wfy3cKUaA1s5WfV2KwEIgQhCCfVZXSpjr2xWylTBvQG31TsVEpKheMwr6DduUP6Xf/UL3rZHuVfKfwc6av5ih+v99Q0o4eUkeHIARZ7Gtifzxo+HqOwsdJxrNo+1r+U7T6kczfCSRn2ILZeTl7/ip52St5KwxMxt7LVH/hhQzFkq2GZhl7YcEGIAiBv1Nd7Obthk+mKxN9ZcwnibG/H/jnTXGLo7wj5DsLWMqi0XM9nNxXFq+R9Sy/GCveNkq4ZrO+Vy/recCqIQiBs8Y+mrfJ8HyGKN/u80ZflK3TGXS3Jdwg61nAUgQSHp/ywDcVm0paZJxWSERPpysS/YRbdxj0GC50VAhC4EnbT7M36u+NE+8aI++tdaq9anXpuieyHhYF3MN2K8DV7+HM3/7vvlf7dDJuGyEQvZetYMRu34UJFQ4Kv0SAm5YBmrNJf8so8bEUee+rrqHuJ/e+cP/EX6vdg2Q9EVhcTsSUiSGp/zj4unxT7IlIJdJXc5S9eoZ2oWNCEAIf9b1s5gb99VHCkxPkvakkJv3vvlenRUyZHT1N1hOBlXhg4r0dA51yDxY6ifTlbGXHELt9l2EIWehgEITAQWUXm7HBcEOM8Fy6jOuoGb197GOJSVhNzXEoRcXT2X9af3LLvrpDsp7IRUHr5ykZ0YJN+i6drKcC64IghKu1u5Flr9f/dYL4dJrsKbim7JuD9Ueeyn4UQ4MOxc/V97lpf34p742i5jJZT+SsoFUzFeN8hJkb9LW9GDB0FPhtAlflrVLplh36T2cq75T57Rgi2ly1Y235+pdnPe3p5CH3ucDaxPmPeTLr0b/ueb687ZSsJ1II9HaW4tZR4uRvDbkaZKFDQBDCFerW0Z27DO+dkPYtVs4OlX0++w+nd71fsPLV2c8F4wUZR5WmTn4483d/3vVcZUe13Od6JEn893TFjdv1rxRJCEO7hyCEK5GnZRPW6V2VdGCxMtZT9hRce2L9BwWfvjrrWewv4eByIqY8MPHeh7c/Wagtkftc88KEvOuV66qlBZv0jX1ynw0sCUEII9Orp0fyDEu26l/KFN/NVrgp5T2dxNi7x/7z3anNb877B1aQASKaHjn1yaxHntz7jz21B+Q+V5SHsGuRcmqwmLpO92E5moZ2C0EII7CuWkr6Sq/tp6IbVEujZb95uga7H9v1TFlrxZtz/xHkFiD36cBWpKmTX5r1zJtHP3z32H8kJu9cB6VIT6WJWxcqPyiXZm7QF7QhDe0QghBMcrRNnLHB8HS+9EGO4pMZigAX2c9Y1Fx6z+aHYr2jXpn9rJezp+znA5syxjf2/YWvVrRXPrzjyabeZrlPl+wn7LtOeUusuGCT/t79Qj16Su0LghAuY38Tu2aL/o59qjtGC/lLlbPkfy9mQD/wryPvP5370v0T7/lt2v8oBNlnZYAt8nb2emnm0xPVqfdseui7k5tlXXqGiESBfjNeLL9ZFexKE7+j3+4zVHejdWgnEIRwYUMSfXZKmvyd/s7dhqVR4vFFg78aJ8q6syARSYxtrtpxx/rf9ep6P170xtSwDHnPBzZOFMRlCTe+Pvf5TVXbf7fl0aLmUrnP6KWiZyewoiXk70wTv9HfuN2wqxFxaPNkftUBbNDhZrbylLS6UkrxF/6SIl4bKYoC9fTIe1KJSbtr9n9a8qWr0uWZnMfiA8bJez6wI9HeESvmv7itevdz+14d7Ru9LOHGhIA4Wc8Y4EL/O1HxeIrik5PS7/cZDIzuGCPePlqI8sC+mDZpZEGo0+lUKtXVHwPWZkiifU3s+xrp62qmEun20WLe9coY+edFEFH7QOfW0zvXVWz0d/X7VcrtaAXCFRBImBs9Y3rE1I2V2/5336sBbv7Xj1mQEzHFWSHjxtAeKvpdvPi7ePGgln1yUpr4jSHaQ7ghRlwYIST7CYhEGyIwZlK7vqen56677tq6datCoXjssccef/zx84958cUX//73vxsMhtmzZ3/yySeenhd+wWHVqlUbNmz47LPPTLzE7u7ui5WCqzFooPxWtlfD9mqkvRoW5yNcEyEujRKS/C7wCPf09Li7uwv8nu7Owa799Yf31h4o0JZkh09aPGaB3H/F24H+/n4nJyeFAoOmlyIxaU/tgQ2ntp5oPZkTMTknYkq6OtmJXyIODg6Konj+n/sGRrsb2TdnpC11rFvHZoeKOWohWy3EeQsiUpE3nU5nMBhcXPi8tmdqEP7lL385evTo+vXrGxoaJk2atG7duqlTp559QF5e3qJFi/Ly8iIiIpYuXZqYmPjCCy9csBSC0FJqe1lFJxW3seJ2lt/KTnSw8T5CtlqYphZmhIh+zpf6v1yCsH2g80TryQJt8fGm4pquuoyQCVnhk7LDJ7mpXK+mrONAEI5Ic1/Lrpr9uXV5J9sq4wPGTQhOSgocP8ZvlKvyqn57XiwIz3a6m+1qZHs0bF8Ta+pjEwKEVH8h0VdI9BXGeguXftbAFJYJwrCwsI8//nju3LlE9OCDD/b19b333ntnH/Db3/5WoVC8+eabRLRz585bbrmlqanpgqUQhPLRS9Q6SC0DrKmfGvpYQx/V97LT3VTdw6q6mKeK4nyEBF8hyU9I9RNS/AUXk3+jjjQIhwxDzX2tml5tY4+mrrvxdEfN6c4zfbr+OP8xiYFxacHJ4wPGqUQMUY8MgvDK9Az1FmhLjjUVlrZUVHZUq90Do70jo7wjwj1DQjzUIe5Bvq4+pr+cbEoQnq1jiI40s4I2VtzOStpZRSdTChTrJcR4ChHuFOEuhLpTiKsQ6EqBLshIU/ENQpN+E/X19TU0NMTHxxs/jY+PX7169TnHnDp1aunSpcMHaLXarq4uLy+vCxbU6XTt7e3Gj0VR9Pb2vtTZBwfb+21v2o6eUa/+8n9k9OrpnI1Au3Vk/ONEYtStZ0RkkKhHT0TUpSNi1Kljeol6ddSrpz4D9Q6xLh116ahjiPXqyc+JfF2EQGcKdBWCXSnWVcgJp3B3IcJd8PjpkzvYQ4NERNSn6zdcZFZyn77fwAxENNA/IKgEPTMQUc9QHyPWrx/QM0O3rndAP9CvH+zW9fbo+rqGujsHu9sHOwcMQwEuPkGuAaHuQeGeIYujZkR5hoeevUzowAA2fRspNjAg6YcEBOEIuRFN8Yuf4hdPRAZmONPdcKarrrq7/mDNIU2ftqmvpXOo28vJw9vJ08vJ09vZ00Pl5qlyd1W6uCid3ZSuzgqVk+ikElUuSiciEiVBpVAplUoiclI4OYsX7nR1Vjg5KVREJBBleFCGB1Hkj//UOkC1PayulzX0U20LO1xDzQOsbZCaB1ifnnxU5OUseCrJ04k8VYKbkpxF8nYmBZGHSlAK5KYiIvJ2IiISiTxVP/55qlKQ609vDRcFOZtws3iqbGkzF3cXZ38vH741TQpCY2h5ePy45L+Xl1dbW9v5xwwfYGzAtbe3XywIt2/fHhsba/w4IiJi//79lzj7l/u+3dD5nSnXaa8EIrrkHCkPIg+iUMX/tXf3MU1dbQDAn34wKFJl09YiipQ4hfmBCOo0KlIVJZNoFOWdqXtBJNONaAK6mfgRE030jb4a4xa3GdxclhmUuWU1Jmb4BcRMZENI2FoLJVKgIIW2Ugq0vT3vH2e76UsR69YP1vv8/jo9fXrvo/b69J577j0A9Es/BMwQdD6HToB6L7YvcvH4L9i8yAUC8sfuwwlP6CK0k08ggvAELhjHQDiBCAZed/HGMTwxA1FOmMDwxE4AMAGYALTs1gze/GkR8jMRQCKA+xVpF49nEfb3CazPBR39Ql6/gPTzwSIgXXzeII/Y+eDgEScfhngAAIN8cP05MmLngf0FNcTOJw6vryQIACQAEvYQdgA4wDnwxyE0IgIcvfCYYJf951//pWeETqfzpfGRkZF8/ksKvVeFcNKkSTwez2Kx0FM3k8k0efLkYTESicRisdC22WymPS/a4Lp167wfGv23YmuRuMDLYOQnPp8sg/4CHBr1H++fY/uqQ6PIH3w7NOrVCXF4eHhCQkJdXR19+fjx48TE4RP8kpKS3APi4+MjIyN9kiJCCCHkP96ODBcWFh4/flyv11dWVl69enXnzp0A0NPTs3bt2u7ubgAoKCj47rvv7t6929bWduzYscLCQj9mjRBCCPmIt9P2SkpKjEbj8uXLo6OjL1y4MHfuXAAghAwNDdF5p7Nnz/7iiy+Ki4tNJtPmzZv379/vx6wRQgghH/H29gkfetXbJ9RqdUJCwmuv+fEJEeilmpubZTLZuHHjgp0Ip+n1erFYHB3t4ylz6JUYDAahUDjKHAgUAEaj0W63T5nim5W6/wGTZnNzc3/7ze/P0kWj2717d2VlZbCz4LqDBw9ev3492Flw3ZkzZy5evBjsLLjuq6++OnnypK+29g8ohAghhJD/YCFECCHEaVgIEUIIcVoQJsvU1NTk5eUNDQ15GW8ymcRiMX2gEQoWi8UiEolwylJw9fX1hYWF+eomYvTXWK1WgUAgEuGT4oNpYGCAYRj2cWajuHHjRlJS0ugxQSiEAKDX6x0OR+D3ixBCiFOmTp360l/wwSmECCGE0BiB1wgRQghxGhZChBBCnIaFECGEEKdhIUQIIcRpY+ueBI1G88033zidzm3bttHneg9jtVovXryo1+uXLVu2adOmwGcY8np7e1UqVWNjo1gs3rRp0+zZsz1jSktLGYah7VmzZqWnpwc2x9Cn0+kqKirYl+vXr/d8pqLD4SgtLX3y5ElycvL27dtfuvQoelUNDQ0///yze49SqRy2utzVq1fp8qsAEBMTk52dHbj8QlpLS0ttba3JZMrNzaXr4FK//PJLWVmZSCTKz8+Pj4/3/GB3d3dpaWl3d/c777yjUCi83N0YOniampoWL15MCImKilq2bFl9/fDF1QkhmZmZ9+7dmzFjxoEDB06dOhWUPEPbvn37fvzxR6lUarFYFi1a9NNPP3nGfPjhhw0NDTqdTqfT0UW4kG/V1taeOHFC96eBgQHPGKVSeeXKlTfffPP8+fN79+4NfJIhz2w2s/8EP/zww+HDhz0X4z169GhVVRWNMRgMQckz9HR1daWmpn7++efvv//+s2fP2P4HDx4oFIpJkybZbLa0tLT29vZhH7TZbEuWLNFoNNOnT9+2bVtZWZm3uyRjxp49e3bu3EnbH3300XvvvTcsoKKiIjY21m63E0Kqq6ulUildBAr50MDAANvet29fTk6OZ0x4eHhHR0cAk+KcsrKy1atXjxKg0WhEIpHJZCKEPH36NCIioqurK1DZcVFOTk5JSYlnf1JSUnV1deDzCW0MwxBCnE4nADx58oTt37Bhw7Fjx2g7Nzf30KFDwz546dKlhQsXulwuQsi33347b948L/c4hs4I79+/n5mZSdtr1qy5f/++Z8DKlSvpj7IlS5bYbLbff/890FmGOvenlgwODr7owQ1ffvnluXPnHj58GKi8OMdgMJw+fbq0tLSzs9Pz3aqqqrS0NLoeU1xcnFwuHzaIh3yop6dHpVLt2LFjxHe///77M2fOuA9lo7/pReP8lZWVa9asoe0RawQN4PF4NKChoaG3t9erPf6NbH3MYDCwS3xJpdLOzk7y/zf7d3Z2sgF8Pl8ikXR0dAQ6S86or6+/fPlycXGx51vp6el9fX1arTYrK+vIkSOBzy3kicXi5ORks9lMnw5VW1s7LMD9WAAAqVSKx4L/fP311ykpKW+99ZbnWykpKQDQ0dGRl5enVCoDnhqHDA4Omkwm9xrhORbtXkQmTpwoFAq9HK8eQ5NlwsLC6LkwADidTqFQSAs7SygUsnM0AMDhcOCjL/2ktbV148aNZ8+eHXHK0q1bt2ijoKAgLS2tqKhIKpUGNsEQl5WVlZWVRdslJSVHjhy5efOmewAeC4F0+fLloqKiEd9iFxgvLi6eOXNmTU3NokWLApgahwgEAj6f714jPL/zQqGQDWAYhmEYL4+LMXRGGBsby/6qbW9vj42N9Qxgr47a7Xaj0eir5YmRu9bW1pUrV3788ccFBQWjR6akpERERDx9+jQwiXHT0qVLdTrdsE73YwEA2tvb8Vjwk5qamqampq1bt44eNmXKFLlc3tLSEpisOCgsLEwikbBf+xG/8+5FhDZiYmK82fgYKoTZ2dnXrl2j7WvXrrETkauqqnp6emjAnTt36JivSqWaNm1aYmJisLINVW1tbQqF4oMPPti1a5d7f11dXWtrKwAMDg6ynbdv32YYZsaMGYHOMtSxf8mEkBs3bsyZM4e+fPToEf2PIDMzs7Gxsbm5mXZaLJYVK1YEK9vQdunSpS1btowfP57tUavVGo0GAOx2u8vlYju1Wu2ItxshX8nOzi4vLwcAQkh5eTmtES6X686dO/39/TRApVLRw6e8vFyhUHizPAXAWJo1ajQaZ86cuXbt2g0bNsTFxbW1tdH+N954Q6VS0XZ+fn5iYmJeXp5EIrl+/Xrwkg1ZW7ZsiYiISP2TUqmk/enp6SdOnCCElJWVzZo16913312/fn1UVNRnn30W1HxD08aNGzMyMpRK5fz58+VyuVarpf0pKSmffvopbR89ejQuLm7Hjh0ymeyTTz4JXrKhzGazRUdHV1VVuXfm5+cXFhYSQh49ehQXF5eTk7N58+bx48cfPHgwSGmGoFWrVi1YsAAA5syZk5qaarVaCSFNTU0ymSwnJycjIyM5Ofn58+eEEKvVCgANDQ2EEKfTmZmZmZqaun379okTJz548MDL3Y2t1SdsNltFRQXDMKtXrxaLxbSzrq5OLpfTCXIAUF1drdfr3377bblcHrxMQ1ZzczN7gzAAREZG0qW81Gr1hAkTYmJinE7n48ePtVptVFRUWlqalyMP6JWYzeaHDx/29vbGxMQsXbqUvc7R2NgokUjYK7K//vqrRqOZN28enoj4SX9/v1qtXrBggft8hZaWFh6PFx8f73K5Ghsb1Wq1UChMTk5OSEgIYqohpr6+nr3aBwDz588XCAQAYDabKyoqIiMjV61aFR4eDgAul6u2tnbu3Ll0hUiGYe7du2c0GtPT02UymZe7G1uFECGEEAqwMXSNECGEEAo8LIQIIYQ4DQshQgghTsNCiBBCiNOwECKEEOI0LIQIIYQ4DQshQgghTsNCiBBCiNOwECKEEOI0LIQIIYQ4DQshQgghTvsfEBP1NNMeHtsAAAAASUVORK5CYII=",
"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 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.321459555244377e-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.00022342332015932995"
},
"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.7.3"
},
"kernelspec": {
"name": "julia-1.7",
"display_name": "Julia 1.7.3",
"language": "julia"
}
},
"nbformat": 4
}