{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Gross-Pitaevskii equation in one dimension\n",
"In this example we will use DFTK to solve\n",
"the Gross-Pitaevskii equation, and use this opportunity to explore a few internals."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## The model\n",
"The [Gross-Pitaevskii equation](https://en.wikipedia.org/wiki/Gross%E2%80%93Pitaevskii_equation) (GPE)\n",
"is a simple non-linear equation used to model bosonic systems\n",
"in a mean-field approach. Denoting by $ψ$ the effective one-particle bosonic\n",
"wave function, the time-independent GPE reads in atomic units:\n",
"$$\n",
" H ψ = \\left(-\\frac12 Δ + V + 2 C |ψ|^2\\right) ψ = μ ψ \\qquad \\|ψ\\|_{L^2} = 1\n",
"$$\n",
"where $C$ provides the strength of the boson-boson coupling.\n",
"It's in particular a favorite model of applied mathematicians because it\n",
"has a structure simpler than but similar to that of DFT, and displays\n",
"interesting behavior (especially in higher dimensions with magnetic fields, see\n",
"Gross-Pitaevskii equation with external magnetic field)."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We wish to model this equation in 1D using DFTK.\n",
"First we set up the lattice. For a 1D case we supply two zero lattice vectors,"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"a = 10\n",
"lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"which is special cased in DFTK to support 1D models.\n",
"\n",
"For the potential term `V` we pick a harmonic\n",
"potential. We use the function `ExternalFromReal` which uses\n",
"Cartesian coordinates ( see Lattices and lattice vectors)."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"pot(x) = (x - a/2)^2;"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"We setup each energy term in sequence: kinetic, potential and nonlinear term.\n",
"For the non-linearity we use the `LocalNonlinearity(f)` term of DFTK, with f(ρ) = C ρ^α.\n",
"This object introduces an energy term $C ∫ ρ(r)^α dr$\n",
"to the total energy functional, thus a potential term $α C ρ^{α-1}$.\n",
"In our case we thus need the parameters"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"C = 1.0\n",
"α = 2;"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"… and with this build the model"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using DFTK\n",
"using LinearAlgebra\n",
"\n",
"n_electrons = 1 # Increase this for fun\n",
"terms = [Kinetic(),\n",
" ExternalFromReal(r -> pot(r[1])),\n",
" LocalNonlinearity(ρ -> C * ρ^α),\n",
"]\n",
"model = Model(lattice; n_electrons, terms, spin_polarization=:spinless); # spinless electrons"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"We discretize using a moderate Ecut (For 1D values up to `5000` are completely fine)\n",
"and run a direct minimization algorithm:"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy log10(ΔE) log10(Δρ) Δtime\n",
"--- --------------- --------- --------- ------\n",
" 1 +169.2606726480 -1.52 361ms\n",
" 2 +126.0086591069 1.64 -0.82 1.50ms\n",
" 3 +72.03783082690 1.73 -0.54 1.53ms\n",
" 4 +13.64379652598 1.77 -0.36 1.50ms\n",
" 5 +11.28133065189 0.37 -0.31 996μs\n",
" 6 +11.05645401970 -0.65 -0.64 829μs\n",
" 7 +10.55002381403 -0.30 -1.06 500μs\n",
" 8 +9.450164406687 0.04 -0.54 818μs\n",
" 9 +7.354132105128 0.32 -0.37 1.00ms\n",
" 10 +4.731470225247 0.42 -0.21 996μs\n",
" 11 +1.741687359323 0.48 -0.04 983μs\n",
" 12 +1.492955275769 -0.60 -0.22 817μs\n",
" 13 +1.323463775139 -0.77 -0.44 673μs\n",
" 14 +1.269640343253 -1.27 -0.54 664μs\n",
" 15 +1.227237297763 -1.37 -0.54 637μs\n",
" 16 +1.176962610400 -1.30 -0.96 654μs\n",
" 17 +1.150536737251 -1.58 -1.16 666μs\n",
" 18 +1.146711574614 -2.42 -1.51 633μs\n",
" 19 +1.145380388805 -2.88 -1.53 654μs\n",
" 20 +1.144422364102 -3.02 -1.86 658μs\n",
" 21 +1.144251253364 -3.77 -2.11 636μs\n",
" 22 +1.144099904777 -3.82 -2.20 660μs\n",
" 23 +1.144080967071 -4.72 -2.59 676μs\n",
" 24 +1.144053337500 -4.56 -2.74 677μs\n",
" 25 +1.144047865766 -5.26 -2.91 699μs\n",
" 26 +1.144042505909 -5.27 -2.43 679μs\n",
" 27 +1.144038676988 -5.42 -2.65 639μs\n",
" 28 +1.144037765972 -6.04 -3.63 664μs\n",
" 29 +1.144037387767 -6.42 -4.14 483μs\n",
" 30 +1.144037052483 -6.47 -4.24 499μs\n",
" 31 +1.144036955532 -7.01 -4.37 658μs\n",
" 32 +1.144036882112 -7.13 -4.22 640μs\n",
" 33 +1.144036868946 -7.88 -4.33 661μs\n",
" 34 +1.144036862487 -8.19 -4.49 639μs\n",
" 35 +1.144036857093 -8.27 -4.54 679μs\n",
" 36 +1.144036854705 -8.62 -4.82 664μs\n",
" 37 +1.144036854215 -9.31 -4.56 635μs\n",
" 38 +1.144036853822 -9.41 -5.13 662μs\n",
" 39 +1.144036853032 -9.10 -4.92 660μs\n",
" 40 +1.144036852876 -9.81 -5.49 641μs\n",
" 41 +1.144036852797 -10.10 -5.41 653μs\n",
" 42 +1.144036852774 -10.63 -6.41 487μs\n",
" 43 +1.144036852766 -11.10 -6.21 636μs\n",
" 44 +1.144036852759 -11.14 -6.10 660μs\n",
" 45 +1.144036852757 -11.68 -6.01 657μs\n",
" 46 +1.144036852756 -12.39 -6.35 638μs\n",
" 47 +1.144036852756 -12.31 -6.54 651μs\n",
" 48 +1.144036852755 -12.76 -7.05 638μs\n",
" 49 +1.144036852755 -13.34 -6.78 656μs\n",
" 50 +1.144036852755 -13.42 -7.27 664μs\n",
" 51 +1.144036852755 -13.45 -7.27 637μs\n",
" 52 +1.144036852755 -14.18 -7.48 660μs\n",
" 53 +1.144036852755 -14.29 -7.67 667μs\n",
" 54 +1.144036852755 -15.35 -8.02 627μs\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown (in Ha):\n Kinetic 0.2682057 \n ExternalFromReal 0.4707475 \n LocalNonlinearity 0.4050836 \n\n total 1.144036852755 "
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"basis = PlaneWaveBasis(model, Ecut=500, kgrid=(1, 1, 1))\n",
"scfres = direct_minimization(basis; tol=1e-8) # This is a constrained preconditioned LBFGS\n",
"scfres.energies"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"## Internals\n",
"We use the opportunity to explore some of DFTK internals.\n",
"\n",
"Extract the converged density and the obtained wave function:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ρ = real(scfres.ρ)[:, 1, 1, 1] # converged density, first spin component\n",
"ψ_fourier = scfres.ψ[1][:, 1]; # first k-point, all G components, first eigenvector"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"Transform the wave function to real space and fix the phase:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ψ = ifft(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n",
"ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"Check whether $ψ$ is normalised:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"x = a * vec(first.(DFTK.r_vectors(basis)))\n",
"N = length(x)\n",
"dx = a / N # real-space grid spacing\n",
"@assert sum(abs2.(ψ)) * dx ≈ 1.0"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"The density is simply built from ψ:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "8.31875022839291e-16"
},
"metadata": {},
"execution_count": 9
}
],
"cell_type": "code",
"source": [
"norm(scfres.ρ - abs2.(ψ))"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"We summarize the ground state in a nice plot:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=3}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3xV9f0/8Pc5997kZpN9sxMIEEI2hJEBhK0oQh1VcVW01lFbrdZBXbU/tVptraNabEtFxT1AlmySsAlkkJCELLJubvYe997z+f1x/abITOBz7nw9/0rC4X1OHjknr3w+n/P5fATGGAEAADgq0dIXAAAAYEkIQgAAcGgIQgAAcGgIQgAAcGgIQgAAcGgIQgAAcGgIQgAAcGgIQgAAcGgIQgAAcGgIQgAAcGgWCMKCgoI1a9aM/Hi9Xi/btcBIGQwGS18CkMFgwJqIFmc0GiVJsvRVODpJkoxGI69qFgjCoqKibdu2jfz4gYEB+S4GRmhgYAC/gi1Or9fjV7DFGQwGjr+C4fIYjUaObSR0jQIAgENDEAIAgENDEAIAgENDEAIAgENDEAIAgENDEAIAgENDEAIAgENDEAIAgENTWvoCAOASGFFBGyvrZGVtopsTi/Zi8T4U4S5Y+roA7ASCEMB6tQ3SB6XSv0slRhTnLYS6kGGAtjUYj7awmDHCPRPFm8eKSnTrAFwZBCGAlfq8UvrtAePiUHHNbMWMAIGI+vuHnJyUCoVCL9H3p6W3TkhvFEr/mqVI9kXr0Cbp9fq77757cHDQ0hdiGx5//PHU1FQ5KiMIAaxOr4Hu3GMs7WDfzFdODzhPyKlEWh4pLo8UPyyXrtpieHiy4ukkNAxtT3d397fffvvvf//b0hdiA957772CggIEIYBDaB2kJVsNk8cIR5crnS6VbneMFxeGitdsNZzuYe+kKxRoGdoaJyenG2+80dJXYQO2bt0qX3H8FQlgRRr6WMYGw9wg4YNZikumoInGhXYtUVZ0s5/vNBqwNQXA6CEIAaxFj56u2Wq8PVp8KXV0TTsPFW1cpBwwsAf3YXsggFFDEAJYBYnRbbuNKX7C5Y32OYn0+Tzl8Vb2agFahQCjgyAEsAqPHTT2G9h76YrLruCqpK/nK94plr6rQRYCjAKCEMDyNteyr6vZZ/OUVzgpMMRN+Hyu4r4cY30v43RpAPYPQQhgYS0DdG+O8T+zFWOcOFSbHiA8EKv4xV4jkhBghBCEABa2Mtt4W7SQFcRt6sPTiWKvgf5WhA5SMKtnnnnm1VdfNX2s1+vT0tKampoucvySJUsKCwvNcmmXgCAEsKTPKqWabvbHKZc/NHgupUhr5yheOm6s7kazEMynp6ent7fX9PEHH3yQkJAQGBh4kePvu+++VatWmeXSLgFBCGAxPXp67KD097SRThkcubEewm/jFI8cQKMQRm39+vWVlZUffvjhk08+aTQajUbj559//txzz61Zs2ZoaMh0TElJyZtvvrlq1aozv3imd99994477iCinp4e09I5JSUlO3bsIKKtW7eWlZUR0dVXX33o0KHq6mqzfWsXgiAEsJg/HjPODxFmaWRZD+bxBLGkg22sRaMQRue111675pprjh8/7u/vbzQalyxZ8tlnn4WHh+/YsWPx4sWMMSL64IMPBgcHo6Oj169ff911151Voba2tqqqatq0aUTU3t7+0EMPEVF2dvbq1auJ6J133tm3bx8RKZXKtLS0LVu2mPs7PAeWWAOwjJIOtqZMKrxeJVN9J5HenKl4aJ9xXrBSzbPnFeR1+26jtt98f778fKx4z8SzW0TXXHONabTv66+/bmtr27x5syAIK1euTElJ2bNnz5w5c15//XXTkXfeeWdYWFh1dXVkZOTwf8/Pzx87dqxSeel8mThxYkFBAbdv5nIhCAEs44lD0qokRaCLjKdYFCrEeQvvFkuPxqPvx2b8Jk7sMON2FBPHnOeLw2tb5+XlNTQ0LFy40PRpQ0NDaWnpnDlz1q5d+/LLL0uS5Obm1t7eXldXd2YQ9vT0uLm5jeTsbm5udXV1V/gtXDkEIYAFHGpmhe3si3myt9RemSbO+d5wT4zoKVfLEzib6mf5pdOdnZ1NH7i4uMyaNeudd94Z/idXV9empqaHHnqooKAgIiKCiCIjI43Gn6zt5+fn197ebvpYrVbr9fozD+jv73dx+fEPwLa2Nn9/f1m/l5HA34kAFvDUYeNzyaKz/D2WE72ERaEiplLA5Vm4cOG2bdt6enq8vb29vb2dnJyMRmNLS4uzs3NQUBAR7dq1q6am5qz/lZqaWlVV1dPTQ0T+/v6BgYHDe0dotdqjR48mJCSYPs3PzzcNJVoWWoQA5vZDPavrpduizfRn6AtTxCnfGH4VIwbI2Q0Ldik1NfX555+fNm1aSkrK0NBQaWnpzp07J02alJKSMnXq1LCwsIGBgfHjx5/1v7y8vObMmbN9+/Zly5YR0Zo1a1auXCmKYk9PT0pKytNPPz1p0iQi6urqys/PX7RokQW+sZ9CEAKY2zNHjP9vqniFq6mNXIS7cPM48S+Fxlen4Z0ZuLRNmzap1erhTx988ME777zz5MmTarV6/Pjxpl7TTZs2mebCx8fH9/T0uLq6EtGLL74oCD/26z766KNvvvmmKQjnz59fWVn58ssv5+bmfv3118P9omvXrl2xYoW3t7eZv8FzIQgBzGp7Pes10PVRZh2VeDJRTPra8FSiwtvZnKcFm+Th4XHWV9zd3adOnXrmV0RRTExMNH3s6ek5fNjwAYsWLdqxY4dOpwsICCAihUIREBDg5eU1nIJEVFVV9cwzz8jxLYwWghDArF7JNz6RKJr5dYhQN2FphPhuibTqsvZ4ArgMw8utmSQlJZ2ZlET0l7/8xbxXdEF4KgDM53grK+ukn4+1wHP3RKL49gljv8H8ZwYgIpo2bdqtt95q6as4PwQhgPn8v+PS7+JF7guqjcREL2F6gLimHK+PwiXMnz9/8+bNclRmjC1atKikpOQixzz00EPr16+X4+wXgSAEMJNTXWxPo7TynFU8zOaJRPH1Qgn7M8HFPfLII3FxcXJU3rBhg5OTk+mV0Qt54IEHnnrqKdNCbmaDMUIAM3m7WLo3RnS33MT2mQFCgJq+Py1dF4G/gOGC+vr6DAYDEWVnZ+v1+t7e3h07diQmJt51111Hjx5dt25dWFjY/fffb3p9NDc3d8eOHe3t7YmJibfddtvwsmpff/11bm5ubGxscnJybW2taT3S995776677jId8Oqrrz744IM6nS4nJ+f222/ftm2bi4tLRkZGbGysq6vr7t27s7KyzPYt43kAMIcePa0tl34ZY+En7sFY8Z1i9I7Cxbz55psnTpwgoh9++GHlypWbN29OSEj405/+dO+997788ssJCQnffvvtY489Zjp43bp1Go0mNTX1iy++GA6511577ZlnnklMTKypqbn++us/+eQTIhoYGNi1a9fs2bNNxzz99NM9PT3l5eWmZWvWr1+/fft20z/Nnj1bpr7ZC0GLEMAcPjolzQkSI9wtvHrWjWPFxw4aT3awmDGWX8cLzqv1X380djSb7XRuMxa7pS+50L+OGzfu3XffJaKhoaEXXnihpqbG1Ld58803v/XWW0T09ttvE1F/f39WVlZUVNQHH3zg7Oz80ksv7dixIyUlhYiqq6sHBweJqLy83NnZ+eI7FJqMHz/+22+/5fUNjgSCEMAc/lEi/XWG5eezO4l0T4z4jxLpzZmWvxg4L69r75YG+812OqX3xZb6HB4sDAgIGD9+vJOTk+nj1tZW09dfeuml1atXe3l5KZVKo9HY0NDg5ubW1dU1PMswJSVl//79RNTb23vmJMKLcHNzMy3PZjYIQgDZ7Wlkeomygq2iEXb/JDH+K8OLUxVYhts6KQNCLX0J/6NQKM77scnRo0fff//9wsJCT09Pg8Hg6uoqSZKXlxcRdXZ2+vj4ENHw6tsBAQHt7e2SJImiSEQuLi79/f/L+76+vrCwMNPHLS0tI2k4coQxQgDZvVsiPRhr7kn0FxLsKmQFi5+cwkghXKnW1lY3NzfTNPk1a9bo9XoiUqvV8+bNM02Wb2pqWrt2rengqKgoLy+v0tJS06cpKSnffPON6eOOjo7hrlQiys/Pnz59ujm/EbQIAeTVOkg/1EnvZ1hR++uXMeIfjhh/NQl/B8MVmT17to+PT2Jioo+PT0hIyPCqoatXr16xYkVQUFBgYOCiRYtMLT9BEG644YaNGzeapk+8++67K1asaGlp6ejoiI2NXbly5fz584lIkqQdO3Y8++yz5vxGEIQA8lpbLl0TLo5xsvR1nGF+sPDLfspvY4k+VtJMBSuyfft2lUpFRE8//fTwF6+99toFCxaYPg4PDzdtveTs7Lx3797y8nK1Wh0REdHR0WFadzQiIiInJ8fUC3rXXXdNmDDB9B9/85vfXH/99Y888ohCoZg8efLx48c/++yzl1566eDBg8PLfK9fv3769OlRUVHm/JYRhADyWlMm/dXK3kwRBbpzvLCmzCre3wFrM5xJZ77bolKpTOlIRKIojhkzZvjjiRMnmj4e/uLWrVu3bdsWHR197Nix7du3Dy86OmHChAceeKC8vDwmJsb0FW9vbxcXlzM3u6itrX3llVfk+t4uAH0jADI60sK69DRbY3UNr19MED8+JQ0aL30kwGhNnTp18uTJbW1tM2bMKCwsNG1AYXL//fcPpyARRUVF3X777Wf+31//+tfjxo0z37USEVqEALL6T5n0iwmitbwnc4ZIDyHeR9hwWrrBvBtCgSPw9fX9xS9+MZIjx48ff+6+vuaHZwBALgNG+qxCunO89cUgERGtnCj+uwzvjgIgCAFks+G0lOwnhFt6NZkLWR4hHtAxrfmmbgNYKQQhgFw+OcVWjLPeR8xFSdeGi59VoFEIZ6uoqHj00Ud//etfb9q0ydLXYg7W+5QC2LT2QdrVKC2PtOpHbEW0+AmCEH6qoqJi1qxZkZGRaWlpDzzwwEcffWTpK5IdXpYBkMWXVdLCENHLmqYPnmtesHBnDyvrZBO8rLT/1gE9tfvF1v52s51uYVTWDTHXnvmVt95667bbbnv44YeJyNXVddWqVbfddpvZrsciEIQAsvi4QvptnFU3B4lIIdBNY8V1Fey5FAShtfj11Hu7h8y35HSoR/BZXykuLl65cqXp4ylTppSVlQ0vEGqvEIQA/DX0saI2dlWoDfzuWBEtrthlfC7FBi7VQQS7ayx7Af39/d3d3aaPu7q63N3d7TsFCWOEAHL4pIItjxSdbWHZlmn+giDQ4WZm6QsBK/Lhhx8ODQ0R0erVqxcuXGjpy5EdWoQA/H1WIf15mi3EIBER/Xys8EWVlOpvMxcMchs3blxycrJKpTIajRs3brT05cgOLUIAzqq7WU0Pm2V9y6pdyA1R4hdVDE1CGLZ8+fIDBw6sX7++sLAwPDzc0pcjOwQhAGdfVLHlkaLSdp6tRB/BWaS8FkQh/I+Hh4cjRKCJ7TysADbiyyrbW8DzZ5HCl1WYUAhERL/61a9iY2MtfRVmZWOPK4CVq+1lld0sK8hm+kVNbogSP69EixCIiFasWBEdHW3pqzArBCEAT19UsmURttQvapLiJwgCHW9FFoIjsrXnFcC6fVUtXW9r/aIm10cKX6B3FBySTT6xANapsY9OdrB5wTbWL2pyfZT4bTVahOCIEIQA3Kw/LV0dJqps86lK9Re69FTWiSwEh4MJ9QDcfFcj3T3BNmOQSCC6Nlz4roY9nmCTLVpbpFQqu7q6xo0bZ+kLsQHNzc2zZ8+WqTiCEICPHj3ta2KfzbXVICSi6yLEF48ZH0+w4W/Btnh6elZXVw8ODlr6QmyDfPMaEYQAfGyqldIDBQ+Vpa/jCswNFm7dxbT9pHGx9KU4jJCQEEtfAmCMEICT72rYdRG2/UCpRFoYKm48jXdHwbHY9nMLYCX0Em2tk64Nt/kH6rpw4bsavC8DjsXmn1sAa7BXy8Z7CUGulr6OK3ZVmLhXK/UaLH0dAGaEIATgYMNpe2gOEpGXE031E3bUo3cUHIg9PLoAFrepli0Jt5NZB0vCxY216B0FB4IgBLhSpZ2sV08JPnYShEvDhQ2nJSQhOI5RTJ+orKz84YcfvL29ly5d6uJy/ter8/LyDhw44OHhkZWVFRoayukiAazaxtPsmnDBTmKQaJyn4K4S8ltZkq/dfE8AFzPSFmFOTs6UKVNOnDjxr3/9KzMz87wzQB977LGlS5cWFBRs27btn//8J9frBLBeG2ulJWF2lRlLwgT0joLjGGmL8IUXXli1atVjjz1mNBqnTp365Zdfrlix4swDtm7d+tFHHxUVFfn5+clwnQBWqktPh5vZ3GC7GmVYEiY+e9S4KsmuvimACxnRjT44OLhjx45ly5YRkUKhuPbaazdv3nzWMV988cUdd9zR2tq6fv366upq7hcKYJ221UnpgYK7LS8oc65ZQUJxB9P1W/o6AMxiRC1CrVbLGAsODjZ9GhwcnJOTc9YxlZWVRUVFe/fujYmJufvuu99444077rjjvNXa29tPnDjx0ksvDX/llltuuciAol6v1+v1I7lOkI/ppyDYz0AYNxtP08JgMs8tavoRSJLscxsEoiwNbTqtXzFW7lPZHr1eL4poK1uYXq83Go0KheKSRyqVykv+4hpREDLGiGi4liiK5z6KQ0NDQ0NDR44cEUVx06ZNt95664oVK857lQaDYWhoqL29ffgrAwMDF3m2JUkyw5MPF2f6KSAIz8KIttaLj0820x1qzp/C4hBha51wSyQevbOZftj4pWRZ0v/hUm1EQajRaARBaGpqioyMJKLGxsbh1uGw4OBgX19f0x9K6enpnZ2dWq32vOvJ+vv7Jycnv/baayO8xKGhIWdn5xEeDDLR6/XOzs4IwrMUtDE3lXGSn5nuT0mSnJycRvJX8JVbEsmeOWZQOTmL+JmfQxRFlcq+esNtjSiKRqORVzSMqIGvVqvT09M3bdpERIyxLVu2zJs3j4gMBsPp06dNmbxw4cKSkhLT8cXFxS4uLgEBAVwuEcBqbalji0PtMyhC3YQAFyGvFe+Ogv0b6Vujf/jDH2655RadTldSUtLR0XHzzTcTUVVV1YQJE5qbm/38/G699dY33njj9ttvj4+Pf++9955//nn8xQR2b2ud9Gi8OdpnFrE4VNhSy6b62WfSAwwb6ZDvokWLdu3apVKpsrKy9u/f7+bmRkQajebjjz/28PAgIldX14MHD2ZmZkqStHbt2t///vcyXjWAFeg10JFmNltjtzmxKFTcikVHwQGMYmWZxMTExMTEM7/i4eFx6623nvnpL3/5S26XBmDddjWwVH97mzhxplkaoaiNdQzRGCdLXwqAnPASMMBl2lonLQq15yfIWUEzA4WdDWgUgp2z58cYQFZ2/KbMsEUh4tY6vC8Ddg5BCHA5KrpYn4Hi7GXHiQtZHCZsQRCCvUMQAlyOH+rZghD7n1Y50UsQBSrtRBaCPUMQAlyO7fVsQYjd5yAR0dwgYVs9ghDsGYIQYNSMjHY3Sna248SFzA8RtiMIwa45xJMMwNeRZhbiJgS5Wvo6zGJBiLi7UdLj1VGwXwhCgFHb5jD9okTkp6YoD+FwMxqFYLcQhACjtr1BWhDiQM/OghAME4I9c6CHGYCLXgPltbBM+11Z7VzzQ8TtmFYP9gtBCDA6exrZFD/BbRSrE9q8WRohv5V1YXtssFMIQoDR2V4vzXekflEiUitoeoCwpxGNQrBPjvRnLQAPOxvY+xlmDcIh49CO6r3FrWWlLadcnVwm+kQnBcbPDJlqzmuYGyzuamDXhpvznABmgiAEGIXWQaruYVPMuEVfdu3+d/L+HeUVMTUoaVbQDKMoneqo/iB/7afFX/966r3R3lHmuYy5QcIvc9AiBPuEIAQYhR310iyNqDRLg1Bi0p8P/L20reL303+dokkgov7+ficnpxkhU2+N/dn3p354bOdzKxNXXBu9yAwXM9VfqO1lTf0U6GKGswGYFYIQYBR2NbK5weZoDkpMenn/3zoGu95f/Lqz4uz9AEVBXDp+8bTglEe2/2HAMHhjzFK5r0chUEaguKdRummsY42PgiPAPQ0wCjsazBGEEmPPZf+5V9/30qxV56bgMI1bwN/m/79vyzZ9eXKD3JdERHODhZ0NmE0IdghBCDBS9b2sfZDFecsehJ+WfN0+0PnHzKdUCtXFjwx08//r/D99XPxlYXOx3Fc1N1jY2YggBDuEIAQYqe0NbG6wKMqcg8UtZZ+XfPdM+qNKUTGS4wNc/Z6c8fAfc/7SOdgl64XF+whdQ6ymB1kI9gZBCDBSO+XvF+0a6n4+59UnZjwc6BYw8v81PXjK3IjMl/f/Tb4LIyKBaE6QuBuNQrA7CEKAkdrVwOYGyRuEHxz/KD009TLmCN6bdEdLX9vOmmw5rmrY3GBhBxYdBbuDIAQYkVNdjBGN95IxCE+1V+6t3feL+Fsv4/8qRcVvU3/1j7z/DBgGuF/YsLnBAlqEYH8QhAAjsquBZcncHHzryAf3JN7u6exxef89zj8mMSDuk+Kv+F7VmaI9BYmoogtZCHYFQQgwInu0bLacQbijem+/YeDqcQuupMivku/8tmxzY08Tr6s615wgNArB3iAIAUZkdyPLku1NGSMzfpD/0UNTVorCFZ3Cz9X3+onXrClcx+vCzoUgBPuDIAS4tPJOJhKN9ZArCLdX7dG4ByQETL7yUjfELN1ff6S+u/HKS53XnCBMqwd7gyAEuDRZm4MSY+uKv74j7udcqrmpXK8bf9WnJd9wqXauaE9BKdIpDBOCHUEQAlza7kYZBwh3n85xUbkkB8bzKnjjpKW7a3KbenW8Cp5ltga9o2BXEIQAl7ZHK9cro4zYR0Vf3BV/M8eank4eS6IXfFryLceaZ5qNYUKwLwhCgEso62QiUZQ8A4RHGo8Lgjg9eArfsjfFXLe9ak/PUC/fsiZZwcLOBuxNCPYDQQhwCXvkHCD8pmzjzyYu4V7Wx8V7evCUzZU7uFcmorEegkoUyjrRKAQ7gSAEuIQ9WjZHnn7Rpt7mouaT8yIy5Si+fOLV35RtlJgscTVbI+zVIgjBTiAIAS4hW8tmaWQJwu/KNy8aO1etVMtRfLJfjLvK7aj2uBzFZwUJezFMCPYCQQhwMZXdzMhonCf/INQb9ZsrdyyNXsS98rDrJlz1TdlGOSrjxVGwJwhCgIvZ08hmy9Mc3HU6N3pMVJhniBzFTeZFzCpqPinHPIrxXoKRUXU3shDsAYIQ4GL2yrbE6ObK7deOl7E5SERqpfO8yFlbKnfJUXxWEIYJwU4gCAEuZo88U+mbepsr2qtnBo9638HRWhw1d2vVTkb8E2uWRtiD3lGwCwhCgAuq62W9BjZBhj0If6jaNTciU6VQca98lom+0U4KpxPNpdwrz0aLEOwFghDggnY3stkaUY6O0R+qdi8amyVD4fNYGDVniwwTCieNEbr0rL4XWQg2D0EIcEEy9YueaDnJiE3yncC98nktjMrac3rfoHGIb1mBKCNQRKMQ7ACCEOCCZJpBuLVy1+Kx87iXvRA/F5+JvtH76g5xrzw7CMOEYA8QhADn19RPzQNssjfnIDRIxt2ncxdGzeFb9uIWjc36oWo397KZGiGnCUEINg9BCHB+e7VShkbkPkJ4VHs83DM0wNWPc92Lygidka8r6h7q4Vs20Udo6GO6fr5VAcwNQQhwftlalilDv+iu07lZEency16ci1KdoknYX3+Yb1lRoJkBwj4ddqIA24YgBDg/OQYIDZJxX92h2WFpfMuOxJzw9F01udzLZmrEbLwvAzYOQQhwHp1DVNnFkn05B6GpX9TP1Zdv2ZFID51e0HyiV9/Ht2ymBqtvg81DEAKcR04Tmx4gqHg/HxbpFzVxUaqTA+Nz6w7yLZvqL5R1sm4936oAZoUgBDiPbK2UqeH8dFiwX9REjt5RJ5FS/IR9eHcUbBmCEOA89jbyf1Mmryk/3DPEIv2iJmkh0/J1Rdx7R2dphGwt3pcBG4YgBDhbv4EK29k0f85BmFN7MDNsJt+ao+KqckkIiD3UkMe3LN6XAVuHIAQ424FmluAjuCp51mTE9tUfTgtJ5Vl09NJDp++r57zEzMxAIa+VDRj5VgUwHwQhwNnk6Bcta61wU7nIug3vSKSFTNtff8Qg8UwtNyVNGiMcbkajEGwVghDgbDlNUkYg50cjt/5geuh0vjUvg6+Ld6hHcGFzMd+yGYECekfBdiEIAX7CINEhHUsL5D9AmB46jW/Ny5MeOo37JIpMjZDThPdlwFYhCAF+4lgri/AQfJx51mzq1bUNtE/ynciz6OVKD52eXXuAb81MjbiviRnRJgTbhCAE+ImcJpbBvTlYdzAtZJooyLHF76iNHRMhCmJVRw3Hmn5qCnIVCtuQhGCTEIQAPyHHWtv76g+nWUe/qElaaGou7wW4MUwItgtBCPA/jCi3ScrgGoQDhoGSlrKUwASONa/Q9OApBxuO8q2JvQnBdiEIAf6nrJO5KIQwN55BmNdUEOM73lXlwrHmFUoKiDvVXtk11M2xZqZG2NuI92XAJiEIAf5Hjn7RA/VHpwdP4VvzCjkpnBL8J+dpCzjWjPIQlKJQ0YVGIdgeBCHA/+Q2Mb79okR0uPGYtQUhydM7mqHBMCHYJAQhwP9kazm/Mnq6q84gGSK9wjjW5GJmyNSDDUcZ8cyt9EAhF8OEYIMQhAA/0vZT+yCbNIZnEB5oODojZCrHgrwEuQe6KF0q2qs51sT7MmCjEIQAP8rWShkaUeTaM3qoIW9aUArPivxw7x2N9xa0fUzXz7EkgDkgCAF+lKNl6Vz7RQcMg8UtpVM0iRxrcjQjZMoBrkEoCjQzUMjFWmtgaxCEAD/ivqZMvq5ovPdYq5o4cabEgLhT7ZV89+lNDxQxTAg2B0EIQETUo6eyTjbFj2cQHmk8PjUoiWNBvpwVTpN8J+TrijjWxDAh2CIEIQAR0T4dm+InOCt41jysteogJKKpQUmHG49zLDjdXzjRzvoMHEsCyA5BCEBElKOV+PaLtvW3t/S1TvSJ5liTu6mapCNcg9BZQQk+wkFs0gs2BUEIQPTjm1iz8i8AACAASURBVDI8H4cj2uMpmgRRsOpHLNp7bNdgt66vhWNNrL4NNseqn1IA89BLdKSF82a8R7T5UzVW3S9KRKIgpGgSjmrzOdbM0Ag5Wrw4CrYEQQhAR1vYOE/By4lnzTxtgZUPEJpM0STy7R3NCBQP6pgBUQi2A0EIQLm8J05UdZ5WCGKwu4ZjTZmkBiUd1R7nuNaatzOFuQsF2KQXbAeCEIBymzhPpT/amJ8alMyxoHwC3QLcndwr2qs41swIxCQKsCUIQnB0ps1407luOnFUm2+1C8qca4omke+WTOkarL4NtgRBCI6urJO5KnluxisxqbC5OCkwjldBuSUHxh9rKuRYMD0Qm/SCLUEQgqPL4b31UmnbqQBXP2/1GI41ZZUSmJCvO2FkRl4Fx3oISlGo7EajEGwDghAcHfcBwjxtQYomgWNBuXk6e2jcAkpbKzjWTAsUcjCbEGwEghAcXQ7vXemPNRUmB9pSEBJRiibhWBPXYUJs0gu2A0EIDq2pn1oGWCy/zXj1kqG4pTQhIJZXQfNIDkzgO0yYgRYh2I5RBOGnn356zz33rFq1SqvVXuSwzz///K233rriCwMwhxytlB4ocNyMt6SlNMwzxMPJnVtFs0gKjDvRclJv1PMqmOgrNPSx1kFe9QBkNNIg/Pvf/75q1arZs2e3trZmZGQMDp7/Bj948OBDDz304osv8rtCABnlNrE0rkuM5jUVJAfGcyxoHm4q13DP0JLWMl4FFQJN8xf2YZNesAUj+hVgNBpff/31d9555/bbb3/vvfc8PDy++uqrcw8bHBx84IEHXnjhBd4XCSAX7pvxHtMWptjaAKFJcmB8HtdhwjRs0gs2YkRBWFdXV1tbm5WVZfp0zpw5ubm55x724osvLlu2LDbWxkZHwGH1Gqikg6X6cwvCIeNQadupeFsbIDRJ0SQca+K9SS+GCcEWKEdyUGNjo7u7u7Ozs+lTf3//I0eOnHVMfn7++vXrDx8+fODAgYtXq6ury87Ovv7664e/8vvf/z4+/oK9Sf39/QoF1/1SYfT6+voEQRAEno0ni9vTJCZ4i9JgXx+nggUtxZGe4WxI6hviVfIn+vv7DQaDTI/DOLfI0tbyju4OJwWf1ccT3Cm/TdXW3ae2r8d3cHBQFEWVSmXpC3Foer3eaDRK0qX73tVqtSheosk3oiBUq9VDQ0PDnw4NDbm4uJx5gMFguOeee957773hsLwIHx+fiIiIn//858NfiY6OVqvVFzper9df5F/BPAwGg1qttrMgPNTGMjWM491V3F6WoomX73ZljDk5OckUhGpSR40Jr+qtTQyYzKkgxY6RCrudM7nOTrE4QRAQhBanUCiMRuNInrVLpiCNMAhDQkIGBwebm5v9/f2JqLa2NiQk5MwDKisr8/Pzb7/9diIaGBhoa2sbN27c9u3bo6Kizq3m6uoaHh5+0003jeTURCSK4ki+E5CV6adgZ0GYqzP8ZrJC5PfOaH7ziRWTr5fvdhX/j0z1EwPiCppPJGu4veyToWH7moXZwXb1/Mr9U4CREEWRMcbrpzCiKv7+/pmZmR9//DERdXR0bNq0afny5UTU1tb2zTffEFFUVNTJkye3bdu2bdu2119/3cvLa9u2baGhoVwuEUAOBokO6Xhuxqs36ktby+P8JvEqaH5JgXHHdTyHCdMDsUkv2IARtQiJ6JVXXlm2bNmuXbuKi4sXLlw4c+ZMIjp58uTPfvYzxphKpRo7dqzpyNraWoVCMfwpgHU63sbC3QWfS/flj1Rxa1nUmHBXlculD7VWCf6TX8h5TW/UqxR8+v0yNeLde41GRgq76koAezPSIExLSzt58uSBAweCgoKSk3/caC0lJaW0tPSsI6dPn3748GGe1wgggxwt55XVjjcVJgbYzI4T5+WqcgnzDDnZVh7vz+fFV381BboIJ9pZgg+SEKzXKDpYfXx8rr766uEUJCK1Wj1hwoSzDlOr1REREXyuDkA23NfaPq4rsqGtly4kOSD+ONdJFFhrDawfxnvBQeU2SRyn0uslw8lWbg0pC0rkPUyYocFu9WDtEITgiE51MVEQIj24BeHJ1rIwzxA3lSuvgpaSGDC5pKXMIHHbmzAjUMhGixCsG4IQHFG2ls3iOkCYrytOsvEBQhM3lWuIR1BpWzmvguO9BCNjNT3IQrBeCEJwRLm8lxjN1xXZ+psywxIDJhfoijkWTA8UMUwI1gxBCI4om+sroxKTiltK4/1teAbhmRICJhfoTnAsmBGIYUKwaghCcDjNA9TUzyZ7cwvC8vbKAFc/T2cPXgUtKyEgtqC5WGLcJsJnaDBMCFYNQQgOJ1srpQcKHKd45+tOJHBan9MajHH28nXxqeyo4VUwyVeo7cEmvWC9EITgcHK0LEPD884v0BXbUxASUWLA5Hx+vaMKgaYHYJNesF4IQnA4OU0sk9+bMoxYUXMxrx0brERCQCznYUIN3pcB64UgBMdi2ox3Kr/NeGs661yULn4uPrwKWoOkgLh83QlG3KIrE8OEYMUQhOBYDuhYkq/AcavYAt2JRNtfWe0s/q5+zgqnuq4GXgWn+wuF7azfwKseAE8IQnAs2VqJ71T6Al1xou2vrHauxMA4jsOErkqK9xYONaNRCNYIQQiOJVvLMrm+KZOvK4oPsMMgTPCfxDEICb2jYMUQhOBA9BIdbmYzA7i1CLW9Or1kCPUI5lXQesT7xxY281xfJlMjZmOTXrBKCEJwIHktLNpT8HLiVrBAZ2/viw4L9wodMAzo+lp4FUwPFA7omAFRCNYHQQgOZK+WZXIeIDyRYI/9okQkkBDnP6mwuYRXQW9nivAQjrehdxSsDoIQHEg27yAsbC62gz0ILyTeP7aQ6+rbGCYE64QgBEfBiPY1SemB3O75zsEuXV/LOO9IXgWtTUIA72FC7E0IVglBCI7iRDvzVQtB/LbOLWwuifObpBD4zUm0MhN8xjX2NPUM9fIqODtIzNFKSEKwNghCcBT8+0V1xfEBdrL10nkpBMVE3+iiFm7DhEGu5KESSjsQhWBdEITgKLK1nDfjLWi2q00nzivBfzLfYcJZQcJe9I6ClUEQgqPI0bLZQdyCcMAwWNlxOsZnPK+C1ikhIDafaxBmYJgQrA+CEBxCRRdjRFEe3IKwpLVs3JhItdKZV0HrFOs38VR7pd6o51VwlkbY3YggBOuCIASHsJdrc5DsegbhmVyU6nCv0JNt5bwKjvcSJEbV3chCsCIIQnAI3KfSFzaXOEIQElGCf2wBhgnBriEIwSHsbWQcN52QmFTSWhbnZ8+vjA6L948t4re+DGE2IVgfBCHYv/pe1q1nMWO4BeGp9ip/F19PZw9eBa1ZQkBsYXOJxLhFF1qEYG0QhGD/9mjZrCCRY8doYXOxXW69dF7e6jFezp41nad5FZzsLbQNsoY+ZCFYCwQh2L9sLc9+USIq0BXH+ztEv6hJfEBsAb+11gSijEAxB41CsBoIQrB/e7gOEBJRUXOJHa+1fa54/0mFOq7DhBr0joIVQRCCnWseIG0/S/DhFoQNPVoShCD3QF4FrV+CP88WIRHNDhL2YDYhWA0EIdi5PY1SRqDIcYSwQFec6EjNQSIK8wzRS3qOm/Qm+Qr1faxlgFc9gCuCIAQ7t6eR81T6wubiOAcLQiKK84vhuOioQqCZAcJeLbarB6uAIAQ7t5t3EBboihPsetOJ84r3n8S7d1RE7yhYCQQh2LPWQarrZcm+3IKwc7Crtb9t7JhIXgVtRXwA593q52CYEKwGghDs2Z5GKT1QUHAdIJzsHyMKDvfgTPAZp+3VdQ118yqY4itU97DWQV71AC6fwz3P4FD2NLLZQTxv8qLmkgTHGyCk/9ukt7illFdBpUgzA4RsDBOCFUAQgj3jP0DoSGvKnCXBf3Ih10VHMUwIVgJBCHarbZCqu1kKvwHCQeNQZUeN3W/GeyHx/pP4bkOBYUKwEghCsFt7tVK6RlDyu8eLW0odYTPeC5nsH1PWVjFkHOJVcKqfUNnN2jBMCJaGIAS7tbuRzdbwvMMLdMUOsgfhebko1RGeoaVtFbwKKkWagWFCsAIIQrBbuxpYVjDnJUbjHGmt7XPFB0wq5DqbcE6QuBu9o2BpCEKwT62DVNPDc4BQYtKJlpMOtenEueL9Oc8mzAoSdjYgCMHCEIRgn3Y1SBmBPAcIKzqq/V19vZw9uVW0QQn+nDfpneInnO7BoqNgYQhCsE+7GllWMN8BwhMOtfXSefm4eHs6e3DcpFcpUnqgsKcRw4RgSQhCsE+7GlgW/yVGHT0IiSghYDLfRUezgsVdGCYEi0IQgh3S9ZO2nyXyGyAk05oyAZM5FrRRCf6x+boTHAtmBQm7MEwIFoUgBDu0s0GarRE5LjFa391IgqBxC+BW0WYlBHAOwiRfoamfafs5lgQYHQQh2KFdjZwnThQ0FyeiOUhERKEewYxJTb06XgVFgTI04u4GDBOCxSAIwQ7tbmRzOA8QnnDMtbbPK94/Np/rJIq5wQKGCcGCEIRgb+p7Wdsgi/Pmuys9Bgj/J94/lu+0+rnBwg4ME4LlIAjB3myrZ/OCRZFfDnYMdnYMdEZ6hXOraOMSAmILuA4TTvYWevWsuhtZCJaBIAR7s6OBzeM7QKgrjvOPEQWeNW1atHdUc19r52AXr4ICUVawiEYhWAqCEOzNrkY2L4T3ACH6Rc8gCuJk/xi+exPOQ+8oWA6CEOxKcQdTCjTWg2cQ5utO4JXRsyT6T+bbOzo/RNjRICEJwSIQhGBXdtSzBVybg736vrruhgk+0Rxr2oHEwMl8N+mNcBc8VMKJdkQhWACCEOwK9wHCouaSGN/xKlHJsaYdiPEZX915uk/Pcxr8vGBhRz2CECwAQQj2w8goWyvN4b3WdoI/+kXPplKoJvpEF7eUcqyJYUKwFAQh2I8jzSzETdC48KyZr8OaMueXEDC5oJnnMOHcYHGvVtJjhRkwOwQh2I/tDWw+137RIeNQRUdVrN9EjjXtRkJAbH4TzyD0U9NYD+GgDo1CMDcEIdiPH+qkhaE8b+niltIorwi10pljTbsR5z+ptO3UkHGIY80FIcJ2LDoKZocgBDvRrae8Vpap4TtxAv2iF+SiVEd4hZW2VXCsuSBE/KEOLUIwNwQh2IndjdJ0f8GN69udBc0n4rEZ74UlBkzmuyVThkYoamftgxxLAlwaghDsxLZ6tiCE5/1skIwlLWXYdOIiEgIm5zcVcSyoVlBaoLC7Eb2jYFYIQrATP9Rxnkpf2lYe4hHk7uTGsaadSQyYfKLlpEEycqy5IETchtmEYF4IQrAHdb2sdZAl+fIMwuNNRUkBcRwL2h8PJ/cg98Dydr7DhMIPCEIwLwQh2IMtdWxBCM+tl4jouK4oMRBBeAmJAXHHufaOxvsIfQZWiS2ZwIwQhGAPtvNeYtTIjCeaT2KA8JKSAjkHoUA0PxjvjoJZIQjB5hkZba+XFoXyDMKytoog90BPZw+ONe1SUmBcUUuJxHi+3rIoFL2jYFYIQrB5h5pZqJsQ7Mp5gDARA4Qj4Onk4e/qV95eybHmolBxV4M0hFdHwVwQhGDzttRKV4Vx3j4+X1eUhAHCkUniPUzop6YJXsK+JjQKwUwQhGDzNtexxVxXVpOYVNR8Mh4DhCPDfZiQiBaHClvq0CQEM0EQgm1rGaDyTpYWyLNFWN5e6efi46324ljTjiUGxBU2F/MdJlwcJm6pRYsQzARBCLZta52UFSyquN7IedqCZE08z4p2zVvt5evizXeYcJq/UN/H6nuRhWAOCEKwbVvq2FVc3xclomNNhcmBCXxr2rcUTcKxpkKOBRUCzQ8Rt+LdUTALBCHYMInRNt4TJ4zMWNRcgk0nRiU5MOGYlmcQkmmYEL2jYBYIQrBhR1qYn1oId+e6xGjrqSB3jZezJ8eadi8pMK6wuZjvoqOLQ8XtDdiwHswBQQg27PvT0jXhnPtF87QFKYEYIBwdTyePIHdNWdspjjUDXSjaU8jBJAqQH4IQbNj3p9mSMM73cF4T3pS5HCmB8XnaAr41l4SJG0+jSQiyQxCCrWroY6d72MwAni1CvWQ42Vqe4I8BwlFL1sTnNXEOwmvChe9Po0UIshtFEObm5i5cuDApKemJJ54YHDx7D2mdTvfUU0/NmjUrNTX1oYceampq4nqdAGf7/jRbFCoquf4tV9xSGuYZgj0IL0NiQFxJa5neqOdYM8VP6DFQeSeyEOQ10t8iLS0tS5Ysufnmmz/55JP9+/c///zzZx1QVlY2MDDwwgsvrF69ur6+fvny5ZyvFOCnNtayJbwHCI83FSZjgPCyuKlcIzzDilvLONYUiK4OEzbi3VGQ2UiDcO3atampqXfffXdsbOyf//zn1atX6/U/+dMvIyPjr3/9a1ZWVlJS0uuvv75///7u7m4ZLhiAiGjQSHsapUVcV1YjoqPa/BQNZhBephRNwlFtPt+aS8KEjbUYJgR5jfT3SGFhYWpqqunjqVOntrW11dfXX+jgvLy84OBgDw9sYQNy2dnAEn0EX2eeNQcMA+XtldiD8LJN0STm8Q7C+SHiQR3rHOJbFeAnlCM8TqfTxcTEmD5WqVTu7u5NTU2RkZHnHtnQ0PDwww+/+eabFyp16tSp7777Lioqavgr77//flpa2oWO7+3tFQTOPWAwWr29vYwx6/lBfFupnB9IPT0DHGsebjoW7RVlGDD0UA/Hshz19/c7OTkpFApLX8j5RbmEnWqv0nU0uypdOJad4afaUNG3LMxa2oWDg4OiKKpUKktfiEPT6/VGo9FgMFzySFdXV1G8RJNvpEHo5eXV29tr+liSpL6+vjFjxpx7mE6nmz9//sMPP3zjjTdeqNS4cePmzZv317/+dfgrERERF3m2GWPu7u4jvE6Qj5ubm5UEISPa3GDYdrXC3V3NsWxxadm0kBRrvtkUCoU1ByERxfpNPNVblRYyjWPNn42VtjYpb5tkLd+1SqVCEFqcKQjVaj6/AUbaNRoVFVVW9uMweEVFhUKhCA0NPeuY1tbWBQsW3HTTTU8++eRFSgmC4O7uPvYM1vxggxU62sLcVTTRi3MqH9Een6JJ4lvT0UzRJB5p5Nw7ujRC2FiLJWZARiMNwhUrVmzatKmiooKI3nzzzWXLlrm5uRHR2rVrt2/fTkRtbW0LFixIT0//zW9+097e3t7ebjTyXG8JYNh3NdJ1EZxTsH2go7mvNcY3mm9ZRzNVk3RUe5xvzWBXYYKXsFeLd0dBLiMNwkmTJj377LNTpkwJDQ09cODAX/7yF9PXN2zYsH//fiLau3dvdXX1p59+Ou7/VFdXy3TR4OC+rWbXRXB+X/SI9nhSYLwoYImJKzLeZ1z7QGdzXwvfstdFiN/VoEkIchnpGCER/e53v3vwwQe7u7v9/f2Hv/j555+bPli2bNmyZcs4Xx3AOSq6WMsAm+bPuUV4VFswBRMnrpgoCMmB8XnagkVj53IsuyxCWLxFenMmWcUYNdid0f39q1arz0xBAPP7toZdFyGKvH8j5mkLMEDIRYom4SjvtdYmjRFcFHS8Fb2jIAt0BIGN+a5G4t4verqrnoiFe4bwLeuYUoOSjzQeY8Q5tJZGCOgdBZkgCMGWNPVTUTubG8y5PXioIW9aUArfmg4r2F3jonSp7KjhW3ZZhPhVFVqEIAsEIdiSr6ulq8NEZ97TbQ415k0LRhByMy04+VBDHt+aMwOFziEq6UAWAn8IQrAlX1ZJ10dybg4OGYeKmkumaBL5lnVkqUEphxuP8a0pEC2LFNAoBDkgCMFmtAxQXgtbzHuh7XzdiXHeUW4qV75lHVlKYHxJa1m/gecCeER0Q5T4VTWGCYE/BCHYjK+rpUWhossopvyMyKHGY9OCkjkXdWxqpXqS74RjTYV8y2YECrp+OtWFRiFwhiAEm/FVlXRDFP+JZIfxpowMUoOSDzdyHiYUBVoWKXyJ3lHgDUEItqF9kA42s6vCON+xzX0t7QOd433G8S0L04JTDjVwHiYkousjxa+q0DsKnCEIwTZ8UyPNDxHdePeLHmzImxqUJFrHrhr2ZOyYiAHDQEOPlm/Z2UHC6V5W1Y1GIfCEIATbsK5CumUs/7g60HBkRsgU7mVBIGF68JT99Yf5llUIdEOU+GklghB4QhCCDdD109EWdjXvflG9UX+sqXB6EIJQFjNDpu6vP8K97C1jxY/K0TsKPCEIwQZ8WildG87/fdFjTYVRXhGezh6c6wIREaUGJRe3lPbp+/mWTdcIfUYqakejELhBEIINWFch3TKO/726v+HwzJCp3MuCiVqpjvWbyH17QoHopihhXQUahcANghCsXU0Pq+xm83mvL0pEB+qPpoWkci8Lw+TqHR0nflrB0CQEXhCEYO0+qWA3RIlK3rdqVedpIzNGjYngXBfOkBYybV/9IYl3ZiX5CmoFHdIhCoEPBCFYu3UV0i1j+d+oB+qPpIVM414WzhTkHujp7FneVsG98s3jxI/ROwqcIAjBquW1sG49pWvk6BfFxAlzSAtJ3cd7EgUR3R4trKuQhhCFwAOCEKzaf8ulu8Zz346eOge7TnVUpQQm8C4MZ0sLSc2pO8i9bKSHEDtG2FyLJAQOEIRgvQwSfV4p3RbNvzmYW3coNSjZSeHEvTKcJc4/tm2gnfsSM0R05wTxv+UYJgQOEIRgvTbWShO8hHGe/IMwu/ZAZugM7mXhXKIgyNQovClK3N0oNXPe6wkcEYIQrNeH5ezOCfxv0QHDYL6uaHowBgjNJDNsRk7tAe5l3VV0dZj4KV6ZgSuGIAQr1TpIOxukG6L436IHG45O9otxd3LjXhnOa0pgYkVHdVt/O/fKd40X/4vl1uCKIQjBSq0tl64NFz1V/Ctn1x7IDEO/qPmoFKrUoOT9Dfxn1s8NFloG6HgrRgrhiiAIwUp9UCrdG8P//jRIxoONR9NDp3OvDBeRGTYzu3Y/97KiQHdPFFeXolEIVwRBCNYoR8uMjDJkmD54vKkwzCPE18Wbe2W4iBnBUwp0xb36Pu6VV04QPq2Q+gzcC4MDQRCCNVpdKv0yhv/0QSLaWZOdFZEhQ2G4GDeVa1JgXG7dIe6VQ9yE9EDx80o0CuHyIQjB6nQO0foa6bZoWfpFc+oOzgqbyb0yXNLciMxdNTlyVL43RkDvKFwJBCFYnbWnpKvCRH81/8pHtMcivMIC3fz5l4ZLSQ+dnq8r6hrq5l756jDxdA92KITLhyAE68KI3i2WfjVJljtzZ01OVjj6RS3DRameokmUo3dUIdA9E8V3i9EohMuEIATrsr2eKQSaJcNrMnqjfl/doVnh6Be1mLkRmTtrsuWofN8kcV2F1D4oR22wfwhCsC5vnZB+GyfLbXmwMS/aO8rPxUeO4jASM0NSi1tKOwY7uVfWuNDVYeIaTK6Hy4IgBCtS3c3266RbxsnUL4r3RS1MrXSeFpSSLcNya0T068niO8WShIFCGD0EIViRt4qluyeIrkr+lXv1fQcbjs4JT+dfGkZjQdTsH6p2y1F5RoDg40yb65CEMGoIQrAWvQb6sFy6X57XZPac3pcSmODl7ClHcRi56cFT6rob5NiViYgeihX/XmSUozLYNwQhWIvVJ6V5wWKkhxzT6Glr5c6FUVlyVIZRUQiKuREZP1TtkqP4LePEk510DEuPwighCMEqGCT6W5H0iDyvyTT16qo6T8/AvkvWYWFU1pbKnYz4x5VKpIdixdcL8coMjA6CEKzCZ5VStCdND5ClObilcte8yEyVQoadLGD0JvpEuyjVRc0n5Sj+q0ni1jqppgeNQhgFBCFYhdcLpccTFDIV3169e1HUXJmKw2VYEDVHpt5RDxWtnCj+rQiNQhgFBCFY3g/1zMhoYagszcEC3QlREGN8x8tRHC7Pwsg5u2tyBwwDchR/eLL4YbnUisn1MGIIQrC8P+YZn0iUZa8JIlpfvvXa6MXy1IbL5OfqGx8waac8a3AHuwo3RIl/LcTrozBSCEKwsG31TDdAPx8ry63YNdh9oOHIwrFz5CgOV2Lp+MXry7fIVHxVkviPEqlFlgYn2CEEIVjYi8eMz6eICnnag1sqd6SFTvN08pClOlyBaUFT2gc6ytoq5Cge7i5cHyW+eQKNQhgRBCFY0rZ6pu2XqzlIRBsrti1Fv6hVEgVhSfSC70/9IFP9PySJ75VIbRgphBFAEIIlvZAnY3PweFMhI4rzj5GlOlyxa8Yt3FmT3afvl6N4uLuwPFJ8HSOFMAIIQrCY72qkbj3dLFtz8JuyTcsnXC1TcbhyPi7eKZqErVU7Zar/XLL4folU14s5hXAJCEKwDCOjVUekV6cpZHpbVNurO9ZUuHjsPFmqAyc3xlz35ckNEpMlq0LchLsnii8ew5xCuAQEIVjGv0slfzUtkmfuIBF9eXL9knELXJRqmeoDF/H+kzyd3ffXH5ap/tNJim9rpOIONArhYhCEYAF9BvrjMem16XItJdOn799atWv5xCUy1QeOboi57vOT38lUfIwTPZ6gWHUYjUK4GAQhWMAr+cZZGmGqn1zNwe9PbZ0WlBLg6idTfeBoTnhaY4/2ZGu5TPUfihXz29jOBjQK4YIQhGBuld3sHyXSq9PkuvcMkvGr0u9vjFkqU33gSyEolk9Y8uXJDTLVVyvojRnir/cZ9WgWwgUgCMHcfrtfejxBEeImV3Pwh6pdIR5BWFzUhiwdv/hw47HarnqZ6i+LECM96O1iJCGcH4IQzGpTLSvrZL+VZ99BIpKY9EnxV3fG3yxTfZCDm8p1+cSrPy7+Sr5T/HWG4uXjRq0sUxbB5iEIwXx6DfTQPuNbaQon2e67bVW7fV18EgMmy3UCkMeNMdftqztU190gU/0JXsI9E8Xf7Mf8ejgPBCGYz9OHjXOChAUhcnWKSkz66MQXd8b9XKb6IB83leuyCVetK/5avlM8l6IoamNfV6ODFM6GIAQzOaBjX1Wzv8g2ZYKIQShQuAAAFolJREFUtlfv8XL2StEkyHcKkM8NMUv31u5v6NHKVN9ZQR/MUjy8X2rHAqTwUwhCMIcBI9291/jmDNHHWa5TDBmHPsj/+L7kO+Q6AcjM08njppjr/nn8Q/lOMTNAWB4hPHoQHaTwEwhCMIfHDxoTfYXro2S8374s3RDjGx3vHyvfKUBuN01aVtxSWtR8Ur5TvJyqyG1iX1ahgxT+B0EIsttcy9afZu+mydgp2j3U83nJd79MQnPQtjkrnH4Rf8u7ef9mJNf8d3cVrctSPLjPWNODKfbwIwQhyEvbTyuzDR/PUXjL1ilKRP8pWJcVkR7qESzjOcAsFo2dO2AYyKk9IN8ppvgJj8Qp7txjNCIKgYgQhCArg0Q37zT8apIiQyPXm6JEVN5eubMm+674W+Q7BZiNKIgPTbnnraP/GjAMyHeW3yeISoGeO4rBQiBCEIKsfnfQ6CzSqiQZbzOJsTcO/eO+pDu8nD3lOwuYU4omITEg9j8F6+Q7hSjQ5/OU6yowWAhECEKQzycV0uY69tk8pUwb0Jt8W7ZJKSoXj8O+g3bloSn3/FC9u7y9Ur5T+DjTV/MVD+4znmhHD6mjQxCCLHKb2CMHjF/PV4xxkvEszX0t/y389HfTHhBIzrAFs/Ny9lyZsOKNQ+9KTMYWW5Kv8OdUxbJtxmYZe2HBBiAIgb9TXeymHcYPZyvjvGXMJ4mxl/b/7YaYayO9wuQ7C1jKkugFbiq3tUVfyHqWuyaIt44Trt5i6DXIeh6waghC4KyxjxZuNr6UKsq3+7zJ5ye/1Rv1KybfIOtZwFIEElalPbK+fHNxS6msJ3p+iiLOR7hlp9GA4UJHhSAEnnT9NG+T4b4Y8c7x8t5ap9or1534+g/pvxMF3MN2y1s95rep9/1p3xv9cr5BKhD9M0PBiN22GxMqHBR+iQA3LQM0f7Ph5nHiE4ny3lddQ93PZv/54an3atwCZD0RWFxm2MwpmsSX9/9Nvin2RKQS6av5yl4DQ7vQMSEIgY/6Xpa10XBdhPBssrw3lcSkP+W+MSts5rzIWbKeCKzEb6be1zHQubboc1nP4iTSF/OUHUPstt3GIWShg0EQAgcVXWzORuP1UcKLU2RcR83kH8fWSEzCamqOQykqns/4/Ybyrbl1h2Q9kVpBGxYqGdHizYYuvaynAuuCIIQrtaeRZWwwPJMsPp8iewp+XvLtgfojz2U8jqFBh+Lj4v3irKdeO/hWYXOJrCdyVtAnWYqJY4SsjYbaXgwYOgr8NoEr8k6xdPNOw0dZyjtkfjuGiLZU7vyydMNf5j7v4eQu97nA2sT4jn82/fFn9r5U2nZK1hMpBPpHuuKWceKM74w5WmShQ0AQwmXq1tMdu43/PCnlLlXOC5Z9PvsPVbtX5699Y96LgXhBxlGlaBJ+N+2Bp3a/WNFRLfe5HosX/zNbccMOw+uFEsLQ7iEI4XIc1LHkbwwuStq/VDnWQ/YU/PLkhg/yP3pj7h+xv4SDywyb+Zup9/1ux7MFuhNyn2thiHBgqfKbamnxZkNjn9xnA0tCEMLo9BrosYPGZdsMr00T389QuCrlPZ3E2PvH/rv+1Ja3F74SgRVkgGh2eNqz6Y89m/3K3tr9cp8r0kPYvUSZFigmfaP/VymahnYLQQij8E21FP+VQddPhderlkfKfvN0DXY/sfuFktaytxe8EuDqJ/fpwFakaBJem/vC20f/9f6x/8q6GCkRKUV6LkXcdpXyg1Ipa6Mhvw1paIcQhDAiR9vEORuNz+dJH2QqPpyj8FPLfsbC5uJfbnl0rFfE6/P+6OnsIfv5wKaM9x67+qo3ytorfrfz2abeZrlPl+Aj5F6rvHmsuHiz4b59Qj16Su0LghAuYV8Tu3qr4fZc1e3RQt5y5Vz534sZMAz8/cjq53Nee3jqvfen/EIhyD4rA2yRl7Pna1nPT9Ek/XLzo+vLt8i69AwRiQL9apJYepMq0IWmrqf7c43V3Wgd2gkEIZzfkEQfn5JmrDfcsce4PEI8vmRw5URR1p0FiUhibEvlzts3PNCr712z5K20kGnyng9snCiIt02+4c0FL22u3PHA1scLm4vlPqOniv6YzAqXka8zTf3WcMMO4+5GxKHNk/lVB7BBh5vZ2lPSpxVSoq/wdKJ4TbgoCtTTI+9JJSbtOb3voxNfuCjVL2Q+Ees3Ud7zgR2J9Ap7d9Gr26v3vJj7RrR35IrJN0z2i5H1jH5q+tNUxZOJig/LpQdzjUZGt48Xb4sWItyxL6ZNGl0Q6vV6lUp15ceAtRmSKLeJfX9a+rqaqUS6LVo8eJ0ySv55EUTUPtC5rWrXN2WbfF18VibelhaSaoaTgp0RSFgQOWd2WNqmiu1/yn3Dz9X3uvGLM8NmOitk3BjaXUUPxIoPxIoHdOzDcmnqt8ZId+H6KPGqMCHBR0Ak2hCBsRG163t6eu68885t27YpFIonnnjiySefPPeYV1999eWXXzYajfPmzfvwww89PM7/gsMnn3yycePGjz/+eISX2N3dfaFScCUGjZTXyrK1LFsrZWtZzBjh6jBxeYQQ73OeR7inp8fNzU3g93R3Dnbtqz+cXbs/X3ciI3T60vGL5f4r3g709/c7OTkpFBg0vRiJSXtr9288te1ka3lm2IzMsJlTNAlO/BJxcHBQFMVz/9w3MtrTyL6tkbbWsW49mxcsZmqEDI0Q4yWISEXe9Hq90WhUq/m8tjfSIHz66aePHj26YcOGhoaG6dOnf/PNN2lpaWcecPDgwSVLlhw8eDAsLGz58uVxcXF//vOfz1sKQWgptb2srJOK2lhRO8trZSc72KQxQoZGmKUR5gSJPs4X+79cgrB9oPNka3m+ruh4U9HprrrUoOT00OkZodNdVS5XUtZxIAhHpbmvZffpfTl1B8vbKmL9JiYHxsf7TxrvM85FeUW/PS8UhGeq6ma7G9leLcttYk19LNlPSPIV4ryFOG9hgpdw8WcNRsIyQRgSErJmzZoFCxYQ0W9/+9u+vr5//vOfZx5w//33KxSKt99+m4h27dp18803NzU1nbcUglA+BolaB6llgDX1U0Mfa+ij+l5W1U3VPayyi3moKGaMMNlbiPcRknyERF9BPeLfqKMNwiHjUHNfq7ZX19ijreturOo4XdVZ06fvj/EdH+cfkxKYMMlvokrEEPXoIAgvT89Qb77uxLGmguKWsoqOao2bf6RXeIRXWKhHUJC7JsgtwNtlzMhfTh5JEJ6pY4iONLP8NlbUzk60s7JOphRorKcQ5SGEuVGYmxDsRkEugr8L+auRkSPFNwhH9Juor6+voaEhNjbW9GlsbOynn3561jGnTp1avnz58AE6na6rq8vT0/O8BYeGhtrb24c/HTNmzEV+w3b39bb32960HQOjXsOl/8joNdBZG4F268n0x4nEqNvAiMgoUY+BiKhLT8SoU88MEvXqqddAfUbqHWJdeurSU8cQ6zWQjxN5qwV/Z/J3EQJdaKyLkBlKoW5CmJvg/tMnd7CHBomIqE/fb7zArOQ+Q7+RGYlooH9AUAkGZiSinqE+RqzfMGBgxm5974BhoN8w2K3v7dH3dQ11dw52tw92DhiH/NRjAlz8gt0CQj2ClkbMifAIDT5zmdCBAWz6NlpsYEAyDAkIwlFyJZrpEzvTJ5aIjMxY091Q01VX3V1/4PQhbZ+uqa+lc6jb08ndy8nD08nDy9nDXeXqoXJzUarVSmdXpYuzQuUkOqlElVrpRESiJKgUKqVSSUROCidn8fydrs4KJyeFiogEolR3SnUnCv/xn1oHqLaH1fWyhn6qbWGHT1PzAGsbpOYB1megMSrydBY8lOThRB4qwVVJziJ5OZOCyF0lKAVyVREReTkREYlEHqoff3mqFOTy01tDrSDnEdwsHipb2sxFpVIGefvyrTmiIDSFlrv7j0v+e3p6trW1nXvM8AGmBlx7e/t5g7C8vHz9+vXbt28f/sq6desyMjIudPYX1j952rn9Qv/qCAQiuugcKff/397dxzR1tQEAf25bR4sgiLajkFVq/AC/EGFuGh2KiCGTaLTOd6QuApItm9kScB/JNmPCEv7QbTHO+LHAolm2VHFbVrdsGTq1zE1hQ0iaUYsl8llGoa1Cgfbee94/7nbTlyLW17aX9T6/v46nD/c+Sq9Pe+655wDEAaRIAbg3/Tgw42C/B3aAliCOr2ApyQMOr2BBSv4+fQyhZCzhOiUE5ISSsjCTgRgCcgZms9RMhopnII6GBIaKpwHACeAEsPJH6wvmb4tQmCkA0gH870izFOWWjdyXDt+T9o7IqBEpGZGAW0r6JdQYRbwS8FGElsA4BQAwJgH2n8/tXgq8D6ghXgnxBX0nQQqgBFDyl7APwAf06N+X0KQIiPTGo5SFo/nHkmbGMwxD0/RD42NjYyWShxT6oArh3LlzKYpyu90JCQkA4HQ6n3zyyQkxSqXS7XZzbZfLxfVMerSFCxfqdLrgh0aP/OcYDo0KLuSTZdD/AYdGwyf4dWwfdWgUhUNoh0aD+kIcExMzf/785uZm7o+3bt1KT584wS8jI8M/IC0tLTY2NiQpIoQQQuET7MhweXn5Bx980NXVde3atXPnzu3btw8ABgcHt2zZMjAwAABlZWUXLlz4+eefu7u7q6qqysvLw5g1QgghFCLBTturrKx0OBzr169PTEw8ceLE8uXLAYAQMj4+zs07Xbp06enTpysqKpxO586dO998880wZo0QQgiFSLCPT4TQoz4+8eGHH5aWls6ePTusWaGpnTp1qrCwUKPRPDwUhY3BYFiyZAn3MRQJ5Ycffpg5c+b69euFTkTUfv31V4fDUVRUFJKj/QsmzZ49e/bu3btCZyF2Fy5cMJvDvic4mtp3333X2NgodBZid+nSJZPJJHQWYvfLL7/89NNPoTrav6AQIoQQQuGDhRAhhJCoYSFECCEkagJMlqmurq6urn7Q4/aBenp6VCoVPr4qLLvdnpCQoFDg6thCGhgYUCgU/BJOSBBDQ0NSqZRbXQQJxe12MwyTlJT00Mji4uKqqqqpYwQohCzLWq3W4Avb+Ph4TAyuRCsw/C1MBz6fTyqVPnS9KBRWNE1TFIXr+wiLYRhCCLfi69TUavVDP8ELUAgRQgih6QM/WiKEEBI1LIQIIYREDQshQgghUcNCiBBCSNSCXXQ7MiwWy+eff07TdHFx8aQLKg4PD3/66addXV3r1q3bsWNH5DOMekNDQ0aj0Ww2x8fH79ixY+nSpYExNTU1DMNw7cWLF+fm5kY2x+hns9n8d67eunVrSkrKhBifz1dTU3P79u3MzMw9e/bgVNKQa21t/e233/x79Hr9hN3lzp07x22/CgBqtTpUS1+ijo6OpqYmp9O5e/du/ydVfv/9d4PBoFAoSkpK0tLSAn9wYGCgpqZmYGDg+eefz8vLC/J00+jiaW9vf+aZZwghcXFx69ata2mZuLk6IaSgoODKlSsLFix45513Dh8+LEie0e3AgQPffvutSqVyu92rV6+edDW/1157rbW11Waz2Ww2bhMuFFpNTU3V1dW2f4yOjgbG6PX6L7/8cuHChceOHXvjjTcin2TUc7lc/K/gm2++ef/99wMf+jp06JDJZOJi+vr6BMkz+vT392dnZ586derll1/+66+/+P7r16/n5eXNnTvX4/Hk5OT09PRM+EGPx7NmzRqLxTJv3rzi4mKDwRDsKcm08frrr+/bt49rv/XWWy+99NKEgPr6+tTUVK/XSwhpaGhQqVTcJlAohEZHR/n2gQMHdDpdYExMTExvb28EkxIdg8GQn58/RYDFYlEoFE6nkxBy9+5duVze398fqezESKfTVVZWBvZnZGQ0NDREPp/oxj0jSNM0ANy+fZvv37ZtW1VVFdfevXv3e++9N+EHa2trn376aZZlCSFffPHFihUrgjzjNPpGePXq1YKCAq69efPmq1evBgZs2LCB+1C2Zs0aj8fz559/RjrLaCeXy/n22NjYgxYx+eyzz44ePXrjxo1I5SU6fX19R44cqampsdvtga+aTKacnJzExEQA0Gg0Wq12wiAeCqHBwUGj0VhaWjrpq19//fVHH33kP5SNHtODxvmvXbu2efNmrj1pjeACKIriAlpbW4eGhoI642NkG2J9fX38umsqlcput5P/fdjfbrfzARKJRKlU9vb2RjpL0WhpaTlz5kxFRUXgS7m5uffv37darYWFhQcPHox8blEvPj4+MzPT5XJdvHgxIyOjqalpQoD/tQAAKpUKr4XwOXv2bFZW1pIlSwJfysrKAoDe3t69e/fq9fqIpyYiY2NjTqfTv0YEjkX7F5E5c+bIZLIgx6un0WSZGTNmcN+FAYCmaZlMxhV2nkwm4+doAIDP53viiScimqJodHZ2bt++/eOPP550ytKPP/7INcrKynJycvbv369SqSKbYJQrLCwsLCzk2pWVlQcPHvz+++/9A/BaiKQzZ87s379/0pf4DcYrKioWLVp08+bN1atXRzA1EeEWF/SvEYHveZlMxgcwDMMwTJDXxTT6Rpiamsp/qu3p6UlNTQ0M4O+Oer1eh8MROJUOPb7Ozs4NGza8/fbbZWVlU0dmZWXJ5XLcNjms1q5da7PZJnT6XwsA0NPTg9dCmNy8ebO9vf2FF16YOiwlJUWr1XZ0dEQmKxGaMWOGUqnk3/aTvuf9iwjXUKvVwRx8GhXCoqKi8+fPc+3z58/zE5FNJtPg4CAXcPnyZW7M12g0PvXUU+np6UJlG626u7vz8vJeffXVV155xb+/ubm5s7MTAMbGxvjOS5cuMQyzYMGCSGcZ7fh/ZELIxYsXly1bxv2xsbGR+4+goKDAbDbfuXOH63S73c8995xQ2Ua32traXbt2zZo1i+9pa2uzWCwA4PV6WZblO61W66SPG6FQKSoqqqurAwBCSF1dHVcjWJa9fPnyyMgIF2A0GrnLp66uLi8vL9itWh53fk/oOByORYsWbdmyZdu2bRqNpru7m+tPSkoyGo1cu6SkJD09fe/evUql8quvvhIu2ai1a9cuuVye/Q+9Xs/15+bmVldXE0IMBsPixYtffPHFrVu3xsXFnTx5UtB8o9P27ds3btyo1+tXrlyp1WqtVivXn5WVdfz4ca596NAhjUZTWlqanJz8ySefCJdsNPN4PImJiSaTyb+zpKSkvLycENLY2KjRaHQ63c6dO2fNmvXuu+8KlGYU2rRp06pVqwBg2bJl2dnZw8PDhJD29vbk5GSdTrdx48bMzMx79+4RQoaHhwGgtbWVEELTdEFBQXZ29p49e+bMmXP9+vUgTze9dp/weDz19fUMw+Tn58fHx3Odzc3NWq2WmyAHAA0NDV1dXc8++6xWqxUu06h1584d/gFhAIiNjc3IyACAtra2hIQEtVpN0/StW7esVmtcXFxOTk6QIw/okbhcrhs3bgwNDanV6rVr1/L3Ocxms1Kp5O/I/vHHHxaLZcWKFfhFJExGRkba2tpWrVrlP1+ho6ODoqi0tDSWZc1mc1tbm0wmy8zMnD9/voCpRpmWlhb+bh8ArFy5ktv3yuVy1dfXx8bGbtq0idsYjmXZpqam5cuXc3stMQxz5coVh8ORm5ubnJwc5OmmVyFECCGEImwa3SNECCGEIg8LIUIIIVHDQogQQkjUsBAihBASNSyECCGERA0LIUIIIVHDQogQQkjUsBAihBASNSyECCGERA0LIUIIIVHDQogQQkjU/gsCEhsiZg9QDgAAAABJRU5ErkJggg==",
"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.0563325319244079e-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.00022340550154961012"
},
"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.11.4"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.4",
"language": "julia"
}
},
"nbformat": 4
}