{
"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.045210e+01 6.561257e+00\n",
" * time: 0.004136085510253906\n",
" 1 2.728247e+01 3.891898e+00\n",
" * time: 0.012537956237792969\n",
" 2 2.013589e+01 5.633218e+00\n",
" * time: 0.03269314765930176\n",
" 3 1.277056e+01 1.930491e+00\n",
" * time: 0.05701398849487305\n",
" 4 1.160173e+01 2.645833e+00\n",
" * time: 0.07322001457214355\n",
" 5 1.126027e+01 1.301441e+00\n",
" * time: 0.08948993682861328\n",
" 6 1.039500e+01 1.146777e+00\n",
" * time: 0.10560798645019531\n",
" 7 9.870690e+00 9.528918e-01\n",
" * time: 0.12183904647827148\n",
" 8 9.541662e+00 6.253213e-01\n",
" * time: 0.13797211647033691\n",
" 9 9.377753e+00 6.056077e-01\n",
" * time: 0.15415096282958984\n",
" 10 9.258907e+00 6.264517e-01\n",
" * time: 0.17035698890686035\n",
" 11 9.130028e+00 7.151177e-01\n",
" * time: 0.1827099323272705\n",
" 12 8.995444e+00 5.639402e-01\n",
" * time: 0.19528412818908691\n",
" 13 8.847151e+00 6.994983e-01\n",
" * time: 0.2080371379852295\n",
" 14 8.767840e+00 6.051972e-01\n",
" * time: 0.29551005363464355\n",
" 15 8.733139e+00 5.843778e-01\n",
" * time: 0.3080611228942871\n",
" 16 8.717371e+00 6.848829e-01\n",
" * time: 0.3201310634613037\n",
" 17 8.649984e+00 4.001283e-01\n",
" * time: 0.331820011138916\n",
" 18 8.588047e+00 3.767248e-01\n",
" * time: 0.34352707862854004\n",
" 19 8.559839e+00 2.977207e-01\n",
" * time: 0.3553769588470459\n",
" 20 8.538579e+00 2.853830e-01\n",
" * time: 0.36740899085998535\n",
" 21 8.528747e+00 2.978936e-01\n",
" * time: 0.379425048828125\n",
" 22 8.511501e+00 2.314970e-01\n",
" * time: 0.3913991451263428\n",
" 23 8.492743e+00 1.527322e-01\n",
" * time: 0.40293312072753906\n",
" 24 8.459099e+00 1.919458e-01\n",
" * time: 0.4149949550628662\n",
" 25 8.437249e+00 2.144477e-01\n",
" * time: 0.4269731044769287\n",
" 26 8.418725e+00 1.585196e-01\n",
" * time: 0.43874406814575195\n",
" 27 8.401146e+00 1.777631e-01\n",
" * time: 0.4505789279937744\n",
" 28 8.387393e+00 1.623046e-01\n",
" * time: 0.4624049663543701\n",
" 29 8.372171e+00 2.212899e-01\n",
" * time: 0.47414302825927734\n",
" 30 8.358771e+00 1.307231e-01\n",
" * time: 0.48568296432495117\n",
" 31 8.345731e+00 1.526229e-01\n",
" * time: 0.4974820613861084\n",
" 32 8.341681e+00 1.684810e-01\n",
" * time: 0.5092079639434814\n",
" 33 8.341530e+00 1.436405e-01\n",
" * time: 0.5287909507751465\n",
" 34 8.340986e+00 1.214301e-01\n",
" * time: 0.5441031455993652\n",
" 35 8.338770e+00 1.601651e-01\n",
" * time: 0.5557460784912109\n",
" 36 8.338770e+00 1.601651e-01\n",
" * time: 0.6651711463928223\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1xU1doH8LVnGEDAFEHkookoKKhc8pamaIimJqKmdIrUojxa2ZulZSc/HbNjHbPsZie17GL5ZpKGx8zyLlBmiZoakje8gaBy84IDc33/2J155+h6dNa4h80wv+/HP2Cx3Ky5wGLt/dvPkqxWKwMAAPBUGrUHAAAAoCZMhAAA4NEwEQIAgEfDRAgAAB4NEyEAAHg0TIQAAODRMBECAECjtnLlym7durVt23batGn19fXXdygqKkpJSQkNDR06dOixY8fkxqVLlw4YMCA8PDwpKemTTz6xdR48eHDP//jb3/7GGJNwHyEAADRaRUVFd95557p167p06XLfffcNHTr073//u30Hq9XatWvXzMzMJ598cuHChRs2bNizZw9jbPbs2QMHDkxISPjtt98yMjJWrVo1bNgwxlhwcPCyZcvatWvHGAsMDIyKisJECAAAjdfMmTMrKys//fRTxtimTZuysrJKSkrsO+Tl5d13333l5eVardZgMLRu3Xrbtm09evSw75ORkdG5c+d//OMfjLHg4OCff/45Ojra9lWcGgUAgMbr0KFDiYmJ8seJiYmlpaWXLl2y71BUVBQfH6/Vahlj3t7ecXFxhw4dsu+g1+t3796dlJRka3nwwQd79eo1bdq0c+fOMca8XP4gGGOMWSyWvQd+j46NbZhv18AsFotG0/T/pMDDbErwMJuSBnuYLXx0N+1z+HCVXm8SOuzlyxdbtza3adPGvjEgIECn0zHGKisrb7vttj8H0KIFY6yiosLWIndo3ry57dOWLVtWVFTYPrVardOmTYuJiRkzZozcsnDhwsTERL1eP3/+/GHDhv36668NNBHW19d/sCpn5MOtG+bbAQCA4sZ2Drtpn/vv/3b//guCB/5dp/syICDAvumxxx5bsGABY6xVq1aXL1+WG+W1YFBQkH3PwMDAK1eu2D69ePGifYfnnnvu4MGDW7ZskSRJbpk0aZL8wVdffdWqVasDBw400EQIAABAuffee3Nycrhf6tSpk+1U56FDh1q3bi2vC+07FBUVWa1WSZLMZvORI0c6deokf2n27Nlbt27dunWr/QrSxsfHx9vbu76+vumfNAAAgAYjSc78u4FJkyZlZ2cfPny4vr5+wYIFtvXcP//5z61btzLGUlJSvLy8Pv74Y8bYBx98EBwc3LdvX8bY3Llzly9fvnTp0pqamuLi4gsXLjDGiouLd+/ebTQaa2trX3jhhYCAgMTEREyEAADQePXs2XPOnDn9+/cPDg728/Oz3TtRUFBw5swZxphWq83Ozn777bf9/PyWLVu2cuVK+Szotm3bfHx8HnjggSFDhgwZMuStt95ijFVXV0+YMCEgICA0NHT37t3r16/38/NroNsn9Hr9k68sGPnwXxvgewEAgCs4co3wjjs+F71GaLX+np5eQ50abQC4RggAAIq64anORginRgEAwKNhRQgAAMq5WfilEcJECAAAipEkSRKdCa0qz5w4NQoAAB4NEyEAAHg0nBoFAADF3PQG+eupfWYUEyEAACjHiWuEqodrcGoUAAA8GlaEAACgHMn9bqjHRAgAAIqRxE91qn5qFBMhAAAoxpn7CNWeCXGNEAAAPBpWhAAAoBwnrhHi1CgAADQl7jYPCk6Ev/3229GjR1u3bn3nnXf6+vrKjfv37z9y5EhcXFzXrl1dMEIAAAAXcvQaodlsnjBhwqhRo1avXv3yyy+vWLFCbv/nP/957733fvfdd6mpqYsWLXLZOAEAwA1I0p95GQFqrwkdXREuXrz44MGDhYWFzZs3Z4zJ+9pXVlbOmzevoKAgNja2oKAgNTX1kUceCQgIcOF4AQCgMXPDa4SOrghXrFjx9NNPV1ZW7tq1q7a2Vk7Hbt68OTo6OjY2ljHWs2fPoKCg3NxcFw4WAAAaNycWhMK3WyjN0RXhsWPHvvrqqyVLljRv3vyPP/5Yv359YmJiaWlpu3btbH3atWtXUlJCHUFeRAIAgJuyWq2qT1qu4OiK8OrVq82bN9+1a9eWLVsmTZo0c+ZMxpjRaNRqtbY+Xl5eRqPRJcMEAAC1ObKekXefEPqnOkcnwvDw8MGDB8t/CwwZMmT//v2MsbCwsAsXLtj6nD9/PiwsjDpCk/w7AgDAc2g0DkwZklP/VOXoRDho0KATJ07IH584cSI8PJwx1r9//71791ZXVzPGysrKjhw50rdvXxcNFAAAwBUcvUb47LPPDhw4sGXLlrfddturr7769ttvM8Y6duyYnp4+evTozMzMTz/9dMKECfIECQAAnkk+3yn4X1Tm6IowLi4uPz//6tWrZ8+ezcnJ+ctf/iK3f/HFFw888EBhYeHkyZOXLFnisnECAIA7EL9GqPp1M4HKMl26dJk3b941jTqdburUqYoOCQAA3FYjuOYnCrtPAACAR0PRbQAAUIx8Q73g/3GTG+oBAABuyolrfmrPgzg1CgAAng0rQgAAUI4bhmUwEQIAgGKcKKKtet0xTIQAAKAYN1wQ4hohAAB4NqwIAcQpsqWY2/3ZDOAIN9yYFxMhAAAox4lrhC4aicNwahQAADwaVoQAAKAYSfwGebVDo5gIAQBAOZIkfjuE2jMhJkJwH0RERTC5InYURWIxXJLoocnfFQK/RJQ4BsANuWFYBtcIAQDAo2FFCAAAinHDBSEmQgAAUI4zJdbUngpxahQAADwaVoQAAKAcNzw3iokQGoQSUU0rFbPktVuJ3mRUkxqK4GFE8H/6ybNKRDu3mTqIVeQg5Aj5ndX/dQaNgRP3Ear+zsGpUQAA8GhYEQIAgHLcsNYoJkIAAFAOrhECAIAnkyT3qzWKa4QAAODRsCIEZ3GzmmRngSAoFfi0WIhjWzj9LdR35HVmjFmpg/OOIxo9FfqDl7q+IhF/tWo0nP4Sr5ExphE5uGh+lfsFREw9jcQk0RvkVb+hHhMhAAAox4lrhGrDqVEAAPBoWBECAIBinAjLqL6CxEQIAACKkZj4fYRqx0YxEcJN0IXNeCkSojMZXTFz2s28RsaYhWg3mzlBF7OJStyIHZwbrqHiPCTejzkZiiHaNVp+u5bXTnfmXw3h9ucemTEmEe38gI5gWkbtX4lwy3CNEAAAwL1gRQgAAMqRhE91qn4aABMhAAAoxg3PjOLUKAAAeDZMhAAAoBzJqX83U1dXd/HixRt0MJvNVVVVlusKUNXU1FzfaLFYqqqqTCaT/ClOjXoesT1y6R1uuYXNyGAnv91k5AQ+TSZ+uTNuZ8aYmXsQqjNxcDKqyk2ZilSMY8QlEDIdSlVH8yKCoF6cP2e9eI2MMS+dQDvVmfsdmWB+laoYR+FX4XK7c3AeQHJmG6ab9J8zZ867776r0+m6deu2Zs2aVq1aXdNhy5YtEydO1Gg0kiR9+eWXAwYMYIw9//zzy5YtkyTJYDD89a9/feONNzQaDWPs119/vf/++w0Gg9ls/uSTT0aMGIEVIQAAKEb6zz31Qv9u4KefflqyZElhYeH58+dbt2798ssvX9PBaDROnDjxvffeKykpee211yZOnCgvAbt06XLw4MHKysrff//9q6++WrlyJWPMarU+/PDDs2bNKi0t/fTTTydMmKDX6zERAgBA47VixYqMjIyIiAhJkqZPn/7FF19cc5pqy5YtOp1u3LhxjLHMzEy9Xp+fn88Yy8rKioiIYIy1b9++X79+hw8fZowVFBSUlpZmZWUxxoYPH96mTZvvvvsOEyEAAChH6WuEx48f79y5s/xx586da2pqqqurqQ4ajSY6Orq4uNi+w4ULF3Jzc1NSUhhjxcXFUVFR3t7e8pe6dOlSXFyMa4QAAKAYJ64RMkmqqanZs2ePfVtERERoaChj7OLFi/7+/nJjQEAAY6ympsb+MuGlS5f8/PxsnwYEBNTU1Ng+ra+vf/DBB8eMGTNo0CCqMyZCAABQ2W+//TZlyhT7lvHjx8+aNYsxFhISYpvY5LVgSEiIfc/WrVvbz3zV1dW2DkajMSMjo0WLFv/6179sne3Tp9XV1T169MBE2KTx984V25mWKs7JjV9SWU2jgWo3cxrrOY1UZ+rg1EjINCmvYCkjcrB0/VWinZ8a5fcVLRMqFPjUeWv57T6c/t5kZ367ItFTKjTLNNx3MwqWNjpO7D4hMTZo0KCcnBzuV7t161ZQUCB/XFBQ0KFDB3ldaN/hwIEDRqNRp9PV1dX9/vvv3bt3Z4yZzeaJEycajcbs7Gwvrz8nu65dux47duzixYstWrSwWCx79uyZPn06rhECAEDj9eijj3777bc5OTlFRUVz5syZOnWq3D558uTs7GzGWN++fSMjI1988cXi4uLnn38+MTExPj5e/o8///zzk08+mZ+fv2XLlqKiIsZYdHR0cnLyjBkziouLX3755dtuu+3uu+/GihAAAJTjxH2EN+zeqVOnVatWzZ8/v7q6esyYMTNmzJDbW7Zs2axZM/njnJycmTNnpqWlxcfHy7MjY8zf3793797Lly+XP01JSYmNjWWMrVixYsaMGWlpaZ07d/7uu+80Go0CE6HBYCgtLW3btq1Op7v1owEAANgbPnz48OHDr2l84403bB+3b9/+66+/vqaD7brgNUJCQr744gv7FkdPjc6cOVOyc/XqVbl93bp1ERERw4cPb9eu3ZYtWxw8GgAANEnO3E2v9jVdgRXh7Nmz582bZ99SX1//6KOPLl++fMSIEdnZ2VlZWSdOnNBq+RfSwYWoqmm8UAeRfWEWqrAZ0c6NqBjq+HkWQ52J217P6091NpAhGl5YhojnUA/HokRYhtywl3fqhzodRG2HqxEJy+iosAyRc/HmtXv7Ep19+b83fHj9uUdmjHl5858rKlyj4f3VrpGI5Bd2/VWXuz3PYmEZvV5v/+mmTZtuu+22ESNGMMbGjRtnMBjy8vKUHB0AAICLCUyEb775ZuvWrUNDQxcuXCi3nDhxIiYm5s8DaTSdOnU6efKk4kMEAAB3ITlF3TE7OhFOnTq1urr6ypUra9euffXVV9esWcMYu3z5si20wxjz8/O7dOkSdYTrN8IAAAA3YjbzL0/Yc+IaodrzoMMTYadOneQ5784775wwYcIPP/zAGAsJCbGv+VZdXd2mTRvyO2lwzyIAgBtzJALixDyo+kVFZyan6upqufJbQkLCvn37DAYDY6y2trawsDAhIUHhAQIAALiSo6nROXPmDBgwoHnz5jt27Pj6669//PFHxljv3r1jYmKeeeaZKVOmvPPOO3379pVvVwRXEQwrcs9GU7XEqNpjVFbToOdkO+t4jYyxuqtEapTXv14vFj018FOj/IOQG/CSqVFOI7lZsSIb85Il1oh2F6ZG+b8ffJvxn1uTkdPfbCIiphZ+O/UcenHvUqbqtIkVuyP6gnPUX+AJc3RFaLVa58+fP3369AMHDuzYsaNHjx5y+9q1a/V6/WOPPabT6VatWuWycQIAgDtww2uEjq4IX3nlFW57eHj4J598otx4AAAAGhRqjQIAgGIk+rQ/+V/cZUUIAABwc254jRAToTuhthKkqqZxczFU7bF6qgqaSP5FX0t1NvLbef2pxA1dek1gP0KhfQcZY1bek0vt0SiE2nhPItqpsAyxC6DAvoOMCMv4+FI17fi/N4y859xsopJcAk84Y4xZOd+U/H1LhWh4Jdkkt/u13bg5sx+h2q8A7u0DAACPhhUhAAAoxpmSaWovCbEiBAAAj4aJEAAAPBpOjQIAgGKcCcu4ZiSOw0TYSHELTVmJDTzIPXV5KT4qHUpVQdPX8gOf+iucdio1Sh2E+02pnCrVzt+Y10iUWCNCjFTpNdelRukSa/z+WmJjXi0vIanzJkqseVN77XLajbw4LiPSoYwxM+99KBTHZXSJNT6J/xtMR/xmlbjnwIjzYmpft3JXEhO+Rug22zABAAA0SVgRAgCAcpy4oV7txTcmQgAAUE4jKKItChMhAAAoxon7CHGNEAAAQE1YEaqMyshxt36lwoomIjXK3VOXyl5yU6CMsdrL/ParVwycg1C1RqnUKG8wBmJjXirvKpYaNVKpUYEapBbyZeM3c1E1RakapFovql2BjXmNBk67cL1W3vuTijqL4m9uTNVrpVKjvC8Qx2DUlSu1Vy+gPEyEAACgGKdOjbpoLI7CqVEAAPBoWBECAIBinKgsg9snAACgaVF7YhOFibChUDEKInbBLUDFLWHFiLQIY8zAi6JQeZarRFiGG4phjF3lhWiESqkxIrljqOPnXLjZH8aYkddOlgFTYsNeZUqsUaEYMixDlFjj5WLMRFiGqiQn9DC5SS5GvJepX4l0noXfruHtSyy6uTG3vyTxnyuiuRFc0WrcJEn4dgjcPgEAAKAmrAgBAEA5Tuw+ofYaGxMhAAAox4lao2rDqVEAAPBoWBECAIBiJCZJ7rb9BCbCBiIYGuWn+KiSV9zYJCMKmFHVzsjAp0jVNLLE2lX+wet51dQMRCk1MjXKLbFm4Hemy4YJxCmFYpNMsDyYlhePZHRq1MvEabfwGhkdBOVXRxOtJMdNjQqmQyXi4Wt4+xILbVbMiOeWjJ5SQxR5mB7ImR3q1X72MBECAICi3GxBiGuEAADg2bAiBAAAxThTdNtFQ3EYJkIAAFCMO9YaxalRAADwaFgRugA3U0cE8KgUHzfESNUUpfba5Vb41BNlP6lyoFTgk3scMpIqVmtUMDVaz92Yl/9cUZsYW6h2fpyS25eOU/L+4KXCityimuwGqVFeWVGTErVGuQVvGb3XLr+3yB657AYPn9fOfeyi7dQTK2n5D4hoph4nv7kpc+KGerWfJUyEAACgHGeuEeI+QgAAaCrccEGIa4QAAODZsCIEAADFOHH7hOqlZTARKk8kK0OGZbhZD6qUWj2xky03ilIvHJZRIolDhGj4JdbqFSmxRpRSI0IxQhvzkmEZkXJiZFiGLL3GP39j0nFTJFpuZ+7DodotgpXkuMiNdsn9hwWyQjpv/sPUeRNhGV47FZahYjsafsU8bl/1T/qpw93OjeLUKAAAeDSsCAEAQDHOFN12zUgch4kQAAAUI0n0xh30f1EXJkIAAFCOG94/gWuEAADg0bAidAFeqI7a45RK8XFTowZis1kDlRrltdcTBcyodio1yk2fUsXeqFwrNyBKPRyx1KhRmY15rUqkRjW8PzipjXnp1ChR8cvM6U+9r8iqaSJJUKoOCHfkGiIF6kVWjOO/cF7enHZv4v3m7cNPk3r7cA5CRU+pOm1W3gtBL2nUXuw0ODdcEGIiBAAABbnhfYQ4NQoAAB4NK0IAAFCU6uc6BYmtCI1GY1pa2oMPPmhrKSoqSklJCQsLGzZsWHFxsdLDAwAAdyLfRyj6T11iK8LXXnvt9OnTmv9c+rdarWPHjn3ooYfWrFnz5ptv3n///bt373bBIBspKl7Abab3HSS2zeOWWCPCMmQ7L3VCRlEEIyrcdqoInFA7ve+gQFiGek7o/QgFNoykXnrq55mbi+EmaBhjGqKUmoUIy1gsnP5CVdAYMXLqSg8Z8+HVJKNKphmIKIqOejOLvFWEclVU1T0LVXWP94RLxKtJ/oZX+1e/6zhRa1T4mqLSBFaEBw8eXLdu3bPPPmtrycvLq6ysfOGFFwIDA//+978fPnx47969LhgkAACAqzg6EZpMpsmTJy9atEin09kai4qKEhIStFotY8zHxyc2NvbQoUMuGSYAALgFSfyf2hydCN98882+ffv269fPvrGysrJ58+a2T1u2bFlRUUEdwWLhn2cAAAC3YDbzzznbc+YCodpzoUPXCE+ePPn+++9v2LChuLj4/Pnz9fX1xcXFHTp0CAwMvHLliq3bxYsXg4KCqINoqEsiAADgDuTzf6ooLCysqanp0aOHr68vt0N5efnhw4djY2NDQkLs2y9cuGAwGCIiImwtp06dss3o/v7+bdq0cWgirKioCA0NzcrKYoxVVVWdO3cuIyMjPz+/U6dORUVFVqtVkiSTyXT06NGOHTs6+SgBAMD9ycs8wf9yIxaLJSMjY9++fW3btj158uS2bduun2iWLVv2wgsv9OjRY8+ePe+++25mZiZjbPXq1c8++2xJSUnPnj1//fVXW+cePXr4+fnJl/lGjBixaNEihybCnj17FhQUyB9/+eWXb7zxhvxpSkqKJEmffvppVlbWBx980Lp16759+zr+4JssXmDPSpwYNhNhRW6FMGqzWaHsHNmZSI3SWU1eapQYIZXhNNQL7D8suDGvWGqULLGmRGqU+wUiesk0VCk1qlYZf4RisVFubE/S8J9Dasda7gipQmVe1FtFJDMsmqPm/1gJVt3jPuHMyn9OiGbVzwW6ktI11jZs2LBv3779+/cHBAQ888wzc+bMWbFihX2HK1euzJgxY9OmTX369NmxY8f48ePHjRvn4+OTmJi4fv36PXv2LF68+Jpjbt26NTo62vap8OnKgIAA2xrTy8tr1apVb7zxhr+//7Jly1auXKl6ChYAAJqS1atXjx8/PiAggDE2adKkb7755prEycaNG8PDw/v06cMYGzRokL+///bt2xljnTp1io+P9/LirPfOnz9fXFxsO0EqXFlm1KhRo0aNsn3at29f29lR0UMBAECTI34f4Q2XhKdPn5YnOcZYhw4d9Hp9RUWF/YXAM2fOREZG2j6NjIw8c+bMDQ6o1WofffRRvV5vMBiWLVt27733KlNiDbMgAAAwp3aoZxI7efLkhx9+aN+WlJTUq1cvxpher/f29pYbfXx8GGNXr16171lXV2frIPe5psM1Dh06JOc6P/744wcffPDkyZOoNQoAAMpx6hphTU3Nnj177NsCAgLkiTA0NLSqqkpurKiokCQpNDTUvmebNm0qKyttn1ZUVISFhd3gu9nubnj00UdnzZp14MABTIQAAKCyxMTEpUuXcr/Up0+fvLy85557jjGWl5cXHx9/zR0UvXv3njZtWm1trb+/f01NTWFhoTyD3lRlZeWlS5eCg4MxEd4KgT1OqVqj1Aaq3DQplRqlM5kCQTtqJ1sjkanjZu2orCaZJlWk1ih3JFQUkCwsSb2a3EbBUp481AUFKjVq4W3Ay6hqqMT7jRwMtxoqmQ4l3hLcoq/ERrs6Az+pZzTy71QTTCkLtIu+JSwiO2974Ia9ztQaveGz8cgjj7zxxhvz5s3r3LnzrFmzXnvtNbl98ODB48aNe/zxx7t27Tpw4MDMzMxHH310yZIlI0eO7NChA2PsxIkT2dnZ+/btKysre/311zt27Dhu3LidO3euXbv2jjvuqKurW7RoUUpKSlxcHG5yBwAAxThRYe3G02abNm1+/PHH0tLStWvXvvXWWxMmTJDbx40bFx8fL3+8atWqpKSkL774ol+/fp9//rncaDKZqqurIyMjMzMzq6ur5fIvkZGRvr6+69aty83NnTx58rp16yRJwooQAAAatdjY2OvvBXz88cdtHzdv3nzOnDnXdIiOjp4/f/41jeHh4a+88so1jZgIAQBAOUrfUN8AMBECAIBiJCduqMNE6MaojXn5JdaoBIRApoMqD0YFQ/h5FrIzFdsR2TpYZCSMeETkwxQJFpEpJCoBIVRijduVxv0xp35XULELK1XEi9udLAJHJXQ4URQtEZbx0gm8mmQRQarWHZly4sXHRN9v3J8IV1bdI98rav/qdx0n7iNU/clAWAYAADwaVoQAAKAcXCMEAACP5sR9hGoX6cSpUQAA8GhYEQIAgGKcK7qtLkyEDqDSoVR3flEufmcLP7DGDzGStaCIYB73IGThKDK/KpAmpatVCVSxohO2AiOnnyvBiCA/kym2My+3WSJOx2iodCj1VpF4bxXiCTeZiKppRs5ozEQpNeotwX/Tirz0jH4huMch3xLkqyxwEOp9qER9PWhcMBECAIBiJKX3I2wAuEYIAAAeDStCAABQjDM31OMaIQAANBnObMOk9kyIU6MAAODRsCK8OfGMmMC+nWRCktdO7u6rxEGE0npUfypQRyUe+ZsYi+Za+QlbZWqNcndhpSt5Eruz8v7gJfoK4/4xbbUI7OLLiFfZbOb/oUwlcsWCnVQ7EaQW+4kgf9yERiKw8zb5zuc306++2icJFeBEZRm1YSIEAADl4D5CAADwZBKTRG+HUHsexDVCAADwbFgRAgCAcrD7BFDIa+nEF7gX6skiWyJX9cnOou38SnJi4QV+AkKJPXXp7W25zSTuT6hQdT1qMIqFZXindTSCD5/7KoumSARL+inQTr9puc1UzkXgx4eRz6FY+qWJZ2WwMS8AAIAbwYoQAAAUI0niN8irfUM9JkIAAFAOrhECAIAnc6bWqGtG4jhcIwQAAI+GFaHKBJODYt35wU4yUCc2FKEwHL0vsUC1KurPRn5sUiv2VyZRkoxZedvkCkdSue3kdxQ8ODV0ge8pdpFG7GEK9qWTt9xGsQSn2A+Q0MPEbr3/gaLbAAAAbgYTIQAAeDScGgUAACWh6DYAAHguZ64RumgoDsNEqDI6vMBt5HcXe9eRmROiWeibUgcnQg3cKvUSkXPRePHP5Gv54QV+KIhM4gjWwOMfhGrnDpx6YonrFRrqZVbitwj32BqN2Ai5I6FGTRybSdQ3FXmYYj9WyjyDChyjiXDD+whxjRAAADwaVoQAAKAYd7x9AhMhAAAoxonKMqrDqVEAAPBomAgBAMCj4dRoAyETf8QXuBFBMsUnEuwkU3xkRFCgnXw4xEG4hdC0RGrUaua3W4T+nBPehVVgP1gylch7+ELPCbvBC8frL/SqUYMRza8SB1FgJMIHF/lxox6mIvlVirudOxTgjtcIsSIEAACPhhUhAAAoyd325cVECAAACnLDG+odnQg3bdq0du3a8vLy4ODgiRMn9u/fX26/ePHi/PnzDx8+3L179+eff97f399lQwUAgMbOiWuEqnP0GmFxcXF8fPwjjzwSExMzbNiw/Px8uX38+PFHjx7961//umfPnkmTJrlsnAAAAC7h6Ipw6q0QPYkAACAASURBVNSp8gdpaWm7d+/evHnzgAEDDh48uHPnzvPnz/v5+fXp0ycsLKy4uDgqKsplo1WH8Cqfm0wjC0sKRAeFAnVUuyIHYdR2uNRBqDKhXpyKoGaTYFFNkT1R6VyrwEFEN2EViu9S7Votf4gaL17ylnrCqTqu3LcK1ZnKtfKDxNy+4qlRbsBYMHkr9BNBh7EFfsA9kBM31Kv+5ImlRs1m84EDB3799dfk5GTG2O7du3v06OHn58cYCwwMjIuLKygocMkwAQDALUhO/VOVwES4dOnSli1bJiQkpKenDx48mDFWXl7eqlUrW4fg4ODy8nLqv1ss/N0AAADALTTVX+MCE+GUKVMuX75cXFy8Y8eOt956izHm7+9fX19v66DX628QltEInXsCAIBGxpFf4xJj0p/nRwX+NcDgb0B4curQocNDDz20ZcsWxljbtm1Pnjxp+9KpU6fatWun4OAAAMC9SMzdpkHHwzJnzpyRJ7n6+vpNmzZ1796dMXbPPfdkZWXt2rXrzjvv3LJlS11d3cCBA104WLUoktGgqlWJ5BG4sYgbtHMTE15kbkUsd8DNblAH0en439Rs4rRLGv7pF6EaY16Slt+bLKXG727hVVOjSqxZqKMoEuchXgjuC+pFPOFUu+BbRaBdQ2V8yLAMt5n6iaAyQQJvcvFMEH+E/M7CX3B/zt1HKJpAU5SjE+GgQYN8fX2DgoL++OOP+Pj4l156iTEWEBDw9ttvjxw5MiEhYf/+/f/61798fHxcOVoAAACFOToRHjly5MiRIzU1Ne3atWvbtq2tPSsra+TIkcePH4+JiQkKCnLNIAEAwE04cUO9JLnHilCr1cbGxnK/FBISEhISotyQAADAXTlzH6FrRuI4JDkBAMCjoeg2AAAoxqn9CF00FkdhIrwFIiWyqHcGVfJKKMVHtvMigl7eYnlCoeigzpuf1bSY+VcAtEZe9FQkHsmoXCIR4OSmQBljZmKEVl47FQ6lcLvTu/vykfFdkdQoHfjkxinF3io63luLCgwLRU+pdvLNKfJNqe8oVHqN/FWu9q94cAQmQgAAUIwT1whV/3MB1wgBAMCjYUUIAACKkZj4NUK1l4SYCAEAQDluuEM9To0CAIBHw4rwVlChMk6j0HajTLSAJBEE5af4iGCnF9ludvzgFrNATVHGmBcvrEg9TOp8i5UX4qRSoMZ6/sMha5Dye/NRm816cV9l4q9gMk1KbinMCzGKbKdMDUaouCsTzK9y3z/sBm9mRSLQ3IOIVt/l7z+M2Oif3HFjXkyEAACgHCdKrKn95wImQgAAUI4LrhHq9fp169ZdunRp6NCh7du35/bZsWPH4cOHu3bt2r9/f1vjxYsXDx486O3t3bt3b/vOu3bt2r9/f8eOHVNTUxmuEQIAQGNWV1fXv3//ZcuW7du3LykpadeuXdf3+dvf/vbYY48dPnx4woQJ8+bNkxs//PDDkJCQjIyMadOm2XdeuHDh+PHj//jjj6eeeurpp59mWBECAICCFF8QZmdnW63WH374QavVtm/f/h//+Md3331n3+H8+fPvvPNOYWFhVFTU5MmTe/XqNW3atJYtW44bN27ixIlff/31okWLbJ1ra2vnzZu3adOmXr16Pffcc506dZo5cyYmQhfg5g6I8IJYiTUiAuBN5Fy8fTjt3MYbHMRAtJu8OdvnUhEVLyN/r10vHW+Evty+ZJ02i5lzcAMRirEQyQgqQ8RNRhgN/IMbiYfJffWpLAa5uS8RouGX9OMfg2moxA23bBiV8BKpF0hVO6NDNPwXQsd735KdRdrJcnREhTn+06J23KPxkG4UHSL+yw27f//996NGjdJqtYyx0aNHz54922QyeXn9/+S1devW2NjYqKgoxlhsbGy7du1yc3PT09NbtWp1/dF27twZEBDQq1cvxlh4eHjPnj03btyIU6MAANB4lZaWhoeHyx9HRESYzeZz585RHRhj4eHhpaWl1NHOnj17fWesCAEAQDlOnRstLCx84YUX7Nv69u2bnp7OGLNYLLYlpvyBxfJfp16sVqv9GlSj0VzT4aadMRECAIBinCm6zZhOpwsMDLRvadGihfxBeHi4bQlYXl6u0WhCQ0Pte4aFhdmvEcvLy+3XfNe4vvPAgQMxEQIAgGKc2o9QiomJmTVrFverKSkpK1asmD17tiRJP/zwQ3Jysk6nY4xduHChWbNmAQEBAwcOfOyxx8rLy0NDQ0+fPn3s2DH7Oyiu0adPnwsXLhQVFcXGxlZXV//yyy9Lly7FNUIAAGi8HnrooXPnzmVmZs6dO/ell16aPXu23D5y5Mhly5Yxxtq1azdp0qQRI0a8/vrrI0eOnDp1akhICGPs0KFDU6ZM+fTTT0+ePDllypT333+fMdayZcunn3567Nixr7/++vDhw8eOHRsdHY0VofPIP3p44T4N8SeH0N6nVBbO25dq57y+3kSc0rueOAjRbjZxzsKT294ScUqzidOfKgNGxSaNBoH6YC2C+O/50Nubc9v9b/O+vrGirJbb+ezJy9z2ulojZ4BEJpNChWa5uVbRtwq/Gh8RMNb5iLxpBdOhZKqZG4Em3/kCD5PKtVKV5Lg/y9RvA9WLh6lA6fsnAgICfvnll5UrV1ZXV2/fvj0hIUFunzt3bmRkpPzx4sWL16xZc+jQoblz544ePVpuDAoKSk1NlW+ZZ4yFhYXJH7z66qt33nnn3r17n3766YyMDIb7CAEAQFGS6LZKN+3fqlWrJ5988prGYcOG2T7WaDTjx4+/pkObNm2ub5SlpaWlpaX9/38XGCwAAECTgxUhAAAoRnLDE8KYCAEAQDluuDEvJkIX4FarEtzajZs7IHMEvFAMY8ynGScXYzTwOxsN/DwL1c7NxZBhGV4VNEZsJUiFF6i0iNaL8zCb+fOf2LYdb+O23zEwgtveOtz/+saj+yu4nQu28+tZlJ/mhGiocnTc54QxZiVuEeY+Xc38+K9yM38dt52bLqF2+yNL+jXjfFMfXiNjzKeZAskvH6Iz1c6N/1D7EVI/m8S9AWr/Lm80nLt9wkWDcRCuEQIAgEfDihAAABTjzA71rhmJ4zARAgCActzwGiFOjQIAgEfDihAAABTjRFhG9RUhJkLlEVEzfmdufSxGVasiUqNURs7EC4iaiGpnVDu3lBojwo2iiUduRQkT9R3JSCp3y1b+2Y6gUE4KlDHWpm0At52bGq0+r+d2bhHM31P4ck399Y3U1sHUwyT26+XXKvNrzk+H+omkRrm7QzP6ueUGO32J1KivH38kvkTeldtORVLpEmsCG/NSP5tCOyF7IDc8M4pTowAA4NmwIgQAAEWpvsQThIkQAACU44Y31GMiBAAAxThzH6HaK0hcIwQAAI+GFaELiMRGqQ17tbw6h1Raj8rOceOXJt5euIxOh4qWxBTCffhGIr9KfkPeF6goIFX1lFsOlBGBT2pjXurvWr8ATkKSGiFVUpU6OjdLTGUvqeAxt66m1ouINBO1RrnpZWokVDVUoTQpFUn1Iarv8jfmJTfgFdlsV+01TSPiRGxUbZgIAQBAMZIkfs1P7XOjODUKAAAeDStCAABQjBMb86q9IMRECAAAypGYxK0YdcP/ojJMhA2EfKWJk9Pc8AIVr6ASEL4mzutLZTEsRIrEwi8ExkddG6ByB9x2L2IrYKr2GPfJpb6joY7/eEqOX+S2a7Sc57z2koHbmUrzNOOGZYjok5lIM1HPLfc41FuFijiZTZxGjYYaCbeZX5KNqnbmQ24drECJNeongvu00BvwcpvV/63d2LlhjTVcIwQAAI+GFSEAAChH/IZ61WEiBAAAxTixDZMkSQrclXwLHJ0Iq6qqNm/efPr06dtvvz09Pd3X98/tZqxW67fffltUVJSQkDBs2DCXjRMAAMAlHL1G2KNHj5UrV1ZWVn7wwQeJiYk1NTVy+1NPPTV79uyrV69Onz79xRdfdNk4AQDAHUhO/VOVoyvC3bt3BwcHM8bMZnO3bt2++eabrKys0tLSjz/++NixYxERERMmTIiPj58xY0ZQUJArB+y2qFfaSiTWeH+i0BulUtvhcl5fqlCZaDv35AeV1dTyspeMqG5VTwQ7qepoXNTJGaq97iovN8mYhfc9TUb+CHVUVpO3HS71alqIcCz53PKOQ52aogrp8b8j8aoJpZe5u/WyG23AK1BizacZtQEv9X7jPVdCpdSY+r+1Gzlnim7TieuG4eiKUJ4FGWNardbb21ur1TLGtm/f3r1794iICMZYp06dIiMj8/PzXTRQAABo/KQ/p0KBf6r/cSF8+0R2dnZFRUV6ejpj7OzZs2FhYbYvhYaGnj17lvqPitRoBgAAtTTVX+NiqdHc3Nxp06b9+9//btmyJZOjPnbPi9VqVX1/RQAAUJm73VAvMBHu3LkzIyPjq6++6tu3r9wSFhZWVlZm61BeXh4eHk79d8yRAABuzZFf4015Y969e/eOHTv2k08+SUlJsTWmpKT8/vvvJSUljLEjR46cPn06OTnZJcMEAABwDUdXhCNGjPD391++fPny5csZY2PGjHnggQfCw8OnTJkydOjQMWPGZGdnz5gxIzAw0JWjbYLoP4U4X6B28aVSfNwklmg6lMLN2mmorVyJdq2OV2vUm5/JNBr47eROtkqQ+BlO/hNOxgx5zxW1HyyFSplyj0O9mlTylhtJpYKdfs29ue3+vHa/5vwUKHezYnaDWqO8sqJUTVHqueKWFSVfNbWXKW7KmRvqXTQUhzk6EX7wwQdm8///GoqNjZU/eOeddzZu3FhYWLh06VL7xSIAAHgiNyy67ehEOHbsWOpL99xzzz333KPQeAAAwJ05UWtU7YkQu08AAIBHQ9FtAABQjCR+j4DoRr6Kw0TYSPHfSFSRLWplz48jiHxHwXYqAEIW5fLmFeXy4Vc7o/bUNRo5ARALUUuM2tyXulOYuxswuc8w8fC52Q0rUeyMW12POsgNBsNFPUzuC8HdT5gx5serGMeInEszMllDHYTf7s0Ny+iIsIzIXruq/xYG1eHUKAAAeDSsCAEAQDHO7UfoosE4CBMhAAAoxonKMqqfnMZECAAAylF/MwlhuEYIAAAeDStCd0KecCD+ntHwcoZEsTPyjzjqm3LDivTWwfx2A6+amrcvPwpoqCdSo7x2bpSU0TvTUnXauNvkCnWmUNdFqOgplQ7lpkypJCR1EO5zTqVGqSpoPrxgJ7kBL68zdRDGmBcv16ol3s1kkpaIjYKCnLpG6KKxOAoTIQAAKKcJ7z4BAADQJGEiBAAAj4ZTowAAoBiJObENE+4jhFtGJyM42Q2JPA1AVaWiMhrcsAy/MxWW0fnwwjL1/DwLtR+hiZeLEerMGDObiPwLbwM/C1Edjb8DJONfAyHzHFQohnwhOI1UWoRKM3HDMqI5F29fTrs3sWUg1U5V4+O+tajnirropPq1KE/gjvcR4tQoAAB4NKwIAQBAOU14Y14AAICbcsdaozg1CgAAHg0rQgAAUIwTYRm1z4xiImwaqOiclfcFXpSUMaYlU6MC+8FqtfzOQhvzGnlRUsaYycjPGXKDoFQpNRORDhUqvcbdrZeRmVE+esdjscQjtyQb9UJoyReC085NgTI68Mk9CPXSe1H7DFN76vLeb6LPFTQQd7tGiFOjAADg0bAiBAAAJal+g7woTIQAAKAYZ64Rqj1vYiIEAADFOHP7hNorSFwjBAAAj4YVYZPG34WUSh8SzUJbnGr4eUKqBqmXjpPV5EZJ2Y2CoLzUKLV3LnEQsj8vIEqmRkVio8KpUeIF4gZEtUT2kqo1yk2TkoFPshyowEjIfYbJjaO5qVGiM6gIlWUAAMCTuaLodklJyUcffVRTUzN69Oi77777+g4Gg2HZsmWFhYUJCQlZWVleXn9ObYWFhZ9//rnFYnnooYcSEhLkxvfee0+v18sfd+nSJT09HadGAQCg8aqpqenTp09VVVVMTExGRsa///3v6/tMmjRp1apVSUlJy5cvnzJlitz4xx9/9OvXLyAgICgoKDk5ef/+/XL7K6+8UlxcXF1dXV1dXVtby7AiBAAABSm+H+Hy5cu7dOmyaNEixpivr+/8+fPT09PtOxw/fnzt2rWlpaWtWrUaOXJkZGTkK6+8EhER8d5772VmZr700kuMsaqqqrfffvuzzz6T/8vMmTOjo6NtR8CKEAAAGq/8/PzU1FT548GDB//yyy8Gg8G+w08//ZSYmNiqVSvGWGhoaJcuXX7++edr/mNqamp+fr7tvyxfvnz+/Pnbt2+XP8WKEP5E/w1HhBp4iRYq0KHRUPW0ePEKXoKGMWYx8w/CrYJGhl+Idir/wg/LEKkYwbCMWCiGSpdw28nOVHSF10511hKvpsTrL5wJ4ndXP08BDlL8PsKysrLWrVvLH4eEhFit1rKysvbt23M7yH3KysoYY+Xl5fb/8ezZs/LHycnJkiRVVlY+8MADGRkZ7733HiZCAABQ2U8//TRkyBD7llGjRj311FOMMZ1OZzKZ5Eaj0cgY8/Hxse/p7e1t6yD38fb2vuY/GgwG2//65ptv5A+mTJnSpUuX6dOnYyIEAADliN9QzyQWHR09a9Ys+7aoqCj5g4iIiNLSUvnjkpISnU5nv/5jjIWHh5eUlNg+LSkpiYiIkP+jrb20tDQ8PPyab9upU6egoKBTp07hGiEAAChHEv/HWEhISOp/s02E6enpOTk58nXB7OzskSNHarVaxtgvv/xy6tQpxtjQoUOPHz9+6NAhxtjevXvPnz8/aNAgxtioUaOys7Plg6xatUqO2Oj1eovlz4svP/74Y3V1dZcuXbAiBACAxmvMmDGLFy++6667IiMj8/LyNm/eLLf/z//8zwMPPDB9+vTAwMA5c+YMGTIkJSVly5Ytr776akBAAGPsiSee+N///d/U1FSdTnf06NF33nmHMbZ9+/Zp06YlJibW1dXl5+cvXLgwLCxMoi77K0uv1z/5yoKRD/+1Ab4XKIt+g3BTJERXMorCbRTYGpBqR1iG3xlhGbgFYzuH3bTP59uPX7hYJ3TY33/ZWvPH9pycHKqD2WzOy8u7ePFicnKynA5ljB05ciQwMNB2mvTQoUOFhYXx8fGdO3e2/ce6urrt27dbrdZBgwb5+fkxxiwWS2Fh4bFjx3x9fRMTE8PCwhhSo3BT5Nl+7q6/1G89ftE0/i99jZXf26qlph9uFTT+d+RObNRBGDG3kX8+UhMh96kS3FRWqL/oLMvtz90Ll/qOdDv2zvUszhTdvll/rVZ7fUGZmJgY+0/j4uLi4uKu6ePr6zt8+HD7Fo1G07179+7du9s3YiIEAABFudtfOQjLAACAR8OKEAAAFOPE5hOqw0QIAACKccU1QlfDqVEAAPBoWBGCswRCozf8yvVdya1Z+c1W7heoXCeVDqVGw/2CaGhUib93hQ5C/okt1Iw7HMAJ2JgXAAA8miR8qlP1U6OYCAEAQDHO7D7hmpE4TuwaYXl5eXl5+TWNZWVlO3bsOHfunHKjAgAAaCCOToQffvhhWFhYeHj4Y489dk179+7d58+f37Vr1y+//NIFIwQAALfiVN1tFTl6arR///7btm1bv359bm6urfHy5cszZ87csmVL7969t23bdv/99993333X7BQFcAOilxJuvS95CCVq7godQ/hn33XPFQqegXIkJkmqz2yCHF0RxsXFxcbGXrPP+MaNG9u1a9e7d2/GWEpKip+f344dOxQfIgAAgOvcUljmzJkzkZGRtk/bt29/+vTpWx0RAAC4LWfCMmovIG9pIqyrq/P29rZ96uPjo9frqc7U3joAAOAWLBaLhtiE6/81gmt+om5pIgwNDa2srLR9WlFRIe/txHXzpw8AABoxR36Ne1yJtd69e+/Zs6e2tpYxVlVVVVRU1KtXL4UGBgAA0BAcXREWFhauX78+Pz//+PHjr7/+enx8/PDhw7t27ZqcnJyZmZmVlbV48eJRo0bZXzIEUJhoEFSRowhtOi82EkHudroJPJMT1whVf287uiI0GAzV1dXdunVLT0+vrq6WV4GMsezs7KSkpBUrViQnJy9fvtxl4wQAAHfgbjcRMsdXhElJSUlJSde3N2/efM6cOYoOCQAAoOGg1igAACjI/cIymAgBAEAxbrgLEyZCgBtT/WcUwL244UyIe/sAAMCjYUUIAACKceKGetVhIgQAAMW4Y61RnBoFAACPhhUhAAAoRpKEb4dQ/VQqJkIAAFCU2qc6ReHUKAAAeDSsCAEAQEnutiDERAgAAMpxx/0IMRECAIByUFkGAADAvWBFCAAAinHHG+oxEQIAgGLcscQaTo0CAIBHw4oQAACU44ZhGUyEAACgGEn8mp/qZ1IxEQIAgGIkJkmCSzzR/orDNUIAAPBoWBECAIBynLhGqDZMhAAAoBzx+whVnzhxahQAADwaVoQAAKAYFN0GAADPhmuEAADgySTmxIrQRWNxFK4RAgCAR8OKEAAAFOOGFdYwEQIAgILccCbEqVEAAPBoWBECAIBiJMmJ2yFw+wQAADQZ2KEeAAA8Gq4RAgAAuBesCAEAQDHO3FCv9pIQEyEAAChGwu4TAAAA7gUTIQAAeDScGgUAAMU4tQ2Ti8biKEyEAACgGCeuEao9D+LUKAAAeDasCAEAQDlueEM9JkIAAFCO+DVC1S8S4tQoAAB4NKwIAQBAMZL4Ak/tBSEmQgAAUI4zt0+ofZGwgSZCLy+vP3I35mZ/3jDfroFVVlYGBgZqNE38PDMeZlNSUVERFBQkvm+cm8HDVFbs+vWxsbE37jM4Mlj0sFH6pJ81dc4OSgGS1WptmO9UVVVVU1PTMN8LAAAU17ZtW29vb7VHobyGmwgBAAAaoSZ+/gcAAODGMBECAIBHw0QIAAAeDRMhAAB4NNxHKOzcuXMFBQWlpaWDBw/u2LGjrf3UqVOfffZZbW3t+PHje/XqpeIIFbF3795NmzZduHAhNjY2MzOzWbNmcvulS5c++uij0tLSQYMGjRo1St1B3rotW7bs3Lmzpqbm9ttvnzBhQlBQkNxeU1Pz0UcflZWVpaamjhgxQt1BKig7O9vHxyc9PV3+tL6+/uOPPz527FhSUlJmZqa73zSybt268vJy+eNWrVqNGzdO/riqqmrZsmXl5eXDhg0bOnSoegNUzPnz55cvX3727NkOHTpMmjSpRYsWzO5NO3jw4HvvvVftMboT937fq2LAgAGvvfbarFmzCgoKbI3l5eW9evWqqakJCQkZMmRIXl6eiiO8dTU1NWlpaRcuXLj99ttXrFgxYMCA+vp6xpjFYklJSdm5c2fHjh2nT5/+7rvvqj3SW7Vq1SqLxRIVFfXjjz8mJiZWVVUxxkwmU3JyckFBQVRU1OOPP75kyRK1h6mMdevWTZ48ecGCBbaWv/zlL6tXr46Ojn7rrbdmzpyp4tgUsWDBgk2bNhUXFxcXF5eUlMiNBoOhf//+Bw4c6NChQ1ZW1meffabqGBVw9OjR+Pj4wsLCyMjII0eOyI/UZDINHDiwoKCgY8eOTz755OLFi9UepluxgiCz2Wy1WhMSEr766itb49y5c0ePHi1/vGDBghEjRqgzOIWYzeb6+nr5Y71e36JFi7y8PKvVumHDhsjISKPRaLVat27dGh4eLn/cBFgslg4dOuTk5Fit1pycnJiYGPmFlh+y/LFbq6mp6dq169y5c/v16ye3/P777/7+/pcuXbJarcePH2/WrFllZaWqY7xVd91117p1665pXLlyZbdu3SwWi9VqzcnJiY6Olj92X8OGDXvxxRevaZQfmvxG/f7779u3b28ymdQYnVvCilAY9/RRXl6e7ZTLkCFDcnNzG3ZQCtNoNLbbZi0Wi8FgaN68OWMsNzf37rvv9vLyYowNHDiwsrLyyJEjag5UOUeOHKmpqenSpQtjLDc3d/DgwfILPXjw4NOnT588eVLl8d2yZ5555plnngkPD7e15OXl9enTR35lo6KiIiIifv31V/UGqIzvv/9+4cKFGzZssP7nDum8vLzU1FS56srQoUOPHj169uxZVcd4S4xG4+bNm9PT0z/55JMlS5bYFr7XvGlLSkqawJu2wWAiVEZZWVnr1q3lj0NCQmpray9duqTukJTy3HPPJScnJyYmMsbKy8ttD1Or1QYFBZWVlak6OgU8//zzERER8fHxb775pjwR2r+a3t7egYGB7v4wt27deuLEiaysLPtG+1eTMRYSEuLWMwRjLC4uzsfH5/z580899VRaWprFYmH//Wr6+fkFBAS49at55swZi8XyxBNPnDx58uDBgwkJCYcOHWL//WrqdLom8KZtSAjLKMPLy8tkMskfyx/odDpVR6SMd955Z/PmzbZLnl5eXmaz2fZVo9HYBOotzZkz59lnn83Pz586dWr37t179eql0+ma0sOsra196qmnvvnmm2tqUTa9V/PDDz+UP3jhhRdiYmI2bdo0bNgw+59NxpjJZHLrh6nRaKxW6xNPPCH/WWM0GhcuXPjxxx83vVezIWFFqIyIiAjbX9OlpaWtWrWyxSzd13vvvff+++9v27YtNDRUbomIiCgtLZU/1uv11dXV9qfa3JS/v39oaOj48eOHDRu2du1a9t8P8/Lly5cvX3brh5mbm1taWvrQQw/17Nlz3rx5Bw4c6Nmzp8VisX+YjLHS0lK3fpj2AgMD4+LiTpw4wf77Z7OioqKurs6tH2ZYWJhGo4mLi5M/7dq166lTp9h1b9pLly659cNsYJgIlZGWlrZmzRr5VMzXX3+dlpam9ohu1bJly956663Nmze3bdvW1piWlrZ582a5eHpOTk7nzp3tbyBxOyaTyWg0yh8bDIYDBw7cfvvtjLG0tLQffvjh8uXLjLHVq1cnJSVFRESoOdBbc9ddd23btm3p0qVLly6dMGFCVFTU0qVLNRrNsGHD9u3bJ/8a/emnn+rr6/v166f2YJ1nNBptK78zZ8789ttvXbt2ZYylpaVtUBWcRgAAAkRJREFU2LDh6tWrjLHVq1f37ds3OFh4e4TGw8fHZ/jw4bt27ZI/3bVrlzwppqWlbdy4UX7TrlmzJjEx0f4nF25C7bSO+3nyySd79OjRrFmzqKioHj16FBQUWK3WK1eu3HHHHcnJyRkZGW3atDl8+LDaw7wlpaWlkiTdfvvtPf7ju+++k7+UmZkZFxf38MMPBwcHf/vtt+qO8xadPn26TZs2o0ePzszMbN++fWpqql6vl780bty4bt26TZo0KTg4eOPGjeqOU0EfffSRLTVqtVpffPHFyMjIrKys0NDQpUuXqjiwW3f8+PGwsLCxY8dmZGQEBgY+8cQTcrvFYklLS0tMTJw4cWJQUND27dtVHaYC9u7dGxISMnHixBEjRnTq1Ons2bNy+/jx421v2h9++EHdQboX7D4h7OjRo/ZBmJiYGDl3V19fv23btitXrqSmpgYGBqo3QAUYDIaDBw/at0RGRso3m1ut1vz8/NLS0n79+rVv316lASrm9OnTv/32W11dXXR0dFJSkq3darXm5uaWl5ffdddd7dq1U3GEyqqoqKisrOzcubOtpaCg4OjRo4mJiTfdZ67xKyoqKioqslgs8fHxMTExtnaLxbJjx47z588PGDDArRf3NpWVldu2bWvZsmX//v1tV2GsVmteXl5ZWVkTe9M2AEyEAADg0XCNEAAAPBomQgAA8GiYCAEAwKNhIgQAAI+GiRAAADwaJkIAAPBomAgBAMCjYSIEAACPhokQAAA8GiZCAADwaJgIAQDAo/0fA/cSuHufxIMAAAAASUVORK5CYII=",
"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.8.1"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.1",
"language": "julia"
}
},
"nbformat": 4
}