{
"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.060136e+01 7.595435e+00\n",
" * time: 0.004382133483886719\n",
" 1 2.880997e+01 5.423184e+00\n",
" * time: 0.013163089752197266\n",
" 2 1.674728e+01 3.684839e+00\n",
" * time: 0.038455963134765625\n",
" 3 1.298976e+01 1.842488e+00\n",
" * time: 0.05840802192687988\n",
" 4 1.157225e+01 2.066631e+00\n",
" * time: 0.07877898216247559\n",
" 5 1.067315e+01 1.395294e+00\n",
" * time: 0.09610199928283691\n",
" 6 1.022136e+01 1.297891e+00\n",
" * time: 0.11374402046203613\n",
" 7 9.733602e+00 7.906879e-01\n",
" * time: 0.1299431324005127\n",
" 8 9.498805e+00 9.571391e-01\n",
" * time: 0.1461939811706543\n",
" 9 9.376034e+00 6.080652e-01\n",
" * time: 0.16211605072021484\n",
" 10 9.256797e+00 4.508729e-01\n",
" * time: 0.1778249740600586\n",
" 11 9.159711e+00 6.484503e-01\n",
" * time: 0.18931198120117188\n",
" 12 9.101580e+00 6.137670e-01\n",
" * time: 0.20073914527893066\n",
" 13 9.062695e+00 4.728745e-01\n",
" * time: 0.21219301223754883\n",
" 14 9.003882e+00 4.749752e-01\n",
" * time: 0.22366809844970703\n",
" 15 8.961750e+00 2.891135e-01\n",
" * time: 0.23530101776123047\n",
" 16 8.940465e+00 3.054697e-01\n",
" * time: 0.2867159843444824\n",
" 17 8.902345e+00 2.605803e-01\n",
" * time: 0.29808616638183594\n",
" 18 8.882018e+00 1.942450e-01\n",
" * time: 0.3094971179962158\n",
" 19 8.869787e+00 2.098498e-01\n",
" * time: 0.3207430839538574\n",
" 20 8.860407e+00 2.653612e-01\n",
" * time: 0.33221006393432617\n",
" 21 8.834069e+00 2.061801e-01\n",
" * time: 0.3436880111694336\n",
" 22 8.792125e+00 2.805051e-01\n",
" * time: 0.355147123336792\n",
" 23 8.756360e+00 2.548108e-01\n",
" * time: 0.36669206619262695\n",
" 24 8.729760e+00 2.118028e-01\n",
" * time: 0.3781859874725342\n",
" 25 8.674094e+00 2.337216e-01\n",
" * time: 0.389707088470459\n",
" 26 8.639844e+00 2.488812e-01\n",
" * time: 0.40123605728149414\n",
" 27 8.610054e+00 2.194887e-01\n",
" * time: 0.4128580093383789\n",
" 28 8.571872e+00 1.824702e-01\n",
" * time: 0.424389123916626\n",
" 29 8.551868e+00 1.659673e-01\n",
" * time: 0.43578314781188965\n",
" 30 8.541628e+00 1.793289e-01\n",
" * time: 0.44706201553344727\n",
" 31 8.535471e+00 1.925776e-01\n",
" * time: 0.45818209648132324\n",
" 32 8.530846e+00 1.393299e-01\n",
" * time: 0.4694051742553711\n",
" 33 8.526631e+00 1.107636e-01\n",
" * time: 0.48075199127197266\n",
" 34 8.524832e+00 8.344340e-02\n",
" * time: 0.49204111099243164\n",
" 35 8.522520e+00 6.546595e-02\n",
" * time: 0.5034000873565674\n",
" 36 8.519161e+00 5.673858e-02\n",
" * time: 0.5149180889129639\n",
" 37 8.517701e+00 5.615832e-02\n",
" * time: 0.5263099670410156\n",
" 38 8.516916e+00 5.937991e-02\n",
" * time: 0.5379111766815186\n",
" 39 8.515851e+00 3.801164e-02\n",
" * time: 0.5496761798858643\n",
" 40 8.515053e+00 5.264755e-02\n",
" * time: 0.5936429500579834\n",
" 41 8.514499e+00 3.867559e-02\n",
" * time: 0.6050460338592529\n",
" 42 8.514251e+00 5.469878e-02\n",
" * time: 0.6164431571960449\n",
" 43 8.513826e+00 3.737006e-02\n",
" * time: 0.6277720928192139\n",
" 44 8.513366e+00 4.034809e-02\n",
" * time: 0.6390640735626221\n",
" 45 8.513194e+00 2.816864e-02\n",
" * time: 0.6503660678863525\n",
" 46 8.513045e+00 3.915307e-02\n",
" * time: 0.6616780757904053\n",
" 47 8.512761e+00 2.934585e-02\n",
" * time: 0.673194169998169\n",
" 48 8.512589e+00 2.754424e-02\n",
" * time: 0.6845731735229492\n",
" 49 8.512434e+00 1.678330e-02\n",
" * time: 0.6960961818695068\n",
" 50 8.512323e+00 1.595735e-02\n",
" * time: 0.7073981761932373\n",
" 51 8.512243e+00 1.906035e-02\n",
" * time: 0.7187399864196777\n",
" 52 8.512171e+00 1.109253e-02\n",
" * time: 0.730100154876709\n",
" 53 8.512089e+00 1.465525e-02\n",
" * time: 0.7413949966430664\n",
" 54 8.512037e+00 9.566615e-03\n",
" * time: 0.7527220249176025\n",
" 55 8.512024e+00 1.233504e-02\n",
" * time: 0.763901948928833\n",
" 56 8.511974e+00 9.108683e-03\n",
" * time: 0.7753121852874756\n",
" 57 8.511942e+00 7.699279e-03\n",
" * time: 0.7867450714111328\n",
" 58 8.511924e+00 9.541852e-03\n",
" * time: 0.7980160713195801\n",
" 59 8.511905e+00 7.157047e-03\n",
" * time: 0.809149980545044\n",
" 60 8.511891e+00 6.120977e-03\n",
" * time: 0.820289134979248\n",
" 61 8.511875e+00 4.611738e-03\n",
" * time: 0.8314640522003174\n",
" 62 8.511867e+00 4.805000e-03\n",
" * time: 0.8427109718322754\n",
" 63 8.511863e+00 3.933121e-03\n",
" * time: 0.8542051315307617\n",
" 64 8.511858e+00 2.674432e-03\n",
" * time: 0.8937420845031738\n",
" 65 8.511855e+00 2.280951e-03\n",
" * time: 0.90513014793396\n",
" 66 8.511852e+00 1.941317e-03\n",
" * time: 0.9165301322937012\n",
" 67 8.511851e+00 2.158678e-03\n",
" * time: 0.9278850555419922\n",
" 68 8.511850e+00 1.499713e-03\n",
" * time: 0.9392681121826172\n",
" 69 8.511849e+00 1.661873e-03\n",
" * time: 0.950592041015625\n",
" 70 8.511848e+00 1.873279e-03\n",
" * time: 0.9620001316070557\n",
" 71 8.511848e+00 1.682758e-03\n",
" * time: 0.9737610816955566\n",
" 72 8.511847e+00 1.516782e-03\n",
" * time: 0.9852659702301025\n",
" 73 8.511846e+00 7.989321e-04\n",
" * time: 0.9965569972991943\n",
" 74 8.511846e+00 7.280981e-04\n",
" * time: 1.0077309608459473\n",
" 75 8.511846e+00 8.588832e-04\n",
" * time: 1.019024133682251\n",
" 76 8.511846e+00 5.718639e-04\n",
" * time: 1.0302910804748535\n",
" 77 8.511846e+00 4.213009e-04\n",
" * time: 1.0414800643920898\n",
" 78 8.511845e+00 4.633857e-04\n",
" * time: 1.0527970790863037\n",
" 79 8.511845e+00 4.753988e-04\n",
" * time: 1.063992977142334\n",
" 80 8.511845e+00 3.306843e-04\n",
" * time: 1.0752429962158203\n",
" 81 8.511845e+00 2.141433e-04\n",
" * time: 1.086557149887085\n",
" 82 8.511845e+00 2.072007e-04\n",
" * time: 1.0977730751037598\n",
" 83 8.511845e+00 2.727225e-04\n",
" * time: 1.108877182006836\n",
" 84 8.511845e+00 2.249107e-04\n",
" * time: 1.1201529502868652\n",
" 85 8.511845e+00 1.443985e-04\n",
" * time: 1.1314470767974854\n",
" 86 8.511845e+00 1.298849e-04\n",
" * time: 1.1428041458129883\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3daXxTxd4H8DlpulG2LtCNQtkKbYEWi+wWLGURKCBLXQoCVS4oKC6IiPiAuyCuXBFEVBRcAC0iirJTvMpDkR0qFMpaWpbuQJekyfPi3JsnF+ZfM2HS0zS/76cv0un0ZJKcZDLn/M6MYjabGQAAgKvSad0AAAAALaEjBAAAl4aOEAAAXBo6QgAAcGnoCAEAwKWhIwQAAJeGjhAAAFwaOkIAAHBp6AgBAMCl6bVuAAAAQHUyMzMXLVpUXFw8YsSIMWPG3FqhtLT03XffPXLkSGxs7JNPPlmvXj3G2J9//pmWlnbq1KmgoKBHHnkkOjparfzyyy/fuHFDvd2pU6cHH3ywhjpCk8m079CRtpGRNXN3NcxkMul0dX9sjYdZl+Bh1iU19jAbebr/bZ3jxwvKyoxCmy0tLQ4J0bVu3Zr718uXL/fu3Xv69Ol9+vR56qmnKioqxo4de1OdBx54QKfTPfzww0uWLDlw4MDq1asZY3PmzImLixs+fPihQ4e6du26e/fujh07MsY++OCDiRMnBgQEMMZ8fHwYY0rNzDVaVlY29eUFQyf8owbuCwAAHGFku+C/rRMbu+LgwSuCGz4yYkRxWloa92/z58/ftWvXhg0bGGOrVq1auHDh/v37rSscO3asS5culy9frl+/fmFhYXBwcGZmZsuWLc1ms6Ioap1hw4bFxcXNnTuXMRYQEPDHH3+0bdvWsoW6/10JAACc1//+7//Gx8ert+Pj4w8ePFheXm5dYc+ePXfccUf9+vUZY76+vtHR0RkZGYwxSy/IGMvNzQ0MDLT8+vbbb8+YMWP16tUmk4mhIwQAAIkUxZ6fauTm5vr7+6u3AwICzGZzXl6edYW8vDw/Pz/LrwEBAbm5udYVFi9eXFhY+NBDD6m/jhw5MioqKjAw8IUXXlCPsiIsAwAAGtuxY0eXLl2sS0aNGvX8888zxurVq1dRUaEWqjfULIyFt7d3ZWWl5dfy8nLrCmvXrn3llVe2bt1qKfz444/VGw8++GB4ePjcuXPREQIAgDSKoijVD/FuYTYrsbGxCxcutC5s1qyZ5ca5c+fU22fPnvXy8lJzLtY1LRUYY+fOnQsLC1Nvr1+/ftq0ab/88ktUVNSt9xsaGurr65ubm4tDowAAIJUi/NO4ceO4/2Y5pTd69Oi1a9dev36dMfbFF1+MHDlSjchu3Ljxr7/+YowNGDAgJydnz549jLGdO3eWlpb27duXMbZp06ZJkyb9+OOPsbGxlqYVFRVZTjH++OOP165d69ChA0aEAABQew0ePPjzzz+PiYkJCwvLzs7eunWrWj5v3rwHHnigffv2DRo0WLhw4ZAhQ7p06ZKRkfHee+95eXkxxh599NHy8vL7779frZ+cnPzGG29kZGQkJyd36NChoqIiKyvr448/DggIQEcIAADy/F34RZSbm9t333135MiRoqKiLl26qJ0cY2zDhg3e3t7q7UmTJg0dOvT48eORkZGWoeT27duNxv+/orFBgwaMsf79+2dmZp48edLLy6t9+/Zq1hQdIQAASGPHOUJm/vv6HTp0uKmkSZMm1r8GBwcHB//XZY7NmzfnbiooKCgoKMi6BOcIAQDApaEjBAAAl4ZDowAAIM3fXiB/KxuOjDoWOkIAAJDGjnOEcsM1dsChUQAAcGkYEQIAgDzqNfJOBR0hAABIo4gf6tT80Cg6QgAAkMae6wi17glxjhAAAFwaRoQAACCPHecIcWgUAADqEmfrBwU7wgMHDmRlZTVp0qR79+6WmU8PHjx44sSJqKio6OhoB7QQAADAgWw9R1hVVTVu3Lhhw4atXbt23rx5K1euVMvfeOONIUOG/PTTT4mJiYsWLXJYOwEAwAkoyr/zMgK0HhPaOiL86KOPDh8+fPToUXUlC7PZzBjLz89/9dVX9+7dGxkZuXfv3sTExIkTJ6qrWgAAgCtywnOEto4IV65cOX369Pz8/N27d1+/fl1Nx27evLlt27aRkZGMsS5duvj7++/cudOBjQUAgNrNjgGh8OUWstk6Ijx58uQ333yzZMmSBg0a/PXXXxs2bIiNjc3JyQkLC7PUCQsLu3DhArUFdRAJAABOymw2a95pOYKtI8IbN240aNBg9+7dW7ZsGT9+/IwZMxhjBoPBzc3NUkev1xsMBoc0EwAAtGbLeEZdfULoR3O2doQhISH9+vVTvwv079//4MGDjLHg4OArV65Y6ly+fPmmBYKt1cnvEQAArkOns6HLUOz60ZStHWHfvn1Pnz6t3j59+nRISAhjrHfv3vv27SssLGSM5ebmnjhxokePHg5qKAAAgCPYeo7w6aef7tOnT+PGjRs2bPjaa6+9++67jLHWrVsPHz58xIgRKSkpn3322bhx49QOEgAAXJN6vFPwXzRm64gwKipq165dN27cuHjxYlpa2v3336+Wf/nllw888MDRo0cnTZq0ZMkSh7UTAACcgfg5Qs3PmwnMLNO+fftXX331pkJ3d/cpU6ZIbRIAADitWnDOTxRWnwAAAJeGSbcBAEAa9YJ6wf9xkgvqAQAA/pYd5/y07gdxaBQAAFwbRoQAACCPE4Zl0BECAIA0dkyirfm8Y+gIAQBAGiccEOIcIQAAuDaMCAHESVlSzOm+NgPYwgkX5kVHCAAA8thxjtBBLbEZDo0CAIBLw4gQAACkUcQvkNc6NIqOEAAA5FEU8cshtO4J0RGC8yAiKoLJFQlbkZOVEd0K+Vkh8CEiYxsA1XLCsAzOEQIAgEvDiBAAAKRxwgEhOkIAAJDHninWtO4KcWgUAABcGkaEAAAgjxMeG0VHCFoyCyUnidrcUnLL1EbEWiJUTLzNZSXMeeXicXSR6KnWH1tQm9lxHaHmHSEOjQIAgEvDiBAAAORxwrlG0RECAIA8OEcIAACuTFGcb65RnCMEAACXhhEhyERnLwWymtRGzCZiI7xyoi69EaE0qUh+lZGpUaFiciJjRccp55WRlRljCu8rsZToqeZf9qGGKUwRvUBe8wvq0RECAIA8dpwj1BoOjQIAgEvDiBAAAKSxIyyj+QgSHSEAAEijMPHrCLU+k4yOEP6GUP5FNOdi4uZcqojKRHkVr5zciMlEbJxbzA/RiCVriMAIGX4hynVu/I3rdJwTHDo3/kbciHJufWojOpHEDU0sKwROA+cIAQAAnAtGhAAAII8ifKhT88MA6AgBAEAaJzwyikOjAADg2jAiBAAAeTDpNjgvMvBJJSR56UtuCpQxZqriZzWrjJxyo4G/ESOvMmOsysDbCFWZKKdTpgLhWKHYKHVehMpkkkFQPee4jp5XyBhzc+eXc+vr3QXukTGmcxPIryqKlAnpoNZR7FmGCZdPAABAXWHHCvWaf8tBRwgAALVacXHx119/XVJSMmTIkOjo6FsrmM3mdevWHT9+PDo6OikpyfJfmzZtys7ODgwMHDlyZMOGDS31N27cePDgwdatW48aNUrHvRQXAADATopdP7Tr169369Zty5YtxcXFvXv33r59+611pk2bNnfuXIPB8Nxzz82cOVMtTExM/OKLL4qLi9euXRsZGZmTk6OWz5s374knnqisrJw/f/7DDz/MMCIEAACJ7DhHWP2x0a+++srX13fNmjWKojRt2vS11167++67rSvk5OQsX748Ozs7JCRk7Nix0dHRM2fODAgI+PXXX/38/NQ68fHxX3755axZs0pKSt5+++0//vijQ4cOU6dODQsLmzNnDkaEAABQe23evHnw4MFq5zpkyJDt27cbDAbrCtu3b+/YsWNISAhjrGXLlq1atdq1axdjzNILMsY8PT31ej1j7I8//vD39+/QoQNjzN/fv1u3blu3bsWI0OWIpkOFZv6ksprGSv5UnoZKTn2DSGWqvpEXJWVExJQRE5Yy4mE6dK5R0WlCuUFQPZEOdffgz1jq7sGpL1SZMabn1afyq+REpsSUqlj111nYsfpE9dUvXrzYv39/9XZwcLDJZMrLywsLC7OuEBQUZPk1KCjIchRU9csvv+zfv3/lypXcyhcvXkRHCAAAGjtw4MDkyZOtS+Lj41NSUhhjOp3O8tVTnTTfze2/vi5ZV2CMmc1m6wr79u0bP378qlWrAgMDuZV1Oh06QgAAkMeO6wgV1rhx47i4OOvCdu3aqTdCQkLy8vLU27m5uW5ubk2bNrWuGRISkpuba/k1Nzc3ODhYvX3o0KEhQ4YsXbp04MCBVOV+/fpJ6AgrKytzcnKaNWvm7u5++1sDAABXEx4e/o9//IP7p0GDBn344Ydz5szR6XTr169PTExUz/adPHmycePGAQEBCQkJDz/88Llz55o3b37ixImzZ8/26dOHMXb8+PHBgwe/8847I0aMsGytR48excXF+/fv79y586VLlzIyMr744gtbwzIzZsxQrNy4cUMtX79+fWho6D333BMWFrZly5bbeiYAAMDJqecIhX6qP0l43333VVZWDh069Mknn3z99ddffPFFtTwlJUU97RcUFDRt2rQBAwbMmjVryJAhM2bM8PX1ZYyNHj1aUZS0tLTk5OTk5OQlS5YwxurXrz979ux777131qxZiYmJEydObNGihcCI8IUXXnj11VetSyoqKh5++OEVK1YMHjx49erVqampp0+fvunoLdQEKv/CXTuXHxYhQzFk/oWXOqms4OdcKsuNRDmnPr0RKkTDKTeIhmWoqde4U6wRz6HQpGHUMrbUFGvUxGbcsIy7YFjGw4tT7uEpUJkx5uHFeQKojVBxHr2ZCtdwS7l16Sm7EK6pGVKfZ29v799///2HH34oLCzMyMho3bq1Wv7uu++qSVHG2FtvvTVw4MBjx44tW7asb9++auFHH31UXl5u2U5oaKh647nnnuvVq9e+ffvU/2Ki1xGWlZV5e3tbft20aVPDhg0HDx7MGBs9evQTTzyRnp5+0xUeAAAAt8PHx+fBBx+8qbBnz57WvyYmJiYmJlqX9O7dm9pg7969rf8qcB3hwoULmzRpEhQU9Pbbb6slp0+fjoiI+PeGdLo2bdqcOXPG9g0CAEAdo9hF2zbb2hFOmTKlsLDw2rVr69ate+2117777jvGWGlpqfUAsV69eiUlJdQW1NgrAAA4qaoq/ukJa3acI9S6H7S5I2zTpo3a53Xv3n3cuHG//PILY6xp06aFhYWWOoWFheqFGvx7wrymAADOzJYIiB39oOYnb+3pnAoLC318fBhjMTEx+/fvr6ysZIxdv3796NGjMTExkhsIAADgSLaGZebOnXvXXXc1aNBgx44da9as+e233xhjXbt2jYiIeOqppyZPnvzee+/16NEjMjLSka0FPm46lBHhxipijVxqTjJDBb+8ghcErSjjp0Op8vIyXmpUJGLKGDPwUqZUapR6mNTSwdwwLbH8sFhslAiHVjP3GBGzFEqNigRBPbz4nw9elfyNcJ/bKiN/I57Exs2e/CdRz/vW7kZ9ldfxN6L5ArAuQfsBnjBbR4Rms/nNN9988sknDx06tGPHDssUAOvWrSsrK3vkkUfc3d2//fZbh7UTAACcgROeI7R1RPjyyy9zy0NCQj799FN57QEAAKhRmGsUAACkUeiVVch/cZYRIQAAwN9zwnOE6AidCbmUIHGJJjcXQ6VFqCgKlXMpu8EpL79uuLWQMVbOq8wYK+dtvIKXoGF0iIa7TiG1AiI1Y5yU9QiFwjKy1iPkhmX07mJLCXJzMZ7exHNYyf/c4IdliCfWTCSOzGYJH0p0iIZTpvlwpI6xZz1CrV8CXNsHAAAuDSNCAACQxp4p07QeEmJECAAALg0dIQAAuDQcGgUAAGnsCcs4piW2Q0dYS3FjiVRYkVxTl7t2LpEOpYKdZUQQtOwap5ysfF0gNVpJ5FSpBXu5qVHuar2MXoC3ykg8t9yFeR2aGiUX5hVYsJdagJdKjfKnqaPSobwnnBHPLbVz0qlRbrEYaio14qklSjX/eHZOChM+R+g0yzABAADUSRgRAgCAPHZcUK/14BsdIQAAyFMLJtEWhY4QAACkseM6QpwjBAAA0BJGhBoj44e8P5DpUCIJyV1Tl5w7lAh83iglyrmpUV4hY6zsBr+c2xhyrlEqNcorp6ZUpZ4rk4zUKPVqcr/viqZGdURqVM9Ljerd+c8VtTAvf75WMmErsogxMRGuaDpU6DmkyvUK57nS8Z8ShjSp60BHCAAA0th1aNRBbbEVDo0CAIBLw4gQAACksWNmGVw+AQAAdYvWHZsodISaE4hdCE2lxhir4K1ky11Ql9E5F24ohjF2o7SSsxFqijXiTit45RXEJHBSwjIuMsUad7Vexpg7lSHirqkruIgxd9Y08rkiUIMJhfe0cAsZvbgx97mlNqIoIi8n/IeiCF8OgcsnAAAAtIQRIQAAyGPH6hNaj7HREQIAgDx2zDWqNRwaBQAAl4YRIQAASKMwhVoPsrp/0hQ6whpCT6XGL+cGRKkpr6g4JXcCs3Jy7VyBBXip+tQCvNTG+VOsEalRbjqUKje4dmrUnUiN0g+fu6YuMZUasaYu+fB56IfPr88NgroR6VDuc0JthE6NCnw6a36Wq/awZ4V6rZ89dIQAACCVkw0IcY4QAABcG0aEAAAgjT2TbjuoKTZDRwgAANI441yjODQKAAAuDSPCGiMwPSMjAnvGSmISTt6cooxKjVJzjVKBT2JNXe70oVQ6lLpTodQoFY41IjV6iyoqNUpME8otp3ZOsdlDReYOZfQ0oUJTqpLlvOeQip7qiACrwl/IV+tBTe1hxwX1Wj956AgBAEAee84R4jpCAACoK5xwQIhzhAAA4NowIgQAAGnsuHxC86ll0BHKx40SkFOpEXkEbqbDUElMsUakS8rLOOXlvHxKNeXcPAuj1tQV3Aj3TqmHQ5UbeBkiF1+Yt8ooFpYRe5gE7gMiH47g7Gh6D15YxoMfXHH34O8q7rz6enf+LkHFdqiYD7+y5kf9NOFsx0ZxaBQAAFwaRoQAACCNPZNuO6YltkNHCAAA0iiK2MIdrBYcQEZHCAAA8jjh9RM4RwgAAC4NI8IaIrQAL2PMaOCUc+ORjJ57jDv1WgUvSsoYqyQDn8Sqv7wMJzU7GjlrmshGDMRMctwwLZUapRY3NslIjVKvMvfIj2hqVEekRvV6zr1S+xU5axq3lEC1nBun1Lnxv22T6VA9sR/yAp8envzKlV78NKkH7x3EjZIyxtx4TyxjTMerrvnBvdrDCQeE6AgBAEAiJ7yOEIdGAQDApWFECAAAUml+rFOQ2IjQYDAkJSU9+OCDlpLMzMyEhITg4OBBgwZlZ2fLbh4AADgT9TpC0R9tiY0IX3/99XPnzlmW6TKbzSNHjhw7dux33323cOHC++67LyMjwwGNrK2ogAEvMkGvO0iEZWRMscYtpxYvJBM3IuUGojJVzm0hFYqhWiIUlnHiKdaq+OXcXYiauk9o1jTq04maY4w7axoZiiGWDHTnTaXGGHPn7RWV5fycCzkbnzdvVyF2Cfcqfkv472Xyo5xaj5Gq7/TsmGtU+JyibAIjwsOHD69fv/7pp5+2lKSnp+fn58+aNcvX1/d//ud/jh8/vm/fPgc0EgAAwFFs7QiNRuOkSZMWLVrk7u5uKczMzIyJiXFzc2OMeXp6RkZGHjt2zCHNBAAAp6CI/2jN1o5w4cKFPXr06Nmzp3Vhfn5+gwYNLL82btz46tWr1BZMJv7xBwAAcApVVfxjztbsOUGodV9oU0d45syZf/7znxMnTszOzr58+XJFRUV2drbZbPb19b127ZqlWnFxsb+/P3lPOlyqAQDgxNTjf3WPTWGZq1evBgUFpaamMsYKCgouXbqUnJy8a9euNm3aZGZmms1mRVGMRmNWVlbr1q0d3GAAAKi91GGe4L9ozKaOsEuXLnv37lVvf/XVV2+99Zb6a0JCgqIon332WWpq6uLFi5s0adKjRw8HNraWIUOj3EJqAV4iNVrFCz1SU6yR5dxgJxE9JctFgqCiqVFuy4VbWPML8xKvphAqe0mlRk0mqlzKmrqcjSs6/qtGralrEEmNGojUqMGTP+bgvvpC+xVVzn2vMfq9qec+4cRzQtH8o9+BHDDHmslkysjIKCgo6NWrV8OGDbl1Tp48efz48aioqJYtW1qXnz9/XqfThYaGWkrOnj1rOcbr4+MTGBgofLiyfv36li3q9fpvv/32rbfe8vHx+eSTT77++mvNU7AAAFCXGI3GIUOGPPLIIx9++GG7du0yMzNvrfPBBx/06tXrs88+69at27Jly9TCr7/+2tfXt02bNvfee6915bi4uOHDhycnJycnJ7/33nvMjpllhg0bNmzYMMuvPXr0sBwdFd0UAADUOeLXEVY7JFy/fv3p06cPHDjg5eX1/PPPz5s379tvv7WuUFxc/MILL+zatSs2NvaPP/4YMmRISkpKvXr1evbsuX///l27di1atOimbX7//fdt27a1/ConwIJeEAAAmANSo2lpaSNHjvTy8mKMPfjgg+vWrbvpGoRNmza1aNEiNjaWMdajRw9fX9/t27czxlq0aBEeHs7dZnZ29qFDh8rKytRfkeQEAAB57LiOsNqO8MKFC82bN1dvt2jRorKy8vLly1QFxljz5s0vXLhQzQY9PT3nzJmTkpISGhq6Zs0ahkm3AQBAcydOnJg/f751SY8ePeLj4xljFRUVlolc1BsVFRXWNSsqKvT6/+/LPDw8bqpw6335+PgwxtasWTNx4sSEhAR0hPJx83pUzFBorlEqCWkQKTcaxIJ25Aq33FyrYAu5GxG6R1kboV4IodSoyFSj0lKjQvlQoQlOdVQ6VM9/Dt14z61e8IWg90NewFjK/kat1UztErxinciazHWbPXONMsVgMBQWFloXFhQUqDeCg4MtU7VcuXJFUZSgoCDrmtYV1DrBwcHV3J3aCzLGxowZ89hjjx0+fBgdIQAASGPf1RPR0dFvvvkm9689e/bcvHnzc889xxjbvn17XFycp6endYXu3btPmTKlpKSkYcOGV69ePXbsWLdu3Wy530uXLhUVFQUGBqIjBACA2mvChAkLFix47rnn2rVrN3v2bEsEtHv37g888MD06dPbtWt3zz33jB49esKECcuWLRszZox6yjArK2v58uVHjx69cOHCrFmz2rdvP2HChPT09FWrVnXp0qW8vHzJkiVJSUmRkZEIywAAgDyywzL+/v67d+82m827d+/+7LPPxowZo5Y/+uijlumvv/rqq4EDB27btm3EiBHLly9XC/V6va+vb+/evadPn+7r61u/fn3GWLt27dq0aZORkZGVlfXCCy8gLAMAAJIpdlxQ93fVW7ZsuWDBgpsKx48fb7nt5eX1zDPP3Ppf6gFVa4GBgc8+++xNhegIbwe1OqvIwrzEohzcmcDIuaBkpEiouceoKAG/hSKVqXKheyQ3IjrFGpGMMHOXvRWcwIxLR3xYmIiIiptZ4MOF+iCikjhG7uxoRn5loRdCyn7FRHdmkXLqHqn3Jv+9LLxL1NkUjR0rzmv+XODQKAAAuDSMCAEAQB4HTLrtaOgIAQBAHjuuI9T6ckscGgUAAJeGESEAAEhjR1gGh0admUBolFwo1cSfT4o/vRO5ii85RxQvI0duhGghUc6tT21EqFzoHsmNUA+TmjdLaGFe4tWkkoPczwUqBEpN1kVRFIEXQkds3c3NUa+mlP2KES8cdw9ngtlgcio14r3JffXJ0KjQtHugEXSEAAAgjSJ7PcIagHOEAADg0jAiBAAAaey5oF7rA8XoCAEAQBp7lmHSuifEoVEAAHBpGBHWECpURqZJeWFFocVjqfrCGxEpr/0b4c4dKrpxSQvzErWpjVABY14QtPa/EMLlUnZm7kaoV1MwGwz/ZsfMMlpDRwgAAPLgOkIAAHBlClNEL4fQuh/EOUIAAHBtGBECAIA8WH2ibqJyLkL1yRPvAifkycr8+aQYsYCoaEsEyrXZiEiehVpTV2zjojEK7vuceNWoqdeIB8R0QnGe2v9qipRTzwn1jhCaHY1qCfHmJOryixlvXjz1D0S581DErwvU/EHj0CgAALg0jAgBAEAaRRG/QF7rC+rREQIAgDw4RwgAAK7MnrlGHdMS2+EcIQAAuDSMCOXjh0apyiJ/IFNscppie0OkbZxbLrzG6e3Wre5ORSLAghsnvgYLb1xGbf6dyng15exXohuntiGW6L79hrgiTLoNAADgZNARAgCAS8OhUQAAkAmTbgMAgOuy5xyhg5piM3SE8vEXn6Mqi/yBqkzuRmJNsb0h0jbOLRe7Rwl1q7tTfrFgzoW7ceFdgqovozb/TmW8mnL2K9GNU9vgb4R4kUU2rvlHeS3ihNcR4hwhAAC4NIwIAQBAGme8fAIdIQAASGPHzDKaw6FRAABwaegIAQDApeHQqA2o8JjQ0ppkMk0gsUZWJr7P6IQ2IqNcm43wHie3kDGmIzZCLYfLf26JdV+FQoxUC6mHSbVc6OE7waspUk48SvIdQbRQoDL5D1IS3XWCM54jxIgQAABcGkaEAAAgk7Oty4uOEAAAJHLCC+pt7Qg3bdq0bt26vLy8gICAhx56qHfv3mp5cXHxm2++efz48Y4dO86cOdPHx8dhTQUAgNrOjnOEmrP1HGF2dnanTp0mTpwYERExaNCgXbt2qeVjxozJysr6xz/+8eeff44fP95h7QQAAHAIW0eEU6ZMUW8kJSVlZGRs3rz5rrvuOnz48O+//3758uV69ep169YtODg4Ozu7VatWDmutExNNpul4eTidm0Blqr7wRkTKa/9GTNTDF1r1V8pcoyIvPWNMEXnhHPtC1HhLqDuV844QTthyi+Hf7LigXvOnVCw1WlVVdejQoT179sTHxzPGMjIy4uLi6tWrxxjz9fWNiorau3evQ5oJAABOQbHrR1MCHeHSpUsbN24cExMzfPjwfv36Mcby8vL8/PwsFQICAvLy8qh/N5moa68AAMAJ1NWPcYGOcPLkyaWlpdnZ2ak+QKUAAB8oSURBVDt27HjnnXcYYz4+PhUVFZYKZWVl1YRldDpcswgA4MRs+RhXGFP+fXxU4KcGGl8N4c6pZcuWY8eO3bJlC2OsWbNmZ86csfzp7NmzYWFhEhsHAADORWHO1g3aHpY5f/682slVVFRs2rSpY8eOjLGBAwempqbu3r27e/fuW7ZsKS8v79OnjwMbW9sITcBEJSPc+BvhntV3I6IBbnr+FxqdG6ec3gjRQqKcW5/aiFA5eY9VRP7FxCnnFjLG3KigC0Gp4t2jSLKGQk6ZRr1AVDCEu6sIvhDc55zciMguJGW/ou6Uu4cz+h1BbETsvSk0T5v2n/E1z77rCGW8p+xma0fYt29fLy8vf3//v/76q1OnTi+++CJjrH79+u++++7QoUNjYmIOHjz44Ycfenp6OrK1AAAAktnaEZ44ceLEiRNFRUVhYWHNmjWzlKempg4dOvTUqVMRERH+/v6OaSQAADgJOy6oVxTnGBG6ublFRkZy/9S0adOmTZvKaxIAADgre64jdExLbIckJwAAuDRMug0AANLYtR6hg9piK3SEt0NgLU5ynVjiuhxu7M3NnahMlOt55dxC6h4ZY3oqgMdtoUhlqlyv558uMFUR5bwQJzXbGUUh1lk28eZe0xGxUbGlmmXMMcaIF1T0heC+ylJeTSn7FRPdmUXKydA18d7kv3CIjTozdIQAACCNHecINf+2gHOEAADg0jAiBAAAaRQmfo5Q6yEhOkIAAJDHCVeox6FRAABwaRgRysc9KkAkAekoIDdOSWTh3EXK9e78KRTdPfjlenfebJtEY6iWVFHlRk65UDqUMWYWSYhSx2yqjCKRVBmTjYqmRunpQwV2FaFySfubWEvo/ZBTLtQS6k6p/Co5B6lQaNT1OOPCvOgIAQBAHjumWNP62Cg6QgAAkAfnCAEAAKQrLS3Nzc2tpoLBYDh//rzRaLRla0aj8fz585WVleqv6AgBAEAaxa6f6s2cOTMsLKxbt27du3e/fPnyrRV+/vnn0NDQhISEsLCwrVu3qoWbN2/u06dPw4YNu3btal35t99+a9GiRUJCQkhIyPfff89waPR2CI3mqWQEucYpN4pC5AjIck9evsCDCikQ5byNMMbcDaZbC6uMnELGWBWRf+GWU1EUwVCM4ARmvKnUGDl/G9ESkTnWRFsoJSwj9OoL7VdUuVBl4RaSGxEoJ+djI96b3PcyJlizUOh9m/yXaqvv3Lnzyy+/PH78eGBg4NixY+fNm7d48WLrCpWVlRMnTly+fHlSUtLq1asnTpx4+vRpNzc3Pz+/p59+Oisra/Xq1ZbKZrM5NTX1lVdeSU1N3bZt26hRowYOHIgRIQAA1F6rVq1KTk4ODAxkjE2bNm3VqlU3fRPdtGlTvXr1kpKSGGNjxowxGAzp6emMsbi4uOHDh6v/aLFnz54rV66MGzeOMZaQkBAaGvrTTz+hIwQAAHlkHxs9ffp0RESEejsiIqKkpKSgoMC6wpkzZywVFEVp06bN6dOnqa2dOXOmZcuW7u7ulg2ePn0ah0YBAEAaeybdZuzy5ctbtmyxLmnTpk14eDhjrLS01NvbWy308fFhjBUXF/v7+1tqlpSUWCqodYqLi6k7urVySUkJOkIAAJDGrvUIlaysrPnz51sXDh06dPr06Yyxpk2bFhUVqYWFhYWMsZuOdlpXUOvcVKH6yl26dEFHCAAAGuvVq1daWhr3T506ddqzZ496e8+ePW3atFHHhdYVDhw4UFlZ6eHhcePGjSNHjsTExFB31LFjx5MnTxYVFTVu3NhkMu3du3fGjBnoCG8D+aWHFyojzsYKTbFGBeo8vPgZOW65hxf/RTdU8KdSM1Tyg6BGGalRfiaTW5XG/fpJZS+NVEzXTUJqlAqTcr8fC+daidSo2Gx8VCaTt1dQ+5VQENSDqCxcLtJCsuW8hy86xZrYwrwuGhsV/xfaww8/HBMTs2rVqnbt2s2ZM2fq1Klq+bhx4wYNGpSSktK1a9d27do9/fTTU6ZMef/997t27RodHc0Yu3TpUnp6+p49ewoKCtasWRMcHNy7d+9WrVolJiY+/vjjM2fOXLFiRZMmTfr06YOwDAAASKSeJRT7qWZzLVu2/OGHH7788ssnnngiJSXliSeeUMvDw8P9/PzU22lpadeuXZswYYLZbP7222/VwkuXLq1Zs+bSpUt33HHHmjVrfvvtN7X8iy++8PHxSU1NzcvL++mnnxRFwYgQAABqtbvvvvvuu+++qfCVV16x3A4JCfn8889vqtCpUyfrKwgt/P39lyxZYl2CjhAAAKRRasFqEqLQEQIAgDxOOOk2OsIaQn1FIsMy7pxyauIooXyBpzcViuHvDEYq/8IrN1XxK8uZNY0q13EeEfXEuhn55WLrETrxFGtCuSp+ZU+RcrKyN9ESb/5+yK1PJb+odwT34XPfa6yasIzWn9q1nH2XTzioMTZCWAYAAFwaRoQAACCNPSvUO6YltkNHCAAA8jjhOUIcGgUAAJeGESEAAEhjR1hG8xEhOkL5hPYBOiIoYYo1r0pOuZFKhxJTqQnNmsbNWDJ67jEu4TglL91n0As+nLqVGnWnUqMis6PRgU8q2MlLKdcTqMwY8xIp96KipyJTrHGfQEa/EMSEedy6rsgJj4zi0CgAALg2jAgBAEAqzYd4gtARAgCAPE54QT06QgAAkMae6wi1HkHiHCEAALg0jAhrDLXIJ7+2zo03gSQ5USQ/rchdO5dbyOg4pYlYU5c/fajgorr8AJ5IOpQqdxN8mEiN3loomhr19nHnFBKpUe96nMqMMW8fIjXK2w7VEmoOUu47iPteY/R+6HxnwGqYHbFRraEjBAAAaRRF/Jyf1sdGcWgUAABcGkaEAAAgjR0L82o9IERHCAAA8ihMUQRPEmrdD6IjrCmiX3m4iQk9MRcUtQxplZHz+nKnRmN0KMbED5eIralL7ench0mFF6ipsPgr07p2WIZcmJdaw1koLEPmXzjl3AQNY8y7PhWW4Zd78crJsAzxMLnvINEFeDUfvtR2TjjHGs4RAgCAS8OIEAAA5BG/oF5z6AgBAEAaO5ZhUhRF8ApkyWztCAsKCjZv3nzu3LnmzZsPHz7cy8tLLTebzT/++GNmZmZMTMygQYMc1k4AAACHsPUcYVxc3Ndff52fn7948eLY2NiioiK1/PHHH3/hhRdu3Ljx5JNPzp4922HtBAAAZ6DY9aMpW0eEGRkZAQEBjLGqqqoOHTp8//33qampOTk5y5cvP3nyZGho6Lhx4zp16vTMM8/4+/s7ssF1DzX1GudQARVvoyKCnryJpvhTo1VTLmNNXaFZ06h0KPUwDbxyaiY5IzWTnIzUKPVcEUu5iqVGdURqlJuElJMaJZa9JadY482aRqVD61GpUaqcO8UaNZUa8fC57yD6MJ7WH8/OyZ5Jt4UnZ5TM1hGh2gsyxtzc3Dw8PNzc3Bhj27dv79ixY2hoKGOsTZs24eHhu3btclBDAQCg9lP+3RUK/Gj+nUP48onVq1dfvXp1+PDhjLGLFy8GBwdb/hQUFHTx4kXqH8UuOwMAgFqmrn6Mi6VGd+7cOW3atB9++KFx48ZMjfpYPS9ms1nz9RUBAEBjznZBvUBH+PvvvycnJ3/zzTc9evRQS4KDg3Nzcy0V8vLyQkJCqH9HHwkA4NRs+Rivywvz7tu3b+TIkZ9++mlCQoKlMCEh4ciRIxcuXGCMnThx4ty5c/Hx8Q5pJgAAgGPYOiIcPHiwj4/PihUrVqxYwRi79957H3jggZCQkMmTJw8YMODee+9dvXr1M8884+vr68jW1kFCiTUdP8TH9Gb+FxqzJzfxyH/RRQ/+c1uuI9cZFpknU2/kVnb3IFKjvMSjobKKW9lF5hp1J9Zwpp5DbmrUg0iHepEL8zpwrlFuVNXdkwgYE8FjIjUqNtcoVM+eC+od1BSb2doRLl68uKrq/z9ZIiMj1Rvvvffer7/+evTo0aVLl1oPFgEAwBU54aTbtnaEI0eOpP40cODAgQMHSmoPAAA4MzvmGtW6I8TqEwAA4NIw6TYAAEijiF8jILqQr3ToCGsp6vw9t5QM0cgY8ZPLk/LX1BUIdDDG9LzsRgUV9Cjnh2gMlZz8i9HAD8tQU68JrVfs2LAM9RyKTLCndxcMy/DmKqOmWKPCMl68WdC4C+oyYso0Vs1au7xVgoWmUmPEc45QDODQKAAAuDSMCAEAQBr71iN0UGNshI4QAACksWNmGa1PEaIjBAAAibRfTEIYzhECAIBLw4jQmZAHHIjvM25yUqMCa+0KJRsZY3peQNTDk58OrSznhxgNFZyAqIFamJcoN1VR5ZwkKLGGsVhqlJhJjUw86twE1it2F12Yl5fJ5EZJGWNeIgv2UilQak1dctY03iNyI54ThdjxtT4V5RLsOkfooLbYCh0hAADIU4dXnwAAAKiT0BECAIBLw6FRAACQRmF2LMOE6wjhtpG7kY4T3qASNNRGhGYCI6dSo7IbHpycSyUvuMEYqyznz5rGXXqQCstUUeXEOoX89Qj5dYWyMmSgg16PkP8PbkJhGWL6On5YRiRZw4hwDbURMj8ltJQgFYrR+oPVlTnjdYQ4NAoAAC4NI0IAAJCnDi/MCwAA8Lecca5RHBoFAACXhhEhAABIY0dYRusjo+gI6wZq7Vz+1F78ytSMX3oimccNN1LTg+n1/HJuiNGDlwJljBm8+WFNbmqUXIDXkQvzUuv1cj8XZC3My02N0jFdgQV7hSozYsI8oRRoNeVYU9eZONs5QhwaBQAAl4YRIQAAyOR013GiIwQAAGnsOUeodb+JjhAAAKSx5/IJrUeQOEcIAAAuDSNCl0N/V6NSfERtXmqUW8iqW7CXk+GkwopGYjpQbhCUrEzNKUqlRrlzjZIL8wrERqkXgpprlEyT8mKZVFaTGzGl6uvdxWY95S4dTD0ceu1cKjXKrw+1DmaWAQAAV4ZJtwEAACT78ssvIyIiAgMDJ0+eXF5efmuFI0eOxMfH+/n5JSQknDhxwlL++uuvN2/ePCwsbN68eZbLf/v169flP55//nmGESEAAEgkfT3CI0eOTJs2bePGje3atRs5cuQbb7zx0ksvWVcwm82jRo1KTU3duHHjwoULk5OTDxw4wBhLS0tbsmTJ1q1b3d3d+/fv37Zt25SUFMbYwYMHP/nkk7CwMMaYr68vw4gQAABqs08//XTUqFE9e/b09/efM2fOJ598clOF9PT0wsLCGTNm+Pj4zJ49+8yZMxkZGYyx5cuXT506tW3btuHh4dOnT1++fLnlX6Kjo+Pi4uLi4lq1asUwIgQL0RCNonCCIeS0YTqBebbc9PzIiXsVfyPc2dHo8AsVluEW82dTc+gUa/TUa/yNc59b0XnauPXJWdDI/IvQLGgIxdRN0q8j/OuvvwYOHKjejo2NvXjxYklJScOGDS0VMjMzO3bs6Obmxhhzd3ePior666+/7rzzzszMzGnTpln+8fXXX7f8y6hRoxRF6dat27x580JCQtARAgCAxiorKwsLC61LfHx8PDw8GGP5+fmWbk+9cfXqVeuOsKCgoEGDBpZfGzVqdPXqVbXc+h/VQsbY+++/37lz5/Ly8jfeeGPQoEF//vknOkIAAJBH/IJ6prDNmze3bt3aumzSpEnz589njPn5+ZWWlqqFJSUljDF/f3/rmn5+fteuXbP8WlxcrFa46R8t/6WeKWSMrVy50s/P7+DBg+gIAQBAHjuuI2RsyJAhaWlp3D9FREQcPXpUvX306NGmTZs2atTIukLbtm2PHTtmMpl0Op3RaDx+/Hjbtm3V8qNHj6qHVY8ePaoWWnN3d9fr9QaDAWEZAACovSZMmLBmzZojR47cuHHjzTffnDhxolr+8ssvb9q0iTHWt29fLy+vpUuXVlVVLVq0KDAwsHv37uo/Ll68ODc398qVK4sWLZowYQJj7NSpU//617/Ky8uLi4ufffbZRo0axcTEoCMEAABp1LCM6E81Onfu/Nprrw0YMCAkJMTPz+/FF19Uy48dO5aXl8cYc3NzW7169ccff9yoUaOvvvrqm2++UY/N3nfffWPGjOnYsWP79u0HDhyodoTFxcWTJ0/29fUNDw8/duzYTz/9VK9ePYXKv8lVVlY29eUFQyf8owbuCzRE700CMUtqI2bebGdUOVGX3ohQEJSqzL9PsXwk9aFApky5KyRTG5ER+BQKgiIFWpeMbBf8t3W+3JF9pZhzzXs1jvzv1sLMbdSh0RqAc4QAACCVs337waFRAABwaRgRAgCANHaFRjWGjhAAAKSxZ2Ferc8k49AoAAC4NIwIQSbRiSXFNkKUm7nzZAqunSuWnqYiqUR1Iqkpco/VPS0igU/hrd9mXXA9WJgXAABcmiJ8qFPzQ6PoCAEAQBp7Vp9wTEtsJ3aOMC8vT72S31pubu6OHTsuXbokr1UAAAA1xNaO8OOPPw4ODg4JCXnkkUduKu/YseObb74ZHR391VdfOaCFAADgVBTxH03Zemi0d+/e27Zt27Bhw86dOy2FpaWlM2bM2LJlS9euXbdt23bfffeNGjXK09PTMU2FOkj0VILDtkwHXW63Lkn4vS8liCRhGwDVUWcP1boVYmwdEUZFRUVGRt60Fvavv/4aFhbWtWtXxlhCQkK9evV27NghvYkAAACOc1thmfPnz4eHh1t+bdGixblz5263RQAA4LTsCctoPYC8rY6wvLzcw8PD8qunp2dZWRlV2WQy3c59AQCAttTFb/+mUi045yfqtjrCoKCg/Px8y69Xr14NDiYX6fj7pw8AAGoxWz7GXW6Kta5du/7555/Xr19njBUUFGRmZt55552SGgYAAFATbB0RHj16dMOGDbt27Tp16tT8+fM7dep0zz33REdHx8fHp6SkpKamfvTRR8OGDbM+ZQggmei0YQ7biqJNbBTACdhxjlDz94KtI8LKysrCwsIOHToMHz68sLBQHQUyxlavXt25c+eVK1fGx8evWLHCYe0EAABn4GwXETLbR4SdO3fu3LnzreUNGjSYO3eu1CYBAADUHMw1CgAAEjlfWAYdIQAASOOEqzChIwSwg+ZvXIBaywl7QlzbBwAALg0jQgAAkMaOC+o1h44QAACkcca5RnFoFAAAXBpGhAAAII2iCF8OofmhVHSEAAAgldaHOkXh0CgAALg0jAgBAEAmZxsQoiMEAAB5nHE9QnSEAAAgD2aWAQAAcC4YEQIAgDTOeEE9OkIAAJDGGadYw6FRAABwaRgRAgCAPE4YlkFHCAAA0iji5/w0P5KKjhAAAKRRmKIIDvFE60uHc4QAAODSMCIEAAB57DhHqDV0hAAAII/4dYSad5w4NAoAAC4NI0IAAJAGk24DAIBrwzlCAABwZQqzY0TooLbYCucIAQDApWFECAAA0jjhDGvoCAEAQCIn7AlxaBQAAFwaRoQAACCNothxOQQunwAAgDoDK9QDAIBLwzlCAAAA54IRIQAASGPPBfVaDwnREQIAgDQKVp8AAABwLugIAQDApeHQKAAASGPXMkwOaout0BECAIA0dpwj1LofxKFRAABwbRgRAgCAPE54QT06QgAAkEf8HKHmJwlxaBQAAFwaRoQAACCNIj7A03pAiI4QAADksefyCa1PEtZQR6jX6//a+evO1V/UzN3VsPz8fF9fX52ujh9nxsOsS65everv7y++bpyTwcOUK3LDhsjIyOrr9AsPEN1sq7LOf+jK7W2UBIrZbK6ZeyooKCgqKqqZ+wIAAOmaNWvm4eGhdSvkq7mOEAAAoBaq48d/AAAAqoeOEAAAXBo6QgAAcGnoCAEAwKXhOkJhly5d2rt3b05OTr9+/Vq3bm0pP3v27Oeff379+vUxY8bceeedGrZQin379m3atOnKlSuRkZEpKSne3t5qeUlJybJly3Jycvr27Tts2DBtG3n7tmzZ8vvvvxcVFTVv3nzcuHH+/v5qeVFR0bJly3JzcxMTEwcPHqxtIyVavXq1p6fn8OHD1V8rKiqWL19+8uTJzp07p6SkOPtFI+vXr8/Ly1Nv+/n5jR49Wr1dUFDwySef5OXlDRo0aMCAAdo1UJrLly+vWLHi4sWLLVu2HD9+fKNGjZjVTtuvX78hQ4Zo3UZn4tz7vSbuuuuu119//bnnntu7d6+lMC8v78477ywqKmratGn//v3T09M1bOHtKyoqSkpKunLlSvPmzVeuXHnXXXdVVFQwxkwmU0JCwu+//966desnn3zy/fff17qlt+vbb781mUytWrX67bffYmNjCwoKGGNGozE+Pn7v3r2tWrV69NFHlyxZonUz5Vi/fv2kSZMWLFhgKbn//vvXrl3btm3bd955Z8aMGRq2TYoFCxZs2rQpOzs7Ozv7woULamFlZWXv3r0PHTrUsmXL1NTUzz//XNM2SpCVldWpU6ejR4+Gh4efOHFCfaRGo7FPnz579+5t3br11KlTP/roI62b6VTMIKiqqspsNsfExHzzzTeWwpdeemnEiBHq7QULFgwePFibxklSVVVVUVGh3i4rK2vUqFF6errZbP7555/Dw8MNBoPZbN66dWtISIh6uw4wmUwtW7ZMS0szm81paWkRERHqC60+ZPW2UysqKoqOjn7ppZd69uyplhw5csTHx6ekpMRsNp86dcrb2zs/P1/TNt6uXr16rV+//qbCr7/+ukOHDiaTyWw2p6WltW3bVr3tvAYNGjR79uybCtWHpu6oGzdubNGihdFo1KJ1TgkjQmHcw0fp6emWQy79+/ffuXNnzTZKMp1OZ7ls1mQyVVZWNmjQgDG2c+fOu+++W6/XM8b69OmTn59/4sQJLRsqz4kTJ4qKitq3b88Y27lzZ79+/dQXul+/fufOnTtz5ozG7bttTz311FNPPRUSEmIpSU9P79atm/rKtmrVKjQ0dM+ePdo1UI6NGze+/fbbP//8s/k/V0inp6cnJiaqs64MGDAgKyvr4sWLmrbxthgMhs2bNw8fPvzTTz9dsmSJZeB700574cKFOrDT1hh0hHLk5uY2adJEvd20adPr16+XlJRo2yRZnn322fj4+NjYWMZYXl6e5WG6ubn5+/vn5uZq2joJZs6cGRoa2qlTp4ULF6odofWr6eHh4evr6+wPc+vWradPn05NTbUutH41GWNNmzZ16h6CMRYVFeXp6Xn58uXHH388KSnJZDKx/34169WrV79+fad+Nc+fP28ymR577LEzZ84cPnw4Jibm2LFj7L9fTXd39zqw09YkhGXk0Ov1RqNRva3ecHd317RFcrz33nubN2+2nPLU6/VVVVWWvxoMhjow39LcuXOffvrpXbt2TZkypWPHjnfeeae7u3tdepjXr19//PHHv//++5vmoqx7r+bHH3+s3pg1a1ZERMSmTZsGDRpk/d5kjBmNRqd+mDqdzmw2P/bYY+rXGoPB8Pbbby9fvrzuvZo1CSNCOUJDQy3fpnNycvz8/CwxS+f1wQcf/POf/9y2bVtQUJBaEhoampOTo94uKysrLCy0PtTmpHx8fIKCgsaMGTNo0KB169ax/36YpaWlpaWlTv0wd+7cmZOTM3bs2C5durz66quHDh3q0qWLyWSyfpiMsZycHKd+mNZ8fX2joqJOnz7N/vu9efXq1fLycqd+mMHBwTqdLioqSv01Ojr67Nmz7JadtqSkxKkfZg1DRyhHUlLSd999px6KWbNmTVJSktYtul2ffPLJO++8s3nz5mbNmlkKk5KSNm/erE6enpaW1q5dO+sLSJyO0Wg0GAzq7crKykOHDjVv3pwxlpSU9Msvv5SWljLG1q5d27lz59DQUC0bent69eq1bdu2pUuXLl26dNy4ca1atVq6dKlOpxs0aND+/fvVj9F//etfFRUVPXv21Lqx9jMYDJaR3/nz5w8cOBAdHc0YS0pK+vnnn2/cuMEYW7t2bY8ePQIChJdHqD08PT3vueee3bt3q7/u3r1b7RSTkpJ+/fVXdaf97rvvYmNjrd+58De0Tus4n6lTp8bFxXl7e7dq1SouLm7v3r1ms/natWt33HFHfHx8cnJyYGDg8ePHtW7mbcnJyVEUpXnz5nH/8dNPP6l/SklJiYqKmjBhQkBAwI8//qhtO2/TuXPnAgMDR4wYkZKS0qJFi8TExLKyMvVPo0eP7tChw/jx4wMCAn799Vdt2ynRsmXLLKlRs9k8e/bs8PDw1NTUoKCgpUuXatiw23fq1Kng4OCRI0cmJyf7+vo+9thjarnJZEpKSoqNjX3ooYf8/f23b9+uaTMl2LdvX9OmTR966KHBgwe3adPm4sWLavmYMWMsO+0vv/yibSOdC1afEJaVlWUdhImIiFBzdxUVFdu2bbt27VpiYqKvr692DZSgsrLy8OHD1iXh4eHqxeZms3nXrl05OTk9e/Zs0aKFRg2U5ty5cwcOHCgvL2/btm3nzp0t5WazeefOnXl5eb169QoLC9OwhXJdvXo1Pz+/Xbt2lpK9e/dmZWXFxsb+7TpztV9mZmZmZqbJZOrUqVNERISl3GQy7dix4/Lly3fddZdTD+4t8vPzt23b1rhx4969e1vOwpjN5vT09Nzc3Dq209YAdIQAAODScI4QAABcGjpCAABwaegIAQDApaEjBAAAl4aOEAAAXBo6QgAAcGnoCAEAwKWhIwQAAJeGjhAAAFwaOkIAAHBp6AgBAMCl/R8jGHGujQeYugAAAABJRU5ErkJggg==",
"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.2"
},
"kernelspec": {
"name": "julia-1.7",
"display_name": "Julia 1.7.2",
"language": "julia"
}
},
"nbformat": 4
}