{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Gross-Pitaevskii equation with external magnetic field"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We solve the 2D Gross-Pitaevskii equation with a magnetic field.\n",
"This is similar to the\n",
"previous example (Gross-Pitaevskii equation in one dimension),\n",
"but with an extra term for the magnetic field.\n",
"We reproduce here the results of https://arxiv.org/pdf/1611.02045.pdf Fig. 10"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using DFTK\n",
"using StaticArrays\n",
"using Plots"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"Unit cell. Having one of the lattice vectors as zero means a 2D system"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"a = 15\n",
"lattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"Confining scalar potential, and magnetic vector potential"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"pot(x, y, z) = ((x - a/2)^2 + (y - a/2)^2)/2\n",
"ω = .6\n",
"Apot(x, y, z) = ω * @SVector [y - a/2, -(x - a/2), 0]\n",
"Apot(X) = Apot(X...);"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"Parameters"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"Ecut = 20 # Increase this for production\n",
"η = 500\n",
"C = η/2\n",
"α = 2\n",
"n_electrons = 1; # Increase this for fun"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"Collect all the terms, build and run the model"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter Function value Gradient norm \n",
" 0 3.068405e+01 7.473877e+00\n",
" * time: 0.004585981369018555\n",
" 1 2.875985e+01 6.225526e+00\n",
" * time: 0.013888835906982422\n",
" 2 2.485722e+01 9.941364e+00\n",
" * time: 0.07479500770568848\n",
" 3 1.415284e+01 2.280556e+00\n",
" * time: 0.10130596160888672\n",
" 4 1.305206e+01 2.413680e+00\n",
" * time: 0.11786484718322754\n",
" 5 1.204064e+01 1.846784e+00\n",
" * time: 0.1341559886932373\n",
" 6 1.057111e+01 1.433523e+00\n",
" * time: 0.15034890174865723\n",
" 7 1.013562e+01 1.747514e+00\n",
" * time: 0.1666250228881836\n",
" 8 9.733598e+00 9.628315e-01\n",
" * time: 0.18261504173278809\n",
" 9 9.366004e+00 9.722948e-01\n",
" * time: 0.19878888130187988\n",
" 10 9.185844e+00 8.222160e-01\n",
" * time: 0.21507596969604492\n",
" 11 9.049656e+00 6.046532e-01\n",
" * time: 0.23162388801574707\n",
" 12 8.979126e+00 5.376082e-01\n",
" * time: 0.24439287185668945\n",
" 13 8.916894e+00 3.240047e-01\n",
" * time: 0.2572929859161377\n",
" 14 8.865773e+00 3.135228e-01\n",
" * time: 0.2699768543243408\n",
" 15 8.838303e+00 2.740715e-01\n",
" * time: 0.2827279567718506\n",
" 16 8.814572e+00 2.353091e-01\n",
" * time: 0.29560089111328125\n",
" 17 8.797703e+00 2.146356e-01\n",
" * time: 0.3090198040008545\n",
" 18 8.790363e+00 1.632228e-01\n",
" * time: 0.3230249881744385\n",
" 19 8.785971e+00 1.080686e-01\n",
" * time: 0.3415999412536621\n",
" 20 8.779421e+00 1.123960e-01\n",
" * time: 0.35556793212890625\n",
" 21 8.771060e+00 2.108575e-01\n",
" * time: 0.3981599807739258\n",
" 22 8.758655e+00 1.775920e-01\n",
" * time: 0.4103658199310303\n",
" 23 8.747922e+00 1.157848e-01\n",
" * time: 0.4228689670562744\n",
" 24 8.735710e+00 2.110615e-01\n",
" * time: 0.43482089042663574\n",
" 25 8.724327e+00 1.405026e-01\n",
" * time: 0.44718289375305176\n",
" 26 8.716358e+00 1.305312e-01\n",
" * time: 0.45949387550354004\n",
" 27 8.704812e+00 1.971143e-01\n",
" * time: 0.4719099998474121\n",
" 28 8.691239e+00 1.931279e-01\n",
" * time: 0.4840669631958008\n",
" 29 8.679364e+00 1.789381e-01\n",
" * time: 0.49630284309387207\n",
" 30 8.668966e+00 2.267373e-01\n",
" * time: 0.5086309909820557\n",
" 31 8.660125e+00 1.535688e-01\n",
" * time: 0.5208728313446045\n",
" 32 8.646998e+00 1.397988e-01\n",
" * time: 0.5330648422241211\n",
" 33 8.633404e+00 1.240944e-01\n",
" * time: 0.5452280044555664\n",
" 34 8.617096e+00 1.758594e-01\n",
" * time: 0.5576088428497314\n",
" 35 8.605580e+00 2.437573e-01\n",
" * time: 0.5703840255737305\n",
" 36 8.594514e+00 1.471578e-01\n",
" * time: 0.583115816116333\n",
" 37 8.580812e+00 2.197228e-01\n",
" * time: 0.595573902130127\n",
" 38 8.560864e+00 1.877918e-01\n",
" * time: 0.6078798770904541\n",
" 39 8.534669e+00 2.329562e-01\n",
" * time: 0.6201050281524658\n",
" 40 8.511610e+00 1.859549e-01\n",
" * time: 0.6325209140777588\n",
" 41 8.492595e+00 2.406591e-01\n",
" * time: 0.6457858085632324\n",
" 42 8.491859e+00 5.175557e-01\n",
" * time: 0.6588449478149414\n",
" 43 8.480363e+00 4.049774e-01\n",
" * time: 0.6721110343933105\n",
" 44 8.474581e+00 4.059790e-01\n",
" * time: 0.6853690147399902\n",
" 45 8.468200e+00 3.218032e-01\n",
" * time: 0.723822832107544\n",
" 46 8.467268e+00 1.737473e-01\n",
" * time: 0.7403459548950195\n",
" 47 8.466971e+00 1.586764e-01\n",
" * time: 0.7605698108673096\n",
" 48 8.466923e+00 1.188528e-01\n",
" * time: 0.7767958641052246\n",
" 49 8.466923e+00 1.188528e-01\n",
" * time: 0.895845890045166\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVzVVf4/8HPvZVVRUVQWcQNMcEHHJXcN9wVRS8txK9RyXMocpxytb2nWmJmVZmqLy6/FJI2mbXLNpcyU3BUTRUURUtkUhLv//vjM3GH0/UbO9XP5cLmv58M/4Hj4cO4Ch/P5vD7vo7Pb7QIAAMBT6bUeAAAAgJYwEQIAgEfDRAgAAB4NEyEAAHg0TIQAAODRMBECAIBHw0QIAAAeDRMhAAB4NEyEAADg0TARAgCAR/OqmG9js9kOHz8ZFR1dMd+ugtlsNr2+6v9JgYdZleBhViUV9jBr+Xrfs8/vv+cWF1ukDnvrVkFoqD4iIsLZcd2vCpoIjUbje5uShz5er2K+HQAAqG7kAyH37PPoo98cO3Zd8sAnhw8vSE5Odm5U96/q/60EAABQhgpaEQIAgCfQ6YROp/5hMzIyCgoKYmJiDAYD2eHmzZvnzp2LiooKCAhwNBYUFFy6dKlBgwYNGjRwNObn5zu2XfLx8alevTpWhAAAUHnZ7fbExMQHH3xwwoQJMTExGRkZd/fZuHFj06ZNp0+f3rRpU8cp1okTJzZu3HjChAnR0dGjR482m81Ke2RkZNOmTSMiIiIiImbOnClwahQAAFSkkydEWUvIbdu27dy5MzU19ciRI3FxcS+//PIdHW7fvj19+vQtW7b88ssvH3/88bRp00wmkxBi9OjRWVlZR48eTU9PP3z48Lp16xxfcujQodzc3Nzc3LVr1wpMhAAAoDKd/D/epk2bRo0aVbt2bSHE5MmTk5KS7thPftu2bUFBQb179xZCDBo0yMvLa8+ePUKIIUOG+Pv7CyFq164dGxublZXl+JLCwsLr1/+b6MFECAAAldelS5eaNWumfBwREVFUVHTjxo3SHTIyMhwdhBDNmjW7dOlS6Q4XL17ctWvX0KFDHS1Dhw5t0aJFRETErl27BMIyAACgJqfCMpmZmV988UXplujo6FatWgkhbt++7efnpzQqHxQWFtar99+b8YqKinx9fR2f+vv7FxYWOj4tKCh4+OGHn3766fbt2ystR44cCQ8Pt9vt77zzzqhRo9LT0zERAgCAav5z2U+GXXf16tU7JsL+/fsrE2GDBg1yc3OVxpycHCFEcHBw6Z4NGjTIy8tzfJqTk+PocOvWrUGDBnXr1m3BggWODuHh4co4Z82a9corrxw/fhwTIQAAaKxjx45JSUnkf7Vv337//v3Kx/v374+JiVGu/JXu8PTTT5eUlPj5+d26devkyZPK4u/27dvDhg2Ljo5+++23ySMXFBQUFRXVqlULEyEAAFRekyZNiomJWb58+QMPPPD3v//9+eefV9rj4+OHDx8+adKk2NjYTp06TZ48ecqUKStXruzTp09UVJQQ4tFHHz1//vyjjz764YcfCiFatGjRs2fPX3/9devWre3btzcajW+99daDDz7YqlUrTIQAAKAaJ26ot5fZPzQ0dOfOncuWLdu2bdvcuXMnT56stPfq1cuRkdm8efOrr7762muvtW3bdt68eUpj69atQ0NDjxw5onyqrCODg4NzcnJWrVrl6+sbHx8/ffp0vV6vuyOH6iLFxcXTFy4Z+viTFfC9AADAFcpTa7Rjp0+OH5erNWqznRg6JA+1RgEAALSBU6MAAKCee90gXwlhIgQAANXohPQ1QlcU6ZaCiRAAAFTjzH2EWs+EuEYIAAAeDStCAABQjxPXCHFqFAAAqhJ3mwclJ8KjR4+mpaXVq1evc+fOjiqox44dO3v2bExMTMuWLV0wQgAAABcq7zVCq9U6fvz4YcOGbd68+eWXX/7kk0+U9n/84x9Dhgz57rvv+vbtu2LFCpeNEwAA3IBOJ783r9ZrwvKuCFetWnXixIlTp04FBAQIIZR6NDk5OYsWLUpJSYmOjk5JSenbt+8TTzxRo0YNF44XAAAqMze8RljeFeEnn3zyzDPP5OTkHDhwoKioSEnHbt++PSoqKjo6WgjRoUOHunXrKvsCAwCAZ3JiQSh9u4XayrsiPHfu3Oeff7569eqAgIAzZ858++23bdu2zczMVDZ2UoSHh1+5coU7QsUUNQUAABex2+2aT1quUN4V4e3btwMCAg4cOLBjx46JEyfOmTNHCGE2mw0Gg6OPl5eX2Wx2yTABAEBr5VnPKLtPSP3TXHknwtDQ0D59+ih/C/Tr1+/YsWNCiJCQkOvX/1tl/Nq1ayEhbG3yKvl3BACA59DryzFl6Jz6p6nyToS9e/e+cOGC8vGFCxdCQ0OFEN27dz98+HBeXp4QIisr6+zZs126dHHRQAEAAFyhvNcIZ8+e3atXr9q1a9esWfPVV1996623hBAREREJCQnDhw8fO3bsunXrxo8fr0yQAADgmZTznZJforHyrghjYmL27dt3+/btq1evJicnP/bYY0r7xx9/PGbMmFOnTk2ZMmX16tUuGycAALgD+WuEml83k6gs06JFi0WLFt3R6O3tPXXqVFWHBAAAbqsSXPOThd0nAADAo6HoNgAAqEa5oV7ya9zkhnoAAIB7cuKan9bzIE6NAgCAZ8OKEAAA1OOGYRlMhAAAoBonimhrXncMEyEAAKjGDReEuEYIAACeDStCACGEEExV/YrfPIz9a9rt/swGz+SGG/NiIgQAAPU4cY3QRSMpN5waBQAAj4YVIQAAqEYnf4O81qFRTIQAAKAenU7+dgitZ0JMhOA27HLBFbn0S8WHYjjcSHTsf0i0ctQ4BoAQwi3DMrhGCAAAHg0rQgAAUI0bLggxEQIAgHqcKbGm9VSIU6MAAODRsCIEAAD1uOG5UUyEUCHUiHCyqVGq3c705g7CR1LJo3OdGapEOJmD0O1SnSVHwvbW+tcZVAZO3Eeo+USIU6MAAODRsCIEAAD1uGGtUUyEAACgHlwjBAAAT6bTuV+tUVwjBAAAj4YVIThLKk0pk+HkAp82G3MQm1Rn2TQp8R9SCVOOXApUsIWMdXqiXU81CiF0zJ++5EFkRyh1hkvzRQC4iE7oZG+Q1/yGekyEAACgHieuEWoNp0YBAMCjYUUIAACqcSIso/kKEhMhAACoRifk7yPU+ooxJkK4B76wGZUiYTpz0RWy3WalO1utVCpGCKtF4iA25iDMAJlwDVsaTmLzXNlQjN7AtFM5FwPX2YtuNxiIqyRS31EwSRymr7DLbgas9aIBygvXCAEAANwLVoQAAKAenfSpTq3PjGIiBAAA9bjhmVGcGgUAAM+GFSEAAKgHRbfBfUnVGBNMnJIPfNLtFjOR4SQbhRAWs7X8ByGjpIKPnnIjJx8m+1zRzfSPOZsaZXKWXIaTDHwavOiTPV7eEu1cZwMXPaW+qZ17ONzZKK4InJ06jta/QOFuOme2YcLtEwAAUFU4sUO95mEZXCMEAACPhhUhAACoB9cIAQDAkzlxjVDzc6M4NQoAAB4NK8IqjS6TKbG9reDLhFotVOCTahRCWEx0u8lIBEHNVKMQwmRiUqPUwdnoKTNCG5MyJR8+l6Rlw29k4FE2Ncq0k1lNLvDp7cO0+xrubvShGoUQ3j50u5cP8bR4MflVLteqZ55DvY58zpmCpVqfavNkTuw+cc/un3/++WuvvZaXlzdixIilS5f6+Pjc0eHMmTMzZsw4ffp0mzZtVq5cGRERIYR4//33P/300/Pnzzdo0GDmzJmPP/640jkjI2PatGmHDx9u3rz5ihUrWrdujRUhAABUXqmpqU899dQ777xz8ODBw4cPL168+I4Odrv94Ycf7tWr18mTJ9u3b//oo48q7RcvXvz73/9+6NChRYsWPf300z/88IPSPm7cuKioqJMnT8bHxyckJFitVkyEAACgHp0TyjreRx99NGLEiIceeigkJOTFF1/84IMP7ujw008/Xbt2bd68eXXq1HnppZfS0tJ+++03IcRrr702cODAkJCQQYMGDRw48OeffxZCpKamHjp06JVXXqlTp87s2bMtFsv27dtVmAhNJtOFCxfMZvP9HwoAAKC01NTUdu3aKR+3a9fuypUrN2/eLN1BOSNqMBiEED4+PjExMadPny7doaSk5NChQ23btlWOFhUVVaNGDSGETqdr27bt6dOnyzsRzpkzp/T0ffv2baX966+/DgsLGzRoUHh4+I4dO+7v8QIAgHtTrhFK/RM6YTKZ8v6XxWJRDnjjxo2aNWsqH9eqVUtpKf0dc3JyAgICHJ/Wrl27dAe73T5jxoyoqKiRI0dynSXCMvPnz1+0aFHpFqPROGnSpA0bNgwePDgpKSkxMfHChQvKtAwVSmafWC4UI7XtrRDCTEVXuJyLsYRpL7bc3WhiOpPJGq7dwiVrmBANVwSOLrFGdpXeVFYu6MGVWCPTKFKhGCGEjx/xq8DXj+vMtPsSB+E6ezOBIy8789c5uXWwnnnVEKLRlvzzvH37diXh4jB58uQlS5YIIerUqXPr1i2lUVkL1qlTp3TPwMDAwsJCx6cFBQWlOzz33HPHjh3bsWOHcgb27s5t2rSROzVaXFxc+tNt27bVrFlz8ODBQohHHnnEZDLt3btX6oAAAABDhgzJ/V/KLCiEiIyMTE1NVT5OTU0NCgqqXbt26a+NiIhITU1VstxWqzUtLS0yMlL5rxdeeGHHjh1bt25VlpJK5/T0dKPR6DhgZGSkxES4dOnSevXqBQcHv/nmm0rLhQsXmjdvrnys1+sjIyMvXrwo/wwAAEAV4UxWpsyl+sSJE5OSks6ePWsymd54442JEycq7YsXL965c6cQok+fPl5eXuvWrRNCrF69uk6dOl27dhVCvPLKK+vWrVu9enV+fn56evr169eFEO3atYuIiHjnnXfsdvumTZsKCgoGDx5c3olw6tSpeXl5hYWFX3311auvvrplyxYhxK1bt/z9/R19qlWrdsc1zNJsNuaUHAAAuAOrlb7iUJoT1wjLPmXdoUOHF198sWvXrkFBQT4+Pi+99JLS/uuvv2ZkZAghDAbDpk2blixZUqNGjVWrVm3cuFGZWXfs2OHn5/fnP/+5X79+/fr1W7ZsmfKFn3zyycaNGwMCAubPn//FF1/4+PiU9xqhY6XZuXPn8ePH//DDDw8//HD9+vXz8vIcffLy8ho0aMAdQc9uuwIAAG6gPBEQZWqTPPA9+s+aNWvWrFk2m630PJKcnOz4uHPnzmfOnLmjw549e8ijtWrV6siRI6U7OzM55eXlVa9eXQgRGxt75MgRk8kkhCgqKjp16lRsbKwTBwQAACjbPVdTUsut0p3LuyJ86aWXevToERAQsHv37i+++OKnn34SQnTq1Kl58+bPPvvsU0899fbbb3fp0iU6Orr84wBpMulQIQR5NtrGVUFj2qWCoGQKVAhRcptpp/pzB2HTpFQ7GWoVTqRGuU14KdyftTLHYI9i4FKj3sQf6VxqlMtw+vpRz6E//fvB10gfxFKN2iHZSh/ExrTb/chm7leV3K6/ZJoUUVKVObH7hNbKO3/a7fbFixfPmjXr+PHju3fvbt++vdL+1VdfFRcXT5482dvbe9OmTS4bJwAAuAO1rxFWgPKuCBcuXEi2h4aGrl27Vr3xAAAAVCjsPgEAAKrRCSG7H6HbrAgBAADuzQ2vEWIidCfSWwlS+RczkxaRrY5Wcpsos15SROdciovomuxkiIYLy0jVaZMNy3DbLnLPLd1ZJhXD9ZUtseZNhmV8mbAMVQVNCGGkwjJ+XDqpGn0Q8q3Flegja9eJMp5DO/VNmUwQ9yTSmxraJXaRhHtyZj9CrZ9q3NsHAAAeDStCAABQzb03GKS+xjVjKS+sCAEAwKNhIgQAAI+GU6MAAKAaZ8IyrhlJ+WEirKTI6Jzsnrpk1TQ2HSpZHY0Mgt4upNOhxUya1EgdnKy7VsYIjSVUatQoV0nOxjyHNurZ4pONUs007pcIX2KN2piX24DXl371fYzErwIuYcuWqaOeWzaOK/kc0k+Ljv4Nxj2HOvIcGLO7r3zlaBBCKSwjfR8hrhECAABoBytCAABQjxM31Gu99sZECAAA6qkERbRlYSIEAADVOHEfIa4RAgAAaAkrQo2x0TnqP7gAHlfOkQyI8rVD6UwmFwQl24uZzuzGvDK1Rrk0qYlKjZq41KiZfvjcc0gmcqUKkAomxMj+Ecz8B5caNVCpUR+TbOCTer+xSVpmI2g660z2Za8KcYsDnZ7cU5fpzLYTo9G721YJoDpMhAAAoBqnTo26aCzlhVOjAADg0bAiBAAA1ThRWQa3TwAAQNWi9cQmCxOh1pi0DLltKVnCSvCb0JqoXAy5oa4Q4jazdy4blrlFhWVkNuDlBsOHZbiHSW7My6RFJDfsJVMkXHkwuaAHu6esXFiGLLFmY7I/XN6KbOero5HNcvsSsw+feq6EEHqqnWwU/CbGeuqb6pihcCOkf827269+19HppG+HwO0TAAAAWsKKEAAA1OPE7hNar6cxEQIAgHqcqDWqNZwaBQAAj4YVIQAAqEYndPJbOWq8hMREWEH4oB3dTlax4pKN3F67JXRqlNloVyYdyvXnUqPFTFS1uJDamJfpzBWHIx++7KaybGqULDPGlQ2TSTyym8cyB+FSo1YrlRqV3A6X7E9Gl8tqJw/OPEw28EkmO5mHb/BiniumnUyTciNho4wIjZbJmR3qtX76MBECAICq3GxBiGuEAADg2bAiBAAA1ThTdNtFQyk3TIQAAKAad6w1ilOjAADg0bAidAEyUieT1hNCWKnUqJlJNnJxSiMVEC0pYlKjbDsT+KTauYKlxUz0lCxkyqVGTUw4lgx88rVDJWqKCi5OKVFTUwhuY16ZiKlgaooKJmDMbR0s9X6TrTVK/l3PPRyDQSIdKpggKPeccO0GL6KdG4lOTz9OZoASEdMqzokb6rV+ljARAgCAepy5Roj7CAEAoKpwwwUhrhECAIBnw4oQAABU48TtE5qXlsFEqD6ZrAybRyCzHlwpNXIDXsHscFvCbXvLRFTYkmxUf7ZOm0y7kdpoV/D7D1vlwjJq7FhLRUtEGbXHqEY2RUIFOoQQ3j5cWIZoJ8MvQghvi0Tihqskp2d+bZEFzNjcijf9KvP5F+LV9/ah3xJcu5ePgToy/VbhYjt28oVDVqY0dzs3ilOjAADg0bAiBAAA1ThTdNs1Iyk/TIQAAKAanY7fuIP/Em1hIgQAAPW44f0TuEYIAAAeDSvC+8AVmqISouyGqEy6z0JVAuNikyYmZknGL8koaRntUmlSLmLKHdxkJNq5cKzFxOxLbJbYmJd7wtnUqMwOyVz9NpKeqezFl1Kj261WIgnpzeZaic6CeS9zJ7jIdKhgqqCZmBCsF9NuYtKk3tTTYvKjHw5Xjc+H+gmycnFc5oXQUz/L/JJG68VOhXPDBSEmQgAAUJEb3keIU6MAAODRsCIEAABVaX6uU5LcitBsNsfHx//5z392tKSmpsbFxYWEhAwcODA9PV3t4QEAgDtR7iOU/actuRXha6+9lpGRodf/e/q02+0jR44cN27cli1bli5d+uijjx46dMgFg6yk2KwM1chnMeh4BVNije5sYlIkZGSAr8emQjsXUiDzLEIIq5l4WrgqaGSASDAhGjYsI1NKjRuMhUktWZiICvljrjdI1mmTGrmdTpFwJErDCb6cGPV3NZcJ4irJGZiIClk1jYuJmf3ph28m3yrM+40NHFG/NbkaiuxveK1/9buOE7VGpa8pqk1iRXjixImvv/569uzZjpa9e/fm5OTMnTs3MDDw//7v/37//ffDhw+7YJAAAACuUt6J0GKxTJkyZcWKFd7e3o7G1NTU2NhYg8EghPD19Y2Ojj59+rRLhgkAAG5BJ/9Pa+WdCJcuXdqlS5euXbuWbszJyQkICHB8Wrt27Rs3bnBHsNkkbrECAIDKxmqlLwqU5swFQq3nwnJNhBcvXnz33XefeOKJ9PT0a9euGY3G9PR0u90eGBhYWFjo6FZQUFC3bl32O+lxqwYAgBtTzv9VPeUKy9y4cSM4ODgxMVEIkZub+8cff4wePXrfvn2RkZGpqal2u12n01kslrS0tIiICBcPGAAAKi9lmSf5JRor10TYoUOHlJQU5ePPPvvsjTfeUD6Ni4vT6XTr1q1LTEx877336tWr16VLFxcO1k2Q+TE7c2KYS0hKbTbLhhipjBwX7OTSpGwQlKxWJbkdrpU6YW6VSdIKJk3K7UwrVUpNMFXTuPwqd3CpUJxF0AfnqvSRv0W4XX+5kTBvWi42SiO/qZcXUzKNKWxmNnGBT+L9RqZAy2gn30JcYTzu1aSfFuYJt3vghr1uWGNN+nRljRo1wsLClI+9vLw2bdr0xhtvVK9e/cMPP9y4caPmKVgAAAAp0pVlhg0bNmzYMMenXbp0cZwdVXVgAADgjuTvI9R6SahOiTXMggAAIJzaof6e8+CFCxdWr15dUFAwfPjwgQMH3t2hpKRk5cqVv//+e8uWLf/yl7/4+PgIIQoLCw8dOnT8+PHAwMAJEyY4Oi9fvry4uFj5uEWLFgkJCUhyAgCAepy4j7DMiTA3N7dz5842m61jx44TJ07cvHnz3X3GjRv3ww8/9OrVKzk5efLkyUrj+vXrZ8+enZSU9O6775buvHDhwvT09Ly8vLy8vKKiIoGi2wAAUJmtX7++bdu2b7zxhhDCYDAsWbLkkUceKd3h7Nmz3333XXZ2dq1atfr37x8eHv7qq6+Gh4fPmDFjxowZH3/88YoVK+445pw5c6KiohyfYiK8D1yKj2pnS1nKbMzLhRW52BuZkWOzl5LtZECUDeBxGU4qTcpFT/ntcKknnCsUybxqUlFV7uBS5EKJQnDRNgPVn00pMw+TxP2ZbvBiUspmYoQWM50C5bZZZl998q3C/UTI7Jwsu1cz93MPCmdqjZa5JPz555/j4uKUj+Pi4hITE41Go6+vr6PDL7/80q5du1q1agkh6tWrFx0dfeDAgfDw8DKOuXr16po1a3bq1GnQoEEC+xECAICK1D4zKrKysoKCgpSP69WrZ7fbs7OzS3fIzs52dBBC1K9fPysrq4wDDhgwIDAw0GKxPPnkk1OmTBFYEQIAgOZ+/vnnfv36lW4ZMmTIrFmzhBC+vr4mk0lpVD7w8/Mr3dPX19dsNjs+NRqNd3S4w6effqp8MGnSpMjIyOeffx4TIQAAqMepG+qjoqKef/750m2Oa3gNGza8cuWK8vHly5d9fHzq1atXumdYWNjly5cdn165csVxs3vZmjRpUrdu3cuXL2MiBAAA1eicuKFOJ+rXr9+3b1/yP4cPH/7CCy+8+OKLfn5+n332WUJCglK5eu/evWFhYREREf37909MTDx27FhsbOyvv/6ak5Pz0EMPcd+qsLDQ399fKZq6a9eu/Pz8mJgYTIT3xmZiZPrLbsxLXtVnQzFcATPq4Nx35LYH4bIbZMxHKvsjJOM8bDKCekR2rmQa83C47AZZ2YuLouiYisTkpXjuINxf03rmjcjUUmMOLkVyV1nyEbEvvcyrKZj3oSo/m2ztOplQjOQexlWZE/cRlt09ISHhgw8+6NSpU+PGjVNSUnbu3Km0/+1vfxszZsysWbNq1aq1aNGiAQMG9OjRY+/eva+//nq1atWEEAcOHJgxY0ZOTs7169c7dOjQo0ePt956a9++fYmJiX/6059KSkoOHTq0YsWKBg0aYCIEAIDKy8vL6/vvvz948GBubm737t1r1qyptG/atElJigohnnnmmSFDhpw5c2bJkiVNmzZVGmNjY5OSkhzHUWbHQYMG7d27Ny0tzd/fv02bNsqOSZgIAQBAPS4ouq3X6zt37nxHY5MmTUp/GhkZGRkZWbrF39+/WbNmdx8tKiqq9E2EAhMhAACoyYn7CLUu0on7CAEAwKNhRQgAAKpxRdFtV8NEqD4yhCaVvRRc7THJEmt07THJwlFcpo58RGwUUCYIyj0nUvW0ZA8itZUr93Nu0NGnWOxUM/fEcuWmuNNHOmZLWKmD6A1Eu8HAnDFiviH56nNPLPdCcOgRetFDMVCdhRB6+gFp/WsYtIaJEAAAVKNzw/0IcY0QAAA8GlaEAACgGmduqNf65DQmQgAAUI0z2zBpPRPi1CgAAHg0rAjvg8S+vGxqVGrDXjIFKmSLcKq0rSg5cm6PU66AJFOwVCIEK5hgJ9uZTY3Sm82Sj4gLJfJvCYk6mUzyVO70EZv15aKqMgfnMsDkc+XjSxdg9atG//IJCPQl22vUItr9a3iTnX386IN7eRODIfOoQsiFSbU+t1eZOFFZRmuYCAEAQD24jxAAADyZTuhkb4fQeh7ENUIAAPBsWBECAIB6XLD7hKthInQBMizD9WXqTJFRFC5Fwm3MS9Zp42I77Gl9mfcov/8w3Z8cDJ+4kdoimNt/2IW7sPLl6GSOwpVMY/qTT5dVR/fW6yVCW+xzyLwQXt7EGabqNX3IzuGRtcn2kMYBZLs3Fbrh3rTePvS5LjK5w9Vp0zMvhNZR/8pOJ/8Uaf6M4tQoAAB4NKwIAQBANTqd/A3yWq+yMRECAIB6cI0QAAA8mTO1Rl0zkvLDNUIAAPBoWBFWFDaUKFEKi0vr2WQSkvwGvHS73L2xkrv72mRqj3EHIXHb1TKpSWHjIoJ0xTzmIMwLRLZyfzXrueQt9xzS35R+S3BJSKmDcPFdX3/i90lQcDWyc0SrOmQ7lxo1GYnkcXGhmezMvVPIgCi3/zC74zH5ymm+qKk0UHQbAADAzWAiBAAAj4ZTowAAoCYU3QYAAM/lzDVCFw2l3DARVgVcNIDeMlCNGmOCSwzIvqNlvikX9CATEHob3dnKBEBsTH871c6Nmn1uSXKV1PggEtkoGTiim5mDcC8EWdiMTNAIIXz86H0KDV709Ro9tQkit5UgV7mQzL/IbgCpdbCj0nPD+whxjRAAADwaVoQAAKAad7x9AhMhAACoxonKMprDqVEAAPBomAgBAMCj4dSoC5CnBdgEGrf5JxVv4zpzIyGrg3G5QeYYbHaODuDJ7WVKPiI+HcqUwqK6cwFOu50+iOGOiQ8AACAASURBVIGrjkblEtlydFwQlOrPhRW5h88lJPXUcfhNZSVeIP47SrQX3TSRna+m3yTbjcX0Js7kYLgnnNwiWAjhrSOiqnbut6BkjhoU7niNECtCAADwaFgRAgCAmtxtX15MhAAAoCI3vKG+vBPhtm3bvvrqq+zs7KCgoAkTJnTv3l1pLygoWLx48e+//966devnnnuuevXqLhsqAABUdk5cI9Rcea8Rpqent2nT5oknnmjevPnAgQP37duntI8aNSotLe3JJ5/87bffJk6c6LJxAgAAuER5V4RTp05VPoiPjz906ND27dt79Ohx4sSJ/fv3X7t2rVq1ag8++GBISEh6enqzZs1cNlr3QIdGuRSfTETQwKX4mHYyl8httMv9CafSCJkCkgaiKCRZO1QIIex0u4X+e47ZmZaucCkMTM6QfF5kdghW+hNfwP3VzL+aEllN2fqZzEjoo3CZTHKEtwqMZOcrTGq0iNlrt3pNn7sbq9XwJjv7V6fbyTe/F5NGRmjUOU7cUK/5AlIuNWq1Wo8fP37w4MGePXsKIQ4dOtS+fftq1aoJIQIDA2NiYlJSUlwyTAAAcAs6p/5pSmIiXLNmTe3atWNjYxMSEvr06SOEyM7OrlOnjqNDUFBQdnY29+U2G1MQHgAA3EFV/TUuMRE+9dRTt27dSk9P371797Jly4QQ1atXNxr/e96juLi4jLCMnrzvFwAA3ER5fo3rhND9+/yoxL8KGHwZpCenpk2bjhs3bseOHUKIhg0bXrx40fFfly5dCg8PV3FwAADgXnTC3abB8odlLl++rExyRqNx27ZtrVu3FkIMGDAgMTHxwIEDnTt33rFjR0lJSa9evVw42MpGqmwYczmYS5eQgREu0MHVHjNQYQe2bpZkPS3m4Ey8gsm/WKlHJLd5rBB66nSNnUnWGJiwjI6Jl5AP32pldvdl6rSR7xXZfV/ZSLpMNT4psgXMyHbuOeFKr3EHJ7Mr3Ai5Hyt6hOwuxnQzUjT34Nx9hJo+q+WdCHv37u3n51e3bt0zZ860adPmxRdfFELUqFHjrbfeGjp0aGxs7LFjx1auXOnr6+vK0QIAAKisvBPh2bNnz549m5+fHx4e3rBhQ0d7YmLi0KFDz58/37x587p167pmkAAA4CacuKFep3OPFaHBYIiOjib/q379+vXr11dvSAAA4K6cuY/QNSMpPyQ5AQDAo6HoNgAAqMap/QhdNJbywkR4b9KvEbnHqeRuq2TZJ3a7UR8mxUe1cwexmLlIqsTIvbyZh+NDhzXJ7XO51CiXP7SSO9MyqVF2Z1puI18bsU+snbmrmBshve0tl/VlM8ASsV7udi+uThvdWTLqTL5puXc+F9W0mOkn12wiXgizke7MHYTerlki6gtVEyZCAABQjRPXCDX/mwPXCAEAwKNhRQgAAKrRCflrhFovCTERAgCAetxwh3qcGgUAAI+GFWE5yN4cSm6Uyu0Hy0QEyWynN5O95Np9/IjX12KSC9qxI6RzrcxIfLnyoUQ7t9OL1UIfxEBlNXU6ie8ohBBmJqpKNXO1Rjnk0+Lrz7yavnQ7F/elN+aV3PWXPLhUYJhr5zpzpDLDdAqUj6SSzwq7aTZ9DO2XL5WcO27Mi4kQAADU40SJNa3/uMBECAAA6sE1QgAAAPeCFSEAAKjGDReEmAjvC3eNnbhSz29vy4UXiMSEty+9gidDMUIIXyNRlcpiorMYUqWtBFO/zYcJetjZUAOVgGAKldnIWmqSB+ESN1w7+fKQmxILPs/iV514gfz86VeN3ZlWBhcA4aJPPtSrycV2uHAWeXC2uCBXBI6rgUdXkpP7sSI3jlanHJ1kexWmk98X+p7drVbrvn37bt682aNHj8DAQLLPmTNnTp482aZNm+bNm5f+wrS0NCFEixYtSnc+f/780aNHW7Ro0bJlS4FTowAAUJmZzeb+/fvPnj17w4YNLVq0OHHixN19li5dGhcXl5yc3LNnz5UrVyqNH3/8ca1atTp06DBhwoTSndeuXdulS5ctW7YMGDBg4cKFAhMhAACoSefUP95XX32VlZX1yy+/bNmy5cknn1ywYMEdHfLy8l5++eVt27Z9+umn33zzzfz584uKioQQcXFxaWlpq1atKt3ZaDTOnTt306ZNn3322e7duxcvXvzHH39gIgQAANUo9xFK/SvbP//5z5EjR/r6+gohHnvssW+++cZq/Z/rNdu3b2/atGmrVq2EEB07dgwKCvrxxx+FEGFhYSEhIXccbf/+/QaDoXfv3kKIyMjI2NjY77//HhMhAACoRueUMg545cqV8PBw5ePw8HCTyXT9+vU7OjRq1MjxaXh4eGZmZtlHc3xHpTPCMgAAoLGzZ8++/vrrpVs6duwYFxcnhDCbzQbDv/NZ3t7eQgij0Vi6Z+kOSp87OpTd2WQyYSJ0nlQFJnaPU66AGZnJZNJ6vn50u5l6fbl0KJ8alUiZcllNO7P5qVxqlImekudWrMzDEYIOwRos9AtEbnrMxSm5F8KvGvFCcOlQrjwY+/ApXK6V28PZmxq5L/MwuWww2c5FTLlqfFK7T3M5am4jaC+qaJyeiY1yqxS6HbFRB6funzCbzXl5eaXbCgsLlQ+Cg4Nv3LihfHzt2jWdThccHFy6Z3BwcOk14vXr10NDQ7lvdXfn7t27YyIEAAAV6WS3VdIJXcuWLRcvXkz+b48ePb777ru5c+cKIXbu3NmpUyfleqHFYtHr9Xq9vlu3blOnTs3Pz69du/a1a9dSU1M7d+7Mfa+OHTtmZmZevHixSZMmRUVFv/zyy5tvvolrhAAAUHlNnDjx1KlTzzzzzKpVq+bOnfv8888r7d26dVu+fLkQIjIyctiwYSNGjFi7du3IkSPHjBnTsGFDIURaWtrcuXOTkpKuXLkyd+7c9evXCyHq1auXmJg4evRopXPv3r1bt26NFSEAAKhGp/ZuEoGBgQcPHnz//ffPnDnz+eefKxcOhRCzZs1y3Cb/8ccff/jhhykpKePHj580aZLS6OXlFRgY2L179+7duwshatSoobQvX758w4YNKSkpQ4cOnTJlikBlGQAAUJMLaqw1bNhQufO9tDFjxjg+9vHxmTZt2h0dmjZt6lg+lmYwGBITExMTEx0tmAjvg0xahivjJLU/HJfR8GF2tvM1S4RlzJL7FFqpSAsX9GCyMswmc1zihjmI3kBVkuO2wZP8ESVTS3yhMq72GPFdufCLbFaIPLjBm8v+SORc2NJ9THE4XyoTxHbm3rRcf2owXGyHywSRm2jqmW0X2X0K5dIyHueet0OQX+KiwZQTrhECAIBHw4oQAABU48wO9a4ZSflhIgQAAPW44T5MODUKAAAeDStCAABQjRNhGc1XhJgI1cdt10tiU6NUvI2t7EWlQwVTZszCdLYwG/BaLczBydQotwEv0261ESO0M8XR2M1mTUS72UAfhdvKld0Olxo496pxI6TDsdxzRY+DHSGZnOSyl1zMksxwkpXhymgnU6Pc/sNkZ1FWypSKpDK5Vu5hkqXX2F182dQo1Uh29UhueGYUp0YBAMCzYUUIAACq0nyJJwkTIQAAqMcNb6jHRAgAAKpx5j5CrVeQuEYIAAAeDStCF5CJjTLBNDrJ5s3kBq3MfrBk4NNiYQKcMjVFBRN65NKhbF1Nup15rph0H1n50+BFh2C57KXVQj98usApVzqVKYdq0xPtbFiRaWe3w6VefXY7XKYIp48vFfiUTY2SwU6mpqh/NW+6vbrEN+Wip+RzIpgapFzlWC41itjoPTgRG9UaJkIAAFCNTid/zU/rc6M4NQoAAB4NK0IAAFCNExvzar0gxEQIAADq0QmdTvIiodbzICbCisL+ySOzYS9bZMtGRwPIyl5c+IVLi7CVwMiwDJMW4fbUJZ8Wvuwc3U5uKSwbimEfvlR1NO5hUrEL7uGwORfmEZFPl1TpPiGED5Vz8ZcppSaYgmdcyTTu4NVqMCEaqp2r38ZtKUw+h1xqSfNlirtywxpruEYIAAAeDStCAABQj/wN9ZrDRAgAAKpxYhsmnU7HbbpSMco7Eebm5m7fvj0jI6NRo0YJCQl+fn5Ku91u/+abb1JTU2NjYwcOHOiycQIAALhEea8Rtm/ffuPGjTk5Oe+9917btm3z8/OV9pkzZ86fP//27duzZs2aN2+ey8YJAADuQOfUP02Vd0V46NChoKAgIYTVam3VqtWXX36ZmJiYmZn50UcfnTt3LiwsbPz48W3atPnrX/9at25dVw646uG2eCUauVpQ3lzFLyrcSMYgy2jnkpBcQJTExanJcKPBy0J2NjOxSTI1auYqxnHpUKbynIXqz1WS4zbm9fYlRu7DpEP1zKvMvRL0C8f8ZuEiqWQBMzKoyXUWTECUC3b6MaXUqlWXKL3GRVLJUmpCCL2BaOdeNe5Kl9tdAKtgzhTd5vejrhjlXREqs6AQwmAw+Pj4GAwGIcSPP/7YunXrsLAwIURkZGSTJk327dvnooECAEDlp/v3VCjxT/MlofTtE0lJSTdu3EhISBBCXL16NSQkxPFfwcHBV69e5b5QagEBAACVTVX9NS6XGt2zZ8+MGTP++c9/1q5dWyhRn1LPi91u13x/RQAA0Ji73VAvMRHu379/9OjRn3/+eZcuXZSWkJCQrKwsR4fs7OzQ0FDuyzFHAgC4tfL8Gq/KG/MePnx45MiRa9eujYuLczTGxcWdPHnyypUrQoizZ89mZGT07NnTJcMEAABwjfKuCAcPHly9evUNGzZs2LBBCDFixIgxY8aEhoY+9dRT/fv3HzFiRFJS0l//+tfAwEBXjrYK4v8UogpIMn+3GJiqlWQu0U6nJqXP/pPfkvuDkS19SQVBTUzkz2Sk99ol2y2SqVGrGqlRLtZL7hPLlf3kCpmamYdvox4QVz+T27GWzHZye+dygU/yIOzuvjLpUO443GbF3HNLvg+5n0GtVynuypkb6l00lHIr70T43nvvWa3//TmMjo5WPnj77be3bt166tSpNWvWlF4sAgCAJ3LDotvlnQhHjhzJ/deAAQMGDBig0ngAAMCdOVFrVOuJELtPAACAR0PRbQAAUI1O/h4B2Y18VYeJsJKi30hcWIb7DyqOIHtDLJt/ob6nVChGCOHtQ1RTMzGBDlOJRFjGbKI7W5kQDRmKEUyIhgsWsTXwyIfP/OxzMR8uRGOgfoi5UAy37W21AJ+7G/2ZPAubf6HauWQNl8ThqqaRuRgvtpQa86almtnfwlqfr4MKg1OjAADg0bAiBAAA1Ti3H6GLBlNOmAgBAEA1TlSW0fwsNCZCAABQj/abSUjDNUIAAPBoWBG6E/ZMup7OE5JbvNJZPf5vOHKLYCGEnkrgkRufCj41SkYB2VJqbGqUiJ6ajUw61EwfRDI1SvalnxPBhBj56Ck9Em5PXTI5yW1vyxU286VSplyAk28nDsLu4uvH7KnLVU2j3kJcJTnuhaB/gtxt+VLJOXWN0EVjKS9MhAAAoJ4qvPsEAABAlYSJEAAAPBpOjQIAgGp0woltmHAfIdw3vgoaEcfQMWXA2DpTXACEaud2RvTmthj0JaIrXFjG7MeFZYh4BVdizWxSYZ9CLufCvRBkM5e44UqpcS+EHxVR4aqj+XARFSqKwhUw4+q3kfkXdstA5uBcrkpP1fTjduhkfxFrfS3KE7jjfYQ4NQoAAB4NK0IAAFBPFd6YFwAA4J7csdYoTo0CAIBHw4oQAABU40RYRuszo5gIqzT6hANTj43b3Ven40KM1EGYEmsGdmNeIlLoIxn4NNMb88qVWCPToUIIq5U4Dhf45JA/59wxuF8K3K6/ZIaTq4Lm7ctkNamDcwFOL6bYGx09ZQ7CPRy2Ohq9tzHSoZWSu10jxKlRAADwaFgRAgCAmjS/QV4WJkIAAFCNM9cItZ43MRECAIBqnLl9QusVJK4RAgCAR8OK0OPwu/syzVwuj4rxkdVNhRB6pgapF7UJLVeY1GJmgqBUbJLdaNcskQ4VQtisErVG2TSpzHaw3AvExSzJop3sLr5cEJR6gbjvyLZT2+SSmxILvnQqWyVU89NnUE6oLAMAAJ4MRbcBAADUZ7fbjUZj2X1KSkrubrRYLGazuewvxEQIAACq0f07LiPjXkvCJUuWBAYG1q9ff+jQoQUFBXd3+Omnn5o1axYSEtK8efODBw8qjTabbfr06XXq1AkKCpoyZYrFYlHaIyMj6/xHYmKiwEQIAACV2cGDB5csWXL48OGcnByDwbBgwYI7OlgsljFjxixYsCAvL2/OnDljx4612WxCiE8//XTXrl2XLl26cuVKSkrKRx99pPTPz8/fvn37+fPnz58/v2LFCoFrhODAn9ZnQg1UGoPfxZep60YlKWxWLhVCH4TMs7Ab7VKdBb8drp1qJxsFXzWNxIZlZHZCFkx0hdshWS7nwlY7kxgh+76S2cQY3Ijq9xFu2LBh9OjRzZo1E0LMmTNnxIgRy5YtK91h586ddrt93LhxQohJkybNnz9///793bt337Bhw9SpUwMDA4UQ06dPX7t27VNPPaV8Sc2aNZV2BVaEAABQeZ07dy46Olr5OCYmJicnJy8v7+4OSq7YYDA88MAD586dE0KcP3++9BeeP3/e8SU9evSoV6/ekCFDzpw5I7AiBAAANcnfUC90Ij8//7fffivd1rhx46CgICFEfn5+jRo1lEblg7y8vNLruYKCgurVqzs+DQgIUGbKO77QMX0mJSW1b9++pKRk4cKFgwcPPnXqFCZCAABQjxP3EQpx9OhRx3lLxciRI+fNmyeEqFevniMgk5+fL4SoX79+6Z5BQUE3b950fJqXl6d0CAoKKv2Fjq+Ki4sTQtSqVWv58uWBgYHHjh3DRAgAABrr3bt3cnIy+V8xMTFHjhxRPj5y5EijRo0c6zxHh+PHj1ssFi8vL6PRePr06ZiYGKX98OHDgwYNUr5QaSzNarXabDYDt3scAACAE5SwjOy/MkyaNCk5OXnr1q0XL15csGDBk08+qbTPnDnzyy+/FEJ069YtJCRk4cKF2dnZL774YosWLdq1ayeEePLJJ1euXHn48OGTJ0++9dZbyheeOnUqKSkpPT399OnTTzzxRKNGjWJjY7EihHuQTZPSXemCX3T+kI2YMu9WMvBpt9F/5PHpUPrgZDU12Y15pbA1xmSymlzgU27bW9lgJ/Ufarx9wJ04U3S7zP4PPPDA+vXrX3jhhYKCgoSEhOeee+7uL09OTp41a1bXrl1bt279xRdfKO1DhgyZN2/exIkTbTbbzJkzH3nkEaXzunXr5s+f7+fn9+CDD27dutXHx0fHVU1UV3Fx8fSFS4Y+/mQFfC/QEP9ukijOyd6cwBycmQhlb5OgD46JkGonmzERVn0jHwi5Z5+Pd6dfv0lUeCnDyQM781J3cadGKwBOjQIAgEfDqVEAAFCNU6FRjWEiBAAA1ah+jbAC4NQoAAB4NKwIQU0ujZhyERUyACK9dy7Tbif/Q5WwjGxxV3bHWrKR6y1xEHYk9DHc8IwYuAI25gUAAI+mkz7VqfmpUUyEAACgGmd2n3DNSMpP7hphdnZ2dnb2HY1ZWVm7d+/+448/1BsVAABABSnvRPj++++HhISEhoZOnjz5jvbWrVsvXry4ZcuWn332mQtGCAAAbkUn/09T5T012r179127dn377bd79uxxNN66dWvOnDk7duzo1KnTrl27Hn300YcfftjX19c1Q4UqSJWMhp1Mi3Dfkg3FcLT+Gf0P2WyN3MEry6MEt6dUD9V6FHLKuyKMiYmJjo7W6/+n/9atW8PDwzt16iSEiIuLq1at2u7du1UfIgAAgOvcV1jm8uXLTZo0cXzauHHjjIyM+x0RAAC4LWfCMlovIO9rIiwpKfHx8XF86uvrW1xczHW22ZiqxgAA4A5sNtsd5wUJleCan6z7mgiDg4NzcnIcn964cSMkhK1Nfu+nDwAAKrHy/Br3uBJrnTp1+u2334qKioQQubm5qampHTt2VGlgAAAAFaG8K8JTp059++23+/btO3/+/Ouvv96mTZtBgwa1bNmyZ8+eY8eOTUxMXLVq1bBhw0pfMgRQmWTVMLlmmapprq2wptoXAGjAiWuEmr+3y7siNJlMeXl5rVq1SkhIyMvLU1aBQoikpKR27dp98sknPXv23LBhg8vGCQAA7sDdbiIU5V8RtmvXrl27dne3BwQEvPTSS6oOCQAAoOKg1igAAKjI/cIymAgBAEA1brgLEyZCAIXMz6LmP7cAlZcbzoS4tw8AADwaVoQAAKAaJ26o1xwmQgAAUI071hrFqVEAAPBoWBECAIBqdDrp2yE0P5WKiRAAAFSl9alOWTg1CgAAHg0rQgAAUJO7LQgxEQIAgHrccT9CTIQAAKAeVJYBAABwL1gRAgCAatzxhnpMhAAAoBp3LLGGU6MAAODRsCIEAAD1uGFYBhMhAACoRid/zU/zM6mYCAEAQDU6odNJLvFk+6sO1wgBAMCjYUUIAADqceIaodYwEQIAgHrk7yPUfOLEqVEAAPBoWBECAIBqUHQbAAA8G64RAgCAJ9MJJ1aELhpLeeEaIQAAeDSsCAEAQDVuWGENEyEAAKjIDWdCnBoFAACPhhUhAACoRqdz4nYI3D4BAABVBnaoBwAAj4ZrhAAAAO4FK0IAAFCNMzfUa70kxEQIAACq0WH3CQAAAPeCiRAAADwaTo0CAIBqnNqGyUVjKS9MhAAAoBonrhFqPQ/i1CgAAHg2rAgBAEA9bnhDPSZCAABQj/w1Qs0vEuLUKAAAeDSsCAEAQDU6+QWe1gtCTIQAAKAeZ26f0PoiYQVNhF5eXmf2bN2T9P8q5ttVsJycnMDAQL2+ip9nxsOsSm7cuFG3bl35fePcDB6muqK//TY6OrrsPn2aBMketllxu1/0Jc4OSgU6u91eMd8pNzc3Pz+/Yr4XAACormHDhj4+PlqPQn0VNxECAABUQlX8/A8AAEDZMBECAIBHw0QIAAAeDRMhAAB4NNxHKO2PP/5ISUnJzMzs06dPRESEo/3SpUvr168vKioaNWpUx44dNRyhKg4fPrxt27br169HR0ePHTvW399fab958+YHH3yQmZnZu3fvYcOGaTvI+7djx479+/fn5+c3atRo/PjxdevWVdrz8/M/+OCDrKysvn37Dh48WNtBqigpKcnX1zchIUH51Gg0fvTRR+fOnWvXrt3YsWPd/aaRr7/+Ojs7W/m4Tp06jzzyiPJxbm7uhx9+mJ2dPXDgwP79+2s3QNVcu3Ztw4YNV69ebdq06cSJE2vVqiVKvWn79OkzZMgQrcfoTtz7fa+JHj16vPbaa88//3xKSoqjMTs7u2PHjvn5+fXr1+/Xr9/evXs1HOH9y8/Pj4+Pv379eqNGjT755JMePXoYjUYhhM1mi4uL279/f0RExKxZs9555x2tR3q/Nm3aZLPZmjVr9tNPP7Vt2zY3N1cIYbFYevbsmZKS0qxZs7/85S+rV6/Wepjq+Prrr6dMmbJkyRJHy2OPPbZ58+aoqKhly5bNmTNHw7GpYsmSJdu2bUtPT09PT79y5YrSaDKZunfvfvz48aZNmyYmJq5fv17TMaogLS2tTZs2p06datKkydmzZ5VHarFYevXqlZKSEhERMX369FWrVmk9TLdiB0lWq9Vut8fGxn7++eeOxgULFgwfPlz5eMmSJYMHD9ZmcCqxWq1Go1H5uLi4uFatWnv37rXb7d9//32TJk3MZrPdbt+5c2doaKjycRVgs9maNm2anJxst9uTk5ObN2+uvNDKQ1Y+dmv5+fktW7ZcsGBB165dlZaTJ09Wr1795s2bdrv9/Pnz/v7+OTk5mo7xfnXr1u3rr7++o3Hjxo2tWrWy2Wx2uz05OTkqKkr52H0NHDhw3rx5dzQqD015o/7rX/9q3LixxWLRYnRuCStCaeTpo7179zpOufTr12/Pnj0VOyiV6fV6x22zNpvNZDIFBAQIIfbs2fPQQw95eXkJIXr16pWTk3P27FktB6qes2fP5ufnt2jRQgixZ8+ePn36KC90nz59MjIyLl68qPH47tuzzz777LPPhoaGOlr27t374IMPKq9ss2bNwsLCDh48qN0A1fGvf/3rzTff/P777+3/uUN67969ffv2Vaqu9O/fPy0t7erVq5qO8b6Yzebt27cnJCSsXbt29erVjoXvHW/aK1euVIE3bYXBRKiOrKysevXqKR/Xr1+/qKjo5s2b2g5JLX/729969uzZtm1bIUR2drbjYRoMhrp162ZlZWk6OhU899xzYWFhbdq0Wbp0qTIRln41fXx8AgMD3f1h7ty588KFC4mJiaUbS7+aQoj69eu79QwhhIiJifH19b127drMmTPj4+NtNpv431ezWrVqNWrUcOtX8/Llyzabbdq0aRcvXjxx4kRsbOzp06fF/76a3t7eVeBNW5EQllGHl5eXxWJRPlY+8Pb21nRE6nj77be3b9/uuOTp5eVltVod/2s2m6tAvaWXXnpp9uzZ+/btmzp1auvWrTt27Ojt7V2VHmZRUdHMmTO//PLLO2pRVr1X8/3331c+mDt3bvPmzbdt2zZw4MDSP5tCCIvF4tYPU6/X2+32adOmKX/WmM3mN99886OPPqp6r2ZFwopQHWFhYY6/pjMzM+vUqeOIWbqv5cuXv/vuu7t27QoODlZawsLCMjMzlY+Li4vz8vJKn2pzU9WrVw8ODh41atTAgQO/+uor8b8P89atW7du3XLrh7lnz57MzMxx48Z16NBh0aJFx48f79Chg81mK/0whRCZmZlu/TBLCwwMjImJuXDhgvjfn80bN26UlJS49cMMCQnR6/UxMTHKpy1btrx06ZK460178+ZNt36YFQwToTri4+O3bNminIr54osv4uPjtR7R/frwww+XLVu2ffv2hg0bOhrj4+O3b9+uFE9PTk5+4IEHSt9A4nYsFovZbFY+NplMx48fb9SokRAiPj7+hx9+8Q1dmwAAAqVJREFUuHXrlhBi8+bN7dq1CwsL03Kg96dbt267du1as2bNmjVrxo8f36xZszVr1uj1+oEDBx45ckT5Nfrzzz8bjcauXbtqPVjnmc1mx8rv8uXLR48ebdmypRAiPj7++++/v337thBi8+bNXbp0CQqS3h6h8vD19R00aNCBAweUTw8cOKBMivHx8Vu3blXetFu2bGnbtm3pn1y4B63TOu5n+vTp7du39/f3b9asWfv27VNSUux2e2Fh4Z/+9KeePXuOHj26QYMGv//+u9bDvC+ZmZk6na5Ro0bt/+O7775T/mvs2LExMTGPP/54UFDQN998o+0471NGRkaDBg2GDx8+duzYxo0b9+3bt7i4WPmvRx55pFWrVhMnTgwKCtq6dau241TRBx984EiN2u32efPmNWnSJDExMTg4eM2aNRoO7P6dP38+JCRk5MiRo0ePDgwMnDZtmtJus9ni4+Pbtm07YcKEunXr/vjjj5oOUwWHDx+uX7/+hAkTBg8eHBkZefXqVaV91KhRjjftDz/8oO0g3Qt2n5CWlpZWOgjTvHlzJXdnNBp37dpVWFjYt2/fwMBA7QaoApPJdOLEidItTZo0UW42t9vt+/bty8zM7Nq1a+PGjTUaoGoyMjKOHj1aUlISFRXVrl07R7vdbt+zZ092dna3bt3Cw8M1HKG6bty4kZOT88ADDzhaUlJS0tLS2rZte8995iq/1NTU1NRUm83Wpk2b5s2bO9ptNtvu3buvXbvWo0cPt17cO+Tk5Ozatat27drdu3d3XIWx2+179+7NysqqYm/aCoCJEAAAPBquEQIAgEfDRAgAAB4NEyEAAHg0TIQAAODRMBECAIBHw0QIAAAeDRMhAAB4NEyEAADg0TARAgCAR8NECAAAHg0TIQAAeLT/DwtcXstToy8JAAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"terms = [Kinetic(),\n",
" ExternalFromReal(X -> pot(X...)),\n",
" LocalNonlinearity(ρ -> C * ρ^α),\n",
" Magnetic(Apot),\n",
"]\n",
"model = Model(lattice; n_electrons, terms, spin_polarization=:spinless) # spinless electrons\n",
"basis = PlaneWaveBasis(model; Ecut, kgrid=(1, 1, 1))\n",
"scfres = direct_minimization(basis, tol=1e-5) # Reduce tol for production\n",
"heatmap(scfres.ρ[:, :, 1, 1], c=:blues)"
],
"metadata": {},
"execution_count": 5
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.3"
},
"kernelspec": {
"name": "julia-1.7",
"display_name": "Julia 1.7.3",
"language": "julia"
}
},
"nbformat": 4
}