{
"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.778722e+02 1.896014e+02\n",
" * time: 0.0008480548858642578\n",
" 1 1.697476e+02 1.356059e+02\n",
" * time: 0.004106998443603516\n",
" 2 1.309580e+02 1.410881e+02\n",
" * time: 0.007585048675537109\n",
" 3 5.074785e+01 9.512703e+01\n",
" * time: 0.01159811019897461\n",
" 4 2.256826e+01 5.890020e+01\n",
" * time: 0.015218019485473633\n",
" 5 2.005575e+01 6.553625e+01\n",
" * time: 0.017641067504882812\n",
" 6 4.395756e+00 7.551468e+00\n",
" * time: 0.020124197006225586\n",
" 7 3.555850e+00 4.648838e+00\n",
" * time: 0.022117137908935547\n",
" 8 2.701007e+00 1.240495e+01\n",
" * time: 0.02421116828918457\n",
" 9 1.979938e+00 4.201658e+00\n",
" * time: 0.026283979415893555\n",
" 10 1.679678e+00 3.302526e+00\n",
" * time: 0.027875185012817383\n",
" 11 1.505921e+00 2.009918e+00\n",
" * time: 0.029410123825073242\n",
" 12 1.375611e+00 1.633936e+00\n",
" * time: 0.03139901161193848\n",
" 13 1.195248e+00 6.281966e-01\n",
" * time: 0.03289508819580078\n",
" 14 1.163141e+00 7.036873e-01\n",
" * time: 0.03447699546813965\n",
" 15 1.151819e+00 4.701986e-01\n",
" * time: 0.03608107566833496\n",
" 16 1.146140e+00 2.700090e-01\n",
" * time: 0.03762006759643555\n",
" 17 1.144899e+00 2.216117e-01\n",
" * time: 0.039214134216308594\n",
" 18 1.144454e+00 9.266961e-02\n",
" * time: 0.04081916809082031\n",
" 19 1.144372e+00 9.631160e-02\n",
" * time: 0.04247403144836426\n",
" 20 1.144130e+00 2.901727e-02\n",
" * time: 0.04359912872314453\n",
" 21 1.144075e+00 2.163725e-02\n",
" * time: 0.045310020446777344\n",
" 22 1.144058e+00 2.442621e-02\n",
" * time: 0.04686713218688965\n",
" 23 1.144049e+00 1.874257e-02\n",
" * time: 0.04839301109313965\n",
" 24 1.144041e+00 6.993482e-03\n",
" * time: 0.0499269962310791\n",
" 25 1.144040e+00 3.926619e-03\n",
" * time: 0.0510561466217041\n",
" 26 1.144038e+00 3.672024e-03\n",
" * time: 0.052146196365356445\n",
" 27 1.144037e+00 2.248227e-03\n",
" * time: 0.0536952018737793\n",
" 28 1.144037e+00 1.058749e-03\n",
" * time: 0.05547904968261719\n",
" 29 1.144037e+00 8.329317e-04\n",
" * time: 0.05704617500305176\n",
" 30 1.144037e+00 8.115895e-04\n",
" * time: 0.058730125427246094\n",
" 31 1.144037e+00 7.980699e-04\n",
" * time: 0.0602719783782959\n",
" 32 1.144037e+00 3.373970e-04\n",
" * time: 0.06185102462768555\n",
" 33 1.144037e+00 1.595301e-04\n",
" * time: 0.06347298622131348\n",
" 34 1.144037e+00 1.599723e-04\n",
" * time: 0.0650780200958252\n",
" 35 1.144037e+00 1.024221e-04\n",
" * time: 0.06667208671569824\n",
" 36 1.144037e+00 1.069512e-04\n",
" * time: 0.06821513175964355\n",
" 37 1.144037e+00 3.880261e-05\n",
" * time: 0.06931400299072266\n",
" 38 1.144037e+00 1.725712e-05\n",
" * time: 0.07048702239990234\n",
" 39 1.144037e+00 7.350636e-06\n",
" * time: 0.07209014892578125\n",
" 40 1.144037e+00 8.843218e-06\n",
" * time: 0.07326698303222656\n",
" 41 1.144037e+00 5.881860e-06\n",
" * time: 0.07446813583374023\n",
" 42 1.144037e+00 3.302919e-06\n",
" * time: 0.07616400718688965\n",
" 43 1.144037e+00 2.048117e-06\n",
" * time: 0.07770919799804688\n",
" 44 1.144037e+00 1.012953e-06\n",
" * time: 0.07932615280151367\n",
" 45 1.144037e+00 6.288406e-07\n",
" * time: 0.08087801933288574\n",
" 46 1.144037e+00 5.185970e-07\n",
" * time: 0.08242106437683105\n",
" 47 1.144037e+00 2.330390e-07\n",
" * time: 0.08398604393005371\n",
" 48 1.144037e+00 2.145736e-07\n",
" * time: 0.08563518524169922\n",
" 49 1.144037e+00 1.615454e-07\n",
" * time: 0.08724117279052734\n",
" 50 1.144037e+00 1.333685e-07\n",
" * time: 0.08934617042541504\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown (in Ha):\n Kinetic 0.2682057 \n ExternalFromReal 0.4707475 \n LocalNonlinearity 0.4050836 \n\n total 1.144036852755 "
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"basis = PlaneWaveBasis(model, Ecut=500, kgrid=(1, 1, 1))\n",
"scfres = direct_minimization(basis, tol=1e-8) # This is a constrained preconditioned LBFGS\n",
"scfres.energies"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"## Internals\n",
"We use the opportunity to explore some of DFTK internals.\n",
"\n",
"Extract the converged density and the obtained wave function:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ρ = real(scfres.ρ)[:, 1, 1, 1] # converged density, first spin component\n",
"ψ_fourier = scfres.ψ[1][:, 1]; # first k-point, all G components, first eigenvector"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"Transform the wave function to real space and fix the phase:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ψ = ifft(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n",
"ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"Check whether $ψ$ is normalised:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"x = a * vec(first.(DFTK.r_vectors(basis)))\n",
"N = length(x)\n",
"dx = a / N # real-space grid spacing\n",
"@assert sum(abs2.(ψ)) * dx ≈ 1.0"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"The density is simply built from ψ:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "1.0922720651564026e-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+gvaeTAAAgAElEQVR4nOzdd3xV9f0/8Pc5997sQfbNTiBAyE5IWEmAsBVFKGpVXBVr626t1kFdtT+1Wm2to1psS8W9UJC9ScImkEFCErLIurnZe9x7z+f3x/WbIjOBz7nz9fwrCYf3OXnknLzy+XzO5/MRGGMEAABgr0RzXwAAAIA5IQgBAMCuIQgBAMCuIQgBAMCuIQgBAMCuIQgBAMCuIQgBAMCuIQgBAMCuIQgBAMCuIQgBAMCumSEICwoK1qxZM/LjdTqdbNcCI6XX6819CUB6vR5rIpqdwWCQJMncV2HvJEkyGAy8qpkhCIuKirZv3z7y4wcGBuS7GBihgYEB/Ao2O51Oh1/BZqfX6zn+CoYrYzAYOLaR0DUKAAB2DUEIAAB2DUEIAAB2DUEIAAB2DUEIAAB2DUEIAAB2DUEIAAB2DUEIAAB2TWnuCwCAy2BEBW2srJOVtYmuDizKk8V7U7ibYO7rArARCEIAy9U2SB+WSv8ulRhRnJcQ4kz6AdreYDjWwqLHCPdOFG8ZKyrRrQNwdRCEABbqy0rpNwcNi0LENbMU0/wFIurvH3JwUCoUCp1EP5yR3j4pvVko/WumItkHrUOrpNPp7rnnnsHBQXNfiHV44okn0tLS5KiMIASwOL16umuvobSDrZunnOp/gZBTibQsQlwWIX5ULl2zRf9IrOKZJDQMrU93d/d3333373//29wXYgXef//9goICBCGAXWgdpMVb9bFjhGPLlA6XS7c7x4sLQsTrturP9LB30xUKtAytjYODw0033WTuq7ACW7dula84/ooEsCANfSxjg35OoPDhTMVlU9BI7Uy7FysrutnPdxn02JoCYPQQhACWokdH12013BElvpw2uqadu4o2LlQO6NmD+7E9EMCoIQgBLILE6PY9hhRf4cpG+xxE+nKu8kQre60ArUKA0UEQAliExw8Z+vXs/XTFFVdwUdK38xTvFkvf1yALAUYBQQhgfptr2bfV7Iu5yqucFBjsKnw5R/GrHEN9L+N0aQC2D0EIYGYtA/TLHMN/ZinGOHCoNtVfeCBG8Yt9BiQhwAghCAHMbGW24fYoISuQ29SHZxLFXj39rQgdpGBSzz777GuvvWb8WKfTzZgxo6mp6RLHL168uLCw0CSXdhkIQgBz+qJSqulmf5x85UOD51OKtHa24uUThupuNAvBdHp6enp7e40ff/jhhwkJCQEBAZc4/le/+tWqVatMcmmXgSAEMJseHT1+SPr7jJFOGRy5se7Cb+IUvz2IRiGM2vr16ysrKz/66KOnnnrKYDAYDIYvv/zy+eefX7NmzdDQkPGYkpKSt956a9WqVWd/8WzvvffenXfeSUQ9PT3GpXNKSkp27txJRFu3bi0rKyOia6+99vDhw9XV1Sb71i4GQQhgNn88bpgXLMxUy7IezBMJYkkH21iLRiGMzuuvv37dddedOHHCz8/PYDAsXrz4iy++CAsL27lz56JFixhjRPThhx8ODg5GRUWtX7/+hhtuOKdCbW1tVVXVlClTiKi9vf2hhx4iouzs7NWrVxPRu+++u3//fiJSKpUzZszYsmWLqb/D82CJNQDzKOlga8qkwuUqmeo7iPTWdMVD+w1zg5ROPHteQV537DFo+k3358vPx4r3Tjy3RXTdddcZR/u+/fbbtra2zZs3C4KwcuXKlJSUvXv3zp49+4033jAeedddd4WGhlZXV0dERAz/9/z8/LFjxyqVl8+XiRMnFhQUcPtmrhSCEMA8njwsrUpSBDjLeIqFIUKcl/BesfRYPPp+rMajcWKHCbejmDjmAl8cXts6Ly+voaFhwYIFxk8bGhpKS0tnz569du3aV155RZIkV1fX9vb2urq6s4Owp6fH1dV1JGd3dXWtq6u7ym/h6iEIAczgcDMrbGdfzZW9pfbqFHH2D/p7o0UPuVqewFmqr/mXTnd0dDR+4OzsPHPmzHfffXf4n1xcXJqamh566KGCgoLw8HAiioiIMBh+srafr69ve3u78WMnJyedTnf2Af39/c7OP/4B2NbW5ufnJ+v3MhL4OxHADJ4+Yng+WXSUv8dyoqewMETEVAq4MgsWLNi+fXtPT4+Xl5eXl5eDg4PBYGhpaXF0dAwMDCSi3bt319TUnPO/0tLSqqqqenp6iMjPzy8gIGB47wiNRnPs2LGEhATjp/n5+cahRPNCixDA1LbVs7peuj3KRH+GvjhZnLxO/+to0V/ObliwSWlpaS+88MKUKVNSUlKGhoZKS0t37do1adKklJSU1NTU0NDQgYGB8ePHn/O/PD09Z8+evWPHjqVLlxLRmjVrVq5cKYpiT09PSkrKM888M2nSJCLq6urKz89fuHChGb6xn0IQApjas0cN/y9VvMrV1EYu3E24ZZz4l0LDa1Pwzgxc3qZNm5ycnIY/ffDBB++6665Tp045OTmNHz/e2Gu6adMm41z4+Pj4np4eFxcXInrppZcE4cd+3ccee+ytt94yBuG8efMqKytfeeWV3Nzcb7/9drhfdO3atStWrPDy8jLxN3g+BCGASe2oZ716Wh5p0lGJpxLFpG/1TycqvBxNeVqwSu7u7ud8xc3NLTU19eyviKKYmJho/NjDw2P4sOEDFi5cuHPnTq1W6+/vT0QKhcLf39/T03M4BYmoqqrq2WefleNbGC0EIYBJvZpveDJRNPHrECGuwpJw8b0SadUV7fEEcAWGl1szSkpKOjspiegvf/mLaa/oovBUAJjOiVZW1kk/H2uG5+7JRPGdk4Z+venPDEBENGXKlNtuu83cV3FhCEIA0/l/J6TfxYvcF1QbiYmewlR/cU05Xh+Fy5g3b97mzZvlqMwYW7hwYUlJySWOeeihh9avXy/H2S8BQQhgIqe72N5GaeV5q3iYzJOJ4huFEvZngkv77W9/GxcXJ0flDRs2ODg4GF8ZvZgHHnjg6aefNi7kZjIYIwQwkXeKpV9Gi27mm9g+3V/wd6Ifzkg3hOMvYLiovr4+vV5PRNnZ2Tqdrre3d+fOnYmJiXffffexY8c+++yz0NDQ+++/3/j6aG5u7s6dO9vb2xMTE2+//fbhZdW+/fbb3NzcmJiY5OTk2tpa43qk77///t1332084LXXXnvwwQe1Wm1OTs4dd9yxfft2Z2fnjIyMmJgYFxeXPXv2ZGVlmexbxvMAYAo9OlpbLt0XbeYn7sEY8d1i9I7Cpbz11lsnT54kom3btq1cuXLz5s0JCQl/+tOffvnLX77yyisJCQnffffd448/bjz4s88+U6vVaWlpX3311XDIvf76688++2xiYmJNTc3y5cs//fRTIhoYGNi9e/esWbOMxzzzzDM9PT3l5eXGZWvWr1+/Y8cO4z/NmjVLpr7Zi0GLEMAUPj4tzQ4Uw93MvHrWTWPFxw8ZTnWw6DHmX8cLLqj1X380dDSb7HSu0xa5pi++2L+OGzfuvffeI6KhoaEXX3yxpqbG2Ld5yy23vP3220T0zjvvEFF/f39WVlZkZOSHH37o6Oj48ssv79y5MyUlhYiqq6sHBweJqLy83NHR8dI7FBqNHz/+u+++4/UNjgSCEMAU/lEi/XWa+eezO4h0b7T4jxLprenmvxi4IM/r75EG+012OqXXpZb6HB4s9Pf3Hz9+vIODg/Hj1tZW49dffvnl1atXe3p6KpVKg8HQ0NDg6ura1dU1PMswJSXlwIEDRNTb23v2JMJLcHV1NS7PZjIIQgDZ7W1kOomygiyiEXb/JDH+G/1LqQosw22ZlP4h5r6E/1EoFBf82OjYsWMffPBBYWGhh4eHXq93cXGRJMnT05OIOjs7vb29iWh49W1/f//29nZJkkRRJCJnZ+f+/v/lfV9fX2hoqPHjlpaWkTQcOcIYIYDs3iuRHowx9ST6iwlyEbKCxE9PY6QQrlZra6urq6txmvyaNWt0Oh0ROTk5zZ071zhZvqmpae3atcaDIyMjPT09S0tLjZ+mpKSsW7fO+HFHR8dwVyoR5efnT5061ZTfCFqEAPJqHaRtddIHGRbU/rovWvzDUcOvJ+HvYLgqs2bN8vb2TkxM9Pb2Dg4OHl41dPXq1StWrAgMDAwICFi4cKGx5ScIwo033rhx40bj9In33ntvxYoVLS0tHR0dMTExK1eunDdvHhFJkrRz587nnnvOlN8IghBAXmvLpevCxDEO5r6Os8wLEu7rp/w2luhtIc1UsCA7duxQqVRE9Mwzzwx/8frrr58/f77x47CwMOPWS46Ojvv27SsvL3dycgoPD+/o6DCuOxoeHp6Tk2PsBb377rsnTJhg/I+PPvro8uXLf/vb3yoUitjY2BMnTnzxxRcvv/zyoUOHhpf5Xr9+/dSpUyMjI035LSMIAeS1pkz6q4W9mSIKdNd4YU2ZRby/A5ZmOJPOfrdFpVIZ05GIRFEcM2bM8McTJ040fjz8xa1bt27fvj0qKur48eM7duwYXnR0woQJDzzwQHl5eXR0tPErXl5ezs7OZ292UVtb++qrr8r1vV0E+kYAZHS0hXXpaJba4hpev5ggfnJaGjRc/kiA0UpNTY2NjW1ra5s2bVphYaFxAwqj+++/fzgFiSgyMvKOO+44+/8+/PDD48aNM921EhFahACy+k+Z9IsJoqW8J3OWCHch3lvYcEa60bQbQoE98PHx+cUvfjGSI8ePH3/+vr6mh2cAQC4DBvqiQrprvOXFIBERrZwo/rsM744CIAgBZLPhjJTsK4SZezWZi1kWLh7UMo3ppm4DWCgEIYBcPj3NVoyz3EfMWUnXh4lfVKBRCOeqqKh47LHHHn744U2bNpn7WkzBcp9SAKvWPki7G6VlERb9iK2IEj9FEMJPVVRUzJw5MyIiYsaMGQ888MDHH39s7iuSHV6WAZDF11XSgmDR05KmD55vbpBwVw8r62QTPC20/9YOPb3npdb+dpOdbkFk1o3R15/9lbfffvv2229/5JFHiMjFxWXVqlW33367ya7HLBCEALL4pEL6TZxFNweJSCHQzWPFzyrY8ykIQkvxcOovu4dMt+R0iHvQOV8pLi5euXKl8ePJkyeXlZUNLxBqqxCEAPw19LGiNnZNiBX87lgRJa7YbXg+xQou1U4EuanNewH9/f3d3d3Gj7u6utzc3Gw7BQljhABy+LSCLYsQHa1h2ZYpfoIg0JFmZu4LAQvy0UcfDQ0NEdHq1asXLFhg7suRHVqEAPx9USH9eYo1xCAREf18rPBVlZTmZzUXDHIbN25ccnKySqUyGAwbN2409+XIDi1CAM6qu1lND5tpecuqXcyNkeJXVQxNQhi2bNmygwcPrl+/vrCwMCwszNyXIzsEIQBnX1WxZRGi0nqerURvwVGkvBZEIfyPu7u7PUSgkfU8rABW4usq61vA82cRwtdVmFAIRES//vWvY2JizH0VJmVljyuAhavtZZXdLCvQavpFjW6MFL+sRIsQiIhWrFgRFRVl7qswKQQhAE9fVbKl4dbUL2qU4isIAp1oRRaCPbK25xXAsn1TLS23tn5Ro+URwlfoHQW7ZJVPLIBlauyjUx1sbpCV9YsaLY8Uv6tGixDsEYIQgJv1Z6RrQ0WVdT5VaX5Cl47KOpGFYHcwoR6Am+9rpHsmWGcMEglE14cJ39ewJxKsskVrjZRKZVdX17hx48x9IVagubl51qxZMhVHEALw0aOj/U3siznWGoREdEO4+NJxwxMJVvwtWBcPD4/q6urBwUFzX4h1kG9eI4IQgI9NtVJ6gOCuMvd1XIU5QcJtu5mmn9TO5r4UuxEcHGzuSwCMEQJw8n0NuyHcuh8olUgLQsSNZ/DuKNgX635uASyETqKtddL1YVb/QN0QJnxfg/dlwL5Y/XMLYAn2adh4TyHQxdzXcdWuCRX3aaRevbmvA8CEEIQAHGw4YwvNQSLydKBUX2FnPXpHwY7YwqMLYHabatniMBuZdbA4TNxYi95RsCMIQoCrVdrJenWU4G0jQbgkTNhwRkISgv0YxfSJysrKbdu2eXl5LVmyxNn5wq9X5+XlHTx40N3dPSsrKyQkhNNFAli0jWfYdWGCjcQg0TgPwU0l5LeyJB+b+Z4ALmWkLcKcnJzJkyefPHnyX//6V2Zm5gVngD7++ONLliwpKCjYvn37P//5T67XCWC5NtZKi0NtKjMWhwroHQX7MdIW4Ysvvrhq1arHH3/cYDCkpqZ+/fXXK1asOPuArVu3fvzxx0VFRb6+vjJcJ4CF6tLRkWY2J8imRhkWh4rPHTOsSrKpbwrgYkZ0ow8ODu7cuXPp0qVEpFAorr/++s2bN59zzFdffXXnnXe2trauX7++urqa+4UCWKbtdVJ6gOBmzQvKnG9moFDcwbT95r4OAJMYUYtQo9EwxoKCgoyfBgUF5eTknHNMZWVlUVHRvn37oqOj77nnnjfffPPOO++8YLX29vaTJ0++/PLLw1+59dZbLzGgqNPpdDrdSK4T5GP8KQi2MxDGzcYztCCITHOLGn8EkiT73AaBKEtNm87oVoyV+1TWR6fTiSLaymam0+kMBoNCobjskUql8rK/uEYUhIwxIhquJYri+Y/i0NDQ0NDQ0aNHRVHctGnTbbfdtmLFigtepV6vHxoaam9vH/7KwMDAJZ5tSZJM8OTDpRl/CgjCczCirfXiE7EmukNN+VNYFCxsrRNujcCjdy7jDxu/lMxL+j9cqo0oCNVqtSAITU1NERERRNTY2DjcOhwWFBTk4+Nj/EMpPT29s7NTo9FccD1ZPz+/5OTk119/fYSXODQ05OjoOMKDQSY6nc7R0RFBeI6CNuaqMkzyNdH9KUmSg4PDSP4KvnqLI9izx/UqB0cRP/PziKKoUtlWb7i1EUXRYDDwioYRNfCdnJzS09M3bdpERIyxLVu2zJ07l4j0ev2ZM2eMmbxgwYKSkhLj8cXFxc7Ozv7+/lwuEcBibalji0JsMyhCXAV/ZyGvFe+Ogu0b6Vujf/jDH2699VatVltSUtLR0XHLLbcQUVVV1YQJE5qbm319fW+77bY333zzjjvuiI+Pf//991944QX8xQQ2b2ud9Fi8KdpnZrEoRNhSy1J9bTPpAYaNdMh34cKFu3fvVqlUWVlZBw4ccHV1JSK1Wv3JJ5+4u7sTkYuLy6FDhzIzMyVJWrt27e9//3sZrxrAAvTq6Wgzm6W22ZxYGCJuxaKjYAdGsbJMYmJiYmLi2V9xd3e/7bbbzv70vvvu43ZpAJZtdwNL87O1iRNnm6kWitpYxxCNcTD3pQDICS8BA1yhrXXSwhBbfoIcFTQ9QNjVgEYh2DhbfowBZGXDb8oMWxgsbq3D+zJg4xCEAFeioov16SnOVnacuJhFocIWBCHYOgQhwJXYVs/mB9v+tMqJnoIoUGknshBsGYIQ4ErsqGfzg20+B4mI5gQK2+sRhGDLEIQAo2ZgtKdRsrEdJy5mXrCwA0EINs0unmQAvo42s2BXIdDF3NdhEvODxT2Nkg6vjoLtQhACjNp2u+kXJSJfJ4p0F440o1EINgtBCDBqOxqk+cF29OzMD8YwIdgyO3qYAbjo1VNeC8u03ZXVzjcvWNyBafVguxCEAKOzt5FN9hVcR7E6odWbqRbyW1kXtscGG4UgBBidHfXSPHvqFyUiJwVN9Rf2NqJRCLbJnv6sBeBhVwP7IMOkQThkGNpZva+4tay05bSLg/NE76ikgPjpwammvIY5QeLuBnZ9mCnPCWAiCEKAUWgdpOoeNtmEW/Rl1x54N+/fkZ7hqYFJMwOnGUTpdEf1h/lrPy/+9uHUX0Z5RZrmMuYECvfloEUItglBCDAKO+ulmWpRaZIGocSkPx/8e2lbxe+nPpyiTiCi/v5+BweHacGpt8X87IfT2x7f9fzKxBXXRy00wcWk+gm1vaypnwKcTXA2AJNCEAKMwu5GNifIFM1BiUmvHPhbx2DXB4vecFScux+gKIhLxi+aEpTy2x1/GNAP3hS9RO7rUQiUESDubZRuHmtf46NgD3BPA4zCzgZTBKHE2PPZf+7V9b08c9X5KThM7er/t3n/77uyTV+f2iD3JRHRnCBhVwNmE4INQhACjFR9L2sfZHFesgfh5yXftg90/jHzaZVCdekjA1z9/jrvT58Uf13YXCz3Vc0JEnY1IgjBBiEIAUZqRwObEySKMudgcUvZlyXfP5v+mFJUjOR4fxffp6Y98secv3QOdsl6YfHeQtcQq+lBFoKtQRACjNQu+ftFu4a6X8h57clpjwS4+o/8f00NmjwnPPOVA3+T78KISCCaHSjuQaMQbA6CEGCkdjewOYHyBuGHJz5OD0m7gjmCv0y6s6WvbVdNthxXNWxOkLATi46CzUEQAozI6S7GiMZ7yhiEp9sr99Xu/0X8bVfwf5Wi4jdpv/5H3n8G9APcL2zYnCABLUKwPQhCgBHZ3cCyZG4Ovn30w3sT7/BwdL+y/x7nF53oH/dp8Td8r+psUR6CRFTRhSwEm4IgBBiRvRo2S84g3Fm9r18/cO24+VdT5NfJd31Xtrmxp4nXVZ1vdiAahWBrEIQAI7KnkWXJ9qaMgRk+zP/4ockrReGqTuHr4rN84nVrCj/jdWHnQxCC7UEQAlxeeScTica6yxWEO6r2qt38E/xjr77UjdFLDtQfre9uvPpSFzQ7ENPqwdYgCAEuT9bmoMTYZ8Xf3hn3cy7VXFUuN4y/5vOSdVyqnS/KQ1CKdBrDhGBDEIQAl7enUcYBwj1ncpxVzskB8bwK3jRpyZ6a3KZeLa+C55ilRu8o2BQEIcDl7dXI9cooI/Zx0Vd3x9/CsaaHg/viqPmfl3zHsebZZmGYEGwLghDgMso6mUgUKc8A4dHGE4IgTg2azLfszdE37Kja2zPUy7esUVaQsKsBexOC7UAQAlzGXjkHCNeVbfzZxMXcy3o7e00Nmry5cif3ykQ01l1QiUJZJxqFYCMQhACXsVfDZsvTL9rU21zUfGpueKYcxZdNvHZd2UaJyRJXs9TCPg2CEGwEghDgMrI1bKZaliD8vnzzwrFznJROchSP9Y12U7ke05yQo/jMQGEfhgnBViAIAS6lspsZGI3z4B+EOoNuc+XOJVELuVcedsOEa9aVbZSjMl4cBVuCIAS4lL2NbJY8zcHdZ3KjxkSGegTLUdxobvjMouZTcsyjGO8pGBhVdyMLwRYgCAEuZZ9sS4xurtxx/XgZm4NE5KR0nBsxc0vlbjmKzwzEMCHYCAQhwKXslWcqfVNvc0V79fSgUe87OFqLIudsrdrFiH9izVQLe9E7CjYBQQhwUXW9rFfPJsiwB+G2qt1zwjNVChX3yueY6BPloHA42VzKvfIstAjBViAIAS5qTyObpRbl6BjdVrVn4dgsGQpfwILI2VtkmFA4aYzQpWP1vchCsHoIQoCLkqlf9GTLKUZsks8E7pUvaEFk1t4z+wcNQ3zLCkQZASIahWADEIQAFyXTDMKtlbsXjZ3LvezF+Dp7T/SJ2l93mHvlWYEYJgRbgCAEuLCmfmoeYLFenINQLxn2nMldEDmbb9lLWzg2a1vVHu5lM9VCThOCEKweghDgwvZppAy1yH2E8JjmRJhHiL+LL+e6l5QRMi1fW9Q91MO3bKK30NDHtP18qwKYGoIQ4MKyNSxThn7R3Wdys8LTuZe9NGelU4o64UD9Eb5lRYGm+wv7tdiJAqwbghDgwuQYINRLhv11h2eFzuBbdiRmh6XvrsnlXjZTLWbjfRmwcghCgAvoHKLKLpbswzkIjf2ivi4+fMuORHrI1ILmk726Pr5lM9VYfRusHoIQ4AJymthUf0HF+/kwS7+okbPSKTkgPrfuEN+yaX5CWSfr1vGtCmBSCEKAC8jWSJlqzk+HGftFjeToHXUQKcVX2I93R8GaIQgBLmBfI/83ZfKa8sM8gs3SL2o0I3hKvraIe+/oTLWQrcH7MmDFEIQA5+rXU2E7m+LHOQhzag9lhk7nW3NUXFTOCf4xhxvy+JbF+zJg7RCEAOc62MwSvAUXJc+ajNj++iMzgtN4Fh299JCp++s5LzEzPUDIa2UDBr5VAUwHQQhwLjn6RctaK1xVzrJuwzsSM4KnHKg/qpd4pparkiaNEY40o1EI1gpBCHCunCYpI4Dzo5Fbfyg9ZCrfmlfAx9krxD2osLmYb9mMAAG9o2C9EIQAP6GX6LCWzQjgP0CYHjKFb80rkx4yhfskiky1kNOE92XAWiEIAX7ieCsLdxe8HXnWbOrVtg20T/KZyLPolUoPmZpde5BvzUy1uL+JGdAmBOuEIAT4iZwmlsG9OVh3aEbwFFGQY4vfURs7JlwUxKqOGo41fZ0o0EUobEMSglVCEAL8hBxrbe+vPzLDMvpFjWaEpOXyXoAbw4RgvRCEAP/DiHKbpAyuQTigHyhpKUsJSOBY8ypNDZp8qOEY35rYmxCsF4IQ4H/KOpmzQgh15RmEeU0F0T7jXVTOHGtepST/uNPtlV1D3RxrZqqFfY14XwasEoIQ4H/k6Bc9WH9satBkvjWvkoPCIcEvNk9TwLFmpLugFIWKLjQKwfogCAH+J7eJ8e0XJaIjjcctLQhJnt7RDDWGCcEqIQgB/idbw/mV0TNddXpJH+EZyrEmF9ODUw81HGPEM7fSA4RcDBOCFUIQAvxI00/tg2zSGJ5BeLDh2LTgVI4FeQl0C3BWOle0V3OsifdlwEohCAF+lK2RMtSiyLVn9HBD3pTAFJ4V+eHeOxrvJWj6mLafY0kAU0AQAvwoR8PSufaLDugHi1tKJ6sTOdbkaFrw5INcg1AUaHqAkIu11sDaIAgBfsR9TZl8bdF4r7EWNXHibIn+cafbK/nu05seIGKYEKwOghCAiKhHR2WdbLIvzyA82ngiNTCJY0G+HBUOk3wm5GuLONbEMCFYIwQhABHRfi2b7Cs4KnjWPKKx6CAkotTApCONJzgWnOonnGxnfXqOJQFkhyAEICLK0Uh8+0Xb+ttb+lonekdxrMldqjrpKNcgdFRQgrdwCJv0glVBEAIQ/fimDM/H4dm2SSgAACAASURBVKjmRIo6QRQs+hGL8hrbNdit7WvhWBOrb4PVseinFMA0dBIdbeG8Ge9RTX6q2qL7RYlIFIQUdcIxTT7HmhlqIUeDF0fBmiAIAehYCxvnIXg68KyZpymw8AFCo8nqRL69oxkB4iEt0yMKwXogCAEol/fEiarOMwpBDHJTc6wpk7TApGOaExzXWvNypFA3oQCb9IL1QBACUG4T56n0xxrz0wKTORaUT4Crv5uDW0V7FceaGQGYRAHWBEEI9s64GW86100njmnyLXZBmfNNVify3ZIpXY3Vt8GaIAjB3pV1Mhclz814JSYVNhcnBcTxKii35ID4402FHAumB2CTXrAmCEKwdzm8t14qbTvt7+Lr5TSGY01ZpQQk5GtPGpiBV8Gx7oJSFCq70SgE64AgBHvHfYAwT1OQok7gWFBuHo7ualf/0tYKjjVnBAg5mE0IVgJBCPYuh/eu9MebCpMDrCkIiShFnXC8ieswITbpBeuBIAS71tRPLQMsht9mvDpJX9xSmuAfw6ugaSQHJPAdJsxAixCsxyiC8PPPP7/33ntXrVql0WgucdiXX3759ttvX/WFAZhCjkZKDxA4bsZb0lIa6hHs7uDGraJJJAXEnWw5pTPoeBVM9BEa+ljrIK96ADIaaRD+/e9/X7Vq1axZs1pbWzMyMgYHL3yDHzp06KGHHnrppZf4XSGAjHKb2AyuS4zmNRUkB8RzLGgariqXMI+QktYyXgUVAk3xE/Zjk16wBiP6FWAwGN5444133333jjvueP/9993d3b/55pvzDxscHHzggQdefPFF3hcJIBfum/Ee1xSmWNsAoVFyQHwe12HCGdikF6zEiIKwrq6utrY2KyvL+Ons2bNzc3PPP+yll15aunRpTIyVjY6A3erVU0kHS/PjFoRDhqHSttPx1jZAaJSiTjjexHuTXgwTgjVQjuSgxsZGNzc3R0dH46d+fn5Hjx4955j8/Pz169cfOXLk4MGDl65WV1eXnZ29fPny4a/8/ve/j4+/aG9Sf3+/QsF1v1QYvb6+PkEQBIFn48ns9jaJCV6iNNjXx6lgQUtxhEcYG5L6hniV/In+/n69Xi/T4zDONaK0tbyju8NBwWf18QQ3ym9TtXX3OdnW4zs4OCiKokqlMveF2DWdTmcwGCTp8n3vTk5OoniZJt+IgtDJyWloaGj406GhIWdn57MP0Ov199577/vvvz8clpfg7e0dHh7+85//fPgrUVFRTk5OFztep9Nd4l/BNPR6vZOTk40F4eE2lqlmHO+u4vayFHW8fLcrY8zBwUGmIHQip8gxYVW9tYn+sZwKUswYqbDbMZPr7BSzEwQBQWh2CoXCYDCM5Fm7bArSCIMwODh4cHCwubnZz8+PiGpra4ODg88+oLKyMj8//4477iCigYGBtra2cePG7dixIzIy8vxqLi4uYWFhN99880hOTUSiKI7kOwFZGX8KNhaEuVr9o7EKkd87o/nNJ1fELpfvdhX/j0z1E/3jCppPJqu5veyToWb7m4VZQTb1/Mr9U4CREEWRMcbrpzCiKn5+fpmZmZ988gkRdXR0bNq0admyZUTU1ta2bt06IoqMjDx16tT27du3b9/+xhtveHp6bt++PSQkhMslAshBL9FhLc/NeHUGXWlreZzvJF4FTS8pIO6ElucwYXoANukFKzCiFiERvfrqq0uXLt29e3dxcfGCBQumT59ORKdOnfrZz37GGFOpVGPHjjUeWVtbq1Aohj8FsEwn2liYm+B9+b78kSpuLYscE+aicr78oZYqwS/2xZzXdQadSsGn3y9TLd6zz2BgpLCprgSwNSMNwhkzZpw6dergwYOBgYHJyT9utJaSklJaWnrOkVOnTj1y5AjPawSQQY6G88pqJ5oKE/2tZseJC3JROYd6BJ9qK4/34/Piq58TBTgLJ9tZgjeSECzXKDpYvb29r7322uEUJCInJ6cJEyacc5iTk1N4eDifqwOQDfe1tk9oi6xo66WLSfaPP8F1EgXWWgPLh/FesFO5TRLHqfQ6SX+qlVtDyowSeQ8TZqixWz1YOgQh2KPTXUwUhAh3bkF4qrUs1CPYVeXCq6C5JPrHlrSU6SVuexNmBAjZaBGCZUMQgj3K1rCZXAcI87XFSVY+QGjkqnIJdg8sbSvnVXC8p2BgrKYHWQiWC0EI9iiX9xKj+doia39TZliif2yBtphjwfQAEcOEYMkQhGCPsrm+MioxqbilNN7PimcQni3BP7ZAe5JjwYwADBOCRUMQgt1pHqCmfhbrxS0Iy9sr/V18PRzdeRU0rwT/mILmYolxmwifocYwIVg0BCHYnWyNlB4gcJzina89mcBpfU5LMMbR08fZu7KjhlfBJB+htgeb9ILlQhCC3cnRsAw1zzu/QFtsS0FIRIn+sfn8ekcVAk31xya9YLkQhGB3cppYJr83ZRixouZiXjs2WIgE/xjOw4RqvC8DlgtBCPbFuBlvKr/NeGs665yVzr7O3rwKWoIk/7h87UlG3KIrE8OEYMEQhGBfDmpZko/AcavYAu3JROtfWe0cfi6+jgqHuq4GXgWn+gmF7axfz6seAE8IQrAv2RqJ71T6Am1xovWvrHa+xIA4jsOELkqK9xION6NRCJYIQQj2JVvDMrm+KZOvLYr3t8EgTPCbxDEICb2jYMEQhGBHdBIdaWbT/bm1CDW9Wp2kD3EP4lXQcsT7xRQ281xfJlMtZmOTXrBICEKwI3ktLMpD8HTgVrBAa2vviw4L8wwZ0A9o+1p4FUwPEA5qmR5RCJYHQQh2ZJ+GZXIeIDyZYIv9okQkkBDnN6mwuYRXQS9HCncXTrShdxQsDoIQ7Eg27yAsbC62gT0ILybeL6aQ6+rbGCYEy4QgBHvBiPY3SekB3O75zsEubV/LOK8IXgUtTYI/72FC7E0IFglBCPbiZDvzcRIC+W2dW9hcEuc7SSHwm5NoYSZ4j2vsaeoZ6uVVcFagmKORkIRgaRCEYC/494tqi+P9bWTrpQtSCIqJPlFFLdyGCQNdyF0llHYgCsGyIAjBXmRrOG/GW9BsU5tOXFCCXyzfYcKZgcI+9I6ChUEQgr3I0bBZgdyCcEA/WNlxJtp7PK+ClinBPyafaxBmYJgQLA+CEOxCRRdjRJHu3IKwpLVs3JgIJ6Ujr4KWKcZ34un2Sp1Bx6vgTLWwpxFBCJYFQQh2YR/X5iDZ9AzCszkrncI8Q061lfMqON5TkBhVdyMLwYIgCMEucJ9KX9hcYg9BSEQJfjEFGCYEm4YgBLuwr5Fx3HRCYlJJa1mcry2/Mjos3i+miN/6MoTZhGB5EIRg++p7WbeORY/hFoSn26v8nH08HN15FbRkCf4xhc0lEuMWXWgRgqVBEILt26thMwNFjh2jhc3FNrn10gV5OY3xdPSo6TzDq2Csl9A2yBr6kIVgKRCEYPuyNTz7RYmoQFsc72cX/aJG8f4xBfzWWhOIMgLEHDQKwWIgCMH27eU6QEhERc0lNrzW9vni/SYVarkOE6rROwoWBEEINq55gDT9LMGbWxA29GhIEALdAngVtHwJfjxbhEQ0K1DYi9mEYDEQhGDj9jZKGQEixxHCAm1xoj01B4ko1CNYJ+k4btKb5CPU97GWAV71AK4KghBs3N5GzlPpC5uL4+wsCIkozjea46KjCoGm+wv7NNiuHiwCghBs3B7eQVigLU6w6U0nLijebxLv3lERvaNgIRCEYMtaB6mulyX7cAvCzsGu1v62sWMieBW0FvH+nHern41hQrAYCEKwZXsbpfQAQcF1gDDWL1oU7O7BmeA9TtOr7Rrq5lUwxUeo7mGtg7zqAVw5u3uewa7sbWSzAnne5EXNJQn2N0BI/7dJb3FLKa+CSpGm+wvZGCYEC4AgBFvGf4DQntaUOUeCX2wh10VHMUwIFgJBCDarbZCqu1kKvwHCQcNQZUeNzW/GezHxfpP4bkOBYUKwEAhCsFn7NFK6WlDyu8eLW0rtYTPei4n1iy5rqxgyDPEqmOorVHazNgwTgrkhCMFm7Wlks9Q87/ACbbGd7EF4Qc5Kp3CPkNK2Cl4FlSJNwzAhWAAEIdis3Q0sK4jzEqNx9rTW9vni/ScVcp1NODtQ3IPeUTA3BCHYptZBqunhOUAoMelkyym72nTifPF+nGcTZgUKuxoQhGBmCEKwTbsbpIwAngOEFR3Vfi4+no4e3CpaoQQ/zpv0TvYVzvRg0VEwMwQh2KbdjSwriO8A4Um72nrpgrydvTwc3Tlu0qsUKT1A2NuIYUIwJwQh2KbdDSyL/xKj9h6ERJTgH8t30dGsIHE3hgnBrBCEYIO0/aTpZ4n8BgjJuKaMfyzHglYqwS8mX3uSY8GsQGE3hgnBrBCEYIN2NUiz1CLHJUbruxtJENSu/twqWq0Ef85BmOQjNPUzTT/HkgCjgyAEG7S7kfPEiYLm4kQ0B4mIKMQ9iDGpqVfLq6AoUIZa3NOAYUIwGwQh2KA9jWw25wHCk/a51vYFxfvF5HOdRDEnSMAwIZgRghBsTX0vaxtkcV58d6XHAOH/xPvF8J1WPydI2IlhQjAfBCHYmu31bG6QKPLLwY7Bzo6BzgjPMG4VrVyCf0wB12HCWC+hV8equ5GFYB4IQrA1OxvYXL4DhNriOL9oUeBZ06pFeUU297V2DnbxKigQZQWJaBSCuSAIwdbsbmRzg3kPEKJf9CyiIMb6RfPdm3AuekfBfBCEYFOKO5hSoLHuPIMwX3sSr4yeI9Evlm/v6LxgYWeDhCQEs0AQgk3ZWc/mc20O9ur66robJnhHcaxpAxIDYvlu0hvuJrirhJPtiEIwAwQh2BTuA4RFzSXRPuNVopJjTRsQ7T2+uvNMn47nNPi5QcLOegQhmAGCEGyHgVG2RprNe63tBD/0i55LpVBN9I4qbinlWBPDhGAuCEKwHUebWbCroHbmWTNfizVlLizBP7agmecw4ZwgcZ9G0mGFGTA5BCHYjh0NbB7XftEhw1BFR1WM70SONW1Ggn9MfhPPIPR1orHuwiEtGoVgaghCsB3b6qQFITxv6eKW0kjPcCelI8eaNiPOb1Jp2+khwxDHmvODhR1YdBRMDkEINqJbR3mtLFPNd+IE+kUvylnpFO4ZWtpWwbHm/GBxWx1ahGBqCEKwEXsapal+givXtzsLmk/GYzPei0v0j+W7JVOGWihqZ+2DHEsCXB6CEGzE9no2P5jn/ayXDCUtZdh04hIS/GPzm4o4FnRS0IwAYU8jekfBpBCEYCO21XGeSl/aVh7sHujm4Mqxpo1J9I892XJKLxk41pwfLG7HbEIwLQQh2IK6XtY6yJJ8eAbhiaaiJP84jgVtj7uDW6BbQHk732FCYRuCEEwLQQi2YEsdmx/Mc+slIjqhLUoMQBBeRqJ/3AmuvaPx3kKfnlViSyYwIQQh2IIdvJcYNTDDyeZTGCC8rKQAzkEoEM0LwrujYFIIQrB6BkY76qWFITyDsKytItAtwMPRnWNNm5QUEFfUUiIxnq+3LAxB7yiYFIIQrN7hZhbiKgS5cB4gTMQA4Qh4OLj7ufiWt1dyrLkwRNzdIA3h1VEwFQQhWL0ttdI1oZy3j8/XFiVhgHBkkngPE/o60QRPYX8TGoVgIghCsHqb69giriurSUwqaj4VjwHCkeE+TEhEi0KELXVoEoKJIAjBurUMUHknmxHAs0VY3l7p6+zt5eTJsaYNS/SPK2wu5jtMuChU3FKLFiGYCIIQrNvWOikrSFRxvZHzNAXJ6nieFW2al5Onj7MX32HCKX5CfR+r70UWgikgCMG6balj13B9X5SIjjcVJgck8K1p21LUCcebCjkWVAg0L1jcindHwSQQhGDFJEbbeU+cMDBDUXMJNp0YleSAhOMankFIxmFC9I6CSSAIwYodbWG+TkKYG9clRltPB7qpPR09ONa0eUkBcYXNxXwXHV0UIu5owIb1YAoIQrBiP5yRrgvj3C+apylICcAA4eh4OLgHuqnL2k5zrBngTFEeQg4mUYD8EIRgxX44wxaHcr6H85rwpsyVSAmIz9MU8K25OFTceAZNQpAdghCsVUMfO9PDpvvzbBHqJP2p1vIEPwwQjlqyOj6viXMQXhcm/HAGLUKQ3SiCMDc3d8GCBUlJSU8++eTg4Ll7SGu12qeffnrmzJlpaWkPPfRQU1MT1+sEONcPZ9jCEFHJ9W+54pbSUI9g7EF4BRL940pay3QGHceaKb5Cj57KO5GFIK+R/hZpaWlZvHjxLbfc8umnnx44cOCFF14454CysrKBgYEXX3xx9erV9fX1y5Yt43ylAD+1sZYt5j1AeKKpMBkDhFfEVeUS7hFa3FrGsaZAdG2osBHvjoLMRhqEa9euTUtLu+eee2JiYv785z+vXr1ap/vJn34ZGRl//etfs7KykpKS3njjjQMHDnR3d8twwQBERIMG2tsoLeS6shoRHdPkp6gxg/AKpagTjmny+dZcHCpsrMUwIchrpL9HCgsL09LSjB+npqa2tbXV19df7OC8vLygoCB3d2xhA3LZ1cASvQUfR541B/QD5e2V2IPwik1WJ+bxDsJ5weIhLesc4lsV4CeUIzxOq9VGR0cbP1apVG5ubk1NTREREecf2dDQ8Mgjj7z11lsXK3X69Onvv/8+MjJy+CsffPDBjBkzLnZ8b2+vIHDuAYPR6u3tZYxZzg/iu0rlvADq6RngWPNI0/Eoz0j9gL6HejiW5ai/v9/BwUGhUJj7Qi4s0jn0dHuVtqPZRenMsew0X9WGir6loZbSLhwcHBRFUaVSmftC7JpOpzMYDHq9/rJHuri4iOJlmnwjDUJPT8/e3l7jx5Ik9fX1jRkz5vzDtFrtvHnzHnnkkZtuuulipcaNGzd37ty//vWvw18JDw+/xLPNGHNzcxvhdYJ8XF1dLSQIGdHmBv32axVubk4cyxaXlk0JTrHkm02hUFhyEBJRjO/E071VM4KncKz5s7HS1ibl7ZMs5btWqVQIQrMzBqGTE5/fACPtGo2MjCwr+3EYvKKiQqFQhISEnHNMa2vr/Pnzb7755qeeeuoSpQRBcHNzG3sWS36wwQIda2FuKproyTmVj2pOTFYn8a1pbyarE482cu4dXRIubKzFEjMgo5EG4YoVKzZt2lRRUUFEb7311tKlS11dXYlo7dq1O3bsIKK2trb58+enp6c/+uij7e3t7e3tBgPP9ZYAhn1fI90QzjkF2wc6mvtao32i+Ja1N6nqpGOaE3xrBrkIEzyFfRq8OwpyGWkQTpo06bnnnps8eXJISMjBgwf/8pe/GL++YcOGAwcOENG+ffuqq6s///zzcf+nurpaposGO/ddNbshnPP7okc1J5IC4kUBS0xclfHe49oHOpv7WviWvSFc/L4GTUKQy0jHCInod7/73YMPPtjd3e3n5zf8xS+//NL4wdKlS5cuXcr56gDOU9HFWgbYFD/OLcJjmoLJmDhx1URBSA6Iz9MULBw7h2PZpeHCoi3SW9PJIsaoweaM7u9fJyens1MQwPS+q2E3hIsi79+IeZoCDBBykaJOOMZ7rbVJYwRnBZ1oRe8oyAIdQWBlvq+RuPeLnumqJ2JhHsF8y9qntMDko43HGXEOrSXhAnpHQSYIQrAmTf1U1M7mBnNuDx5uyJsSmMK3pt0KclM7K50rO2r4ll0WIX5ThRYhyAJBCNbk22rp2lDRgfdte7gxb0oQgpCbKUHJhxvy+Nac5i90DlFJB7IQ+EMQgjX5ukpaHsG5OThkGCpqLpmsTuRb1p6lBaYcaTzOt6ZAtDRCQKMQ5IAgBKvRMkB5LWwR74W287Unx3lFuqpc+Ja1ZykB8SWtZf16ngvgEdGNkeI31RgmBP4QhGA1vq2WFoaIzqOY8jMihxuPTwlM5lzUvjkpnSb5TDjeVMi3bEaAoO2n011oFAJnCEKwGt9USTdG8p9IdgRvysggLTD5SCPnYUJRoKURwtfoHQXeEIRgHdoH6VAzuyaU8x3b3NfSPtA53nsc37IwJSjlcAPnYUIiWh4hflOF3lHgDEEI1mFdjTQvWHTl3S96qCEvNTBJtIxdNWzJ2DHhA/qBhh4N37KzAoUzvayqG41C4AlBCNbhswrp1rH84+pgw9FpwZO5lwWBhKlBkw/UH+FbViHQjZHi55UIQuAJQQhWQNtPx1rYtbz7RXUG3fGmwqmBCEJZTA9OPVB/lHvZW8eKH5ejdxR4QhCCFfi8Uro+jP/7osebCiM9wz0c3TnXBSIiSgtMLm4p7dP18y2brhb6DFTUjkYhcIMgBCvwWYV06zj+9+qBhiPTg1O5lwUjJ6VTjO9E7tsTCkQ3RwqfVaBRCNwgCMHS1fSwym42L0iGAcL6YzOC07iXhWFy9Y6OEz+vYGgSAi8IQrB0n1awGyNFJe9btarzjIEZIseEc64LZ5kRPGV//WGJd2Yl+QhOCjqsRRQCHwhCsHSfVUi3juV/ox6sPzojeAr3snC2QLcAD0eP8rYK7pVvGSd+gt5R4ARBCBYtr4V16yhdLUe/KCZOmMKM4LT9vCdRENEdUcJnFdIQohB4QBCCRftvuXT3eO7b0VPnYNfpjqqUgATeheFcM4LTcuoOcS8b4S7EjBE21yIJgQMEIVguvURfVkq3R/FvDubWHU4LTHZQOHCvDOeI84tpG2jnvsQMEd01QfxvOYYJgQMEIViujbXSBE9hnAf/IMyuPZgZMo17WTifKAgyNQpvjhT3NErNnPd6AnuEIATL9VE5u2sC/1t0QD+Yry2aGoQBQhPJDJ2WU3uQe1k3FV0bKn6OV2bgqiEIwUK1DtKuBunGSP636KGGY7G+0W4OrtwrwwVNDkis6Khu62/nXvnu8eJ/sdwaXDUEIVioteXS9WGih4p/5ezag5mh6Bc1HZVClRaYfKCB/8z6OUFCywCdaMVIIVwVBCFYqA9LpV9G878/9ZLhUOOx9JCp3CvDJWSGTs+uPcC9rCjQPRPF1aVoFMJVQRCCJcrRMAOjDBmmD55oKgzzCPZx9uJeGS5hWtDkwuaSXl0f98orJwifV0h9eu6FwY4gCMESrS6V7ovmP32QiHbVZM8Oy5ChMFyKq8ol0T82t+4w98rBrkJ6gPhlJRqFcOUQhGBxOodofY10e5Qs/aI5dYdmhk7nXhkua0545u6aHDkq/zJaQO8oXA0EIVictaela0JFPyf+lY9qjod7hga4+vEvDZeTHjI1X1vUNdTNvfK1oeKZHuxQCFcOQQiWhRG9Vyz9epIsd+aumpws9IuaibPSabI6UY7eUYVA904U3ytGoxCuEIIQLMuOeqYQaKYMr8noDLr9dYdnhqFf1GzmhGfuqsmWo/KvJomfVUjtg3LUBtuHIATL8vZJ6TdxstyWhxrzorwifZ295SgOIzE9OK24pbRjsJN7ZbUzXRsqrsHkergiCEKwINXd7IBWunWcTP2i2Vnh6Bc1Jyel45TAlGwZllsjoodjxXeLJQkDhTB6CEKwIG8XS/dMEF2U/Cv36voONRybHZbOvzSMxvzIWduq9shReZq/4O1Im+uQhDBqCEKwFL16+qhcul+e12T2ntmfEpDg6eghR3EYualBk+u6G+TYlYmIHooR/15kkKMy2DYEIViK1aekuUFihLsc0+hpa+WuBZFZclSGUVEIijnhGduqdstR/NZx4qlOOo6lR2GUEIRgEfQS/a1I+q08r8k09WqrOs9Mw75LlmFBZNaWyl2M+MeVSqSHYsQ3CvHKDIwOghAswheVUpQHTfWXpTm4pXL33IhMlUKGnSxg9CZ6RzkrnYqaT8lR/NeTxK11Uk0PGoUwCghCsAhvFEpPJChkKr6jes/CyDkyFYcrMD9ytky9o+4qWjlR/FsRGoUwCghCML9t9czAaEGILM3BAu1JURCjfcbLURyuzIKI2Xtqcgf0A3IUfyRW/KhcasXkehgxBCGY3x/zDE8myrLXBBGtL996fdQieWrDFfJ18Yn3n7RLnjW4g1yEGyPFvxbi9VEYKQQhmNn2eqYdoJ+PleVW7BrsPthwdMHY2XIUh6uxZPyi9eVbZCq+Kkn8R4nUIkuDE2wQghDM7KXjhhdSRIU87cEtlTtnhEzxcHCXpTpchSmBk9sHOsraKuQoHuYmLI8U3zqJRiGMCIIQzGl7PdP0y9UcJKKNFduXoF/UIomCsDhq/g+nt8lU/w9J4vslUhtGCmEEEIRgTi/mydgcPNFUyIji/KJlqQ5X7bpxC3bVZPfp+uUoHuYmLIsQ38BIIYwAghDM5vsaqVtHt8jWHFxXtmnZhGtlKg5Xz9vZK0WdsLVql0z1n08WPyiR6noxpxAuA0EI5mFgtOqo9NoUhUxvi2p6tcebCheNnStLdeDkpugbvj61QWKyZFWwq3DPRPGl45hTCJeBIATz+Hep5OdEC+WZO0hEX59av3jcfGelk0z1gYt4v0kejm4H6o/IVP+ZJMV3NVJxBxqFcCkIQjCDPj398bj0+lS5lpLp0/Vvrdq9bOJimeoDRzdG3/Dlqe9lKj7GgZ5IUKw6gkYhXAqCEMzg1XzDTLWQ6itXc/CH01unBKb4u/jKVB84mh02o7FHc6q1XKb6D8WI+W1sVwMahXBRCEIwtcpu9o8S6bUpct17esnwTekPN0Uvkak+8KUQFMsmLP761AaZ6jsp6M1p4sP7DTo0C+EiEIRgar85ID2RoAh2las5uK1qd7B7IBYXtSJLxi860ni8tqtepvpLw8UId3qnGEkIF4YgBJPaVMvKOtlv5Nl3kIgkJn1a/M1d8bfIVB/k4KpyWTbx2k+Kv5HvFH+dpnjlhEEjy5RFsHoIQjCdXj09tN/w9gyFg2z33faqPT7O3on+sXKdAORxU/QN++sO13U3yFR/gqdw70Tx0QOYXw8XgCAE03nmiGF2oDA/WK5OUYlJH5/86q64n8tUH+TjqnJZOuGaz4q/le8Uz6coitrYt9XoIIVzIQjBRA5q2TfV7C+yTZkgoh3VVPl/dQAAFoJJREFUez0dPVPUCfKdAuRzY/SSfbUHGno0MtV3VNCHMxWPHJDasQAp/BSCEExhwED37DO8NU30dpTrFEOGoQ/zP/lV8p1ynQBk5uHgfnP0Df888ZF8p5juLywLFx47hA5S+AkEIZjCE4cMiT7C8kgZ77evSzdE+0TF+8XIdwqQ282Tlha3lBY1n5LvFK+kKXKb2NdV6CCF/0EQguw217L1Z9h7M2TsFO0e6vmy5Pv7ktActG6OCodfxN/6Xt6/Gck1/91NRZ9lKR7cb6jpwRR7+BGCEOSl6aeV2fpPZiu8ZOsUJaL/FHyWFZ4e4h4k4znAJBaOnTOgH8ipPSjfKSb7Cr+NU9y112BAFAIRIQhBVnqJbtml//UkRYZarjdFiai8vXJXTfbd8bfKdwowGVEQH5p879vH/jWgH5DvLL9PEJUCPX8Mg4VAhCAEWf3ukMFRpFVJMt5mEmNvHv7Hr5Lu9HT0kO8sYEop6oRE/9j/FHwm3ylEgb6cq/ysAoOFQIQgBPl8WiFtrmNfzFXKtAG90Xdlm5SictE47DtoUx6efO+26j3l7ZXyncLbkb6Zp3hwv+FkO3pI7R2CEGSR28R+e9Dw7TzFGAcZz6Lta/lv4eePT3lAIDnDFkzOw9F9ZcKKNw69Z2Ay9l4m+Qh/TlMs3W5olrEXFqwAghD4O93Fbt5p+GiWMs5LxnySGHvlwN9uil4S7hkq31nAXBZHzXdzcF1b9KWsZ7l7gnjbOOHaLfpevaznAYuGIATOGvtowWbDy2mifLvPG31Rsk5n0N0Wu1zWs4C5CCQ8Nf3R78o2n2yRcVohEb0wWRHnLdy6y6DHcKG9QhACT9p+mrtJ/6to8a7x8t5ap9srPy9e94f034kC7mGb5evs/bsp9/8p980+nYzbRghE/8xQMGK378GECjuFXyLATcsAzdusv2Wc+GSivPdV11D3c9l/fiT1l2pXf1lPBGaXGTo9NTDp1YNvyTfFnohUIn0zT9mrZ2gX2icEIfBR38uyNupvCBeeS5b3ppKY9KfcN2eGTp8bMVPWE4GFeDT1Vx0DnXIPFjqI9NVcZccQu32PYQhZaGcQhMBBRRebvdGwPFJ4abKM66gZ/eP4GolJWE3NfihFxQsZv99QvjW37rCsJ3JS0IYFSka0aLO+SyfrqcCyIAjhau1tZBkb9M8miy+kyJ6CX5Z8d7D+6PMZT2Bo0K54O3u9NPPp1w+9XdhcIuuJHBX0aZZi4hgha6O+thcDhvYCv03gqrxbLN2yS/9xlvJOmd+OIaItlbu+Lt3wlzkvuDu4yX0usDTRPuOfS3/i2X0vl7adlvVECoH+ka64dZw47XtDjgZZaBcQhHCFunV05x7DP09JuUuUc4Nkn8++rWrP6vy1b859KQAvyNirFHXC76Y88PSelyo6quU+1+Px4n9mKW7cqX+jUEIY2jwEIVyJQ1qWvE7vrKQDS5Rj3WVPwa9Pbfgw/+M35/wR+0vYuczQ6Y+m/up3O58r0J6U+1wLgoVDNyjXVUuLNusb++Q+G5gTghBGp1dPjx8yLN2uf32K+EGGwkUp7+kkxj44/t/1p7e8s+BVrCADRDQrbMZz6Y8/l/3qvtoDcp8r3E3Ys1g5I0BMWqf7VymahjYLQQijsK5aiv9Gr+2nwuWqZRGy3zxdg91P7nmxpLXsnfmv+rv4yn06sBYp6oTX57z4zrF/fXD8vxKTd66DUqTnU8Tt1yg/LJWyNurz25CGNghBCCNyrE2cvdHwQp70Yabio9kKXyfZz1jYXHzflsfGeoa/MfePHo7usp8PrMp4r7Grr3mzrL3id7uea+ptlvt0Cd5C7vXKW8aKizbrf7VfqEdPqW1BEMJl7G9i127V35GruiNKyFumnCP/ezED+oG/H139Qs7rj6Ted3/KLxSC7LMywBp5Onq8nvVCqjrpvs2PrS/fIuvSM0QkCvTrSWLpzaoAZ0pdT/fnGqq70Tq0EQhCuLAhiT45LU1br79zr2FZuHhi8eDKiaKsOwsSkcTYlspdd2x4oFfXu2bx2zOC0+Q9H1g5URBXxN741vyXN1fufGDrE4XNxXKf0UNFf0xmhUvJx5FSv9PfuNOwpxFxaPVkftUBrNCRZrb2tPR5hZToIzyTKF4XJooC9fTIe1KJSXvP7P/45FfOSqcXM5+M8Z0o7/nAhkR4hr638LUd1Xtfyn0zyitiReyNsb7Rsp7R14n+lKp4KlHxUbn0YK7BwOiO8eLtUUK4G/bFtEqjC0KdTqdSqa7+GLA0QxLlNrEfzkjfVjOVSLdHiYduUEbKPy+CiNoHOrdX7V5XtsnH2Xtl4u1oBcIVEEiYHzF7VuiMTRU7/pT7pq+Lzw3jF2WGTndUyLgxtJuKHogRH4gRD2rZR+VS6neGCDdheaR4TaiQ4C0gEq2IwNiI2vU9PT133XXX9u3bFQrFk08++dRTT51/zGuvvfbKK68YDIa5c+d+9NFH7u4XfsHh008/3bhx4yeffDLCS+zu7r5YKbgagwbKa2XZGpatkbI1LHqMcG2ouCxciPe+wCPc09Pj6uoq8Hu6Owe79tcfya49kK89mREydcn4RXL/FW8D+vv7HRwcFAoMml6KxKR9tQc2nt5+qrU8M3RaZuj0yeoEB36JODg4KIri+X/uGxjtbWTf1Uhb61i3js0NEjPVQoZaiPYURKQibzqdzmAwODnxeW1vpEH4zDPPHDt2bMOGDQ0NDVOnTl23bt2MGTPOPuDQoUOLFy8+dOhQaGjosmXL4uLi/vznP1+wFILQXGp7WVknFbWxonaW18pOdbBJY4QMtTBTLcwOFL0dL/V/uQRh+0DnqdbyfG3RiaaiM111aYHJ6SFTM0Kmuqicr6as/UAQjkpzX8ueM/tz6g6Vt1XE+E5MDoiP95s03nucs/KqfnteLAjPVtXN9jSyfRqW28Sa+liyr5DkI8R5CXFewgRP4dLPGoyEeYIwODh4zZo18+fPJ6Lf/OY3fX19//znP88+4P7771coFO+88w4R7d69+5ZbbmlqarpgKQShfPQStQ5SywBr6qeGPtbQR/W9rKqbqntYZRdzV1H0GCHWS4j3FpK8hUQfwWnEv1FHG4RDhqHmvlZNr7axR1PX3VjVcaaqs6ZP1x/tMz7OLzolIGGS70SViCHq0UEQXpmeod587cnjTQXFLWUVHdVqV78Iz7Bwz9AQ98BAN3Wgq7+X85iRv5w8kiA8W8cQHW1m+W2sqJ2dbGdlnUwp0FgPIdJdCHWlUFchyJUCnQU/Z/JzQkaOFN8gHNFvor6+voaGhpiYGOOnMTExn3/++TnHnD59etmyZcMHaLXarq4uDw+PCxYcGhpqb28f/nTMmDGX+A3b3dfb3m9903b0jHr1l/8jo1dP52wE2q0j4x8nEqNuPSMig0Q9eiKiLh0Ro04d00vUq6NePfUZqHeIdemoS0cdQ6xXT94O5OUk+DmSn7MQ4ExjnYXMEApxFUJdBbefPrmDPTRIRER9un7DRWYl9+n7DcxARAP9A4JK0DMDEfUM9TFi/foBPTN063oH9AP9+sFuXW+Prq9rqLtzsLt9sHPAMOTrNMbf2TfI1T/EPXBJ+Oxw95Cgs5cJHRjApm+jxQYGJP2QgCAcJRei6d4x071jiMjADDXdDTVdddXd9QfPHNb0aZv6WjqHuj0c3Dwd3D0c3D0d3d1ULu4qV2elk5PS0UXp7KhQOYgOKlHlpHQgIlESVAqVUqkkIgeFg6N44U5XR4WDg0JFRAJRmhuluRGF/fhPrQNU28PqellDP9W2sCNnqHmAtQ1S8wDr09MYFXk4Cu5Kcncgd5XgoiRHkTwdSUHkphKUArmoiIg8HYiIRCJ31Y+/PFUKcv7preGkIMcR3CzuKmvazEWlUgZ6+fCtOaIgNIaWm9uPS/57eHi0tbWdf8zwAcYGXHt7+wWDsLy8fP369Tt27Bj+ymeffZaRkXGxs7+4/qkzju0X+1d7IBDRJedIuRG5Ef3/9u4+pqmrDQD4c9s6WgRBtB2FrFLjB/iFCHPT6FBEDJlEo3W+I3URkGzZzJaA+0i2GROW8Idui3HGjwUWzbKlituyumXL0KllbgobQtKMWiyRzzIKbRUKtPfe8/5xt5u+FLG+tr2s9/n9dTx9uPdRen3ac889J0UKwL3px4EZB/s9sAO0BHF8BUtJHnB4BQtS8vfpYwglYwnXKSEgJ5SUhZkMxBCQMzCbpWYyVDwDcTQkMFQ8DQBOACeAlT9aXzB/W4TCTAGQDuB/R5qlKLds5L50+J60d0RGjUjJiATcUtIvocYo4pWAjyK0BMYpAIAxCbD/fG73UuB9QA3xSogv6DsJUgAlgJK/hH0APqBH/76EJkVApDcepSwczT+WNDOeYRiaph8aHxsbK5E8pNAHVQjnzp1LUZTb7U5ISAAAp9P55JNPTohRKpVut5tru1wurmfSoy1cuFCn0wU/NHrkP8dwaFRwIZ8sg/4PODQaPsGvY/uoQ6MoHEI7NBrUF+KYmJj58+c3Nzdzf7x161Z6+sQJfhkZGf4BaWlpsbGxIUkRIYQQCp9gR4bLy8s/+OCDrq6ua9eunTt3bt++fQAwODi4ZcuWgYEBACgrK7tw4cLPP//c3d1dVVVVXl4exqwRQgihEAl22l5lZaXD4Vi/fn1iYuKJEyeWL18OAISQ8fFxbt7p0qVLT58+XVFR4XQ6d+7c+eabb4Yxa4QQQihEgn18IoQe9fGJDz/8sLS0dPbs2WHNCk3t1KlThYWFGo3m4aEobAwGw5IlS7iPoUgoP/zww8yZM9evXy90IqL266+/OhyOoqKikBztXzBp9uzZs3fv3hU6C7G7cOGC2Rz2PcHR1L777rvGxkahsxC7S5cumUwmobMQu19++eWnn34K1dH+BYUQIYQQCh8shAghhEQNCyFCCCFRE2CyTHV1dXV19YMetw/U09OjUqnw8VVh2e32hIQEhQJXxxbSwMCAQqHgl3BCghgaGpJKpdzqIkgobrebYZikpKSHRhYXF1dVVU0dI0AhZFnWarUGX9jGx8djYnAlWoHhb2E68Pl8Uqn0oetFobCiaZqiKFzfR1gMwxBCuBVfp6ZWqx/6CV6AQogQQghNH/jREiGEkKhhIUQIISRqWAgRQgiJGhZChBBCohbsotuRYbFYPv/8c5qmi4uLJ11QcXh4+NNPP+3q6lq3bt2OHTsin2HUGxoaMhqNZrM5Pj5+x44dS5cuDYypqalhGIZrL168ODc3N7I5Rj+bzea/c/XWrVtTUlImxPh8vpqamtu3b2dmZu7ZswenkoZca2vrb7/95t+j1+sn7C537tw5bvtVAFCr1aFa+hJ1dHQ0NTU5nc7du3f7P6ny+++/GwwGhUJRUlKSlpYW+IMDAwM1NTUDAwPPP/98Xl5ekKebRhdPe3v7M888QwiJi4tbt25dS8vEzdUJIQUFBVeuXFmwYME777xz+PBhQfKMbgcOHPj2229VKpXb7V69evWkq/m99tprra2tNpvNZrNxm3Ch0Gpqaqqurrb9Y3R0NDBGr9d/+eWXCxcuPHbs2BtvvBH5JKOey+XifwXffPPN+++/H/jQ16FDh0wmExfT19cnSJ7Rp7+/Pzs7+9SpUy+//PJff/3F91+/fj0vL2/u3LkejycnJ6enp2fCD3o8njVr1lgslnnz5hUXFxsMhmBPSaaN119/fd++fVz7rbfeeumllyYE1NfXp6amer1eQkhDQ4NKpeI2gUIhNDo6yrcPHDig0+kCY2JiYnp7eyOYlOgYDIb8/PwpAiwWi0KhcDqdhJC7d+/K5fL+/v5IZSdGOp2usrIysD8jI6OhoSHy+UQ37hlBmqYB4Pbt23z/tm3bqqqquPbu3bvfe++9CT9YW1v79NNPsyxLCPniiy9WrFgR5Bmn0TfCq1evFhQUcO3NmzdfvXo1MGDDhg3ch7I1a9Z4PJ4///wz0llGO7lczrfHxsYetIjJZ599dvTo0Rs3bkQqL9Hp6+s7cuRITU2N3W4PfNVkMuXk5CQmJgKARqPRarUTBvFQCA0ODhqNxtLS0klf/frrrz/66CP/oWz0mB40zn/t2rXNmzdz7UlrBBdAURQX0NraOjQ0FNQZHyPbEOvr6+PXXVOpVHa7nfzvw/52u50PkEgkSqWyt7c30lmKRktLy5kzZyoqKgJfys3NvX//vtVqLSwsPHjwYORzi3rx8fGZmZkul+vixYsZGRlNTU0TAvyvBQBQqVR4LYTP2bNns7KylixZEvhSVlYWAPT29u7du1ev10c8NREZGxtzOp3+NSJwLNq/iMyZM0cmkwU5Xj2NJsvMmDGD+y4MADRNy2QyrrDzZDIZP0cDAHw+3xNPPBHRFEWjs7Nz+/btH3/88aRTln788UeuUVZWlpOTs3//fpVKFdkEo1xhYWFhYSHXrqysPHjw4Pfff+8fgNdCJJ05c2b//v2TvsRvMF5RUbFo0aKbN2+uXr06gqmJCLe4oH+NCHzPy2QyPoBhGIZhgrwuptE3wtTUVP5TbU9PT2pqamAAf3fU6/U6HI7AqXTo8XV2dm7YsOHtt98uKyubOjIrK0sul+O2yWG1du1am802odP/WgCAnp4evBbC5ObNm+3t7S+88MLUYSkpKVqttqOjIzJZidCMGTOUSiX/tp/0Pe9fRLiGWq0O5uDTqBAWFRWdP3+ea58/f56fiGwymQYHB7mAy5cvc2O+RqPxqaeeSk9PFyrbaNXd3Z2Xl/fqq6++8sor/v3Nzc2dnZ0AMDY2xndeunSJYZgFCxZEOstox/8jE0IuXry4bNky7o+NjY3cfwQFBQVms/nOnTtcp9vtfu6554TKNrrV1tbu2rVr1qxZfE9bW5vFYgEAr9fLsizfabVaJ33cCIVKUVFRXV0dABBC6urquBrBsuzly5dHRka4AKPRyF0+dXV1eXl5wW7V8rjze0LH4XAsWrRoy5Yt27Zt02g03d3dXH9SUpLRaOTaJSUl6enpe/fuVSqVX331lXDJRq1du3bJ5fLsf+j1eq4/Nze3urqaEGIwGBYvXvziiy9u3bo1Li7u5MmTguYbnbZv375x40a9Xr9y5UqtVmu1Wrn+rKys48ePc+1Dhw5pNJrS0tLk5ORPPvlEuGSjmcfjSUxMNJlM/p0lJSXl5eWEkMbGRo1Go9Ppdu7cOWvWrHfffVegNKPQpk2bVq1aBQDLli3Lzs4eHh4mhLS3tycnJ+t0uo0bN2ZmZt67d48QMjw8DACtra2EEJqmCwoKsrOz9+zZM2fOnOvXrwd5uum1+4TH46mvr2cYJj8/Pz4+nutsbm7WarXcBDkAaGho6OrqevbZZ7VarXCZRq07d+7wDwgDQGxsbEZGBgC0tbUlJCSo1Wqapm/dumW1WuPi4nJycoIceUCPxOVy3bhxY2hoSK1Wr127lr/PYTablUolf0f2jz/+sFgsK1aswC8iYTIyMtLW1rZq1Sr/+QodHR0URaWlpbEsazab29raZDJZZmbm/PnzBUw1yrS0tPB3+wBg5cqV3L5XLpervr4+NjZ206ZN3MZwLMs2NTUtX76c22uJYZgrV644HI7c3Nzk5OQgTze9CiFCCCEUYdPoHiFCCCEUeVgIEUIIiRoWQoQQQqKGhRAhhJCoYSFECCEkalgIEUIIiRoWQoQQQqKGhRAhhJCoYSFECCEkalgIEUIIiRoWQoQQQqL2X62XGyIR7PcFAAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"using Plots\n",
"\n",
"p = plot(x, real.(ψ), label=\"real(ψ)\")\n",
"plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n",
"plot!(p, x, ρ, label=\"ρ\")"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"The `energy_hamiltonian` function can be used to get the energy and\n",
"effective Hamiltonian (derivative of the energy with respect to the density matrix)\n",
"of a particular state (ψ, occupation).\n",
"The density ρ associated to this state is precomputed\n",
"and passed to the routine as an optimization."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"E, ham = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)\n",
"@assert E.total == scfres.energies.total"
],
"metadata": {},
"execution_count": 11
},
{
"cell_type": "markdown",
"source": [
"Now the Hamiltonian contains all the blocks corresponding to $k$-points. Here, we just\n",
"have one $k$-point:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"H = ham.blocks[1];"
],
"metadata": {},
"execution_count": 12
},
{
"cell_type": "markdown",
"source": [
"`H` can be used as a linear operator (efficiently using FFTs), or converted to a dense matrix:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ψ11 = scfres.ψ[1][:, 1] # first k-point, first eigenvector\n",
"Hmat = Array(H) # This is now just a plain Julia matrix,\n",
"# which we can compute and store in this simple 1D example\n",
"@assert norm(Hmat * ψ11 - H * ψ11) < 1e-10"
],
"metadata": {},
"execution_count": 13
},
{
"cell_type": "markdown",
"source": [
"Let's check that ψ11 is indeed an eigenstate:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "1.7955189615047208e-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.00022344388438416256"
},
"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.8.5"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.5",
"language": "julia"
}
},
"nbformat": 4
}