{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Custom potential"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We solve the 1D Gross-Pitaevskii equation with a custom potential.\n",
"This is similar to Gross-Pitaevskii equation in one dimension and we\n",
"show how to define local potentials attached to atoms, which allows for\n",
"instance to compute forces."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using DFTK\n",
"using LinearAlgebra"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"First, we define a new element which represents a nucleus generating\n",
"a Gaussian potential."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"struct ElementGaussian <: DFTK.Element\n",
" α # Prefactor\n",
" L # Width of the Gaussian nucleus\n",
"end"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"We extend the two methods providing access to the real and Fourier\n",
"representation of the potential to DFTK."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function DFTK.local_potential_real(el::ElementGaussian, r::Real)\n",
" -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)\n",
"end\n",
"function DFTK.local_potential_fourier(el::ElementGaussian, q::Real)\n",
" # = ∫ V(r) exp(-ix⋅q) dx\n",
" -el.α * exp(- (q * el.L)^2 / 2)\n",
"end"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"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": 4
},
{
"cell_type": "markdown",
"source": [
"In this example, we want to generate two Gaussian potentials generated by\n",
"two \"nuclei\" localized at positions ``x_1`` and ``x_2``, that are expressed in\n",
"``[0,1)`` in fractional coordinates. ``|x_1 - x_2|`` should be different from\n",
"``0.5`` to break symmetry and get nonzero forces."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"x1 = 0.2\n",
"x2 = 0.8\n",
"nucleus = ElementGaussian(1.0, 0.5)\n",
"atoms = [nucleus => [[x1, 0, 0], [x2, 0, 0]]];"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"We setup a Gross-Pitaevskii model"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"C = 1.0\n",
"α = 2;\n",
"n_electrons = 1 # Increase this for fun\n",
"terms = [Kinetic(),\n",
" AtomicLocal(),\n",
" PowerNonlinearity(C, α),\n",
"]\n",
"model = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,\n",
" spin_polarization=:spinless); # use \"spinless electrons\""
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"We discretize using a moderate Ecut and run a SCF algorithm to compute forces\n",
"afterwards. As there is no ionic charge associated to `nucleus` we have to specify\n",
"a starting density and we choose to start from a zero density."
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy Eₙ-Eₙ₋₁ ρout-ρin α Diag\n",
"--- --------------- --------- -------- ---- ----\n",
" 1 -0.143496231904 NaN 3.81e-01 0.80 7.0\n",
" 2 -0.156017549875 -1.25e-02 7.94e-02 0.80 1.0\n",
" 3 -0.156767019303 -7.49e-04 2.74e-02 0.80 2.0\n",
" 4 -0.157035650937 -2.69e-04 5.95e-03 0.80 2.0\n",
" 5 -0.157001609561 3.40e-05 8.16e-03 0.80 2.0\n",
" 6 -0.157056392414 -5.48e-05 1.79e-04 0.80 1.0\n",
" 7 -0.157056406877 -1.45e-08 1.17e-05 0.80 2.0\n",
" 8 -0.157056406917 -3.99e-11 1.35e-06 0.80 1.0\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Energy breakdown (in Ha):\n Kinetic 0.0380293 \n AtomicLocal -0.3163462\n PowerNonlinearity 0.1212605 \n\n total -0.157056406917"
},
"metadata": {},
"execution_count": 7
}
],
"cell_type": "code",
"source": [
"basis = PlaneWaveBasis(model; Ecut=500, kgrid=(1, 1, 1))\n",
"ρ = zeros(eltype(basis), basis.fft_size..., 1)\n",
"scfres = self_consistent_field(basis, tol=1e-8, ρ=ρ)\n",
"scfres.energies"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"Computing the forces can then be done as usual:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "2×1 Matrix{StaticArrays.SVector{3, Float64}}:\n [-0.05567808722973941, 0.0, 0.0]\n [0.0556776143063973, 0.0, 0.0]"
},
"metadata": {},
"execution_count": 8
}
],
"cell_type": "code",
"source": [
"hcat(compute_forces(scfres)...)"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"Extract the converged total local potential"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"Extract other quantities before plotting them"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=4}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3xTVd8A8HPvzR7N6N57t5S2lF2gyCqUIVMQUEQeeBEVRAXFgeJWHnlEUQQEFXEgQ6FsEAoyW6ClpXvvkb2TO94/wlN5GKUjTdLkfD/8kYaTm5OTm/O759wzEIqiAARBEAQ5K9TWGYAgCIIgW4KBEIIgCHJqMBBCEARBTg0GQgiCIMipwUAIQRAEOTUYCCEIgiCnBgMhBEEQ5NRgIIQgCIKcGgyEEARBkFODgRCCIAhyanYXCPPy8kwmUycTkyQJl4izIYIgbJ0FpwbL37Zg+duWBcvf7gLhxIkTW1paOpnYYDCQJNmr+YE6oNVqbZ0FpwbL37Zg+duWBcvf7gIhBEEQBFkTDIQQBEGQU4OBEIIgCHJqMBBCEARBTg0GQgiCIMipwUAIQRAEOTUYCCEIgiCnBgMhBEEQ5NRots4A1LtMJChWUCUKqlIFcBIYSSBgAA8WCBcg/cQIHV4IQZDtEBQoklPFCqpBC2QGwEABgoAALogQIFFChAOrZ2uBJe2YmnTgl3LyRD35dxPlzUGihEgIHzAwwEBBuRJcbAZFcrJMSfV3RaYFoTODkCA+YussQ5CzaNWDA1XkvkrycgvlzUFiRYgPB4iYQIMDggTZraBEQZarqERXZKwv+kQIEi6AP8/eBQOhozlZT31RQPzdTD0eiD4bif40ChUxH5xSi4OLzdS+KnLQn0Q/MbIqDkv3R+APDoLu8emnn167du3+5wmCwDCsS4dSGEGNmpIYgBsLeLDAY0yEhgIAQDMAzXcliwAglAJyA9hnoL7SAQ4N+HGBJ9vZf56jRo1avnx5bxwZsbdFq/38/K5cueLr69uZxDqdjsFgdPVcdFRHa6n11wm1CbzSD50VgnI7fZFjIsFvFeS/80kTCT4ZiE3w6+zPTaVS8fn8bmYX6jFY/taRmpo6bty4qKgoW2fEqV2/fj0/P//QoUPtz1jw/IctQkdQJKdWXSaq1OD9Aei0QBTt4nUjHQVPhqFPhqF/VJOrLhOfc8FXw7AwFye/+oSgf6SlpQ0fPtzWuXBqbDY7Pz+/lw4OB0v0bSYSvHuDHHEYH++H5k2nTQ/qchS829RA9NZ0Wro/OuRPfOMtkrCvzgIIgqBeAQNhH5Yvowb+gV9pIW88TlsZh1pkCCgNBSvj0CtTaZk15OhMvF4DgyEEQQ7OMoFw9+7d4eHhHh4eixcv1ul0D0t25MiRAQMGHD582CJv6uS2F5OjM/EXYtHM8TRfroW7MUP4yKmJtHF+6ICD+PE6GAshCHJkFgiEhYWFzz333K5du4qLi6urqzds2PDAZAqFYs2aNTKZTCKR9PxNnZkWB/PPEpsLyKwM2qKI3mrTowhY1x/99THa4vPExltw92MIghyWBarR7777btq0acOGDROJRG+88caOHTseOBJ19erVq1atcnNz6/k7OrMaNZV6GMcQcHkKLUrY6+NZRnghl6ZgP5WRz2QRRhgNIQhyRBYIhEVFRQkJCebHCQkJLS0tMpnsnjRnzpypqKhYtGhRz9/OmV1spgb/ic8PQ78fibGtNeDXn4ucn0yTGUDGcVxlstKbQhAEWY0FalOpVNo+mcPFxQUA0NbWJhaL2xNoNJoXXnjh4MGDSCfmg2o0Gj8/v/Y/f/zxx2nTpj0ssVPNIzxYi67Kpn07GB/rbVCrrf3uuwaDl7JpIw+R+0aY3Fl3WvwajaYz3ynUS2D5WwdJOkhnSEtLS1JSUl1dnfnPbdu2NTY2vvXWWw9Lf+zYsf3793/77bfWyuAjEAShvqvu6+T5z+FwUPQRTT4LBEKxWKxUKs2PzQ9cXV3vTrB+/fqUlBSFQpGTk6PRaKqqqqqrqwMDAx94NC6Xm5+f38kJ9RiGOUkg/Dyf/PwWeWoSliBm2CoPO9LA+uvExLPoqXTMPDyHoigej2er/ECw/K3jkdVoX0GSZFtbm/mxTqd79913c3JyOkg/fvz4V155paCgIDY21ioZfAQMw+4+4S14/lsgEIaHhxcUFJgf5+fnu7q63t0cBACQJHnr1q2lS5cCAKqqqnbt2oVh2BtvvNHzt3YGFACvXyMO1VAXp2B+lh4d2lXrkzAenRyZSZyeiAXyYFsEgmyguro6Nzc3NDR07969I0aMGD16dFFR0eHDh00m05QpU8xBS6PRZGZmFhQUsNnsyZMn3x/J9u3bl5yc7OHhAQA4c+aMh4dHXFzc999/P336dI1Gk5WVNXv2bARBnnzyya+//vrLL7+0wee0Igtc6SxatGj//v03b97UarUfffTRokWLzM3VDRs2HD16FACwcePG7P+KjY1dv349jIKdRFBgyXniXCOVlUGzeRQ0ezkeXRmHjsokqlRwWgUE2UBeXt7y5ctXrFghEolQFM3MzExPT0cQhMPhpKenZ2VlAQCuX7+elZUVEBBAUdRjjz2WnZ19z0EyMzPT0tLMj3fs2HHmzBkAwKpVqyQSSWlp6fr1683/lZaWlpmZab3PZiMWaBHGx8d/+umnGRkZarU6PT397bffNj9fXFzs7+9/T2JfX1+4OmInGUkw/y9CYaROTqR1fuFQK1gRg6IAjD5CZI5CouGXCTmf623UmmuEdd5rsAeyIfneuz9KpfLgwYMCgQAAEBERsWPHjtGjRwMAPD09P/zwwxEjRqSmpqamppoTM5nM7777bsCAAXcf4datWwsXLnzku0dERFRVVSmVSvP4D0dlmfp1yZIlS5YsuefJ3bt3359y//79FnlHh6fDwczTOANF/hxHY9rfPdDlMShBgUl/MbImU3bSVIUgqwniI2v6Weln6cZ6wJORkZHmKKhUKktLSzds2PDhhx8CAFQqVUtLCwCgsbFx2bJlhYWFdDpdp9Pd3zWqUqk4HM4j3918E06lUsFACFmbBgdTTuDeHGTXCIxmr/fpn49FVTpi7FH03CSaB9vWuYEgKxIzwRhfW17/MZl3NlczjxbcunWru7u7+Rnz0J433ngjKirKPFb/q6++OnLkyD1HcHd3b5/nxmKx9Hp9+3/p9XoW6074lUgkKIreM/7R8dhrLevENDiYfBwP5iM/jLTfKGj2fCQ+JwQZexSXGmydFQhySiwWKy0tbc+ePaL/MhqNAIDm5uaQkBAEQYxG488//3z/C1NSUvLy8syPExISTpw4geO4+c9Dhw7179/f/Dg3NzchIaE9Ljoq+65onY/SBMYewaOEyLZUrCf7SFjN+iTsMR9kyglci9s6KxDklLZt23b06NGkpKTHH388Li7uo48+AgAsW7bs9ddff/zxxwcMGPDA2WizZs0yD2YEACxdulQgEERERKjV6jFjxuTk5LSvlHns2LGZM2da7bPYCuwatSMKI5hwDE8QI18Nw/pCELxj42BsyXli6kk8czyNAa+sIKiXTZgwYdSoUe1/BgUFXb58uaqqqrm5OTg42NPTEwCQkZFRWFhYXl4eERHB5/PNzUQPD4/a2lrzq9LS0gwGQ15eXr9+/ZhM5r59+1paWsLDw3fv3j148GBzGq1We+jQoYsXL1r7E1odrLfshcwAxh7FB3sgXw/vS1EQAIAA8M0wjI0hi7PgDoYQ1OvodPo9Y+8RBAkODh48eLA5Cpp5eHgMGTLE1dWVwWCYx7ygKNp+KxFBkM2bN1+5cuXu9BiGeXl5tT9z7dq1119/3TzX0LHBQGgXzFEw1Qv5fHAfi4JmNBT8OhqrUlOvXLHSmHIIgnpo2LBh94z2f/HFF82DUc1Gjhz5zDPPWD1fNgADoe1JDWDMUTzNG9k4yP7mSXQamwb+HEs7VkfBPZsgqI96++23RSKRrXNhAzAQ2pjEAMYcwcf4IJ/25ShoJmKCYxOwLwrIPeUwFkKOw94Ggh09enTs2LG9dPDDhw8vX768gwRFRUXjxo174F57fRcMhLYkM4AJR/ExvsjHA/t8FDTz4yJHJ2AvXSZO1DvU7wRyWmVKqkZtXydzXFzcypUre+PIFEWtWbNmxYoVHaSJiopisViHDh3qjQzYCgyENtOiA6My8fF+yCeOEgXNYoTIb4/RFpzFc6X2VX1AUFc1asG4o4QX275u3OM4rtVqAQBGo/Hjjz+ur69fv3792rVrq6urtVrt559//sorr1y/ft2cuK2tbdu2batWrdqwYUNJSUn7QWpqat59991169aVlZVt2rRJoVAAAM6dO8dms2NiYgAAWVlZ5gVId+3aVVNTo1KpNm7caH7tggULvv76ayt/6l4FA6FtNOlAWib+eBDy3gCHioJmI7yQr4ZiGccJe7uUhqDOU5nApOP44khUyLR1Vv5Xfn7+F198AQAwGAxr165dvHixn5+fQqEYO3asefNzsViclpbW2toKALh06VJdXd3AgQMZDMbw4cPLysoAAK2trYMGDdLr9VFRUcuWLXvttdekUikA4MiRI+0TM44fP25ej+aLL76oqKiQy+Vvvvmm+b9GjRr1119/6XQ6G3z43gHnEdpAvYZ67AixMBx9vb/DXojMDEbrNSD9GHFhMk1kZ/UIBD2SiQSzTuMD3ZF1/dFj9/2vsbpI/vtX1skJIyRO+PjSDhK89957AwYMeOaZZ9zc3MaMGWMeCHrixImzZ8/OmjVr8uTJkydPpihKLpfX1dX9/PPPb7755vbt20eOHPnBBx8AAJKSkuLi4syHys/PnzJlyiOz5O7uzmazy8rK4uPjLfERbQ8GQmurUlFjjhLLotGX4x02Cpq9GIfWaqhpJ/Hj6TSWA7Z7IYdFAbDkPMFAka+GPfjEpXkGCGe/YJ3MoGxuxwnMYQxFUTc3t/aQ5uHhIZFIAAAFBQXPPvusVCrl8/lNTU0ZGRkAgJKSksTERHPKmJiY9pVLtVotm92phYM5HI5KperWB7JHMBBaVamCGnuUeKUf+lyMg0dBs08GYk+eJRacJX4d3TdWjIMgAMCb2USRgjozkfawWb0oi8PwD7duph4Kw/6J1uYVtwEACIKYB3a++OKLTz/9tHlf9Ndee80cHYVCYfuK2xqNxrzuDADA3d3d3EcKAGCz2XK5vP3IOp2uPUaSJCmVSu+evN/XOUV1bCduSam0I8RbSc4SBQEAKAJ2jcDa9NRLcKI91Ed8U0j+VkkdGkfjOEQzoT1iyeXy3377zfzkpEmT9uzZ09DQAABoHwIDABg0aFBubq75cVJS0qlTp9qbfb///ntycrL5cUlJCZ/PDwkJsdqn6G3OUiPb3JUWauxR/N+D0GcinKvMmRg4MJZ2poH6DE60h+zen9XkhhvksQmYu6Nst7B69epnn302PT09NTU1KSnJ/OSYMWOWLl2amJgYFBREkiSbzTbvTTh79uz2bSgmTpw4ZcqUyMjI4uLi+fPnHzx4cPPmzeaXZ2Zmzpo1C0EcqJOHsjO+vr51dXWdTKzVanEc79X8WMSJOtJjt/FIDWnrjFiYUqnsZMo6NRn4s+mnMqJX8+NsOl/+UGf83UR67DZmt977Ox0+fPj58+dtkqUHwnFcp9NRFGXuomx/Xi6Xm0wm82O1Wq3X682PJRJJbm6uTqfT6XQajaY9PUmSJEkWFBQIBAKSvPOpZ82atW/fvvY0Wq02Ojr6wIED7c8QBBEfH19UVNRrn+/BDh06lJGRcfczFjz/HaLxb98OVJHL/yb2PUYb7uVAF1Bd5MtFjkzARmfirkxkvJ/zlgNkt8qU1KzTxA8jaclu9n5+Yhhmvi+IIMjdK6LdvUwol/vPEBuxWCwWi+85yJIlS+Li4kwm0zfffLNu3br25t1HH31092R5NpvNYrGEQmH7M2VlZcuWLYuMjLToZ7IxGAh713/yyY23yJMTaXEie/919bYYIbJvDG36KTxzPG2A3dc1kFOp01BjjhAfD0Sd5yrtySefzMnJQRDkhx9+GDp0aPvzISEhL7744t0pn3766YCAgPY/IyIiIiIirJdRq4CBsLdQALxzndhbQZ2fjAXynOXX1bFhnsiOVNrk4/jZDFqkAJYJZBcURpBxnHg+Fp0f5kT370eNGnX3poYdeOEFK00UsSEn+uKtSYeDOaeJvxqoC5NpMAreLSMAeT8FSz9GNGjhojOQ7WlxMOk4PtYXWe3o83qhDsDv3vKadGBUJs7EwIl0uKjKAzwTgS6NQscfJWQGW2cFcm7m5WPCXJBP+v7eL1BPwK5RC7shoR4/STwTib6ZCGeQP9SaBLRNT006jp+cSOPCcxCyBZICT58jaCiyPbXv7YZNkuS3335748YNf3//FStW3D2YBeoG2CK0pN8qyAnH8M8GoW/BKPgonwzCooXI9JO4AU61h2zh+UtEvZb6dTRG64O14NKlSw8fPpyenl5dXZ2WlkYQ8FfUI/Bq3DJwEqzLJvZWUifTaf3EMAg+GgLAt6nYE2eIeX8RfbQygvqu164R11qpUxO7uQru7baSTde+sXSmHqyfR8yK5Gfvfqaurm7v3r21tbV8Pn/q1KkxMTFHjx41LyIKdQ8MhBbQrANzz+BMDFybRnOFNwU7DUPAT2nYtJP4M1nErpFwMVLISj64SR6uoc5m0Fzo3TxCoMBv9aCOtnG3IB793kW3i4uLQ0ND+Xw+AABBkKSkpOLiYhgIewIGwp463UA9dY5YEom+mYjCqryrGCjY9xht4nF82d/E1uF971YN1Od8nk9+X0qendSja1YunRMpDrNcprpGq9Uqlcr2P5VKpYuLi60y4xhgh1T3mUiwLpt46hzxw0js7SQYBbuJTQOHxtFuy6gXL8H7HFDv2nKb/LKAPD0R8+bYOis9U1ZWdvbsWQBAZWXluXPn0tLSbJ2jvg0Gwm4qklND/8RzJdSNx2mjfWAM7BEeHWSOp11poV66DGMh1Fu+LSI/ySNPT8T8uH3+B9u/f/8NGzYMHTo0JSXl/fffDwuzWfPUMcCu0S4jKPBFAfnhTeK9Adi/ouCVhGUIGOB4Om3sEXz1FWIjnNQFWdq2IvL9m+SZiVgQv89HQQAAl8s9ffp0VVWVu7v73cuKQt0D6/GuKZBRww/hh6rJS1NoMApalpABTk6kZTVSq68QcNUZyIK2FZHv3STPTMRCXRwhCrYLCgqCUdAiYFXeWVocvH6NSMvEF0WgpyfRHOwXZSeEDHAinXahiXrhIoyFkGV8eZv8INehomB0dPSyZctsnQuHAgNhp/xeScbuw2s0IG8G/V9RcFhMLxIxwcmJtBsSaukFgoTBEOqZz26Rm/LJs5McJwoCAMLCwubPn9/tl1MUpdfrLZIT85aHFjmUbcFA+AjXWqmRh/H3b5I7R2C7R2FebFtnyAm40MGxCbQyBbXgLGGC29pD3fV2DrGjmDw3yQG3f8nPz8/Nze04zbVr10pKSu5//ubNm1FRURbJRkxMzI0bN7r9chzH9+7dS5K2/5HDQPhQt+XUjFPE9FPEwnA0ZxptlLej/ZbsGY8OjkygqUxgxilCD0eSQl1EAbDyMnGohsrKoPn2/TGi9/vll19+/PHHjtNs2bLl7i127ZBOp5s9e7Y9rA8HR40+QJ6U+uAmebaRfKUftnsUxoaFZAssDOwbgz19jphwDD84liZk2DpDUB9hIsHiLKJSTZ2Z5JinTWFh4ZkzZ3AcX7t2bWho6JIlS3Ac37p1a3Z2tq+v74oVK7y8vK5cuZKTk1NTU9Pa2pqYmDhnzpwHHspoNH799dc3btwIDAxcsWKFu7u7+fmzZ88ePHhQKpUmJSW9+OKLBEHs2rXr6tWrOI6npqY+9dRTKPrQRtTGjRvT0tL27t3b0tIyc+bM9PR08/MlJSXbt2+XSCSjRo2aP38+giBbtmwBAKxbtw5F0eeee87f39/SRdVZsEX4P07VUxOP4+nHiBR3pHwOfXU8CqOgDdFR8OMorL8rMvIw3qi1dW6gvkCDg6kncbkRnJjgmFEQACAQCLy9vT09PZOTk82bxc+fP//QoUOzZs0iCCI5OVkul7u5ubm6uvr5+SUnJwcHBz/sUDNmzDhz5szs2bNVKlVKSopGowEAfPfddwsXLkxMTJw3b15bWxtJknq9vqSkZMqUKdOnT9+xY8d7773XQfa2bds2b9682NjY8ePHP/vss3/++ScAoLy8fMiQIR4eHlOnTt20adNrr70GADD30CYlJSUnJ9t2+Cus5gEAQG4EP5aSW4tIBIBV8eiBMSgTzmSzDygCNg3GPsolhx3CM8dj0UIH7OaCLKVFByafwONEyNbhvbuMu7JSW7Knrhff4C7CcG7YbN+7n/Hx8YmOjtbr9bNmzQIAlJeX//nnn/X19SKRaOLEidnZ2bt27Vq5cmVQUFBcXJw5zQPdunUrKyurvr6ex+NNnDjx8uXLe/bsWbJkyRtvvPH999+PHTsWADBhwgQAAI/H++STT3Q6XXNz85o1a9auXfvWW291kOdly5aZh/NoNJrPPvtsypQpmzdvnjFjxssvvwwACAsLS0xMfPvtt0ePHg0AmDFjBp3e3VVfLcSpA6GJBCfrqd1l5NFacoI/+uVQbKQ3AitaO7Q2AfXjgrRM/NfRtJHwZi30IMUKatJxYmG4NbYC5fmx4pYF9fKb3IEyHhHSy8rKgoKCRCKR+c+kpKQHjpG5X2lpaUREBI/Hu/uFcrm8sbFx8ODBd6fUaDRz584tKioKDAzU6XSNjY0dHzkhIcH8oH///ubMlJaWti8LHhMTg2FYdXW1r6/vQw9hXc4YCDU4OFVPHqymDlWTkUJkfhj65VC6GO4aYd/mh6G+XGTOGfzTQdiCMNilD/2Ps43U3DP4RwOxp8KtcW6gdJTlauOO1/Z5C0KhUC6Xtz8vk8nab/V1TCgUymSyu1/o4+PD4/EYDIZEIjHvbmG2c+dODMPMIe3atWvjxo3r+Mjt+ZHL5eYILRKJ2p/UarV6vb49ctsDZ6lQDAS40ES9f5McexT3/sn05W0y0RW5OZ3292Ta/0WjMAr2CWneyF+TaO9cJ1+/BqcYQv/YUUzOPYPvGU2zThS0B+7u7jU1NebH/fr1Q1F07969AIDa2tr9+/ebx6fcneaBUlJSlEplZmYmAKCioiIzM3PChAk0Gm3KlCnvvfeeeTBna2srRVE6nQ5BEAAAQRCbNm16ZPa+/fZbHMdJkvzmm2/Mnavp6ek//PCDORZu3rw5KSnJ09OTx+Ox2eyOM2kdjtkiJChQr6HKlOC2nMqXUjlt1G05FStCRnghL8ZiaWMRrmN+bscXLUSuTKVNP4nPOEX8MArj2/jOAmRjOAleuUocqaWyMmjhAifqM58zZ86+ffuCgoKGDx++e/fun3/++emnn37nnXdaWlrWrl07fPhwAMAzzzyzcOHCoKCg6dOn//vf/25/LYqiLBYLAMDn8/fs2bN48WI+n9/a2vree+8lJSUBAL766qtFixaZu1v1en1hYeFTTz31/fffR0dHkyQ5a9YsoVBoPhSLxXrg8NHw8PC4uDij0RgUFPTNN98AAObOnZudnR0ZGSkSiZhM5p49ewAACIJs2LBhzJgxKIru37+/vUPV+hB7WxfAz8/vypUrnew7blDo/1NMQxBEYwJKE1AYQYuOatKBBi3lzkLCXEC0EIkTIYluSKIr0r2tqKEOqFSqu/tPrMZIghcuEheaqYNjsTAHWjGkq2xV/nZCYgBzTuN0FPw8uncHiKampn744Yfm6GLPWltbxWIxhnW5pmtpaXF1db3nheZdDz09PZH/DpxoaGgQi8XmINqBqKionTt3JiYmqtVqNze3u//LZDIplUpXV9eu5hAAcPjw4Rc++mboe394scFngzBg0fO/b7eMMIQS0gGKIr4c4MIAQgZwZ6HeHODDQeCwTwfGQME3w7GtReTwQ/j2VFpGgPPGQqd1vY2adZqYGYx8kILBDZ3NOnlr8H4eHh73P8nhcDic/9m20cfHp/PHZLFY94dMOp3evSho5s0BE/yQ3liGuG8HQhc6eLUfgmHOcmMAutvSKDRBjMw5Q1xpRdYnwdrQiewoJl/PJrYMxWYEw9++3Xn33XdDQkJ648hiJjK/dwbKwdMI6sMGeyDZ02iXW6ixR/AGrX118kO9QW0CC88Sm/LJrAwajIL2afbs2Z6enrbORdfAMwnq29xZ4PgE2mgfdMBB/FgdjIWO7KaEGnAQZ2Dg6lRapDMNjYF6GwyEUJ+HIuCNRPTX0bRlF4hVl+Ei3Q6IAuDft8jxx/D1Sej2VLj8L2RhMBBCDiLVC7nxOK1BCwb+gedKYdPQcdRqqPFH8QPV5JUptCdCYZUFWR48qyDHIWKCX0djr/RDxx3F379J4rbf5gzqqe9LyQEH8VHe6NlJtCA+7A6FegXsYoAczYIwNM0befY8cbCK3DEC6yeGtWefVKehll0g6rXgZDrNtl9icnLynDlz7p8MQFEUAhcnthatVjts2LBeOnjfnlCv0+kYDEY3JpBCFmHPE7opAHaVkGuvEf+KQtf1xxxyOQV7Lv+eICmwtYh8O4d4IRZbk4DSbd1v9bBlpjUajW03D3I2Xl5ed09thBPqIegREAAWRaAT/NAXLhHx+/Cvh2FjfOHFex+QK6WWXSBoKDibQYuxj1232Gz2AyfGOeqFiBOCgRByZN4csPcxLLOW+tcFIsUd2TgI9ePaRd0K3U9hBG/lEL9UkO8PwBZH9vpWShDUztadDhDU+yb5IwUzaFECkHgA/+AmqcNtnSHof5EU2FFMRv9u0hPg9kz6szAKQtYFAyHkFNg08E4ydnUq7aaEiv4d31MO93GyF6cbqOSD+A+l5KFxtK3DMVe4JxpkdbBrFHIiwXzkt8ewC03U6ivExlvkxynwxqEt3ZRQa68R5UrwUQoK10uDbAgGQsjpDPdCLk+l7askn7tI+HLAewOwoZ4wHFpVkZxaf50810i+kYj9K8r240IhJwdPQMgZIQDMDEYLZtDmh6NPniXSj+EXm2FfqTUUyamFZ4mRmXh/V6RsDv25GBgFIduD5yDkvGgoeCYCLZ5Fmx6Ezj9LPHYEP1kPw7scrgUAACAASURBVGFvud5GPXGGGJWJRwmR0tn0tQkoF3ZIQfYBnomQs2OgYEkUuigC3VNOvnSZYGDg5Xh0ZjBsqVgGBcCJOmrjLaJIDlbFo9tT6Ty6rfMEQf8LBkIIAgAAGgoWhqMLwtHMGurzfOLVq+RzMejiSNT93nW1oM7S4GB3GflFPsnAwKo4dG4ovLaA7BQMhBD0DwSAjAAkI4CWK6U2F5CRe00Z/ujSaHQYHE3TFQUy6tsi8qcycqQ3+uUwLM0blh5k12AghKAHSBAj21OxTwdiO0vIJecJBIBnItEnw1Avtq1zZseUJvBbBbmzhKxWg2cikBvTaf5wHR+oL4CBEIIeSsQEL8WjL8Wj55uonSVkzO+moR7IvDB0aiAc6PEPEwlO1FN7ysgjteRoH3RtAjrRH8VgBIT6DvhrhqBHS/VCUr2wzTh2oIr8qYx87m9ivB86KxhJ90c5zvobwknwVyO1t5I8WEVGCJB5oegXQ+lwXRioL3LWHzEEdR2XBuaHofPDUIkBHKgivy0iF58n0rzRKYHIRH/U0zl6TRVGcKKe/LOaOlpLhguQWcFo9jRaAA82AKE+DAZCCOoyVyZ4NhJ9NhKVGcCRWvLPGmr1FVMoHxnvh4z1RYd6Ig42PJKkwA0JdaKeOllHZrdRw72QKQHoxwNpPhwY/yBHYMlAaDAYmMyOekZMJhOdDucQQY5DxARPhqFPhgGcxP5upo7XkS9fIUoU1BBPZJQ3OtwTSXFHmH1zT2CCArkS6kIzdbaROtdIerKRsb7I6n5YmjfitL3BkKOyzBldV1c3d+7cvLw8BoPx+eefz58//54Eb7/99s6dO1tbW7lc7ooVK9avX2+R94UgO0FDwUhvZKQ39kEKkBlAVhN5tpFadZkslFMJrsggdyTFHUl2Q8JcEHveYahGTeW0Udlt1OUWKruV8uMiw72QmcHIV0Pp3pxHvxyC+ijLBMKVK1fGxcWdO3cuJydn9OjRaWlpvr6+dydwd3c/efJkZGRkUVHRiBEjYmNjZ82aZZG3hiB7I2KCqYHo1EAAANDgILuVutJK7aui1mWTEj0VL0bixUg/MRIlRKKFiA3nY0gNoFhBFcioAhl1S0rdlFB0FCS7ISnu6Mvx6CAPRAxHvkDOAaGoni6uKJfL3d3dS0tLg4KCAAAZGRkjRox49dVXH5Z+5syZsbGx77zzzgP/18/P78qVK/fE0YfR6XQMBgPD+mbfU9+nUqn4fL6tc9GXyAwgT0rlSakCGVUopwrllJ4AYS5IMB8J5oMAHhLAAz4cxIcD3FmP7lPtTPnjJGjRU806UK8BtRqqRk1VqUGliipVUAQFIgVIrAiJFiIJrki8CIHNvi6B579tWbD8LdAirK6uZjAY5igIAIiJiSkvL39YYqlUeuHChWXLlj0sAUVRCoWCw7nzi+RyuQwGo+eZhCB7IGKae1D/6R6VG0GZkqpUUVUqUKakzjSARi3ZoAUtOoqJATcWImYCIQMIGAibBrg0wKYB1n8DpNFIYzAI82MTCdQmoCeADgcyI6UwArkRtOkplQm4s4AnG/HjAj8u4sdFJgeAYD4a5oLA1eMgyMwCgVAul3O53PY/+Xx+ZWXlA1OaTKaFCxeOHz9+zJgxDzuaVqsdMmQIit4Zdbd169aJEyc+LDFsEdqWRqNBEDu+5dUX0ACIYoEoFgDu9/6X0oRIDEBuBHIjojQBHQF0BKLDgYG8U+Y0ysgE/32MUUEswEApDg0I6EDAAAI6JWYCMeMhXT44UKt760M5CXj+21Yny5/D4bQHlIexQCB0c3NTKBTtf8pkMk9Pz/uTEQSxYMECAMC2bds6OBqXy+181yiGYTAQ2hBFUTwez9a5cFg8AHw6TKBSmfh82KyzGXj+25YFy98C052CgoJoNFpBQYH5zxs3bkRHR9+ThiTJRYsWSaXS33//HXZ1QhAEQfbDAoGQy+XOmzdv3bp1jY2NP//8840bN+bOnQsAyMvLy8jIMKdZvnz56dOnn3/++QsXLpw6daqwsLDn7wtBEARBPWeZ6RMbN25ctWrV0KFDvb29//jjD7FYDACgKArHcXMCtVodExPzxRdfmP+cOHHi/a1GCIIgCLI+C0yfsCw4faIPgcPHbQuWv23B8rctC5a/Yy2JCEEQBEFdBBcNdFgyvTy/tTC/tahaWduobpboZCYS1+N6FwZfwHLx4LiFioLDRcHJXv1d2SJbZxaCnIvaqLnenFcsKS2XVzWomxV6pdKopKF0JsYQMF28eZ7+Lj4xbpHx7jFeXA9bZ9bxwUDoaNq0kpNV587XXq5R1sW5R8W6RU0Jn+DN9XTliBkonUVjKQxKhUHVpGkul1X9XXd1c852D47bCP+hE0JGe3Lvm8sGQZDlKA2qk1Xn/qo+XyGv7ucRE+0amRE23o/vI2DyBUwXE4kbCaNMr2hSN1cpay/UXtlyfaeA6TLCf/BjgSMCBf62zr7DgvcIHUdOU+6Bkszc5oIRAUNGBQxL9OxHQx9dMiRFFkpKTlaeO1N9PtI1bG709CSvfp18R3iPxLZg+dtWl8q/TFb58+39lxuyh/ikjAsZlegRT8cevRUPSVFFkpJzNRdPVp315/tOjUgfFTAMReAtLQAsev7DQOgILtZf/TH/N61JNytqymNBI9m07kyyNhGmU9VZv9zez6KxliQsGODd/5EvgRWxbcHyt61Oln+prGLbzR/L5VWzoqZMDhvPpXdnRVecJC7UXd5XfFiqk82PnTkuJA1DnL3eg4HwDhgIb7cVb7n+nQ7XPxX/xHC/wWiPF3yiAJVVc2l77o8eXPfnkhaHCAM7SAwrYtuC5W9bjyz/Vm3b1hs/XG/OfSr+iYmhY+moBW5F5bYUfH/rl1atZGniU8P9BvX8gH0XDIR3OHMglOpkW67vzG3JX5wwf1xwWs9D4N0IijhcdmJn3p5xwWmL+s17WBMTVsS2Bcvftjoof4Iifi869FPB79MiJs6Nmd69TpoOXG24/vWNnUKWYFXKsgAXP8sevK+AgfAO5wyEJEUdLDny/a1fMsLHLYidzaL11q5xcoPim+u7rjffWjP4+WSvhPsTwIrYtmD529bDyr9cXvXhxU0ilnBlylJfvncvvTtJkQdKjvxw69eMsLFPxT/BwJxu6UoYCO9wwkBYo6z/5PIXKIK+PGhFgEunSqmHcppyP768ebBP8vKkRaz/vbCFFbFtwfK3rfvLn6TInwp+31d8aFniogkho62QB6lOtjlne5ms4tXBL8S7O9dyXTAQ3mGpQEgYSKPCZNIQhIEkDAQAAKWjGANliuhMER1B7WKnFQpQvxcd2p2/9+l+T0wNn2jZvtCOaUzazdnbCtqK3hz2coQ4tP15WBHbFix/27qn/Js1rRv+3sikMV4b/KIbx9WaOTlfe3lT9tbHAlOfTZhvP01Dg9yklxhJE2WuVDEGijJQOo/G4NFoXAu0XmAgvKMbgdCkxrXNBl2rUd9m0LUa9RKjQWaiCIohoNF5NJSB0lgYAIAwkqSR1EuNJhXO8WGJIniiaL5LsM028G7Rtn14aZOJwF8futKH52WTPJyuPr85+9v5sbNnRGUgAAGwIrY1WP62dXf5n6+9tPHqljnRj8+JftyaF6ntlAbVxqtbqpV1bwxdFSYKsX4GzLSNeultlaxYrarR0Zgoy42BMVGMeadSJQykSY2bVDhJUCwRneXKYLkx2O5MtjuD7c5kCumgKyUHA+EdHQRCEqeMCpNBZtLLjHqJSd9m1LUZ9C1GgAG2O5PjwWS7M1huTJYbgyWi0zgPDaUkTqmqtfIStSRPiWCIT6qrxwAhgln1RD9b8/ema9/Mipo6N2a6bacQNaqb377wsTfX89XBz3PpHFgR2xYsf9sylz9OEltvfp9Vc/Gd1DVRruG2zdKJyrNbru+YFztzVtQUpEtRpYcoIMlXNpyX6NuMbgkuwgieSwgXYz60siIMpEFq1EmM+jajrtWoazXoWgy4nmS7MdjuDJbrnZqZKaIzhHSM8eDjwEB4h0qiaTotRxCEIilCTxJGEtcRuJYwqXDSRDJc6Oa+zTvXHW4Mthuz+01yCshL1fVn2wwyU9hsX+u0DvW4YXPOtpvN+W8OW23z35iZiTBtztl+vTlvw4jX3FARrIhtCAZC21KpVCYa/tb5j3kMzutDV7kw7OK7aFQ3b/h7I5fOeX3oShFLaIV31LUYyvY2EEbSb7Sba7xLt+8lEQZS12rQtxl1bXf66gxyk1FuAgDQ+TQ6F6NxMIyJ0fm00OneAAbCdmq5RlmgQ1EUIAiNjaJ0lMbB6ByMzqfR2L01gkaSp6w42OjazyV4ilev3j6slFevv/BJpGv4qpRlFh9+3UPmC8+lcU+lR46xdV6cFwyEtnW1+von17+cEj5+Qdxsqza/HoWgiF15vxypOPX6kJUPHO9tQQ1ZktpTrQHj3L2HufZSGRAG0qjEcS2B6wjCQFAkcE8UABgI29lq1CiuJ4p/qAUIErXQv4Pmf08cKju+/ebu5UmLxltl7Fk3lMkq1519f3RQ6pL+C+CaTzYBA6EN/Vl6bEfuT68PXTnIJ9nWeXmwG8233r/4+fiQtGf6zeuNZWgokqo40Kis1MYsDmSKHr1cnMXBQHiHDadPUCRVvr9RVa2N/7/gDm4xdoPGpN145atqZd3bw1+1zgSJbmuQNH52cwsNpb017GUeg2vr7DgdGAhtAieJL7K/zWspeH3AygivMFtnpyNyg+L9i5/rccNbw1a7c9wseGSKoAp31VAEFfVUQC81Bh4J7kdoewiKhM30EUbwbn9XQ+IWu5gokZYvObqKx+B9Pf5TO4+CAAA+g/fp6PUBLr7Ljr9co6y3dXYgqNfJDYqXTr8h0Um3jP/Um+tp6+w8gpAp+CTt7cE+yf86tvpSfbbFjkuB0t/qAQAxzwbaKgpaliN8BhsKzvBiielF39dQpAVi4f7izFf/Wr8kYcFLA//PfiYDdQxDsBXJzz4ZM+OFk69dbbhu6+xAUC8ql1ctO/ZygkfchhGvc+hsW2enUxCAPBk7893UtZ9f+/rr6ztxkuj5MauPNmubDZEL/O1kjnXPwUDYMwgIm+1L6MnaU609OYzKqH4j68Njlae3jP80LXC4pXJnNemhY94b8frHl7/4rfCgrfMCQb3ifO2l1aff/Ff/hYsTnrTJTMGeiHeP3p6+qUZZ//zJNY3q5p4cqi1X2XpTEfts4MNmNfRFjvNJbAWlIVFP+TddlKqqdd07wq3WwsVHVnpzPb4a94mtJsv3XJx71NcTPjtZde7DS5tMhMnW2YEgi6EA9f2tXzfnbP8kbf3owFRbZ6ebXJj8D0ateyxw5P8df/mv6gvdO4hRiZfvb4ic70/nOdSm7jAQWgCdRwud4VO8u5YwkF16IUmR39/65a3zH61KWfZc8mKL7NJiQx4ct81jPzIQxhdPrZPoZLbODgRZgB7Xrz//ydXGnG/Gf3b3+oJ9EQKQmVGTP0lbvyPvp48vf6HH9V17PQVKf633GSbmB/SNbuHOg4HQMlzjXVyCOVWZXehzaNK0vHDy9dyWgu3pm4b4Dui9vFkTi8Z8e/grQ3wHLDu2ukhS2ntvROKUttkgL1Gb/2mbDF29CoH6IhKndC0GRZnG/L1r6vWEsRe/9yZNy3Mn1nDo7E2PvS9mi3rvjawpQhy6Lf1ziqKWHF1VLCnr/AsbL0lxLeE3xr338mYrfbsJYldCHvfO+ajUa4iI6/3oye/HK85sub5zbuz02VHT+tz9ho4hAFkQNztEGLT27Lv/l2jheZC6FkNbrlKSp9A2G5gixp3FCSlgVJr0MhONjfED2S5BHGEkrzPfAtQn6KVGeZFaWalV1WgNMhNTRGcI7iyFb1LhulYDU8Rwjee79RPwLNpSudF8692/P5sfO3NG5GQLHtYesGmstUNe/Kv6wpqz78yMnDIvdsYjpwLjWqLmWEv88mCHGSBzNziP0JIazkukt1VxS4M6SCM3KDZe2VKranhj6EthomBrZa1XdDyPp1pRuy7rg4HeSc8lP9Pz+bx6qbH6aIuiRO2WKDC3v+//QeqlRlWVTlmlkRWqKYJyjXdxSxS4BHLsadEPS3LseYTaJkPrDbkkT2nSEqIoniCUyw/kcDyY936bFFA36CV5itbrCpYrIyjDk+dngXC4t+jPPbf3vTXs5UTP+IelcYDyb9W2fXjpP3rc8PrQlX58nw5SVhxsJHEqbGZHaawMTqi/w94CIUVQ1z8tC5nqJYp+8Ndzrubif7K3jg8e/Uy/eXTMBmsxWNYjT0S1UfPexY1ak+6d1DU9Wfaw/lxb7alWn1RX35FunZy3pGsxtN1UtN5UkEbKI0XoMUDIcu0bM1I6zwEq4vuZ1HhLjrzlmhzXEW79Be4JAp4/uzOXMhRBNV+R1ZxocUsQBE/x6vbK+Hrc8OmVL2uUde+NeM2T69FBSscofwpQB4qP7Lr181Pxcx6PyHhgB5Wu1Zj3RXnSmnC7GiMDA+Ed9hYIAQDSAlVVZlPiy2H3tFdkevmma1srFTVrBj8f6xZlq+xZVmdORJKivr/1S2b5yfXDX41z7/IHJ4xk2a/1Ookx+ukAprA7lw6aen3zNVnrdQXXh+U1ROQa52LlzUN6j2NUxHdQQF6qbroklZdoxHF8z4EiQQi3G015Qk+W7KkzaYmop/wZ/C7X2vWqxjezPgwXh7w0cDnzUXN5Han861QNH1/eTFHkq4NfuH8pj8KdNfxAtt9o+7o7CAPhHXYYCAEAeZsrfFJd3foLzH9SgDpSfmrbzR8mhY59Kv6JvjJTvjM6fyJeqs/++PIXT8bOnBk1ufPLExMGMv+bSo4nK3SmD0rrUfQicUqar2y6JNU2GzwHiryGiG2yOqJlOUZFjGuI5muyxotSGhP1GiJ2TxJirJ4N4qNA7anW5iuy+BXBXbp4Mu8puKjfvKnh6Z1J7xjl346kqD9Kj+zM+3lG5OR5sTPaB7FrGvQF26oHrIvo4W/Q4mAgvMM+A6H0tqrmWEv/l0IBAJWKmn9f/RonTasHPtfX7wjer0snYpOm5a3zH3lxPdYMfoFLf/QmViRO3d5ezRLTw2b5WvAmn67F0HhR2pIjdwnieA9zFUXy+u4dxL5eEatqdI1/S6T5KnEc33uoKz/QkkNd6s+1NV+Wxa8IoXdi5zUTiX974/us2ktd2lOwr5f/AzVrWv+TvbVO1fhSyrL+nvEAgJKf6jjeLL/Rllyq1CJgILzDPgMhoMD1T0t9MsT7dAdPVp1d1G/e5LAJDjY01KyrJ6KJMG25sfNyffZbw1+Odo3oKCkFin6oAQCJXODXG6PUSCPZekPR+LcU1xFeQ8WeA0WdqS7tTR+tiAkj2Xpd0XRRiusJ7yFiz4Gi7u8S2qHqI82yYnX8c8Edr4HSqG5+58KnrmzRmiEvdGlPwT5a/p1xvvbylznbY90il4QsrP5amrIusqfN9F4AA+Ed9hkISYo8mfl323VF4/jqJf0XCJguts5Rb+neiXi+9tK/r30zO2rqnOjHH3Z90JAlabupiH8uuLfv56lqdE0XpZJbSlE0z2uIuHv3pWylz1XEmgZ90yVp6w2FIITrNVRsheZ4yZ46lI6GzXroWMfT1ec3Z2+bHztrRlRGV/cU7HPl3yUGwvhTwe+q48ZQt8AxC4aw7GxLVAADYTs7DIQX6q5sv/mjG9N1xt+z45eE8nzt7uyxoG6fiC3atvf//jeKoq8PWXn/7jCaBn3+N1UJK0NYYivdT8V1REu2vOmSlCKA5yCRR4qwG+MsrK+vVMS4nmi7oWi+IjMqcc9BIs9Bou6Ne+oGwkDe2FgWPNnLNf7e61GNSbvp2tZiadmbw1aHi0K6cfC+Uv7dhmuIqx8UXxj911VlztPxc9JDx/TGvobdBgPhHXYVCHOacnfk/mQgDEv6LxjsM6DurzZdiyF8jr1vpdQTPTkRSYrac3vf70V/LE9aPC541D/P49TNf5f7Pebmkdz96RbdpqrWNl2WSfKULsEcjxShONbF3gYI3M3OK2KKpBSlmuZsuey2ShjB9RxkjSbg/VTVusLvqvuvDmO4/HNxc6P51keX/jPIJ3l50jMsGrObR7bv8u+5+rNt2kZD+FzfIknptzd/aNG2Leo3Ly1guJ3c6IGB8A47CYTZjTe/z/9VYVA+FTcnLTDVfJaY1HjOh6Upb0U6xn5dD9TzE7FUVvH+xc/9+T4vDfw/80TD6mMtuhZD1EJ/C+WxOwgjKclTtmTL1XU613gX90SBIIxrhwtq2GlFTAFVja71hrztpoIppLsnCz2ShL10F7CTao63aBv1UU8HAAD0uGFb7g/nai69OmjFQJ+knhzWTsvfcq5/XBo2x9cl6M7Qtpym3J15e1RG9cK4OWmBwx+5GE1vg4HwDtsGQoIismou7bm9z0Ti82Nnjg4ccc+FUtGuGmEU32uwgyxReD+LnIgmwrTr1s+Z5aeWJz0zUjTs5say/qvDrNZ11jGjEjdX6HqJyTXexS3BRRDKtZ9piPZVEVNAVaNty1NKcpUoHXHrL3BPErDdu9nYsiwSp65/VBo+17eCV/HJ5c1x7tHPD3i2S+NiHsi+yt/SlJXast/qk9bcO4Y2u/HmD/m/tmolc6KnTQh5rNvt6Z6DgfAOWwVCpUGVWX7yQEmmN9dzTszjQ3wHPPA2u6xQVXOiJeHFvr1ifQcseCKWSMs/vvzF2KLx8WFRsVO6c8OmV+mlRkmusi1PqWsziKL4rrF8YSSPxrZxV4Q9VMSkkVSUayT5KultFY2NuvUTuPZz4frY3a3x2mtNBccqd8bvWDVw2WAfy6xxbw/l33tKf67n+DB9Rz541kRBW9HPtw/cark9MXTMtIj0jpfg6SUwEN5h5UBIASqv5fah0uOXGq6l+g+ZEZnR8T12iqSy3yuJXRLIcdAFoC1bEcirNLk7Sr+I/8/jMRPnRE+zzyXojEpcelslLVAqyjVcb5YokieM4PEC2DbpOLVZRUwBTZNeXqKWF6uVlVqeP1scy3eNdWG52eNiERSgjlWc2Xbjx6XFS2MfC/UbZLEq24EDIaEnr20oTn7tEWuqNaqb95dkHqs4HesWOTlswmDfZGuOpoGB8A6rBcI6VcOpqnPHK/5iYoxJYePGh6R1sl+l+lgLoSdCpnn3dg5twrIVwa0tlR4pQirG9GXO9kp5zfMDnrXUlXtvIHFKWa6RFavlJWqDzOQSzBGEcl2COTx/ttX6Tq1ZEVMkpW0yKCu0igqNolyDMVFhBE8UwRNG8Oxwhlm7YknZpuytAICVKUt9VL5F39cMWBdhqS/IgQNh0yWpvEQd9VRAZxIbCOOZ6vOZZSfqVI1jg0aOCR4ZKQ7r7RwCGAjb9XYgrFHWX6i7/Ff1BYlOOipg+PiQtK5+wXqpMfc/FQPfjrTDoRY9Z8ETUVWtK/6xNvn1cHNBXW28vjl7mxfXc3nSomBhoEXeoveY1LiiQqss1ygqNPo2I9eHxQ/k8APYPH82S8zovXGSvV0RGxUmdZ1eVaNV1ejU1Tq6C80c7wVhXDu5iduBVm3bttzd2Y03lvRfOCFktPnmRcG3VW79BZ4DLXPb3oEDYd7mCv8x7g/bPOBh6lQNJyr/OlWVhQBkVOCwEX5DIlxDuzo7s/NgILyjNwKhkTDeai280pBzqT5bh+uH+w0aGTA0wSO220OkcjeVB07yFIbzLJhJO2HBE7FwV40wjOs93LX9GZwk/ig98mP+b8P8Bi2Kn+vGce3g5faDMJCqGp3aHDxqdYSB5PqwuD4sjjeL48nkeDEteGfRshUxYSB1LQZts0HbqFc36DUNekABnh+LF8DhB7D5gZy+sviOxqT9+fb+P0qPTgmf8GTMTA79n5Xb5KWaiv0NSa+GW6RydtRAaFSYbnxWNnB9VLebzsXSsnM1F8/XXtbj+kE+yYN8kpO9Eu7+IiwCBsI7LBUI9bihSFqa11Jwszm/UFISIgw0f3kRYgtcztT/1aaTGO1qHy9LsdSJqGs15G2uTHkjAr1vKSy1UbPn9r5DZccnhY59IuZxIVPQ87ezJpOG0DToNA16bZNB26TXNRsRDLDdmSw3BsuVwXJlMEV0pojOcKF3Y8Ji98qfIimTCtdLTQaZySAz6iVGXZtR12rEtQTHk8n2ZHK9mObIbf/Nvnvocf2BkiO/Fh4Y4puyqN88j/vWagAA3Py8PGCchzjWAuetowbChiyJplFvkTnQtcr6Sw3ZVxpybrcVBwkCkjzj4z1i4tyieQxuzw8OA+Ed3Q6EetxQqaguk1WWSMuLJKU1yrpQYXCce1R/z/gEj9jOLAndhfdy3N5RS52IZXsbGC60gPEPHcXQppPuzt97ujprctj4WVFTRaw+Fg7vZlThuhaDXmLUS4wGmUkvNRlkRqMSp7ExBp9GN//jYDQuRmNjNDaGMVGMiWIsDEGBuTWJYIh5cqparebxeAAA0kSRJhIAQOhJiqJwLUGaKEJP4HoS1xEmNY5rCZOGMClNRhVu0hB0LsYUM5hCOkvMYLnSWa4MtjuTKaT3oeXl7qHH9X+UHvul8ECCR+yi+LmBgofOQ227qWg4L+n3vAVGJjtqIMzbXOE/1kMUZclOLCNhLGgrvtF861br7SJJqTvHLdo1PEIcGioKDhMGdy8uWrD8+8A6Uj2kMChbtZIWbWujurlO1VivaqxR1kl1skCBf5goOFwUmh4yJlwU3HtjFFliBktEV5RrheEWuAhyPCY13parSH6tozW43djilSlL58XO+Kng94WHlo8NHjUnepon1752R+skBp/G4NMEofeeDCYVblTjJhVuVOG4hjBpCV2rkdAThIEkDCSuJwAJcB0BAKAIijCQAACKohAEAQCgdASlowAAjIUiCELjYCgdwZgYjYXSOBhTROf5suk8jO5CZ/BpdB7mSNdkSqPqQPGRAyWH+3vG/3v0u4+8o+zaz6Uq2oCepwAAIABJREFUs1lVrbPsZhcOw6gw6VoMFq+sGBgj0TM+0TMeAEBSZKWipkhSWiwpO1N9vkJezaGxAwX+vnxvX763D8/Lg+PmwXETsoRWW8Kmb7cIJSrp0erTKIpqTTqCJPSEXo8btCadyqhWGFQKg1Kul7PpbHe2qyfX3ZPr4cv39uP7BAr8vLge1lwWoe6vNr0j9o5a5Iqs/lybttEQ/kRn+2GkOtlvRX9klp9M8U6cEz3NOuPT7JOjtkg6qUHdtLfoj1OVWcP9B82NmXH/drIPU/dXm77VEDa7p11/Dln+FuwX7bwWbVuNoq5O1VCvamzUNLdo2lq0bSqjSsB0ETBdXJguLgweh85mYkwuncOiMZ+KfwLAFmE7kiLVRg2CIBw6G6Njbpgri8bk0jk8OlfAchEyXYQsYfv2kjbkluCSu6kidLq3I12JW0rzFVnYrC786sRs0bLEpxfEzT5cduLNrI/c2OIZkRkjAobawxcNWQFJUdmNN/aXHC5sK80IG7sr40tXdtdGgXoOEOZ8VBo81duBlz/strZchf9Ya8+ONzcBB3j3v/tJgiLkeqXcoFAaVAqDUofrDbhBi+tYmOWnZfftuoNH5y5JWGDztUYfiSVmMEV0ZaX2/g4xJ6eq1lIkaF/MsPO4dM6c6Gmzoqb8XXf1QEnm5pztE0JGZ4SN8+M7WrMbaifRyY5WnM4sO8FjcB+PmPRO6lom1p0p/HQ+TRDGbbup8BzksMsfdo9RhWubLd8v2j0YgrmyRV29yumevh0I+xBxDF9WqIKB8B7NV2Seg0TdHqOBImiq/+BU/8F1qobDZSeeP/maL89rQshjowKGWWRYGmQPjITxYv21YxWn81uLRgYMXT/81UjXnvaHew4S1Z1uhYHwHrLbKmEEz35W07UaGAitRBTFL/21PijD1vmwJ4SRbMtTJr1676q+3eDH91mW+PSS/guuNOQcqziz5fp3SV4JY4JGDPYZYMNFgaGewEnienPumarzF+quRIhDJ4SMXj98jaW+TVEUr+y3em2zgeMJT49/yIpU4hiH3Ui8AzAQWgk/gG1S4waZiSnqY3Ozek9brkIQwr17l7gewhBsqO/Aob4DNSZtVs3FzPKTn1zenOKdOMJ/yGDfAZadFQP1EiNhzG66eb728t91V/1dfEYFDF/Sf6HF+8cQFPFIEbVclQVN9rLskfsuiqDkJZrQ6c54cwEGQmtBgCiKJytUeQ0V2zor9qL1uqKX9qji0jnpoWPSQ8coDaq/666cqsraeHVLlGv4YN8Bg30GdH54IWQ1rdo283JON5pvRYhDU/0HP2xGvKV4JAkKtlUHZXj13dmTlqWs1LLcGXS+MwYFZ/zMtiKK4rfekMNAaIbrCHW1TrSoU6v6dpsLk2+OiHpcn9OUd6n+2u9FfyIASfFOTPZKSPLqJ2A6Y0eQndDh+pvN+debcq813ZTqZCneiWmBw9cMeaHnOwV2BsebhTFRVQ2cUHiHrEgl7uLiog4DBkLrEUXzyn9vIHGqG4tpOR5JnlIYycXuW1Otl7BorGF+A4f5DQQAVCtqrzXePFF59tMrX3py3ft7xvVzj413j+4ry5n2aUqDKr+tMK/ldm5LQaW8OtotItkzYe3gFyLEYVabPd3ONUHQlqeAgdBMeltt5emD9gMGQuuhsTGON1NZrhFGOuAC3F3Vlqf0TBHa5K0DBf6BAv+ZUZNJiiyRlue2FJysOrvp2lYmjRnnFhntFhntGhEmCoajbCwCJ4kKeVWRpPR2W/FtSUmbVhLtFtHPPWZZ4tPRruGMbs1/sBS3fi6F39UEw95RAAxyk0mN8wOc9JoABkKrEkXzpYUqGAhxHaGq1EYtfOiCkNaBImiUa3iUa/ic6GkAgFplfaGk5HZbyemqrEpFjTfPM1wUHC4KCRUFhwiD+vQCp9akMWkr5FXlsqpSWUWptKJaWefL84p0DYt1j5oZNSVEGGjNRZ06xvVhIRiirtPx/J00ALSTFapEkTynvSCAgdCqhOG8sr31ts6F7UnylYJwrr2t6+Hv4uvv4jsuOA0AgJNEpaK6VFpRKqv4u/5ahawKRdAQUWCAi1+gi3+Ai6+fi48Hx936vXn2pk0nrVM11Crra5T1VfKaKkWN2qQJEgSECoMixKETQ8eECu26be2W4NKWq4SBUF6qcdobhAAGQivj+bMMMpNJjdN5Tl3ykjyle3+7bmDRUCxcFBIu+mePAolOVqWoqVLUVilqLtRdrlXWyw1Kb56nL8/Lh+/lxfX04np4ct09uG59bq+oztCYtC2a1kZNS7OmpUnd0qBualA31asaWTSmH983UODn7+Kb5NkvWBjgyXXvvb1YLc41QVC0qyYow9PWGbEpCijKNMFOPJPEqatj60NQxCWYoyjXuCU4YF3ZSaSJVJRrIp70s3VGusa82lOyV0L7MwbC2KBqrFc3NaibGtXNN5rzmjVtLdpWPW5w57i6scXuHDcxSyhmi8QsoZAlEDIFQpZAwOSzaJZfLLGHjIRRaVDJDUqZXi43KGR6hVQnk+hkEp20TSdt0bSiCOrBcfPkenhy3b15ntFuET48L1++d1+fncnzZQEKOPnMem2THmOizjzFGQZCaxOEcxVlTh0IFWUanh+bxrL3FWIfiYkxgoWB9+/7YyCMrdo2iU7Wqm2T6uUSrbRSUSPTy2U6ucKgVBpVJEXxGVw+g8dn8Lh0LpfO5tA5PAaXTWMxMAaPwWWgdCaNycAYTIyBISj7v1t78xn/c3dZo9WoEM3dz2hNWoIiAQB63ICTOE7iOlxvIk163KAz6QyEUWvSaUxarUmnMWnUJq3aqFYZNSqjiiBJAZMvYLoIWQIxSyRiCcRsUYgw0JUtdmWLPThuFt9e3H4Io3iyQpUzB0J5mcZO1he1FRgIrU0Qxmu6VGvrXNiStFAtinLkuxFMjOHH9+lg+W8DYVQZ1SqjWm1Uq41arUmrxXUqo1qP6+V6Rb2q0UAYjYTRQBiMhImgSJ1JBwCgAKU2/k/YI0kSRf/nPiuHzsEQFADApDHpKA1DMQ6NTUfpLBqTTWMxaAw+k+fF82DT2DwGh0vn8v4bj9n210i1GnE0vyFL4juqFyfv2zknvzQHMBBaH8+HhWtwo8LEEDhpR4SsSBXdy/Po7RwTYzDZYjd2T5dWcMj98KxPEM4t/qmWMJD2NnrLSiigrNCEOtxuqV3ilF+8bSHAJZSrKNc8OqUj0rUYKJziejlv+wOyNxgD5Qdw5KVqW2fk/9u78ygpyrtf4NV7Ld09CzM93dOzsA0gqAyCKOIVBVE4SsS4E5eIiSiS45pzNPiir0mOGnxz1OSqkYAek5uoGDkRvDcaj6IkgAgqSwSGZZi1e6Zn7eqq6q267h9tJuOsPdPVXdv381dPT03XQ9PV36p6fs/zKCPSIticVrshZ1brgyBUQOFUpvekQYOw+1ik6CyXdooKwRCKznJ2HzVoEPae4ApqjD6yGUGogIIaZ88JgwZh11G2aIbRjzpQm+KzXN3HWKVboYzeU1zBVENXyhAIQkXQHocYS8V6Eko3JN/EeIpt4AunIQhBXSiPw2Q28cGY0g3JO4kIn+ELJmt7DEz2EIRKMBGuaopt4JVuR76FT3HOCsqgJQmgbkXTnd3HDXdRyAejNsZi8Pk9CAShUtyTmHC94YKw9yRXaPibMKBO6QG+Srci38L1vHsSDkkEoULck2gDBmHPSXTLg0oVTGHCp3kpJSndkLwK1/OuiUa/L0ogCJXirKSE9pgYSyndkPxJRkWhPYbZjUGdbE6rvcDGNUeVbkhehc/w7kkIQgShQsxWE+MjI02C0g3Jn/Ap3lVNY1FiUK3CGqbHSHdHE2wyyYu0x7hzy/WRLQh5nm9sbEylhr3EEUWxoaEhGjXWCdcIXAa7O9pzImLw+QxB5QqmMr0nDTSasLeed0+iMaiXkCsIX3zxRb/fv3jx4hkzZhw9enTwBvv37588efLSpUvLy8vfeOMNWXaqde6JdLjeQKefvScxXAlUrWAKEz7DS6JRugnZeg73RdNkCMKGhoaf/exnu3fvPnny5KpVq+6///7B26xZs+bhhx+uq6v78MMP165d29nZmf1+tc49mWHPCAbpnE/yYqwr4axAByGol5W2UBPsxumwQKVMHxmC8M0337zsssvOOussgiDWrl378ccft7W19d/g6NGj33zzzY9+9COCIObNmzd79uxt27Zlv1+tszEWm8vCtxliDG/vSc41iTZZcBcGVK1gqtMg3YSpRIoPxlwoXiMIQpbVJ+rr66dNm5Z+7PF4XC5XQ0NDWVlZ/w38fj9Nf3vqUVNTc+bMmeFeLZVKHTp0KBgMpn+cPHlyUVHRcBtLYjLR3CiatVryw5Qmu79qsCW12n6R5+N0RmeUXV8nnSWmeNOJXDfJUDJ//yFDTGEq+JVYNr0nk401/f6zrRJZLCXbTindkDEzmS02/2R5X1OGIGRZtqTkP0t5MQzT29s7YAOKokbYoL9YLPbTn/7UZvt2iaKnnnrq0ksvHW5jvr019ZcXxt1y5fFnd9WXWI7sVLod45RKpYTMzkJ6Om8sdX3WeSKY6yYZSubvP2QoJTki7bd1/Pl/m4jR+yw0/f73cOdaxILOP+9SuiFjZqKczA//iyCISCSjyiaapi2WUZYBlyEIPR5Pd3d334/d3d39LwfTG/T09PTfYObMmcO9GkVRH3zwgd/vz2TXVmuV/eHfjPqPVC1Xo3Bia4vv4RuVbsg4ZbgenhhPndlwrOrR/8LYCXlhPcJcaHv2RMEP/ofxj75SmKbf//Afm7wzXJ55tyjdkKzI9f7LcDoze/bsffv2pR8fOXLEbDZPmTKl/wYzZ85sb29vaWlJ/7hv377Zs2dnv18dYPxkNBQX4zofVs828IyfRAqCJhhk1ie2QXBVoYPwWzIE4Q033NDY2Pjcc88dPHjwgQceuPPOOxmGIQji0UcffeaZZwiCKCsru/7669etW3fo0KEnnnhCkqTly5dnv18dMFlMtNfBteh8bCVbj9krQDNcE+nwGZ0HYZIXk7xIlWIo/bdkCEKGYT766KPdu3evWbNm3rx5v/rVr9LPV1RU+Hy+9OOXX365qqrqxz/+cV1d3QcffGC1Gn2y8z7OKppt1PlRFz6DKm3QDPdEmtV7ELKNgrOCwlD6PvIE0jnnnPPuu+8OeHLdunV9j91u9wsvaLmqJWdcVZTOVwSVCLZRmLYKQQjaQJU6xHgq1pNwFNqUbkuuRBp5J+6L9qPVkifdcFZSbKOeB/BywaiNsWLBM9AMAywXyjYJGEHYH4JQYbTHkeTEJCcq3ZBcYc/wbtwXBU3Rfb1MpEnAFWF/CEKlmQhnBcnqd1ancD3vQqUMaIpb1/Uyse4EIRE6vvE7DghC5TkraR1PbxjGFSFojbOSEoKxlE7HNbG4HBwEQag8ZxWl18LRJCcmIyJdhipt0BKzzUx7HRGdjmuKoINwEASh8pwVpF6HErJNgrMSVdqgPc4qSq/3aSLpoxL6QRAqjyyyi/FUgk0q3RD5sajSBm1yVdJ6LefmWgTGj6PyOxCEKmAinH5Sl/dhIo2Yxgk0Sa8dFrHuBGE22d0YzvQdCEJVYPwU16LD089IM27CgCbpdVxTpFnA+tiDIQhVweknI816uyJElTZomIlg/GSkWW+np1xL1Fkx+sIaRoMgVAWmgoro7oqQbRRc1Rg4AVrlqtJhN2GkRXCig3AQBKEq0B5Hgk0mo7q6DxNpQqUMaJguuwkjzVEGV4SDIAjVwUQwPr0NomBRKQNa5qqiIvq6IkxEkql4iiyyK90Q1UEQqgWjs9GEEsE1R3ETBrTLUWgjTESsJ6F0Q2QTaY46K0iM6x0MQagWTj+lp555vj1mc1msjEXphgCMn7NSVxeFXIvAoGR0KAhCtWAqKD0VjkaaUKUNmuespPQ0IX6kJeosRwfhEBCEasF4HdGueCqhk3l+I8049wTNc1boaoAvrgiHgyBUC5PFRJXY+WBM6YbIA8OVQAecflI3M46K8VS8N0mVolJmCAhCFWHKSa5VF3dHJYJrQaUMaJ69wGYym/RRL8O3Rimvw2RGqcwQEIQqopsgFDpiVtpipVEpA5rHVFD6KOfmWtFBOCwEoYow5SQX0MUh14JBu6ATTr1MtMa1RmkfjsqhIQhVRDdXhBGMIAS9YPQyrolrjTK4IhwGglBFbE6r2aKHDolIs4BKGdAHnaybLRFcEEE4LAShujB+PVwUYuVP0A2yWA/rZke74lbKYqXQbT80BKG6MD7NByFW/gRd0cW62VxrlEEH4fAQhOqig25CrPwJOqODdbPRQTgyBKG60NoPQgylB53RwbrZCMKRIQjVhfY4Yj0JTU+0FkEHIeiLDtbN5gIIwpEgCNXFZDFRpdqeaI1rjTr9OORAP6hSeyKcFKNaPT0VY6lEOEmWYHK1YSEIVUfT9TJJQUzyIlmMQw70w2Q2UV4HH9TqUckHolQZJlcbCYJQdRgfqd1D7tviNBxxoC+aPj3FfdFRIQhVh/aRXECrt0bRJw+6xJSTEc0GIR+M0V6H0q1QNQSh6tBeB6/ZGUe51iiDDkLQHU2Pa+ICGEQ4CgSh6jgKbSlRSkQ0OZMF14IrQtAhppzkA1FCUrod48IHY5hue2QIQjWivaQWC0ellCS0x2gvDjnQGytlsdKWaGdc6YaMWZxNEpJkd2Gmp5EgCNWI8Tq0uB6TEIrbC6wWBz5UoEMa7SbkA1h9aXT4zlIj2qfJK0JUyoCOMX6K12IQBmPoIBwVglCNGJ8m62W41ihTjjllQJ80Wi/DBaIoGR0VglCNaB/JabBnHleEoGNObQYhbo1mAkGoRlbKYiEtmluhl2sREISgV+QEe4JLJgVR6YaMhUTwbTG6DFeEo0AQqhTj01i9TJITUwnJUWhTuiEAuWEiaB+prT6LaFfcSmM93tEhCFVKc4ccl74Dg8nVQL8YH8lpqooNQ+kzhCBUKcarxUMOd2BAzxifQ1uFoxhKnyEEoUrRWisc5YMYSg86R3tJbXVY8EGUjGYEQahStMchhOJSSjOVo1wrKmVA55hykg/ENFTOzQdiDE5PM4AgVCmz3Wx3WzUzpVO6OA3nnqBrVtpitpu0Us4tpSShM055sDjo6BCE6kV7HVqZXybaFbdSKE4D/WPKNXN3NNoRdxRYzTZ8yY8O75F6aWjqbRSngUFoqJybC6DbPlMIQvVivA6tLFXPB1CcBobAaGfdbFTKZA5BqF5auyLEIQf6p6ErQhRyZw5BqF5UmUPo1EbhKOYzBIOg00elqIWjEleEGUMQqpfZanIU2ISQ2gtHU0kp2p2gPTjkQP/MVpOj0Ca0q/1WjSRK0a4EVYqjMiMIQlXTROGo0BYji+0mC2ZXA0PQRDehEIo7imxmK47KjCAIVY3WQr0MF0QHIRgI7dPAUckHowzui2YMQahqmqiXQckoGAqjhYnWOMwyOhYIQlXTxBUh+uTBUGifBjoscFSOCYJQ1WiPI9qVUHmJGkbTg6GQE+wJNinGUko3ZCQYOzEmCEJVM1lMZLGqS9TEWCrJi2Qx5jMEozCZTZTHwbep96hMJaVYd4IqwVGZKQSh2tFlqj7k+ECUKnNgPV4wFNqr6lXShPYYOcGGQu7MIQjVjvKSag5CLoh1XsBwVF7FxrfF6DIclWOAIFQ7ukzVPfPokwcDYnwOvk29V4RYE22sEIRqp/LCUfTJgwHRXlWPqcfp6VghCNVO5YWjPKbbBuNxFNpS8ZQoqLRwlA/i1ujYIAjVzmQxOYpsQocaZxxN8qlUUrIX2JRuCEB+mQi6zBFtV+NRKYlSrDtBlaJkdAwQhBqg2ruj0bY4Zq8AY6J9pNCeULoVQ+DbY2QxSkbHBkGoAaqtl4mFEpjPEIyJ9jpiITUGoYBKmbFDEGoAXabSERRCewKVMmBMtJeMtqkxCFG/Ng4IQg1Q9a1RnHuCITE+h6DKPkK+LUqX4agcG9mCUBCE/fv3NzY2DrdBR0fHgQMHmpub5dqjcVAeR7RTjYWj0VACQQjGZHNaTSZTPJxUuiED8cEYgnCs5AnCAwcOTJky5YEHHjj//PMfe+yxwRvcfPPN06ZNW7t27Zw5c6699tp4XI1nUqr17aLYKiscjfcmzBaTzWlVuiEAyiA9NrXdqkkvTE9iYfoxkicIH3nkkfvvv/8f//jHl19++fLLL3/zzTcDNli9enUwGPz8889Pnz595MiR1157TZb9GocKl6rnAjGyDCXaYFyOUpvajkqhPUYWY2H6MZMhCNva2j799NO77rqLIAi/379s2bK33357wDZXXHGF3W4nCMLlcs2aNautrS37/RqKCrsJ+baoowSXg2BcVJlNbSv08m24LzoeMnyRNTU1uVyukpKS9I+TJk0aoafw5MmTn3766ZNPPjncBqIo7tq1q+/Vzj33XI/Hk30jtY4uIzuPhJVuxXfwgRjpxRUhGBfpsfceVtlRiZLRcck0CJ9++unDhw8PeHLevHkPPfQQz/Ppq700kiQ5jhvyRbq6uq677rpHHnmktrZ2uB3F4/Hf/OY3Dse3JzWPPvroRRddNNzGgiDY7XaLxZLhv0LD3GKkVYhEIkq34z/YVr64hlRVk4yG4ziTCTfBFJNyJrhWIcJG1LMMWbglUjiLMchRmeHnn6Zps3mUe5+ZBuHChQunTJky4Em/308QRFlZWU9PTyqVSu+ss7PT6/UOfoXe3t5ly5Zdfvnl69evH2FHFEW9/fbb6VcelcViMUgQ0hOlE91BhmLUMmGERMRCicJqj9PpVLopxiVJEt5/BUmSZKXCtqTDUaSWWQbjHYHi6gLaaYiLQhk//5kG4SWXXDLcryZNmlRYWLh37970pds///nPBx98cMA2HMd973vfq62tfe6558bdViPrKxxVSQdAtCtupS0WB8ahgqExPgcfjKokCFEyOm4yfJHZ7fa1a9euW7fub3/722OPPRYKha6//nqCIHbv3l1VVZXe5vrrr29sbJw7d+6mTZteffXVzz77LPv9Go2qCkfRFQFAEARdpqL1mFAyOm7yVP1t2LBhwoQJL730kt/v//TTT0mSJAiirKzspptuSm+waNGi2bNn19fXp3/sq4WBzNFl6bVA3Uo3hCAIgg9i9SUAgvY6ek8NXRKRfygZHTeTJKlrvpKKiorPP/88wz5CAxXLEEToy97OI+EZt1cq3RCCIIjj/6e5aLqTmm5xuVxKt8W4WJbF+68glmWJbuupd1prHxpYP6GIxr+1EwRRtcwoZfYyfv7Rx6MZ6ro1GsAS2AAE7XXw7TEppYrLCb4NR+U4IQg1g/I4op1xNcw4KqUkIRSnPDjkwOgsdrPdaY11qWIZCswyOm4IQs0wW9WyVH20I24vsFrs+PAAELTPoYb5ZdIlozg9HR98l2mJSlbo5QIxBgvTAxAEQRC0l1TDUSlgYfosIAi15N+Fowrjg1EaQQhAEARBMF5VXBGiZDQbCEItUcm5Jx+IMuiTByAIgiBoH6mGCfExtDcbCEItUckaFBwOOYB/S6+bnUoqXMWGktFsIAi1RA1L1aeSUqw7QZVi3QkAgiAIs9VEFtuEkMK3argAbo2OH4JQS74tHA0pWTjKB6NUqR198gB9aB/JK9pN+O3pKUpGxwtBqDHpSX4VbAC6IgAGUHyyC6E9Rk7A6en4IQg1hvaSnKKHHGYZBRiA8ZHKFo5ipqcsIQg1hvY6lL0JwwVwRQjwHYqXc3PBGAq5s4Eg1Bjaq3CtNh+I0rgiBOiHKrEnIkkxllKqARjamyUEocZQpfZYbzKVUOaQS0bFpCCSRSgZBejHRFClSnYTouc+SwhCjTGZTVSJnW9T5pDjAzHa6yDQJQ/wXbRyVWxiPBVnk+QEmyJ71wcEofYo2CHBB6OYZRRgMEa5Pgs+GKM9DpMZ56fjhyDUHsanWL0MKmUAhkT7SK5VsdNTlIxmCUGoPbSXVOrWKNcaZcoRhAADMeUk1yoosmt0EGYPQag9tHKz3WO4EsCQ7G4rYTLF2WT+d42hvdlDEGoPWWxP8mIyKuZ5v7HuhNlmsjmted4vgCYwPgffqsAZKjossocg1CATQZU5hLzXy3AB3BcFGBbtI7m8B2FSEMWY6ChEyWhWEISapMiUTlwrSkYBhqXkUYmK0ewgCDVJkUOOD2D2CoBh4fRUuxCEmsSUK1CrjVujACOgfQ4hFM/zcqF8IErjqMwaglCTmHKSb40SeTziUkkp2okFzwCGZbaZHYX5Xi6UC+CKUAYIQk2y0haz3RTrSeRtj3xbjCqxm63oiwAYFpPnehkpPYgQp6fZQhBqFVOe10OOb0UHIcAoGF9ex/hGO+M2xmqlLHnbo14hCLWKKc9rzzwXwKBdgFHQPjKf0x9yrVgTTR4IQq3K86AlVMoAjEqB01MclXJAEGoVk99zT9waBRgVWWxPCmJSyNOsTzwqZWSCINQquswR7UrkZ4XeeDiZSkmYvQJgFCaC8ZFcS57OUCM4PZUJglCrTBYTVZqnFXq5FsHpp/KwIwCtc1ZQkZZ8LEMhxlOJcJIqtedhX7qHINSwvM1kEWmJOv048QQYHVOepytCPhCjsB6vTBCEGkb7yPzMds+1oE8eICOMP09BiKH0MkIQapjTT0bycshFWgSmArdGAUbH+EihM55K5nzaJ65FYHCfRiYIQg1jKiiuJecTrYlRdEUAZMpkMVEl9jxUdEdaoghCuSAINczGWCwOc7Q7t3Mbcq0C7SPRFQGQIaefyvmtGongAwhC2SAItY3xk1xzbg85nHgCjAnjJ7kcF47y7TGby2olMbmaPBCE2ub057xWm2tFySjAGOShXgYjmuSFINS2fBxyzQKDQw4gY04/xQVy23kfaYk6K3B6KhsEobY5K6hIcw6vCCVR4kNxBuu8AGTMQpptTqsQyuFkFzg9lReCUNscRbaUKMXZZI5en2+LkUWGR7OFAAASn0lEQVQ2sx2fE4AxyPXQpkgreu7lhC84zXOWU1zOLgojTYKzEieeAGPjrMzhrZpoV9xsMdld1hy9vgEhCDXPWZHDc0+2UXBVIQgBxsZZSUUacxWEXEvUiQkuZIUg1DzGT+WuVjvSJDgr6Ry9OIBeOauoSIsgpXJSMBNpiWKmJ3khCDXPWUFGcjOUMJWUhPYYU45KGYCxsZIWu8sqhHIy2QXXLGBEk7wQhJpHlToSXDIXa4FyrVGq1G624UMCMGbOSjpHd0cjzQJujcoL33HaZyKcFTnpkIg08s4q3BcFGA9XFcU28bK/bKwnQUiEowirZMsJQagHrmqabZD/kEOlDMC45ahehm0QXNU4PZUZglAPXFUUm4srQoydABgvp5/kgjHZ12OKNPJOnJ7KDUGoB65qWvYgFGOpWE+CxpwyAONitpupEjsflLmQjW0UXNUIQpkhCPXA7raaraZol5wlapEmgSnH6ksA4yf73VEpJaFSJhcQhDrhrKLYBjkPORb3RQGy46qi2CY5j0o+GHMU2qwUVl+SGYJQJ1xVdKRRznqZSCOPShmAbDirZK5iYxsFdBDmAoJQJ1xyXxGG63n3JEbGFwQwGsbniPUkkpxsY3zZBh4lo7mAINQJZxXFBaKSKE+JmhCKEyYTxioBZMNkNrmq6PAZ2S4KIxjRlBsIQp2w2M1ksZ0LyFOiFq7nCqbgxBMgW+7JdLheniAUY6loZ5zxYXI1+SEI9cNVTbEynXuyZ3j3RNwXBciWexITPs3J8lJso8D4SZMFhdzyQxDqh3sK03tKniDsPc27J+OKECBb7mqKC0RlGVYfPsUVTMHpaU4gCPWjYArTe4ojsj7iEpyYCCcxlB4ge2a7mS5zyFLR3YsgzBkEoX44Cm0Wu5lvj2X5OuHTnGsSjaH0ALJwT2Ky7yZMJaVIs+CaiPs0OYEg1JWCKUz4VLYdEuF63j0JxxuAPNyT6N7T2QYh28DTXofFgW/snMDbqivuqXSvDEHIIQgB5OKeTLP1fJar1eO+aE4hCHXl227CLIixFB+MYawSgFxsTqvNbeUDWfVZhE/xCMLcQRDqCllsN1lMQmj8s2/3nuRc1TRWpQeQUdF0Z/fxyLj/XBIltpF34T5NzuD7Tm+yvCjsPs4WzXDK2B4AKJrh7DnOjvvP2UaB8jisJObazhUEod5kG4THIoXTEYQAcnJPYdgmQYylxvfnvae4gsm4L5pDCEK9KZzu7D7Gjq9nPtoVT8VTjBdzOAHIyWI3uyrHX8jW/Q2L09OckicIU6nUnj173n///Z6enhE2i8ViBw4c6OzslGWnMCRHoc1RYBvfShTdRyNFM1wEBhACyK1whrP72Hi6CZOcyAdjBVNxRZhDMgShKIpXX3313Xff/eqrr86YMePw4cPDbblhw4b58+fv2LEj+53CCIpnurq+GU+HRPcxdBAC5ETR9HF2E3YdZQumMWYrzk9zSIYg3LFjx8mTJ/ft2/fXv/717rvvfuKJJ4bcbN++fbt27ZozZ072e4SRFc9yd/0rPNa/kkQpfJrHiSdALjA+Mr18xFj/sOtf4eJZrlw0CfrIEITbtm279tprKYoiCGLVqlXbt29PJpMDtonH4/fcc89LL71ksaDwKedcVVSCE6NdYzvkwvU8VWq3Oa05ahWAoZmIoumusd4dlUSpp44rnoEgzC0ZvvWam5vnzZuXflxVVZVMJoPBYEVFRf9tfvGLX6xYsaK2tnbUV0smk++9915xcXH6xwsvvHDAS/UniqIoyrb6s54UzXB2HO71XVyc+Z+EDvYUznSO6f3E+68svP/KGuv7XziTCfyjy3NhQeZ/0nuCI0vtZtqE/+jBMnz/zWazyTTKjeWMgvDgwYMPPfTQ4Oe3bNlSXV0dj8dttm+XMk8/iMW+M4fCoUOHtm3b9sUXX2SyL1EUt2/fnr6+JAiirKystLR0uI1jsZgkSbjKHMxZ4wh9Hi4+P9P7nFJK6jgYnr6mfMD/3cji8fiYtgd54f1X1ljff3qyjd8aZds5e0GmVyAdh3vd0yn8Lw8pw/efJEl5gnDixIkbNmwY/Hw6orxeb0dHR/qZUChEEITP5+u/2S9/+Uu/3//kk08SBNHU1LR161aXy/X9739/yH05HI5Nmzb5/f5MGmYymex2O4JwMMc5ZOO7HXaTw0pl9OZ0H2WpEkeRfwznqgRBiKJI05jtQjF4/5U1jvd/wjkF3LF44WXujLaWiPDx5rPuqqJpjGgagoyf/4yCsKCgYNGiRcP9duHChe+999769esJgvjkk0/mzJkzoHGrV69uaGhIP3Y4HKWlpSUlJVm0GUZncZgLZzhDX/b6FmZ0dzT0Za/nvLGlIACMVel5hWe2B/2XZfQF2HuKszjMjA8pmHMy9BHefvvtzz777MMPPzxz5szHH3/8+eefTz+/cOHCa6+99pFHHrnyyiv7Nt68efOll156ySWXZL9fGJn3wqL694KZBGEqnur6hp30PW8eWgVgZIVTmTib5NtidNnoC18H93aXXViUh1aBDFWjRUVFe/futdvt+/fv37Jly0033ZR+/t577x0cePfdd9/cuXOz3ymMqrDGKcZTkabRR9Z3HmFd1ZTNhXpRgBwzEaW17o6vekfdMMmJ3UdZz7zCPDQK5Pnuq6qqevrppwc8eeuttw7e8vbbb5dljzA6E+G9oCi4t3tq5ShrKrXv7y49D8cbQD6Uzi089npT5RWlJvNIFRztB3qKZ7ky7OOHLGGuUT3zzC/q+Lp35Kl+I80C3xYrqUUHIUA+OCsoR7Et9OUoF4XBvV1e3BfNFwShntld1oKpTPv+kSaAbfxbe8WSUkzgBJA3VVd6mv7ePsLM+L0nOSlFuCdhmqc8QRDqXNUyT+OH7Qlu6GGnkWaBC0TL5uPEEyB/CqYw9oJhLwqllHT6r4HqZR5Mf583CEKdY3xkyeyChv/XNuRvcTkIoIjq5WWNHw59URj4R5eNsaK3Ip8QhPpXvdzTdSQcaRxYPtq2r5tvj+FyECD/3JNoqsTe8P7AM9Q4m2z6KDTl+74h/wpyBEGof1bKMvFqb92fm/vPfN/1Ddvwf9tm/bgal4MAiph+a2XXN2zLpx19zyQ5se6PTWUXFFGe0UcZgowwdMwQPPMKk4J48IXTk1d6HUX2SJPQ9FFo5o+qqVIcbwDKsNKWWWsmHvrN6aSQKphCExJx4q2WktkF1cs8SjfNcBCERlH+vya4J9En3mwxW02MnzprdZWrapTxhQCQU45C29lrJgb+2dX0YSjWm5hyXXnxTKy4pAAEoYE4K6g5j0xVuhUA8B+UxzH5WvQIKgx9hAAAYGgIQgAAMDQEIQAAGJq2g7Cjo4NlWaVbYVyBQCAajSrdCuNqamoSxaHnDIJcS6VSjY2NSrfCuGKxWGtrq1yvpu0gXL9+/bvvvqt0K4zrjjvuOHDggNKtMK5FixZ1dXUp3QqD6u3tvfjii5VuhXF9/fXXQy5wND7aDkIAAIAsIQgBAMDQEIQAAGBoJkkadk0sRVxzzTWHDh0ymzNKaJZlbTYbSZK5bhUMqaenh2EYm82mdEMMqrOzs6ioKMODBeQlSVJnZ2dJSYnSDTGoZDLJsmxR0ehrBuzYseOss84aeRvVBSHLsqFQSOlWAACAHlRUVNjt9pG3UV0QAgAA5BNuqgAAgKEhCAEAwNAQhAAAYGgIQgAAMDTNrEcYCoU2b94cCoWuuuqqxYsXD95AFMU33njj8OHDM2bMuPPOO1HTL6NQKLR9+/ajR48WFhbeeOONNTU1AzZIpVK///3v+348++yzL7roovy2Uc927txZV1eXfmy1WlevXj14m9bW1i1btvT09KxcuRJTf8lr06ZN/YsKB3+86+rqdu7c2ffjypUrPR6sMp+V9vb2AwcONDU1XXrppdOmTet7vrm5+bXXXguHw9ddd92FF144+A/j8fjmzZtPnDhRW1t76623Zji4SBtXhDzPL1iw4Pjx49XV1atWrXrrrbcGb3Pfffe99NJLNTU1f/zjH3/4wx/mvY16tm7dug8++MDn87W3t9fW1u7evXvABslkcs2aNceOHTt9+vTp06c7OzsVaadevf7662+99Vb6va2vrx+8QU9Pz/z581taWioqKq655pr3338//43Usfr6+tP/9uCDDx4+fHjABnv27Nm4cWPfNrFYTJF26smSJUueeuqpxx57bM+ePX1PdnR0nH/++aFQyOfzLV++/O9///vgP7zlllveeuutmpqaF1544aGHHsp0f5IWbNmy5fzzz0+lUpIk/elPfzr33HMHbBAIBBwOR1NTkyRJXV1dJEmeOnVKgYbqlCAIfY/XrFmzevXqARukj3yWZfPbLqO44447fv3rX4+wwfPPP7948eL041deeWXhwoV5aZfh7N+/n6Ko7u7uAc+//vrrK1asUKRJeiWKoiRJF1xwweuvv9735LPPPrt8+fL04xdeeKHvM9/n6NGjNE339vZKklRfX09RVEdHRya708YV4WeffbZ06VKTyUQQxNKlSw8dOtTd3d1/gz179kydOrWiooIgiKKiovPOO2/Xrl3KtFWP+s/dE41GnU7nkJv97ne/e/HFF7/88st8tctA9u7du3Hjxq1btyYSicG//eyzz6644or046VLl+7Zs2fIzSBLmzdvvuGGGwoLCwf/qrm5eePGjVu2bMF8ILIY8pbmgM/5rl27UqnUgA3mz5/vdrsJgpg4cWJlZeXnn3+e0e6ybnA+BAKB0tLS9OMJEyZYrdZAIDDcBgRBlJWVybhUFfTZs2fPtm3bfvKTnwz+1ZIlSzo7O48ePbp48eKNGzfmv206VlVVVVpa2tXV9fTTT8+fP5/n+QEb9P/8ezyeVCoVDAbz3kydEwThzTffHLKDtqCg4Oyzz+7t7d22bduMGTMG3zsFWQz4nCcSiY6Ojv4bBIPB/kHg8XgyDAJtFMtYrdZkMpl+LIqiKIoDpsyxWq39VyhNJBKjzqkDY3X8+PEbbrhh06ZNU6dOHfAru93+0UcfpR/fcsstl19++dq1axmGyXsb9empp57qe3Deeedt2bJl3bp1/Tfof4CkH+DzL7u//OUvRUVFl1xyyeBfrVy5cuXKlenH995773//93+/8847+W2dIYz6OR93EGjjitDv9/cFe/qBz+cbsEFLS0vfjy0tLeXl5flsoe7V1dUtWbLkmWeeufHGG0fecsGCBclksrm5OT8NMxSbzTZ//vzB9TL9D5CWlhabzYbJoGW3ZcuWu+66K91BM4KLLrro9OnT+WmS0Qz4nNM0PeA29biDQBtBuGLFiu3bt0ejUYIg3nnnncWLF6evNr7++uuGhgaCIBYtWtTR0ZFeLf3EiRPHjh3ru5UM2WtoaLjyyis3bNgwYEnoL774Iv2xEwSh78kdO3bQND1x4sQ8N1LH+t5elmV37tw5a9YsgiDi8fjHH3+cLlNasWLFtm3b0v2CW7duveqqqywWi4IN1p/6+vpdu3bddtttfc/wPP/xxx+nr0v6/oMkSXr//ffPPvtsZVqpdytWrHj33XfT13xbt25dsWJF+vkvvvgiHZBXXnnlwYMH02eKe/fu5Xl+4cKFGb20LBU+uZZMJq+44oq5c+fedtttEyZM2L17d/r5xYsX//znP08/fv75530+3+rVqysrK/ueBFksW7aMYZi5/3bPPfekn58zZ85vf/tbSZI2bdo0a9asH/zgB8uWLXO73X/4wx8Uba/eFBcXr1ixYtWqVT6f7+qrr04kEpIkpY/8M2fOSJIUjUYvvvjiBQsWrFq1qqSk5KuvvlK6yXqzfv36q666qv8zx44dIwiis7NTkqRly5YtWbLk1ltvPeecc2pqahoaGhRqpn48+OCDc+fOZRhm4sSJc+fOTX/n8zx/wQUXXHzxxTfffLPH4zly5Eh643PPPfeVV15JP3788cerq6tXr17t9XpffvnlDHenmdUnRFHcuXNnR0fHokWLvF5v+snjx4+7XK6+i98jR46kB9TPmTNHuZbqUF1dHcuyfT+6XK70ENd//etfpaWl6V7rAwcO1NfXFxQUzJs3D6OJ5XXmzJmvv/46Ho9PmzattrY2/WQikfjqq69qa2vTvSCJROKTTz7p6em57LLL+tcLgCzq6urcbnffNw9BENFo9NChQ3PnzrVYLF1dXfv27evu7vb7/QsWLMBsHtk7depUT09P3481NTXpWtD0jZBwOLxkyZIJEyakf3vkyJGysrK+j/2BAwfq6upmz549c+bMDHenmSAEAADIBW30EQIAAOQIghAAAAwNQQgAAIaGIAQAAENDEAIAgKEhCAEAwNAQhAAAYGgIQgAAMDQEIQAAGBqCEAAADA1BCAAAhvb/AeY7Si7IS/y3AAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"ρ = scfres.ρ[:, 1, 1, 1] # converged density, first spin component\n",
"ψ_fourier = scfres.ψ[1][:, 1]; # first k-point, all G components, first eigenvector\n",
"ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]\n",
"ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));\n",
"\n",
"using Plots\n",
"x = a * vec(first.(DFTK.r_vectors(basis)))\n",
"p = plot(x, real.(ψ), label=\"real(ψ)\")\n",
"plot!(p, x, imag.(ψ), label=\"imag(ψ)\")\n",
"plot!(p, x, ρ, label=\"ρ\")\n",
"plot!(p, x, tot_local_pot, label=\"tot local pot\")"
],
"metadata": {},
"execution_count": 10
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.0"
},
"kernelspec": {
"name": "julia-1.7",
"display_name": "Julia 1.7.0",
"language": "julia"
}
},
"nbformat": 4
}