{
"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.620254e+02 1.447057e+02\n",
" * time: 0.0006029605865478516\n",
" 1 1.387679e+02 1.079574e+02\n",
" * time: 0.001878976821899414\n",
" 2 1.005616e+02 1.217431e+02\n",
" * time: 0.0030519962310791016\n",
" 3 2.097994e+01 4.182211e+01\n",
" * time: 0.004261016845703125\n",
" 4 1.149958e+01 2.748740e+01\n",
" * time: 0.005278825759887695\n",
" 5 9.125888e+00 9.082666e+00\n",
" * time: 0.006073951721191406\n",
" 6 7.043560e+00 1.917138e+01\n",
" * time: 0.0068929195404052734\n",
" 7 3.988754e+00 5.945317e+00\n",
" * time: 0.007701873779296875\n",
" 8 2.726113e+00 3.024515e+00\n",
" * time: 0.008566856384277344\n",
" 9 1.579391e+00 2.833376e+00\n",
" * time: 0.009402036666870117\n",
" 10 1.413821e+00 2.784967e+00\n",
" * time: 0.010223865509033203\n",
" 11 1.316640e+00 2.903979e+00\n",
" * time: 0.010879993438720703\n",
" 12 1.225098e+00 1.191196e+00\n",
" * time: 0.011534929275512695\n",
" 13 1.183430e+00 5.483881e-01\n",
" * time: 0.01220703125\n",
" 14 1.157298e+00 2.956812e-01\n",
" * time: 0.012887954711914062\n",
" 15 1.148557e+00 2.021497e-01\n",
" * time: 0.013510942459106445\n",
" 16 1.146062e+00 1.049308e-01\n",
" * time: 0.014078855514526367\n",
" 17 1.144636e+00 1.088491e-01\n",
" * time: 0.014629840850830078\n",
" 18 1.144163e+00 3.123473e-02\n",
" * time: 0.015193939208984375\n",
" 19 1.144122e+00 2.765210e-02\n",
" * time: 0.01575303077697754\n",
" 20 1.144109e+00 4.863846e-02\n",
" * time: 0.016216039657592773\n",
" 21 1.144064e+00 1.773289e-02\n",
" * time: 0.01669001579284668\n",
" 22 1.144058e+00 1.421626e-02\n",
" * time: 0.017282962799072266\n",
" 23 1.144051e+00 1.498466e-02\n",
" * time: 0.017679929733276367\n",
" 24 1.144044e+00 9.733992e-03\n",
" * time: 0.018204927444458008\n",
" 25 1.144038e+00 3.650658e-03\n",
" * time: 0.018730878829956055\n",
" 26 1.144037e+00 3.983450e-03\n",
" * time: 0.019284963607788086\n",
" 27 1.144037e+00 3.426817e-03\n",
" * time: 0.019841909408569336\n",
" 28 1.144037e+00 1.854724e-03\n",
" * time: 0.020424842834472656\n",
" 29 1.144037e+00 1.374631e-03\n",
" * time: 0.020973920822143555\n",
" 30 1.144037e+00 5.714051e-04\n",
" * time: 0.021512985229492188\n",
" 31 1.144037e+00 3.041005e-04\n",
" * time: 0.02206897735595703\n",
" 32 1.144037e+00 3.243686e-04\n",
" * time: 0.022650957107543945\n",
" 33 1.144037e+00 3.128633e-04\n",
" * time: 0.023231029510498047\n",
" 34 1.144037e+00 2.327409e-04\n",
" * time: 0.023794889450073242\n",
" 35 1.144037e+00 1.546077e-04\n",
" * time: 0.02422499656677246\n",
" 36 1.144037e+00 5.593509e-05\n",
" * time: 0.02480602264404297\n",
" 37 1.144037e+00 2.952055e-05\n",
" * time: 0.025367021560668945\n",
" 38 1.144037e+00 1.658163e-05\n",
" * time: 0.02593088150024414\n",
" 39 1.144037e+00 9.627992e-06\n",
" * time: 0.026556968688964844\n",
" 40 1.144037e+00 7.030475e-06\n",
" * time: 0.027151823043823242\n",
" 41 1.144037e+00 5.072021e-06\n",
" * time: 0.02775287628173828\n",
" 42 1.144037e+00 3.266220e-06\n",
" * time: 0.02842998504638672\n",
" 43 1.144037e+00 1.813561e-06\n",
" * time: 0.029094934463500977\n",
" 44 1.144037e+00 1.296418e-06\n",
" * time: 0.029744863510131836\n",
" 45 1.144037e+00 1.452421e-06\n",
" * time: 0.030373811721801758\n",
" 46 1.144037e+00 1.276097e-06\n",
" * time: 0.03082895278930664\n",
" 47 1.144037e+00 5.119396e-07\n",
" * time: 0.03144383430480957\n",
" 48 1.144037e+00 3.767626e-07\n",
" * time: 0.032013893127441406\n",
" 49 1.144037e+00 1.981879e-07\n",
" * time: 0.032629966735839844\n",
" 50 1.144037e+00 1.486966e-07\n",
" * time: 0.03319883346557617\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown:\n Kinetic 0.2682057 \n ExternalFromReal 0.4707475 \n PowerNonlinearity 0.4050836 \n\n total 1.144036852755 \n"
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"Ecut = 500\n",
"basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))\n",
"scfres = direct_minimization(basis, tol=1e-8) # This is a constrained preconditioned LBFGS\n",
"scfres.energies"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"## Internals\n",
"We use the opportunity to explore some of DFTK internals.\n",
"\n",
"Extract the converged density and the obtained wave function:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ρ = real(scfres.ρ)[:, 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": "1.3777528460094525e-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+gvaeTAAAgAElEQVR4nOzdd3xV5f0H8O85997sQfbNHgTIXmySAGEjglC3Alq31qp1VK3WirZaB1arrbbanxYc1IWIIHsmbALZJCEDsm5u9h733vP8/riWIjOB59z5ef8VwuF7Tl6ck899nuc8zyMwxggAAMBeiea+AAAAAHNCEAIAgF1DEAIAgF1DEAIAgF1DEAIAgF1DEAIAgF1DEAIAgF1DEAIAgF1DEAIAgF1DEAIAgF0zQxC2trZWVVUN/XiDwSDfxcAQSZJk7ksAkiQJayKaHWMM/wtmxxjj+EvJDEG4adOm559/fujH9/b2yncxMES9vb14+M1uYGAAn0jMbnBwUK/Xm/sq7J1erx8cHORVDV2jAABg1xCEAABg1xCEAABg1xCEAABg1xCEAABg1xCEAABg1xCEAABg1xCEAABg15TmvgAAuAxGlN/KyjpYWavo6sCiPVmiN4W7Cea+LgAbgSAEsFytA/RRqfR/pRIjSvASQpxJ309b6w1Hm1nMCOGeMeItUaIS3ToAVwdBCGChvqyUHjtgmBcifjJNMclfIKK+vkEHB6VCodBJ9MNp6d0i6a0C6V9TFak+aB1agZaWlocffhiLJ1+xp556avz48XJURhACWJwePd2x21DaztbOUk70v0DIqURaEiEuiRBXlUvzN+kfiVf8LgUNQ0vX0NCQk5OzcuVKc1+IVfrggw/y8/MRhAB2oWWAFmzWx48Qji5ROlwu3ZaPEueEiNdu1p/uZn9LVyjQMrRsnp6eN954o7mvwipt3rxZvuL4FAlgQep7WcZ6/YxA4aOpisumoJHamXYuUFZ0sZt3GPTYmgJg+BCEAJaiW0fXbjYsixZfGT+8pp27ijbMVfbr2a/2YfwJYNgQhAAWQWK0dJchzVe4stE+B5G+nKk83sJez0erEGB4EIQAFuHJg4Y+PfsgXXHFFVyU9O0sxd+KpXWnkIUAw4AgBDC/H2vYt9XsPzOVVzkpMNhV+HKG4v5sQ10P43RpALYPQQhgZs39dG+24eNpihEOHKpN9BceilP8co8BSQgwRAhCADO7e69habSQFcht6sPvksUePb1diA5SgCFBEAKY038qpVNd7KWxVz40eD6lSKunK145bqjuQrMQ4PIQhABm062jJw9Kf50y1CmDQxflLjyWoPjNATQKQS4vvvjiunXrjF8PDg4uW7asv7//YgdLknTHHXd0dnaa6uqGB0EIYDYvHTPMChamqmVZD+apJLGknW2oQaMQZFFSUtLQ0GD8+oMPPvDx8XFycrrYwaIoxsTEvP7666a6uuFBEAKYR0k7+6RM+vN4np2iZ3MQ6Z3Jisf2G/oxyR4uLicnp7Ozc8+ePWvWrDF+5/Dhw59//nleXt6ZY9rb2zdu3Lh69er8/PzzKzDG3n333bvvvpuIOjs7c3JyiKiysrKsrIyIDh061NLSQkR33HHHhx9+eIlWoxkhCAHM4+lD0nMpigBnGU8xN0RI8BL+XowOUrio5cuXL1q06J133snJyZEk6eabb37ssceOHj26fPnyF154wXjM9ddfv27duuPHj994440vv/zyORUKCgq6u7sTExOJ6MSJE3fccQcRrVq16v333yeihx566OjRo0QUFBQUGBiYnZ1t0h9vaLDoNoAZHGpiBW3sq5lyNQfP+PMEcfoP+ntiRA+V3KeCK/GfSmlHvem6r9+bolCd1/zJzMw0xtvnn39eV1eXnZ0tCEJXV9eoUaPuvffe0NDQ7du3G4985plnIiMjn3nmGZXqf/dTbm5ufHz8UM6ekJBw5MiRWbNm8flh+EEQApjBs4cNf0gVHWXPQRrjKcwNEd8ulF5IRfePJQp1Fcb6mu504oXGo2fPnm38Yu/evb29vQ888MCZvyosLAwNDd24ceN777136tSp/v7+/v7+urq6iIiIM8e0tbV5eHgM5eyenp6tra1Xc/0yQRACmNqWOlbbQ0ujTZRMK8aKY9fqH4gR/eXshoUrMyVAmBJg5t2zXFxcjF8YDIb4+PgzG0XdeOONSUlJFRUVd9555/r169PS0pRKpbu7u06nO/ufe3p6dnV1Gb9WqVTn/K1OpzvTfOzs7AwJCZH3h7ki+JAIYGq/P2L40zjxKldTG7pwN+GWkeKbBXhnBi4jIyOjqKho+vTps2bNmjVr1syZM/38/MrLy8PCwiZOnKhSqXbv3t3T03POv0pJSSkqKjJ+HRUVpdVqa2trjX9samqqqKiIjo42/rGoqCg1NdVkP87QoUUIYFLb6liPnq6PNOln0GeSxZRv9c8mK7wcTXlasDJLly7duHHj+PHjZ82a1dvbu3Xr1n379k2cOFGj0SxfvtzX1/fw4cPn94KmpKQ4ODiUlJTExsZ6enquWLEiPT3d29t7cHDw+++//81vfhMaGkpEWq22pqZm6tSp5vjJLgMtQgCT+nOe4enkC47UyCjEVVgULv69BK+PwrlWrVo1evRo49eiKK5Zs+Zf//pXQkLCjBkzdu3a5evr6+XllZeXN3fu3PT09M2bN69fvz44OJiIXnzxxeuuu874rx5++OGPP/7YWOS3v/3tnj17Ro8eHRkZuW3btjNvma5evfrOO+880w1rUdAiBDCd4y2srINujjLDB9Cnk8XpP+gfTxCd8dDDWdLT08/5TlpaWlpa2tnf8fHxuf32241fn2nSxcbGnjng4Ycfvvnmm/v6+pydnYkoPDw8Nja2o6MjMjLSeIAkSfv27fvwww9l+imuEp4JANP503HpiUSR+4JqQzHGU5joL35SLj0Yi34g4MzJyenMcmtGUVFRZ48miqL4zTffmPy6hgpBCGAiJzvZ7gbp46lmm9D3dLK4bJfhvhhRYea3FMH2LV++3NyXMAz4bAhgIu8VS/fGiG7mm9g+2V/wd6IfTmOkEOBnEIQAptCto9Xl0n0xZn7ifhUn/g0rrgH8HIIQwBQ+PSlNDxTD3czcKXljlFjQyk60Y0sK+Mnjjz9eXV0tU/FXXnlFq9Ve4oA//elPlz7ANBCEAKbwfon0qzjzP24OIt0TI76PeRTwX52dnXq9Xo7Ku3fv3r59u7+//yWO8fX1XbFihRxnHxbzP5kANm93A9NJlBVkEe+oPBgrfnpS6tRd/kiwBy+99FJ4eDgRNTU1tbe3Nzc379mzx7hxUl9fX3Z2dl1d3ZmD9Xr98ePHc3Jy2trazi7S3Ny8d+/e1tbWtrY2478lopUrV957771EZDAYKisriairq6uxsZGIGhsbjZv0Ll269KuvvjL7AqQIQgDZ/b1E+lWcqSfRX0yQi5AVJH5+Eo1CICLKyMgw7jL47LPPLl26dP78+a+++uro0aM3btyYlZX1+uuvJyYmbt682XhwQkLC888/v3Llyri4uO+//974zbVr18bFxb3zzjtz5869+eab//jHPxJRd3f3li1b5s2bR0TNzc3GOfvffPPNgw8+SEQPPfSQcftDV1fXsWPHbtiwwRw/+v9g+gSAvFoGaEut9I8MC9oG6b4Y8fkjhgcwodACDJQdG6wpN9np3LOuJ/Gim560trbu379fqVTef//9d999d1FRkbe39z/+8Y+333577ty5RHT8+HHjNvT79u375S9/uWjRosHBwYceeujLL7+cPn16f39/UlKScUumY8eOBQYGjhgx4rKXlJSUdOjQoWXLlvH7KYcNQQggr9Xl0rVh4ggHc1/HWWYFCff1UV4rS/a2kGaq/ZIGB6TebnNfxU/mz5+vVCqJKD4+XqvVent7G7/+61//ajzg4MGDq1atqq+vHxwcPHnypE6nq6ioGBwcnD59OhE5OTnNnz/feGRTU5Pxn1+Wj4+PfG/rDBGCEEBen5RJf5ks/8aDwyEKdMco4ZMy6S+TLOvC7JBzwiTnhEnmvoqfGFt7RKRQKIyLpRm/NhgMRJSXl3fbbbf985//jI2N7erqSklJ0el0kiSJ4v+6FhQKBWOMiNzc3M6sLGP8piT9rzfeYDAoFD/dez09Pe7u7vL/cJeCvhEAGR1pZp06mqa2uIbXL0eLn52UBrA1EwzZkSNHJk2atGDBgqioqDMbLY0cOdJgMBw9epSIDAbDtm3bjN9PTEysrq4eGBggIh8fH1dX15KSEuNfGQyGoqKiM8uQlpaWJiYmmvqH+Tm0CAFk9HGZ9MvRoqW8J3OWCHch0VtYf1q6wbQbQoH1mjx58uOPP/7qq6+KovjDDz8Yv+nk5PTmm28uXLhw8eLFRUVFbm5uxqZeYGBgXFzcgQMHpk2bJgjCH/7wh8WLF8fGxlZUVMyfPz8sLGzatGlEJElSdnb2Sy+9ZM4fDC1CAPn0G+g/FdIdoywvBomI6O4x4v+V4d1Re/fWW28ZG2cPPPDAokWLjN+cO3fuo48+avx69OjRb7/9NhHFxcXt2rVLr9e7u7uvW7fuyy+/dHR0JKK77rpr9+7dM2fO/Pjjj8eMGXNmD/qHHnpo1apVxq+feOKJr7/+2tPTUxTFxx9/fOvWrca83LZtW1xc3Jl9oMwFLUIAuaw/LaX6CmHmXk3mYpaEiw/vM2j6SO1s7ksB81m8eLHxi3Hjxp35ZnR09Jlt5X18fIyzIIgoNTX1zBbzN954o/GLw4cPi6KYmZm5Z8+etWvXvvDCC8bv33nnnR9//HFDQ0NgYCARJScnz5w5s6en50w1Inr//fdfffVVGX+8oUEQAsjl85Ps9pGW2+nirKSFYeJ/KqRHEyz3IsHy9fb2vvbaa01NTWFhYevXrz8z+KdUKnNycs4+0sfHZ+TIkWd/Z+3ataa70ItDEALIom2AdjZIn0yzoOmD57s9Wvz9EQOCEK7GtGnTjAN+l7Vw4cKFCxfKfT1XAA8AgCy+rpLmBIueljR98Hwzg4TT3aysA2twg11DEALI4rMK6bZoCx0dPEMh0E1R4hcVCEKwawhCAP7qe1lhK5sfYgXP1+3R4qdYdxTsmxU8qABW5/MKtiRCdLSGZVsm+AmCQIeb0CgE+4UgBODvPxXSrRb8vug5bo4SvqpCo9CuNTQ0GDdIsk9W86wCWIvqLnaqm021vGXVLuaGSPGrKoYmoX0aGBi44YYbMjMzMzMzp0+ffs5Gg3YC0ycAOPuqii2JEJXW8yEz2VtwFCm3mY31tZrwthkn2ypruxpMdrqpoZNF4We35ptvvqnX68vKyojIuJvgypUrTXY9FgJBCMDZ11XSH8dZw/DgWX4RIXxdJY31tbLLtgHlrZUH6o+a7HQZIZPOWfn2hx9+eP755407SNx3332PP/64yS7GciAIAXiq6WGVXSwr0MqaVjdEijduN7w63tzXYX/mj5w1f+QsM15AU1OTh4eH8WsPD4+Ojg4zXoy5WE/3DYA1+KqSLQ63pn5RozRfQRDoeAsGCu0OY+z48ePGr48dOzZmzBjzXo9ZoEUIwNM31dLvU62yg/H6COGrKinFxyovHq7GBx984O7uLoriihUrPv74Y3NfjhlY2wdXAAvW0Esn2tnMICvrFzW6PlL8rhotQnv01ltvlZeX79+//7PPPjt7awj7gRYhADffn5auCRVV1vnxcryf0Kmjsg422tMqgxyumI+Pz5/+9CdzX4U5WecjC2CR1p2Srgu31hQRiBaGCetOoVFoXzw9PZVKe28RIQgB+OjW0b5GNtca1he9mOvCxXWnsMSMfcnNzU1JSTH3VZiZFT+0ABZlY42UHiC4W/T+g5cxI0goaWeaPnNfB4BpIQgB+Fh3il0Xbt0PlEqkOSHihtNoFIJ9se7nFsBC6CTaXCstDLP6B+o6DBOC/bH65xbAEuzRsFGeQqCLua/jqs0PFfdopB69ua8DwIQQhAAcrD9tC81BIvJ0oHG+wvY69I6CHbGFRxfA7DbWsAVh1jpx4hwLwsQNNegdBTti79NHAK5eaQfr0VGSt40E4aIw4Y18wweksJGfx5L09fUdPWq6vSZsSXNzs3zFhxGElZWVW7Zs8fLyWrRokbOz8wWPyc3NPXDggLu7e1ZWVkhICKeLBLBoG06za8MEm4mNkR6Cm0rIa2EpPjbzM1mEiIiIzMzM+++/39wXYpWUSuXo0aPlKj7E47KzsxcuXLh06dLS0tI33ngjJyfH0dHxnGOefPLJNWvWXHvttb29veXl5S+99BLvqwWwRBtqpEfjbWqUYUGosKEGQciZm5ubfS5pbfmGGoQrVqx47rnnnnzySYPBMG7cuK+//vr2228/+4DNmzd/+umnhYWFvr6+MlwngIXq1NHhJjYjyMaCUHzhqOG5FJv6oQAuZkg3+sDAwPbt2xcvXkxECoVi4cKFP/744znHfPXVV8uXL29pafn++++rq6u5XyiAZdpaK6UHCG7WvKDM+aYGCsXtTIslZsA+DKlFqNFoGGNBQUHGPwYFBWVnZ59zTGVlZWFh4Z49e2JiYu6666633npr+fLlFyt48uTJV1555ezvhIeH33TTTRc8WKfT6XS6oVwnyMf4vyDYzkAYNxtO05wgMs0tavwvkCTZ5zYIRFlq2nhad3uU3KeyPjqdThTRVjYznU5nMBgUistvn6lUKi/7i2tIQcgYI6IztURRPP9RHBwcHBwcPHLkiCiKGzduvO22226//faLXeXg4GBbW9vZ3/H29r7Y4y1JkgmefLg04/8CgvAcjGhznfhUvInuUFP+L8wLFjbXCrdG4NE7l/E/G7+UzEv6Ly7VhhSEarVaEITGxsaIiAgiamhoONM6PCMoKMjHx8f4QSk9Pb2jo0Oj0QQHB1+wYFxc3BtvvDHESxwcHDz/xRwwMZ1O5+joiCA8R34rc1UZYn1NdH9KkuTg4DCUT8FXb0EE+/0xvcrBUcT/+XlEUVSpbKs33NqIomgwGHhFw5Aa+E5OTunp6Rs3biQixtimTZtmzpxJRHq9/vTp08ZMnjNnTklJifH44uJiZ2dnf39/LpcIYLE21bJ5IbYZFCGugr+zkNuCmfVg+4b61ujzzz9/6623arXakpKS9vb2W265hYiqqqpGjx7d1NTk6+t72223vfXWW8uWLUtMTPzggw9efPFFfGICm7e5Vno80RTtM7OYFyJsqmHjfG0z6QHOGOqQ79y5c3fu3KlSqbKysvbv3+/q6kpEarX6s88+c3d3JyIXF5eDBw9mZmZKkrR69erf/va3Ml41gAXo0dORJjZNbbM5MTdE3IxFR8EODGNlmeTk5OTk5LO/4+7uftttt539x/vuu4/bpQFYtp31bLyfrU2cONtUtVDYytoHaYSDuS8FQE54CRjgCm2uleaG2PIT5KigyQHCjno0CsHG2fJjDCArG35T5oy5weLmWrwvAzYOQQhwJSo6Wa+eEmxlx4mLmRcqbEIQgq1DEAJciS11bHaw7U+rHOMpiAKVdiALwZYhCAGuxLY6NjvY5nOQiGhGoLC1DkEItgxBCDBsBka7GiQb23HiYmYFC9sQhGDT7OJJBuDrSBMLdhUCXcx9HSYxO1jc1SDp8Ooo2C4EIcCwbbWbflEi8nWiSHfhcBMahWCzEIQAw7atXpodbEfPzuxgDBOCLbOjhxmAix495TazTNtdWe18s4LFbZhWD7YLQQgwPLsb2FhfwXUYqxNavalqIa+FdWJ7bLBRCEKA4dlWJ82yp35RInJS0ER/YXcDGoVgm+zpYy0ADzvq2T8yTBqEg4bB7dV7ilvKSptPujg4j/GOTglInBw8zpTXMCNI3FnPFoaZ8pwAJoIgBBiGlgGq7mZjTbhF396a/X/L/b9Iz/BxgSlTAycZROlke/VHeavXFH/763H3RntFmuYyZgQK92WjRQi2CUEIMAzb66SpalFpkgahxKTXDvy1tLXitxN/naZOIqK+vj4HB4dJweNui/vFDye3PLnjD3cn374weq4JLmacn1DTwxr7KMDZBGcDMCkEIcAw7GxgM4JM0RyUmPTq/rfbBzr/MW+lo+Lc/QBFQVw0at6EoLTfbHu+Xz9wY8wiua9HIVBGgLi7Qbopyr7GR8Ee4J4GGIbt9aYIQomxP+x9rUfX+8rU585PwTPUrv5vz/rTd2Ubvz6xXu5LIqIZQcKOeswmBBuEIAQYqroe1jbAErxkD8I1Jd+29Xe8lPmsSqG69JEBrn5/mfXHz4q/LmgqlvuqZgQJOxoQhGCDEIQAQ7Wtns0IEkWZc7C4uezLknW/T39cKSqGcry/i+8zkx55KfvNjoFOWS8s0VvoHGSnupGFYGsQhABDtUP+ftHOwa4Xs19/etIjAa7+Q/9XE4PGzgjPfHX/2/JdGBEJRNMDxV1oFILNQRACDNXOejYjUN4g/Oj4p+kh469gjuC9Kcube1t3nNorx1WdMSNI2I5FR8HmIAgBhuRkJ2NEozxlDMKTbZV7avb9MvG2K/i3SlHx2PgH3s/9uF/fz/3CzpgRJKBFCLYHQQgwJDvrWZbMzcF3j3x0T/IyD0f3K/vnCX4xyf4Jnxd/w/eqzhbtIUhEFZ3IQrApCEKAIdmtYdPkDMLt1Xv69P3XjJx9NUUeSL3ju7IfG7obeV3V+aYHolEItgZBCDAkuxpYlmxvyhiY4aO8Tx8ee7coXNUpfF18rh9z7ScFX/C6sPMhCMH2IAgBLq+8g4lEUe5yBeG2qt1qN/8k//irL3VDzKL9dUfquhquvtQFTQ/EtHqwNQhCgMuTtTkoMfZF8bfLE27mUs1V5XLdqPlrStZyqXa+aA9BKdJJDBOCDUEQAlzergYZBwh3nc52VjmnBiTyKnhj7KJdp3Iae7S8Cp5jmhq9o2BTEIQAl7dbI9cro4zYp4Vf3Zl4C8eaHg7uC6Jnryn5jmPNs03DMCHYFgQhwGWUdTCRKFKeAcIjDccFQZwYNJZv2ZtirttWtbt7sIdvWaOsIGFHPfYmBNuBIAS4jN1yDhCuLdvwizELuJf1dvaaGDT2x8rt3CsTUZS7oBKFsg40CsFGIAgBLmO3hk2Xp1+0saepsOnEzPBMOYovGXPN2rINEpMlrqaphT0aBCHYCAQhwGXs1bCpalmCcF35j3OjZjgpneQoHu8b46ZyPao5LkfxqYHCHgwTgq1AEAJcSmUXMzAa6cE/CHUG3Y+V2xdFz+Ve+YzrRs9fW7ZBjsp4cRRsCYIQ4FJ2N7Bp8jQHd57OiR4RGeoRLEdxo5nhUwubTsgxj2KUp2BgVN2FLARbgCAEuJQ9si0x+mPltoWjZGwOEpGT0nFmxNRNlTvlKD41EMOEYCMQhACXslueqfSNPU0VbdWTg4a97+BwzYucsblqByP+iTVVLexG7yjYBAQhwEXV9rAePRstwx6EW6p2zgjPVClU3CufY4xPtIPCoaiplHvlaWgRgq1AEAJc1K4GNk0tytExuqVq19yoLBkKX8CcyOmbZJhQGDtC6NSxuh5kIVg9BCHARcnUL1rUfIIRi/UZzb3yBc2JzNp9et+AYZBvWYEoI0BEoxBsAIIQ4KJkmkG4uXLnvKiZ3MtejK+z9xif6H21h7hXnhaIYUKwBQhCgAtr7KOmfhbvxTkI9ZJh1+mcOZHT+Za9tLlRWVuqdnEvm6kWshsRhGD1EIQAF7ZHI2WoRe4jhEc1x8M8QvxdfDnXvaSMkEl52sKuwW6+ZZO9hfpepu3jWxXA1BCEABe2V8MyZegX3Xk6Jys8nXvZS3NWOqWpk/bXHeZbVhRosr+wT4udKMC6IQgBLkyOAUK9ZNhXe2ha6BS+ZYdielj6zlM53MtmqsW9eF8GrByCEOACOgapspOl+nAOQmO/qK+LD9+yQ5EeMjG/qahH18u3bKYaq2+D1UMQAlxAdiOb6C+oeD8fZukXNXJWOqUGJObUHuRbdryfUNbBunR8qwKYFIIQ4AL2aqRMNeenw4z9okZy9I46iJTmK+zDu6NgzRCEABewp4H/mzK5jXlhHsFm6Rc1mhI8IU9byL13dKpa2KvB+zJgxRCEAOfq01NBG5vgxzkIs2sOZoZO5ltzWFxUzkn+cYfqc/mWxfsyYO0QhADnOtDEkrwFFyXPmozYvrrDU4LH8yw6fOkhE/fVcV5iZnKAkNvC+g18qwKYDoIQ4Fxy9IuWtVS4qpxl3YZ3KKYET9hfd0Qv8UwtVyXFjhAON6FRCNYKQQhwruxGKSOA86ORU3cwPWQi35pXwMfZK8Q9qKCpmG/ZjAABvaNgvRCEAD+jl+iQlk0J4D9AmB4ygW/NK5MeMoH7JIpMtZDdiPdlwFohCAF+5lgLC3cXvB151mzs0bb2t8X6jOFZ9Eqlh0zcW3OAb81MtbivkRnQJgTrhCAE+JnsRpbBvTlYe3BK8ARRkGOL32GLGhEuCmJV+ymONX2dKNBFKGhFEoJVQhAC/Iwca23vqzs8xTL6RY2mhIzP4b0AN4YJwXohCAH+hxHlNEoZXIOwX99f0lyWFpDEseZVmhg09mD9Ub41sTchWC8EIcD/lHUwZ4UQ6sozCHMb82N8RrmonDnWvEop/gkn2yo7B7s41sxUC3sa8L4MWCUEIcD/yNEveqDu6MSgsXxrXiUHhUOSX3yuJp9jzUh3QSkKFZ1oFIL1QRAC/E9OI+PbL0pEhxuOWVoQkjy9oxlqDBOCVUIQAvzPXg3nV0ZPd9bqJX2EZyjHmlxMDh53sP4oI565lR4g5GCYEKwQghDgJ5o+ahtgsSN4BuGB+qOTgsdxLMhLoFuAs9K5oq2aY028LwNWCkEI8JO9GilDLYpce0YP1edOCEzjWZEf7r2jiV6Cppdp+ziWBDAFBCHAT7I1LJ1rv2i/fqC4uXSsOpljTY4mBY89wDUIRYEmBwg5WGsNrA2CEOAn3NeUydMWjvKKsqiJE2dL9k842VbJd5/e9AARw4RgdRCEAERE3Toq62BjfXkG4ZGG4+MCUzgW5MtR4RDrMzpPW8ixJoYJwRohCAGIiPZp2VhfwVHBs+ZhjUUHIRGNC0w53HCcY8GJfkJRG+vVcywJIDsEIQARUbZG4tsv2trX1tzbMsY7mmNN7sapU45wDUJHBSV5C0Q7w3YAACAASURBVAexSS9YFQQhANFPb8rwfByOaI6nqZNEwaIfsWivqM6BLm1vM8eaWH0brI5FP6UApqGT6Egz5814j2jyxqktul+UiERBSFMnHdXkcayZoRayNXhxFKwJghCAjjazkR6CpwPPmrmafAsfIDQaq07m2zuaESAe1DI9ohCsB4IQgHJ4T5yo6jitEMQgNzXHmjIZH5hyVHOc41prXo4U6ibkY5NesB4IQgDKaeQ8lf5oQ974wFSOBeUT4Orv5uBW0VbFsWZGACZRgDVBEIK9M27Gm85104mjmjyLXVDmfGPVyXy3ZEpXY/VtsCYIQrB3ZR3MRclzM16JSQVNxSkBCbwKyi01IPFYYwHHgukB2KQXrAmCEOxdNu+tl0pbT/q7+Ho5jeBYU1ZpAUl52iIDM/AqGOUuKEWhsguNQrAOCEKwd9wHCHM1+WnqJI4F5ebh6K529S9tqeBYc0qAkI3ZhGAlEIRg77J570p/rLEgNcCagpCI0tRJxxq5DhNik16wHghCsGuNfdTcz+L4bcark/TFzaVJ/nG8CppGakAS32HCDLQIwXoMIwjXrFlzzz33PPfccxqN5hKHffnll+++++5VXxiAKWRrpPQAgeNmvCXNpaEewe4ObtwqmkRKQEJR8wmdQcerYLKPUN/LWgZ41QOQ0VCD8K9//etzzz03bdq0lpaWjIyMgYEL3+AHDx58+OGHX375ZX5XCCCjnEY2hesSo7mN+akBiRwLmoaryiXMI6SkpYxXQYVAE/yEfdikF6zBkH4FGAyGlStX/u1vf1u2bNkHH3zg7u7+zTffnH/YwMDAQw89tGLFCt4XCSAX7pvxHtMUpFnbAKFRakBiLtdhwinYpBesxJCCsLa2tqamJisry/jH6dOn5+TknH/Yyy+/vHjx4rg4KxsdAbvVo6eSdjbej1sQDhoGS1tPJlrbAKFRmjrpWCPvTXoxTAjWQDmUgxoaGtzc3BwdHY1/9PPzO3LkyDnH5OXlff/994cPHz5w4MBlCx45cuT6668/+ztxcXHPPvvsBQ/u6+tTKLjulwrD19vbKwiCIPBsPJnd7kYxyUuUBnp7ORXMby6O8Ahjg1LvIK+SP9PX16fX62V6HEa6RpS2lLd3tTso+Kw+nuRGea2q1q5eJ9t6fAcGBkRRVKlU5r4Qu6bT6QwGgyRdvu/dyclJFC/T5BtSEDo5OQ0ODp754+DgoLOz89kH6PX6e+6554MPPjgTlpcWHBx88803n/0dPz8/JyenCx6s0+ku9ldgMnq93snJycaC8FAry1QzjndXcVtZmjpRvtuVMebg4CBTEDqRU+SIsKqemmT/eE4FKW6EVNDlmMl1dorZCYKAIDQ7hUJhMBiG8qxdNgVpiEEYHBw8MDDQ1NTk5+dHRDU1NcHBwWcfUFlZmZeXt2zZMiLq7+9vbW0dOXLktm3bIiMjL1gwMDDwpptuGsqpiUgUxaH8JCAr4/+CjQVhjlb/aLxC5PfOaF5T0e3x18t3u4r/JVP9ZP+E/KaiVDW3l30y1GxfkzAtyKaeX7n/F2AoRFFkjPH6XxhSFT8/v8zMzM8++4yI2tvbN27cuGTJEiJqbW1du3YtEUVGRp44cWLr1q1bt25duXKlp6fn1q1bQ0JCuFwigBz0Eh3S8tyMV2fQlbaUJ/jG8ipoeikBCce1PIcJ0wOwSS9YgSG1CInoz3/+8+LFi3fu3FlcXDxnzpzJkycT0YkTJ37xi18wxlQqVVRUlPHImpoahUJx5o8Alul4KwtzE7yH1Jc/JMUtZZEjwlxUzpc/1FIl+cWvyH5DZ9CpFHz6/TLV4l17DAZGCpvqSgBbM9QgnDJlyokTJw4cOBAYGJia+tNGa2lpaaWlpeccOXHixMOHD/O8RgAZZGs4r6x2vLEg2d9qdpy4IBeVc6hH8InW8kQ/Pi+++jlRgLNQ1MaSvJGEYLmG0cHq7e19zTXXnElBInJycho9evQ5hzk5OYWHh/O5OgDZcF9r+7i20Iq2XrqYVP/E41wnUWCtNbB8GO8FO5XTKHGcSq+T9CdauDWkzCiZ9zBhhhq71YOlQxCCPTrZyURBiHDnFoQnWspCPYJdVS68CppLsn98SXOZXuK2N2FGgLAXLUKwbAhCsEd7NWwq1wHCPG1xipUPEBq5qlyC3QNLW8t5FRzlKRgYO9WNLATLhSAEe5TDe4nRPG2htb8pc0ayf3y+tphjwfQAEcOEYMkQhGCP9nJ9ZVRiUnFzaaKfFc8gPFuSf3y+tohjwYwADBOCRUMQgt1p6qfGPhbvxS0Iy9sq/V18PRzdeRU0ryT/uPymYolxmwifocYwIVg0BCHYnb0aKT1A4DjFO09blMRpfU5LMMLR08fZu7L9FK+CKT5CTTc26QXLhSAEu5OtYRlqnnd+vrbYloKQiJL94/P49Y4qBJroj016wXIhCMHuZDeyTH5vyjBihU3FvHZssBBJ/nGchwnVeF8GLBeCEOyLcTPecfw24z3VUeusdPZ19uZV0BKk+CfkaYsYcYuuTAwTggVDEIJ9OaBlKT4Cx61i87VFyda/sto5/Fx8HRUOtZ31vApO9BMK2lifnlc9AJ4QhGBf9mokvlPp87XFyda/str5kgMSOA4Tuigp0Us41IRGIVgiBCHYl70alsn1TZk8bWGivw0GYZJfLMcgJPSOggVDEIId0Ul0uIlN9ufWItT0aHWSPsQ9iFdBy5HoF1fQxHN9mUy1uBeb9IJFQhCCHcltZtEegqcDt4L5Wlt7X/SMMM+Qfn2/treZV8H0AOGAlukRhWB5EIRgR/ZoWCbnAcKiJFvsFyUigYQEv9iCphJeBb0cKdxdON6K3lGwOAhCsCN7eQdhQVOxDexBeDGJfnEFXFffxjAhWCYEIdgLRrSvUUoP4HbPdwx0anubR3pF8CpoaZL8eQ8TYm9CsEgIQrAXRW3Mx0kI5Ld1bkFTSYJvrELgNyfRwoz2HtnQ3dg92MOr4LRAMVsjIQnB0iAIwV7w7xfVFif628jWSxekEBRjfKILm7kNEwa6kLtKKG1HFIJlQRCCvdir4bwZb36TTW06cUFJfvF8hwmnBgp70DsKFgZBCPYiW8OmBXILwn79QGX76RjvUbwKWqYk/7g8rkGYgWFCsDwIQrALFZ2MEUW6cwvCkpaykSMinJSOvApapjjfMSfbKnUGHa+CU9XCrgYEIVgWBCHYhT1cm4Nk0zMIz+asdArzDDnRWs6r4ChPQWJU3YUsBAuCIAS7wH0qfUFTiT0EIREl+cXlY5gQbBqCEOzCngbGcdMJiUklLWUJvrb8yugZiX5xhfzWlyHMJgTLgyAE21fXw7p0LGYEtyA82Vbl5+zj4ejOq6AlS/KPK2gqkRi36EKLECwNghBs324NmxoocuwYLWgqtsmtly7Iy2mEp6PHqY7TvArGewmtA6y+F1kIlgJBCLZvr4ZnvygR5WuLE/3sol/UKNE/Lp/fWmsCUUaAmI1GIVgMBCHYvt1cBwiJqLCpxIbX2j5fol9sgZbrMKEavaNgQRCEYOOa+knTx5K8uQVhfbeGBCHQLYBXQcuX5MezRUhE0wKF3ZhNCBYDQQg2bneDlBEgchwhzNcWJ9tTc5CIQj2CdZKO4ya9KT5CXS9r7udVD+CqIAjBxu1u4DyVvqCpOMHOgpCIEnxjOC46qhBosr+wR4Pt6sEiIAjBxu3iHYT52uIkm9504oIS/WJ5946K6B0FC4EgBFvWMkC1PSzVh1sQdgx0tvS1Ro2I4FXQWiT6c96tfjqGCcFiIAjBlu1ukNIDBAXXAcJ4vxhRsLsHZ7T3SE2PtnOwi1fBNB+hupu1DPCqB3Dl7O55Bruyu4FNC+R5kxc2lSTZ3wAh/XeT3uLmUl4FlSJN9hf2YpgQLACCEGwZ/wFCe1pT5hxJfvEFXBcdxTAhWAgEIdis1gGq7mJp/AYIBwyDle2nbH4z3otJ9Ivluw0FhgnBQiAIwWbt0UjpakHJ7x4vbi61h814LybeL6astWLQMMir4DhfobKLtWKYEMwNQQg2a1cDm6bmeYfna4vtZA/CC3JWOoV7hJS2VvAqqBRpEoYJwQIgCMFm7axnWUGclxhNsKe1ts+X6B9bwHU24fRAcRd6R8HcEIRgm1oG6FQ3zwFCiUlFzSfsatOJ8yX6cZ5NmBUo7KhHEIKZIQjBNu2slzICeA4QVrRX+7n4eDp6cKtohZL8OG/SO9ZXON2NRUfBzBCEYJt2NrCsIL4DhEV2tfXSBXk7e3k4unPcpFcpUnqAsLsBw4RgTghCsE0761kW/yVG7T0IiSjJP57voqNZQeJODBOCWSEIwQZp+0jTx5L5DRCScU0Z/3iOBa1Ukl9cnraIY8GsQGEnhgnBrBCEYIN21EvT1CLHJUbruhpIENSu/twqWq0kf85BmOIjNPYxTR/HkgDDgyAEG7SzgfPEifym4mQ0B4mIKMQ9iDGpsUfLq6AoUIZa3FWPYUIwGwQh2KBdDWw65wHCIvtca/uCEv3i8rhOopgRJGCYEMwIQQi2pq6HtQ6wBC++u9JjgPB/Ev3i+E6rnxEkbMcwIZgPghBszdY6NjNIFPnlYPtAR3t/R4RnGLeKVi7JPy6f6zBhvJfQo2PVXchCMA8EIdia7fVsJt8BQm1xgl+MKPCsadWivSKbels6Bjp5FRSIsoJENArBXBCEYGt2NrCZwbwHCNEvehZREOP9YvjuTTgTvaNgPghCsCnF7UwpUJQ7zyDM0xbhldFzJPvF8+0dnRUsbK+XkIRgFghCsCnb69hsrs3BHl1vbVf9aO9ojjVtQHJAPN9NesPdBHeVUNSGKAQzQBCCTeE+QFjYVBLjM0olKjnWtAEx3qOqO0736nhOg58ZJGyvQxCCGSAIwXYYGO3VSNN5r7Wd5Id+0XOpFKox3tHFzaUca2KYEMwFQQi240gTC3YV1M48a+ZpsabMhSX5x+c38RwmnBEk7tFIOqwwAyaHIATbsa2ezeLaLzpoGKxor4rzHcOxps1I8o/La+QZhL5OFOUuHNSiUQimhiAE27GlVpoTwvOWLm4ujfQMd1I6cqxpMxL8YktbTw4aBjnWnB0sbMOio2ByCEKwEV06ym1hmWq+EyfQL3pRzkqncM/Q0tYKjjVnB4tbatEiBFNDEIKN2NUgTfQTXLm+3ZnfVJSIzXgvLtk/nu+WTBlqobCNtQ1wLAlweQhCsBFb69jsYJ73s14ylDSXYdOJS0jyj89rLORY0ElBUwKEXQ3oHQWTQhCCjdhSy3kqfWlrebB7oJuDK8eaNibZP76o+YReMnCsOTtY3IrZhGBaCEKwBbU9rGWApfjwDMLjjYUp/gkcC9oedwe3QLeA8ja+w4TCFgQhmBaCEGzBplo2O5jn1ktEdFxbmByAILyMZP+E41x7RxO9hV49q8SWTGBCCEKwBdt4LzFqYIaiphMYILyslADOQSgQzQrCu6NgUghCsHoGRtvqpLkhPIOwrLUi0C3Aw9GdY02blBKQUNhcIjGer7fMDUHvKJgUghCs3qEmFuIqBLlwHiBMxgDhEHg4uPu5+Ja3VXKsOTdE3FkvDeLVUTAVBCFYvU010vxQztvH52kLUzBAODQpvIcJfZ1otKewrxGNQjARBCFYvR9r2TyuK6tJTCpsOpGIAcKh4T5MSETzQoRNtWgSgokgCMG6NfdTeQebEsCzRVjeVunr7O3l5Mmxpg1L9k8oaCrmO0w4L1TcVIMWIZgIghCs2+ZaKStIVHG9kXM1+anqRJ4VbZqXk6ePsxffYcIJfkJdL6vrQRaCKSAIwbptqmXzub4vSkTHGgtSA5L41rRtaeqkY40FHAsqBJoVLG7Gu6NgEghCsGISo628J04YmKGwqQSbTgxLakDSMQ3PICTjMCF6R8EkEIRgxY40M18nIcyN6xKjLScD3dSejh4ca9q8lICEgqZivouOzgsRt9Vjw3owBQQhWLEfTkvXhnHuF83V5KcFYIBweDwc3APd1GWtJznWDHCmaA8hG5MoQH4IQrBiP5xmC0I538O5jXhT5kqkBSTmavL51lwQKm44jSYhyA5BCNaqvped7maT/Xm2CHWS/kRLeZIfBgiHLVWdmNvIOQivDRN+OI0WIchuGEGYk5MzZ86clJSUp59+emDg3D2ktVrts88+O3Xq1PHjxz/88MONjY1crxPgXD+cZnNDRCXXz3LFzaWhHsHYg/AKJPsnlLSU6Qw6jjXTfIVuPZV3IAtBXkP9LdLc3LxgwYJbbrnl888/379//4svvnjOAWVlZf39/StWrPjwww/r6uqWLFnC+UoBfm5DDVvAe4DweGNBKgYIr4iryiXcI7S4pYxjTYHomlBhA94dBZkNNQhXr149fvz4u+66Ky4u7rXXXvvwww91up999MvIyPjLX/6SlZWVkpKycuXK/fv3d3V1yXDBAEREAwba3SDN5bqyGhEd1eSlqTGD8AqlqZOOavL41lwQKmyowTAhyGuov0cKCgrGjx9v/HrcuHGtra11dXUXOzg3NzcoKMjdHVvYgFx21LNkb8HHkWfNfn1/eVsl9iC8YmPVybm8g3BWsHhQyzoG+VYF+BnlEI/TarUxMTHGr1UqlZubW2NjY0RExPlH1tfXP/LII++8884lqm3fvj0yMvLs70yaNOnDDz+84ME9PT2CwLkHDIarp6eHMWY5/xHfVSpnBVB3dz/Hmocbj0V7Rur79d3UzbEsR319fQ4ODgqFwtwXcmGRzqEn26q07U0uSmeOZSf5qtZX9C4OtZR24cDAgCiKKpXK3Bdi13Q6ncFg0Ov1lz3SxcVFFC/T5BtqEHp6evb09Bi/liSpt7d3xIgR5x+m1WpnzZr1yCOP3HjjjZeoNnny5JUrV579HWdnZzc3twsezBi72F+BKbm6ulpIEDKiH+v1W69RuLk5cSxbXFo2ITjNkm82hUJhyUFIRHG+Y072VE0JnsCx5i+ipM2NyqWxlvJTq1QqBKHZGYPQyYnPb4ChBmFkZGRZ2U/D4BUVFQqFIiQk5JxjWlpaZs+efdNNNz3zzDOXrubi4hIVFTXcawUwOtrM3FQ0xpNzKh/RHH9q4q/51rQ3Y9XJRxry+AbhonDhd0cMOknBd2l1gDOGemfdfvvtGzdurKioIKJ33nln8eLFrq6uRLR69ept27YRUWtr6+zZs9PT0x999NG2tra2tjaDged6SwBnrDslXRfOOQXb+tubeltifKL5lrU349QpRzXH+dYMchFGewp7NHh3FOQy1CCMjY194YUXxo4dGxIScuDAgTfffNP4/fXr1+/fv5+I9uzZU11dvWbNmpH/VV1dLdNFg537rppdF865dXBEczwlIFEU0Oi4KqO8R7b1dzT1NvMte124uO6UpYwRgu0ZatcoET3xxBO/+tWvurq6/Pz8znzzyy+/NH6xePHixYsXc746gPNUdLLmfjbBj3OL8KgmfywmTlw1URBSAxJzNflzo2ZwLLs4XJi3SXpnMlnEGDXYnOF9/nVycjo7BQFM77tT7LpwUeT9GzFXkz9WncK5qF1KUycd5b3WWuwIwVlBx1vQOwqyQEcQWJl1pyTu/aKnO+uIWJhHMN+y9ml8YOqRhmOMOIfWonABvaMgEwQhWJPGPipsYzODObcHD9XnTghM41vTbgW5qZ2VzpXtp/iWXRIhflOFFiHIAkEI1uTbaumaUNGB9217qCF3QhCCkJsJQamH6nP51pzkL3QMUkk7shD4QxCCNfm6Sro+gnNzcNAwWNhUMladzLesPRsfmHa44RjfmgLR4ggBjUKQA4IQrEZzP+U2s3m8F9rO0xaN9Ip0VbnwLWvP0gISS1rK+vQ8F8AjohsixW+qMUwI/CEIwWp8Wy3NDRGdhzHlZ0gONRybEJjKuah9c1I6xfqMPtZYwLdsRoCg7aOTnWgUAmcIQrAa31RJN0Tyn0h2GG/KyGB8YOrhBs7DhKJAiyOEr9E7CrwhCME6tA3QwSY2P5TzHdvU29zW3zHKeyTfsjAhKO1QPedhQiK6PkL8pgq9o8AZghCsw9pT0qxg0ZV3v+jB+txxgSmiZeyqYUuiRoT36/vruzV8y04LFE73sKouNAqBJwQhWIcvKqRbo/jH1YH6I5OCx3IvCwIJE4PG7q87zLesQqAbIsU1lQhC4AlBCFZA20dHm9k1vPtFdQbdscaCiYEIQllMDh63v+4I97K3RomflqN3FHhCEIIVWFMpLQzj/77oscaCSM9wD0d3znWBiIjGB6YWN5f26vr4lk1XC70GKmxDoxC4QRCCFfiiQrp1JP97dX/94cnB47iXBSMnpVOc7xju2xMKRDdFCl9UoFEI3CAIwdKd6maVXWxWkAwDhHVHpwSP514WzpCrd3SkuKaCoUkIvCAIwdJ9XsFuiBSVvG/Vqo7TBmaIHBHOuS6cZUrwhH11hyTemZXiIzgp6JAWUQh8IAjB0n1RId0axf9GPVB3ZErwBO5l4WyBbgEejh7lrRXcK98yUvwMvaPACYIQLFpuM+vSUbpajn5RTJwwhSnB4/fxnkRBRMuihS8qpEFEIfCAIASL9u9y6c5R3Lejp46BzpPtVWkBSbwLw7mmBI/Prj3IvWyEuxA3QvixBkkIHCAIwXLpJfqyUloazb85mFN7aHxgqoPCgXtlOEeCX1xrfxv3JWaI6I7R4r/LMUwIHCAIwXJtqJFGewojPfgH4d6aA5khk7iXhfOJgiBTo/CmSHFXg9TEea8nsEcIQrBcq8rZHaP536L9+oE8beHEIAwQmkhm6KTsmgPcy7qp6JpQcQ1emYGrhiAEC9UyQDvqpRsi+d+iB+uPxvvGuDm4cq8MFzQ2ILmivbq1r4175TtHif/Gcmtw1RCEYKFWl0sLw0QPFf/Ke2sOZIaiX9R0VArV+MDU/fX8Z9bPCBKa++l4C0YK4aogCMFCfVQq3RvD//7US4aDDUfTQyZyrwyXkBk6eW/Nfu5lRYHuGiN+WIpGIVwVBCFYomwNMzDKkGH64PHGgjCPYB9nL+6V4RImBY0taCrp0fVyr3z3aGFNhdSr514Y7AiCECzRh6XSfTH8pw8S0Y5Te6eHZchQGC7FVeWS7B+fU3uIe+VgVyE9QPyyEo1CuHIIQrA4HYP0/SlpabQs/aLZtQenhk7mXhkua0Z45s5T2XJUvjdGQO8oXA0EIVic1Sel+aGinxP/ykc0x8I9QwNc/fiXhstJD5mYpy3sHOziXvmaUPF0N3YohCuHIATLwoj+Xiw9ECvLnbnjVHYW+kXNxFnpNFadLEfvqEKge8aIfy9GoxCuEIIQLMu2OqYQaKoMr8noDLp9tYemhqFf1GxmhGfuOLVXjsr3x4pfVEhtA3LUBtuHIATL8m6R9FiCLLflwYbcaK9IX2dvOYrDUEwOHl/cXNo+0MG9stqZrgkVP8HkergiCEKwINVdbL9WunWkTP2ie7PC0S9qTk5KxwmBaXtlWG6NiH4dL/6tWJIwUAjDhyAEC/JusXTXaNFFyb9yj673YP3R6WHp/EvDcMyOnLalapcclSf5C96O9GMtkhCGDUEIlqJHT6vKpQfleU1m9+l9aQFJno4echSHoZsYNLa2q16OXZmI6OE48a+FBjkqg21DEIKl+PCENDNIjHCXYxo9ba7cMScyS47KMCwKQTEjPGNL1U45it86UjzRQcew9CgME4IQLIJeorcLpd/I85pMY4+2quP0JOy7ZBnmRGZtqtzBiH9cqUR6OE5cWYBXZmB4EIRgEf5TKUV70ER/WZqDmyp3zozIVClk2MkChm+Md7Sz0qmw6YQcxR+IFTfXSqe60SiEYUAQgkVYWSA9laSQqfi26l1zI2fIVByuwOzI6TL1jrqr6O4x4tuFaBTCMCAIwfy21DEDozkhsjQH87VFoiDG+IySozhcmTkR03edyunX98tR/JF4cVW51ILJ9TBkCEIwv5dyDU8ny7LXBBF9X755YfQ8eWrDFfJ18Un0j90hzxrcQS7CDZHiXwrw+igMFYIQzGxrHdP2081RstyKnQNdB+qPzImaLkdxuBqLRs37vnyTTMWfSxHfL5GaZWlwgg1CEIKZvXzM8GKaqJCnPbipcvuUkAkeDu6yVIerMCFwbFt/e1lrhRzFw9yE6yPFd4rQKIQhQRCCOW2tY5o+uZqDRLShYusi9ItaJFEQFkTP/uHkFpnqP58iflAitWKkEIYAQQjmtCJXxubg8cYCRpTgFyNLdbhq146cs+PU3l5dnxzFw9yEJRHiSowUwhAgCMFs1p2SunR0i2zNwbVlG5eMvkam4nD1vJ290tRJm6t2yFT/D6niP0qk2h7MKYTLQBCCeRgYPXdEen2CQqa3RTU92mONBfOiZspSHTi5Mea6r0+sl5gsWRXsKtw1Rnz5GOYUwmUgCME8/q9U8nOiufLMHSSir098v2DkbGelk0z1gYtEv1gPR7f9dYdlqv+7FMV3p6TidjQK4VIQhGAGvXp66Zj0xkS5lpLp1fVtrtq5ZMwCmeoDRzfEXPfliXUyFR/hQE8lKZ47jEYhXAqCEMzgz3mGqWphnK9czcEfTm6eEJjm7+IrU33gaHrYlIZuzYmWcpnqPxwn5rWyHfVoFMJFIQjB1Cq72Psl0usT5Lr39JLhm9IfboxZJFN94EshKJaMXvD1ifUy1XdS0FuTxF/vM+jQLISLQBCCqT22X3oqSRHsKldzcEvVzmD3QCwuakUWjZp3uOFYTWedTPUXh4sR7vReMZIQLgxBCCa1sYaVdbDH5Nl3kIgkJn1e/M0dibfIVB/k4KpyWTLmms+Kv5HvFH+ZpHj1uEEjy5RFsHoIQjCdHj09vM/w7hSFg2z33daqXT7O3sn+8XKdAORxY8x1+2oP1XbVy1R/tKdwzxjx0f2YXw8XgCAE0/ndYcP0QGF2sFydohKTPi366o6Em2WqD/JxVbksHj3/i+Jv5TvFH9IUha3s22p0kMK5EIRgIge07JtqZzTQtAAAFX1JREFU9qZsUyaIaFv1bk9HzzR1knynAPncELNoT83++m6NTPUdFfTRVMUj+6U2LEAKP4cgBFPoN9BdewzvTBK9HeU6xaBh8KO8z+5PXS7XCUBmHg7uN8Vc98/jq+Q7xWR/YUm48PhBdJDCzyAIwRSeOmhI9hGuj5Txfvu6dH2MT3SiX5x8pwC53RS7uLi5tLDphHyneHW8IqeRfV2FDlL4HwQhyO7HGvb9afb3KTJ2inYNdn9Zsu6+FDQHrZujwuGXibf+Pff/GMk1/91NRV9kKX61z3CqG1Ps4ScIQpCXpo/u3qv/bLrCS7ZOUSL6OP+LrPD0EPcgGc8BJjE3aka/vj+75oB8pxjrK/wmQXHHboMBUQhEhCAEWeklumWH/oFYRYZarjdFiai8rXLHqb13Jt4q3ynAZERBfHjsPe8e/Ve/vl++s/w2SVQK9IejGCwEIgQhyOqJgwZHkZ5LkfE2kxh769D796cs93T0kO8sYEpp6qRk//iP87+Q7xSiQF/OVH5RgcFCIEIQgnw+r5B+rGX/mamUaQN6o+/KNipF5byR2HfQpvx67D1bqneVt1XKdwpvR/pmluJX+wxFbeghtXcIQpBFTiP7zQHDt7MUIxxkPIu2t/nfBWuenPCQQHKGLZich6P73Um3rzz4dwOTsfcyxUd4bbxi8VZDk4y9sGAFEITA38lOdtN2w6ppygQvGfNJYuzV/W/fGLMo3DNUvrOAuSyInu3m4Lq68EtZz3LnaPG2kcI1m/Q9elnPAxYNQQicNfTSnB8Nr4wX5dt93ug/JWt1Bt1t8dfLehYwF4GEZyY/+l3Zj0XNMk4rJKIXxyoSvIVbdxj0GC60VwhC4EnbRzM36u+PEe8YJe+tdbKtck3x2ufTnxAF3MM2y9fZ+4kJD/4x561enYzbRghE/8xQMGJLd2FChZ3CLxHgprmfZv2ov2Wk+HSyvPdV52DXC3tfe2TcvWpXf1lPBGaXGTp5XGDKnw+8I98UeyJSifTNLGWPnqFdaJ8QhMBHXQ/L2qC/Llx4IVXem0pi0h9z3poaOnlmxFRZTwQW4tFx97f3d8g9WOgg0lczle2DbOkuwyCy0M4gCIGDik42fYPh+kjh5bEyrqNm9P6xTyQmYTU1+6EUFS9m/HZ9+eac2kOynshJQevnKBnRvB/1nTpZTwWWBUEIV2t3A8tYr/99qvhimuwp+GXJdwfqjvwh4ykMDdoVb2evl6c++8bBdwuaSmQ9kaOCPs9SjBkhZG3Q1/RgwNBe4LcJXJW/FUu37NB/mqVcLvPbMUS0qXLH16Xr35zxoruDm9znAksT4zPqhfSnfr/nldLWk7KeSCHQ++mKW0eKk9YZsjXIQruAIIQr1KWj5bsM/zwh5SxSzgySfT77lqpdH+atfmvmywF4QcZepamTnpjw0LO7Xq5or5b7XE8mih9PU9ywXb+yQEIY2jwEIVyJg1qWulbvrKT9i5RR7rKn4Ncn1n+U9+lbM17C/hJ2LjN08qPj7n9i+wv52iK5zzUnWDh4nXJttTTvR31Dr9xnA3NCEMLw9OjpyYOGxVv1b0wQ/5GhcFHKezqJsX8c+/f3Jze9N+fPWEEGiGha2JQX0p98Ye+f99Tsl/tc4W7CrgXKKQFiylrdv0rRNLRZCEIYhrXVUuI3em0fFVyvWhIh+83TOdD19K4VJS1l783+s7+Lr9ynA2uRpk56Y8aK947+6x/H/i0xeec6KEX6Q5q4db7yo1Ipa4M+rxVpaIMQhDAkR1vF6RsML+ZKH2UqVk1X+DrJfsaCpuL7Nj0e5Rm+cuZLHo7usp8PrMoor6gP579V1lbxxI4XGnua5D5dkreQs1B5S5Q470f9/fuEOvSU2hYEIVzGvkZ2zWb9shzVsmghd4lyhvzvxfTr+/965MMXs994ZNx9D6b9UiHIPisDrJGno8cbWS+OU6fc9+Pj35dvknXpGSISBXogViy9SRXgTOO+pwdzDNVdaB3aCAQhXNigRJ+dlCZ9r1++27AkXDy+YODuMaKsOwsSkcTYpsody9Y/1KPr+WTBu1OCx8t7PrByoiDeHn/DO7Nf+bFy+0ObnypoKpb7jB4qeimVFSwmH0ca953+hu2GXQ2IQ6sn86sOYIUON7HVJ6U1FVKyj/C7ZPHaMFEUqLtb3pNKTNp9et+nRV85K51WZD4d5ztG3vOBDYnwDP373Ne3Ve9+OeetaK+I2+NviPeNkfWMvk70x3GKZ5IVq8qlX+UYDIyWjRKXRgvhbtgX0yoNLwh1Op1Kpbr6Y8DSDEqU08h+OC19W81UIi2NFg9ep4yUf14EEbX1d2yt2rm2bKOPs/fdyUvRCoQrIJAwO2L6tNApGyu2/THnLV8Xn+tGzcsMneyokHFjaDcVPRQnPhQnHtCyVeXSuO8MEW7C9ZHi/FAhyVtAJFoRgbEhteu7u7vvuOOOrVu3KhSKp59++plnnjn/mNdff/3VV181GAwzZ85ctWqVu/uFX3D4/PPPN2zY8Nlnnw3xEru6ui5WCq7GgIFyW9heDdurkfZqWMwI4ZpQcUm4kOh9gUe4u7vb1dVV4Pd0dwx07qs7vLdmf562KCNk4qJR8+T+FG8D+vr6HBwcFAoMml6KxKQ9Nfs3nNx6oqU8M3RSZujkseokB36JODAwIIri+R/3DYx2N7DvTkmba1mXjs0MEjPVQoZaiPEURKQibzqdzmAwODnxeW1vqEH4u9/97ujRo+vXr6+vr584ceLatWunTJly9gEHDx5csGDBwYMHQ0NDlyxZkpCQ8Nprr12wFILQXGp6WFkHFbaywjaW28JOtLPYEUKGWpiqFqYHit6Ol/q3XIKwrb/jREt5nrbweGPh6c7a8YGp6SETM0Imuqicr6as/UAQDktTb/Ou0/uyaw+Wt1bE+Y5JDUhM9Isd5T3SWXlVvz0vFoRnq+piuxrYHg3LaWSNvSzVV0jxERK8hAQvYbSncOlnDYbCPEEYHBz8ySefzJ49m4gee+yx3t7ef/7zn2cf8OCDDyoUivfee4+Idu7cecsttzQ2Nl6wFIJQPnqJWgaouZ819lF9L6vvpboeVtVF1d2sspO5qyhmhBDvJSR6CyneQrKP4DTk36jDDcJBw2BTb4umR9vQrantaqhqP13VcapX1xfjMyrBLyYtICnWd4xKxBD18CAIr0z3YE+etuhYY35xc1lFe7Xa1S/CMyzcMzTEPTDQTR3o6u/lPGLoLycPJQjP1j5IR5pYXisrbGNFbaysgykFivIQIt2FUFcKdRWCXCnQWfBzJj8nZORQ8Q3CIf0m6u3tra+vj4uLM/4xLi5uzZo15xxz8uTJJUuWnDlAq9V2dnZ6eHhcsKBOp2tra/vZdSiVF0u73oGBtj7rm7ajZ9Sjv/yHjB49nbMRaJeOjB9OJEZdekZEBom69UREnToiRh06ppeoR0c9euo1UM8g69RRp47aB1mPnrwdyMtJ8HMkP2chwJminIXMEApxFUJdBbefP7kD3TRARES9uj7DRWYl9+r7DMxARP19/YJK0DMDEXUP9jJiffp+PTN06Xr69f19+oEuXU+3rrdzsKtjoKttoKPfMOjrNMLf2TfI1T/EPXBR+PRw95Cgs5cJ7e/Hpm/Dxfr7Jf2ggCAcJheiyd5xk73jiMjADKe66k911lZ31R04fUjTq23sbe4Y7PJwcPN0cPdwcPd0dHdTubirXJ2VTk5KRxels6NC5SA6qESVk9KBiERJUClUSqWSiBwUDo7ihTtdHRUODgoVEQlE491ovBtR2E9/1dJPNd2stofV91FNMzt8mpr6WesANfWzXj2NUJGHo+CuJHcHclcJLkpyFMnTkRREbipBKZCLiojI04GISCRyV/308VSlIOef3xpOCnIcws3irrKmzVxcnRx9PEbwrTmkIDSGlpvbT0v+e3h4tLa2nn/MmQOMkdbW1naxINy+fXtUVNTZ30lPTz8/XI2+ylm3oeP7oVynrRKI6JJzpNyI3IiCFETGm36ADAOk6SQNUd4Q6jtL/9/encc0sa0BAD9t2QqURSg8tBotilQrxo1ESGN8cakRSFxyfUHjcwET9aloNdGooIIxuK+BsGk0IVAgEkOQWBOXECCAEXAjxZLaCha0LUOtpXSZ98dcCZei1HvbTu18v7+mZ77OfCyTr3PO9Bwa/QeHZ1oRA//z9L44zcuKE410HPnhNIYVBViQL478LCjUSguw0FgWFGhGwRYay4wQ0iKkRahr5Gif7PlpAXAyJkKxCI0ekbbSaJiXXsf4Osjo1XvR9AxcT0cYA++j04Zo+DAdmWi4mY6MNIQQGqIj6/eekWEaGv5BDRmm4ya7RxIYCLERYo9cwiaETMhs+PMSGheOKDrwyB3+V+5/LhF3hGazecJ4f39/On2CQm9XIQwPD6fRaBiGBQcHI4S0Wm1kZOSYGDabjWEYsT0wMEC0/OiAQqHQ/q7R//77j/+xdtoZDJzE4Q/LgL8Bukadx/55bH+1axQ4g2O7Ru26Ifb19eVyuS9fviRetrW1xcaOfcCPx+ONDpg+fbq/v79DUgQAAACcx96e4fT09JycHKVS+fz5c7FYnJaWhhBSq9WrV6/+/PkzQmjnzp1VVVVPnjz5+PFjdnZ2enq6E7MGAAAAHMTex/ZEItGXL18EAkFISEheXt68efMQQjiOG41G4rnTuXPnFhQUHDp0SKvVbtiw4ciRI07MGgAAAHAQe78+4UDw9YnfEYwRugMYI3QHMEboDkgYIyRXRkaGTCYjOwuqO3XqVEtLC9lZUN2VK1ckEgnZWVDd7du3KyoqyM6C6u7fv19YWOioo/0GhbC1tVWn05GdBdV1dHSo1Wqys6C6zs7O3t5esrOgOplMJpfLyc6C6hQKRVdX18Rx9vkNCiEAAADgPFAIAQAAUBoJD8s0Nzdv27bNaDTaGa/ValksFjGhESALhmFMJtPHx4mL2oAJ6XQ6b29vRz0gAP6er1+/MhgMJhNmiieTwWCwWCwj05n9RE1NDY/H+3kMCYUQIaRUKk0mk+vPCwAAgFI4HM6En+DJKYQAAACAm4AxQgAAAJQGhRAAAAClQSEEAABAaVAIAQAAUJobFUIcx48ePTpp0qSQkJD9+/dbLBbbmJqaGi6XGxAQsHz5cqVS6fokPV5ZWdnixYv9/PzYbPaePXuGhoZsY3g8XvR3R48edX2SHq+2tjZ6lKamJtuYxsZGPp/PZDIXL178+vVr1yfp8YqKiqL/ilhpdbRVq1aN7E1NTSUlT8+j1+szMjIEAkF0dPToSXyGh4fT0tKCg4PDw8NzcnLGfW9JSUlUVBSLxVq3bp3t3+uHcLchFou5XG5vb69arZ43b15BQcGYAI1Gw2KxHj58aDKZ9u3bl5SUREqenu327dtPnjwxGAwfPnyIi4vLzMy0jfH19W1sbJTJZDKZrL+/3/VJerzy8vLExETZd9++fRsTYDKZOBxOUVGRxWLJzc2Ni4sjJU/PNjAwMPInOHz4cGJiom0Mj8cTi8VETE9Pj+uT9EgYhmVlZVVXVyOEpFLpSPuFCxfi4+MxDJPL5VOmTKmrqxvzxs7OThaL9eLFC4PBsH79+j179th5RjcqhGvWrDl//jyxXVhYuHTp0jEBeXl5I/+LKpXKy8vr06dPLk2RYnJycpKTk23bfX19e3t7XZ8PdZSXl69YseInAbW1tdOmTbNarTiOG43GoKCglpYWV2VHRXPmzCkpKbFt5/F49fX1rs+HCsxm85hCOGfOnLKyMmL7xIkTmzZtGvOWY8eOpaamEtstLS0sFmt4eNiec7lR16hUKuXz+cQ2n8+3nVC1q6uLWAcRIRQZGTlp0iRYlcJ5rFZrXV3d0qVLx92bkJAwY8aMzZs3Qwe1kzQ3N0+ePHn+/Pm5ubm2wwRdXV18Pp9YFcvHxycmJkYqlZKRJiU0NTUpFIqNGzeOuzc1NZXD4aSkpEAHtVPhOC6TySasEaMDdDqdSqWy5+BuVAi1Wu3IfDksFkuj0eB//bK/VqsNCAgYeRkUFKTRaFyaIpWcOXNmcHDw4MGDtrvEYnFDQ4NEIrFarUlJScQHN+BACxYskEgk7e3tV69ezc/Pv3DhwpgAuBZcqbi4eNOmTeOuinr58uX6+vqGhoaZM2euXLnyFwalwC/S6/VGo3F0jbBdD2d0EfHz8/P29rbzunCjQhgeHo5hGLE9MDDAZrPHLAMbFhY2ODg48pKIcWmKlHHlypXS0tK6urpxp7VMSUmJioqaOXNmSUnJu3fv4F7E4WbNmhUfH89ms5cvX56ZmWm7+h1cCy6j1+vFYvGOHTvG3SsUCqdOnTpt2rRLly75+PjU19e7OD3qCAwMZDKZo2tERETEmJiwsLCRAL1ebzKZ7Lwu3KgQxsbGdnR0ENsdHR2zZ8/+SUBPTw+GYdHR0S5NkRoKCwuvXbv26NGjqKgoe+JxmKXPmcZ8HCTExsa+evWK+M0PDQ1JpVLb6wU4REVFRVRU1I/GCEaj0WDGSueKiYmxv0Z0dHSEhobaFsvx/cPxTAd68OABh8Pp7OyUy+UxMTF3794l2rds2UI8CIBhWEhISGlpKYZh27dv37hxI6n5eqbi4uLAwMDKysrW1tbW1ta3b98S7cePH6+qqsJxvL29va6uTqVSdXd3b968mc/nm0wmUlP2QFVVVW1tbWq1+tmzZ1wu9+zZs0T7wYMHa2trcRy3WCxcLvfixYs6ne748eNLliwhNV9PJhAIRh7iI1y/fv3WrVs4jisUisrKyp6eHqVSeeTIkcjISGJAB/xzbW1tzc3NCKH79++3traazWYcx2/cuBEXF9fT09Pe3h4REfH06VMcx4eGhlJSUuRyOY7j3d3dgYGBjx8/VqvVQqHw0KFDdp7OjdY2Sk5O7uzsFAqFFotlx44dW7ZsIdr7+/uJNZuCgoKqq6tFItGBAwcEAkF+fj6p+XqmN2/ezJ49+9y5c8RLHo937949hJBGo9Hr9Qgho9F4+vTp7u5uf3//hISEmpoaWCHL4d6/f3/y5EmVSsXhcHbt2iUSiYh2jUZjMBgQQnQ6vbq6eu/evefOnVu4cGFZWRmp+Xqsvr4+g8GwdevW0Y2Dg4MMBgMhZLVab968uW/fPi8vr4ULF0okktDQUJIy9TQikWhgYGDRokXE9wWfPXsWEBCwe/duhUIRHx/v5+eXlZW1bNkyhBCO4319fcSTCjNmzLhz586BAwe+fPmydu3a7OxsO08H9/IAAAAozY3GCAEAAADXg0IIAACA0qAQAgAAoDQohAAAACgNCiEAAABKg0IIAACA0qAQAgAAoDQohAAAACgNCiEAAABKg0IIAACA0qAQAgAAoDQohAAAACjt/y+wXkZHusu5AAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"using Plots\n",
"\n",
"p = plot(x, real.(ψ), label=\"real(ψ)\")\n",
"plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n",
"plot!(p, x, ρ, label=\"ρ\")"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"The `energy_hamiltonian` function can be used to get the energy and\n",
"effective Hamiltonian (derivative of the energy with respect to the density matrix)\n",
"of a particular state (ψ, occupation).\n",
"The density ρ associated to this state is precomputed\n",
"and passed to the routine as an optimization."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"E, ham = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)\n",
"@assert E.total == scfres.energies.total"
],
"metadata": {},
"execution_count": 11
},
{
"cell_type": "markdown",
"source": [
"Now the Hamiltonian contains all the blocks corresponding to kpoints. Here, we just have one kpoint:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"H = ham.blocks[1];"
],
"metadata": {},
"execution_count": 12
},
{
"cell_type": "markdown",
"source": [
"`H` can be used as a linear operator (efficiently using FFTs), or converted to a dense matrix:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ψ11 = scfres.ψ[1][:, 1] # first kpoint, first eigenvector\n",
"Hmat = Array(H) # This is now just a plain Julia matrix,\n",
"# which we can compute and store in this simple 1D example\n",
"@assert norm(Hmat * ψ11 - H * ψ11) < 1e-10"
],
"metadata": {},
"execution_count": 13
},
{
"cell_type": "markdown",
"source": [
"Let's check that ψ11 is indeed an eigenstate:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "2.5509757446571617e-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.0002234290790088227"
},
"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.1"
},
"kernelspec": {
"name": "julia-1.6",
"display_name": "Julia 1.6.1",
"language": "julia"
}
},
"nbformat": 4
}