{
"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.107107e+01 7.291843e+00\n",
" * time: 0.004687786102294922\n",
" 1 2.876062e+01 5.700959e+00\n",
" * time: 0.013830900192260742\n",
" 2 2.304645e+01 8.538451e+00\n",
" * time: 0.035635948181152344\n",
" 3 1.369228e+01 3.329003e+00\n",
" * time: 0.06287288665771484\n",
" 4 1.344589e+01 4.949706e+00\n",
" * time: 0.07953095436096191\n",
" 5 1.267577e+01 2.016079e+00\n",
" * time: 0.09606385231018066\n",
" 6 1.030306e+01 1.828672e+00\n",
" * time: 0.11199092864990234\n",
" 7 9.804730e+00 1.001051e+00\n",
" * time: 0.12796688079833984\n",
" 8 9.348828e+00 1.100985e+00\n",
" * time: 0.14391088485717773\n",
" 9 9.085456e+00 9.168852e-01\n",
" * time: 0.16001176834106445\n",
" 10 8.894131e+00 5.946768e-01\n",
" * time: 0.17604899406433105\n",
" 11 8.780160e+00 5.250814e-01\n",
" * time: 0.18827080726623535\n",
" 12 8.722700e+00 5.630382e-01\n",
" * time: 0.20036077499389648\n",
" 13 8.643762e+00 3.687352e-01\n",
" * time: 0.2127208709716797\n",
" 14 8.596611e+00 2.993913e-01\n",
" * time: 0.2248549461364746\n",
" 15 8.579058e+00 1.886540e-01\n",
" * time: 0.23724794387817383\n",
" 16 8.569429e+00 1.507300e-01\n",
" * time: 0.24960684776306152\n",
" 17 8.558628e+00 1.775237e-01\n",
" * time: 0.26216983795166016\n",
" 18 8.552135e+00 1.312452e-01\n",
" * time: 0.27478480339050293\n",
" 19 8.548353e+00 1.458555e-01\n",
" * time: 0.3371129035949707\n",
" 20 8.537814e+00 1.403037e-01\n",
" * time: 0.34923481941223145\n",
" 21 8.536068e+00 1.913872e-01\n",
" * time: 0.36132287979125977\n",
" 22 8.533663e+00 1.213611e-01\n",
" * time: 0.3736109733581543\n",
" 23 8.531155e+00 1.087756e-01\n",
" * time: 0.3859109878540039\n",
" 24 8.528869e+00 8.732103e-02\n",
" * time: 0.39804577827453613\n",
" 25 8.525036e+00 7.347369e-02\n",
" * time: 0.41023993492126465\n",
" 26 8.523419e+00 6.792758e-02\n",
" * time: 0.42258691787719727\n",
" 27 8.522093e+00 8.064870e-02\n",
" * time: 0.43480682373046875\n",
" 28 8.519904e+00 9.278107e-02\n",
" * time: 0.4468247890472412\n",
" 29 8.518735e+00 6.592605e-02\n",
" * time: 0.45891880989074707\n",
" 30 8.517826e+00 5.937242e-02\n",
" * time: 0.4708578586578369\n",
" 31 8.516552e+00 6.000060e-02\n",
" * time: 0.48300695419311523\n",
" 32 8.515120e+00 7.825918e-02\n",
" * time: 0.4949328899383545\n",
" 33 8.514677e+00 5.167538e-02\n",
" * time: 0.5069518089294434\n",
" 34 8.514196e+00 5.071020e-02\n",
" * time: 0.5190749168395996\n",
" 35 8.513556e+00 3.271990e-02\n",
" * time: 0.5310459136962891\n",
" 36 8.512923e+00 2.751282e-02\n",
" * time: 0.5431389808654785\n",
" 37 8.512733e+00 3.217487e-02\n",
" * time: 0.5551638603210449\n",
" 38 8.512515e+00 2.426274e-02\n",
" * time: 0.567166805267334\n",
" 39 8.512408e+00 1.659509e-02\n",
" * time: 0.5791699886322021\n",
" 40 8.512355e+00 1.780821e-02\n",
" * time: 0.5912129878997803\n",
" 41 8.512252e+00 1.697188e-02\n",
" * time: 0.6032118797302246\n",
" 42 8.512147e+00 1.524332e-02\n",
" * time: 0.6154689788818359\n",
" 43 8.512066e+00 1.338087e-02\n",
" * time: 0.6278908252716064\n",
" 44 8.512043e+00 1.209438e-02\n",
" * time: 0.6402428150177002\n",
" 45 8.512022e+00 1.011911e-02\n",
" * time: 0.6899709701538086\n",
" 46 8.511975e+00 8.841898e-03\n",
" * time: 0.7020578384399414\n",
" 47 8.511932e+00 7.106516e-03\n",
" * time: 0.7141478061676025\n",
" 48 8.511911e+00 8.613331e-03\n",
" * time: 0.7266199588775635\n",
" 49 8.511889e+00 8.436513e-03\n",
" * time: 0.7390577793121338\n",
" 50 8.511876e+00 4.374544e-03\n",
" * time: 0.7510759830474854\n",
" 51 8.511871e+00 4.957880e-03\n",
" * time: 0.7632238864898682\n",
" 52 8.511863e+00 4.498580e-03\n",
" * time: 0.7753629684448242\n",
" 53 8.511857e+00 3.170509e-03\n",
" * time: 0.7872359752655029\n",
" 54 8.511854e+00 2.437927e-03\n",
" * time: 0.7992448806762695\n",
" 55 8.511853e+00 2.246858e-03\n",
" * time: 0.8112249374389648\n",
" 56 8.511853e+00 3.645975e-03\n",
" * time: 0.823185920715332\n",
" 57 8.511851e+00 2.312424e-03\n",
" * time: 0.8352530002593994\n",
" 58 8.511849e+00 1.738604e-03\n",
" * time: 0.8471779823303223\n",
" 59 8.511848e+00 1.317958e-03\n",
" * time: 0.8592209815979004\n",
" 60 8.511848e+00 8.898985e-04\n",
" * time: 0.8712368011474609\n",
" 61 8.511847e+00 1.677463e-03\n",
" * time: 0.883171796798706\n",
" 62 8.511847e+00 1.617500e-03\n",
" * time: 0.895258903503418\n",
" 63 8.511846e+00 1.070956e-03\n",
" * time: 0.9072108268737793\n",
" 64 8.511846e+00 9.476841e-04\n",
" * time: 0.9191999435424805\n",
" 65 8.511846e+00 7.572531e-04\n",
" * time: 0.9311749935150146\n",
" 66 8.511846e+00 6.234727e-04\n",
" * time: 0.9431848526000977\n",
" 67 8.511846e+00 6.651777e-04\n",
" * time: 0.9551668167114258\n",
" 68 8.511846e+00 5.189622e-04\n",
" * time: 0.9671077728271484\n",
" 69 8.511845e+00 4.145397e-04\n",
" * time: 0.9792499542236328\n",
" 70 8.511845e+00 3.287611e-04\n",
" * time: 1.0219318866729736\n",
" 71 8.511845e+00 3.203461e-04\n",
" * time: 1.0339758396148682\n",
" 72 8.511845e+00 2.453704e-04\n",
" * time: 1.0459418296813965\n",
" 73 8.511845e+00 2.266407e-04\n",
" * time: 1.0580909252166748\n",
" 74 8.511845e+00 2.465120e-04\n",
" * time: 1.070326805114746\n",
" 75 8.511845e+00 1.801527e-04\n",
" * time: 1.0823619365692139\n",
" 76 8.511845e+00 1.499721e-04\n",
" * time: 1.0944569110870361\n",
" 77 8.511845e+00 1.511243e-04\n",
" * time: 1.1066389083862305\n",
" 78 8.511845e+00 1.294225e-04\n",
" * time: 1.1186778545379639\n",
" 79 8.511845e+00 1.227820e-04\n",
" * time: 1.1305859088897705\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVyTR/4H8HnCJeIBgsqNiqCIonjjXbGiVjzqWY96tLVuT9tatbo93Lq2tbZ1t4dnbd219arV2qpV8QC71lW8xQMVrAriwSGKHAnJ749nN7+snS/NxAkPIZ/3iz/CMDyZJE8ymef5PDOKyWRiAAAAzkqndQMAAAC0hI4QAACcGjpCAABwaugIAQDAqaEjBAAAp4aOEAAAnBo6QgAAcGroCAEAwKmhIwQAAKfmqnUDAAAAKpKenv75558XFBQMHTp0yJAhv69QVFT0t7/97cyZMzExMS+99FKNGjUYYydOnNi8efOlS5caNmz41FNPNW/eXK38l7/85f79++rtmJiYMWPGVFJHaDQaj548HREVVTl3V8mMRqNOV/3H1niY1QkeZnVSaQ+zrofbH9Y5fz6vuNggtNm7d+8EBurCw8O5f719+3bXrl2nTp3asWPH5557rqSkZPTo0Q/UeeKJJ8rLyydNmrRs2bJjx46tWbOGMfb666/HxsYmJCQcP368ffv2//73v6Ojoxljf//73ydNmuTn58cY8/LyYowplTPXaHFx8fN/WTBw4pRKuC8AALCHx5sF/GGdNm1WnThxS3DDp4cMubNp0ybu3z788MO9e/du27aNMfbPf/5z0aJFR44csaxw7ty5tm3b3rhxo3bt2nl5eYGBgefPnw8LC7P8fpCYmNihQ4e33nqLMebn5/frr79GRESYt1D9vysBAIDj+vXXX3v16qXe7tWr17Fjx0pKSiwrHDx4MDY2tnbt2oyxevXqRUdHHzp0iDFmOUq+ceOGOgRUffbZZ7Nmzdq8ebM6FERHCAAA0iiKLT8VyMnJ8fX1VW/7+fmZTKacnByqglonOzvbssLSpUtv3bo1YcIE9ddBgwaFhYXVqlXr1VdfnThxIkNYBgAANLdv37727dtblowYMWLmzJmMsRo1apSVlamFpaWljLGaNWta1vT09NTr9eZfS0pKLCv88MMPb7/9dlJSkno6kDG2cuVK9caECRMaN248Z84cdIQAACCNoihKxUO83zGZlDZt2ixcuNCyMCgoSL0RHBx85coV9faVK1dq1KhheZBTrWmuwBi7evVqcHCwenv79u1TpkzZtm1by5Ytf3+/ISEh9erVy87OxqFRAACQShH+8fb2bve//P391Y0NGzZs48aN6gUPq1evHjJkiHryLykpKT09nTHWt2/fq1evqgmaX3755c6dO+o5xV27dk2cOHHLli3t2rUzN62wsNA8vty+fXthYWF0dDRGhAAAUHUNHDhw5cqVbdu2DQ0NPXPmzO7du9XyOXPmPPHEE5GRkXXq1Pnggw/69+/fuXPnX3/99aOPPvL09GSMTZ06tbi4eMyYMWr9kSNHvvfee7/++uv48eNbtWpVUlKSlpa2ZMmS+vXroyMEAAB5/ij8IsrFxWXz5s3Hjx8vLCzs0KGD+fzf5s2bzaf9/vSnPw0cOPD8+fNLliwJDAxUC/fu3Wsw/P8VjWqsNCEh4ejRoxcvXvTw8GjRokXdunUZwjIAACCRDecImekP6iuKEhsb+0BhQMD/XNQYEhISEhJiWRIaGsrdWnBwsPkkogrnCAEAwKmhIwQAAKeGQ6MAACDNH14g/3t/dGTU7tARAgCANDacI5QbrrEBDo0CAIBTw4gQAADkUa+RdyjoCAEAQBpF/FCn5odG0RECAIA0tlxHqHVPiHOEAADg1DAiBAAAeWw4R4hDowAAUJ04Wj8o2BEeP378woUL9evX79y5c40aNdTCEydOpKent2jRIjo62g4tBAAAsCNrzxGWl5ePHz9+0KBB33333TvvvLN69Wq1/L333nvssce2bt3ap0+fTz/91G7tBAAAB6Ao/8nLCNB6TGjtiHDx4sWnTp1KS0tTV7IwmUyMsdzc3Hnz5qWmpkZFRaWmpvbp02fSpEm1atWyY3sBAKAqc8BzhNaOCFevXv3yyy/n5uYePHiwqKhITcfu2rUrIiIiKiqKMda+fXtfX9/k5GQ7NhYAAKo2GwaEwpdbyGbtiPDixYtr165dsmRJ7dq1z50799NPP7Vp0yYrK8ty/aeQkJBr165RW1AHkQAA4KBMJpPmnZY9WDsivH//fu3atQ8ePJiUlDRhwoTp06czxvR6vYuLi7mOq6urXq+3SzMBAEBr1oxn1NUnhH40Z21HGBgYGB8fr34XePTRR0+cOMEYCwgIuHXrlrnOzZs3H1gy2FK1/B4BAOA8dDorugzFph9NWdsR9urVKzMzU72dmZkZGBjIGOvWrdvRo0fz8/MZY9evX09PT4+Li7NTQwEAAOzB2nOEr776as+ePb29vevUqfPXv/71k08+YYyFh4cPHjx4yJAhY8eO/eqrr8aPH692kAAA4JzU452C/6Ixa0eELVq02L9///3797Ozszdt2jR69Gi1/J///OcTTzyRlpb2zDPPLFmyxG7tBAAARyB+jlDz82YCM8s0b9583rx5DxS6ublNnTpVapMAAMBhVYFzfqKw+gQAADg1TLoNAADSqBfUC/6Pg1xQDwAA8IdsOOendT+IQ6MAAODcMCIEAAB5HDAsg44QAACksWESbc3nHUNHCAAA0jjggBDnCAEAwLlhRAggTsqSYg73tRnAGg64MC86QgAAkMeGc4R2aonVcGgUAACcGkaEAAAgjSJ+gbzWoVF0hAAAII+iiF8OoXVPiI4QHAcRURFMrkjYipysjOhWyM8KgQ8RGdsAqJADhmVwjhAAAJwaRoQAACCNAw4I0RECAIA8tkyxpnVXiEOjAADg1DAiBAAAeRzw2Cg6QtCSSSg5SdTmlpJbpjYi1hKhYuJtLithzisXj6OLRE+1/tiCqsyG6wg17whxaBQAAJwaRoQAACCPA841io4QAADkwTlCAABwZorieHON4hwhAAA4NYwIQSY6eymQ1aQ2YjISG+GVE3XpjQilSUXyq4xMjQoVkxMZKzpOOa+MrMwYU3hfiaVETzX/sg+VTGGK6AXyml9Qj44QAADkseEcodZwaBQAAJwaRoQAACCNDWEZzUeQ6AgBAEAahYlfR6j1mWR0hPAHhPIvojkXIzfnUk5UJsrLeeXkRoxGYuPcYn6IRixZQwRGyPALUa5z4W9cp+Oc4NC58DfiQpRz61Mb0YkkbmhiWSFwGDhHCAAA4FgwIgQAAHkU4UOdmh8GQEcIAADSOOCRURwaBQAA54YRIQAAyINJt8FxkYFPKiHJS19yU6CMMWM5P6tZbuCUG/T8jRh4lRlj5XreRqjKRDmdMhUIxwrFRqnzIlQmkwyCunKO67jyChljLm78cm59VzeBe2SM6VwE8quKImVCOqhyFFuWYcLlEwAAUF3YsEK95t9y0BECAECVdvfu3XXr1hUWFg4YMKB58+a/r2AymbZu3ZqWlta6det+/fqZ/yspKenSpUv+/v5DhgypVauWuf7u3btTU1ObNWs2ePBgRRG8AhYAAKAiik0/tPv373fq1Gnr1q05OTlxcXHJycm/r/Pyyy/PmjWrqKho2rRpb7zxhlrYu3fvZcuW3bhxY/Xq1VFRUdevX1fL33333SlTphQVFc2dO/eZZ55hGBECAIBENpwjrPjY6Lffflu7du3vv/9eUZSAgIB58+b17NnTskJ2dvby5csvXLgQHBw8YcKEVq1aTZ8+3dfXd/v27X5+fmqd7t27r1q1atasWYWFhR9++OEvv/wSExPz4osvhoWFzZ49GyNCAACounbt2jVw4EC1c01MTNyzZ4/BYLCssHfv3ujo6ODgYMZYeHh448aNU1JSGGPmXpAx5unp6eLiwhg7ePCgj49PTEwMY6x+/fodO3ZMSkrCiNDpiKZDhWb+pLKahjL+VJ76Mk59vUhlqr6BFyVlRMSUEROWMuJh2nWuUdFpQrlBUFciHermzp+x1M2dU1+oMmPMlVefyq+SE5kSU6pi1V9HYcPqExVXz8rK6tOnj3o7ICDAaDTm5OSo3Z4qOzs7ICDA/GtAQEB2drblFvbu3ZuamvrVV19RldERAgCAxo4fP/7ss89alvTs2XPMmDGMMZ1OZ/7qqd544EukoiiW301NJpNlhePHjz/xxBOrV68OCgqiKqMjBAAAeWy4jlBh3t7e7dq1syxs2rSpeiMgICAnJ0e9ff36dRcXl4YNG1rWtKzAGMvJyQkMDFRvnz59un///p999tmAAQOoyvHx8RI6wrKysqysrODgYDc3t4ffGgAAOJtGjRpNmTKF+6eEhISlS5f++c9/1ul0P/74Y+/evV1dXRljly9frlu3ro+PT+/evZ9++umrV6+GhIRcvHjx8uXLPXr0YIylp6f369dv4cKFw4cPN28tLi4uPz//xIkTrVu3vnnz5qFDh77++mtrwzLTp09XLNy/f18t37JlS1BQUP/+/UNCQpKSkh7qmQAAAAenniMU+qn4JOHo0aOLiooGDx48ffr0d999980331TLR40atWrVKsZYQEDA1KlTExIS5syZ89hjj02bNq1evXqMMbX/++GHH0aOHDly5MglS5YwxmrXrj1z5syhQ4fOmTPn0UcfHT9+fKNGjQRGhHPmzJk3b55lSWlp6VNPPbVq1aoBAwasX79+8uTJmZmZajIHKhWVf+GuncsPi5ChGDL/wkudlJXycy5lJQainFOf3ggVouGU60XDMtTUa9wp1ojnUGjSMOoiXmqKNWpiM25Yxk0wLONeg1Pu7iFQmTHmXoPzBFAboeI8riYqXMMt5dalp+xCuKZySH2ea9asefDgwY0bN965c+fQoUMRERFq+YcffmiOzHzyySc7duw4ffr0559/bk7WvP/++0VFRebtNG7cWL0xZ86crl27HjlyZP78+eohU7FDo8XFxZ6enuZfd+7cWadOHXVDw4cPf+mll1JSUh555BEbHy4AAMDv1KpVa8KECQ8Uqsc/zRISEhISEixLzOcFf69Xr169evUy/ypwHeHChQvr16/v7+//0UcfqSWZmZmRkZH/2ZBO17Rp08uXL1u/QQAAqGYUm2jbZms7wqlTp+bn59+7d2/z5s1//etfN27cyBi7e/eu5QCxZs2ahYWF1BaMRuJwEgAAOILycv7pCUs2nCPUuh+0uiNs2rSp2ud17tx5/PjxP//8M2OsQYMG+fn55jr5+fkPpFr/5550mMUGAMCBWRMBsaEf1PzkrS2dU35+vpeXF2OsdevWx44dKysrY4wVFRWpM39LbiAAAIA9WRuWefvtt7t37167du19+/Zt2LDhl19+YYx17NgxMjLylVdeefbZZxctWhQXFxcVFWXP1gIfNx3KiHBjObFGLjUnmb6UX17KC4KWFvPToVR5STEvNSoSMWWM6XkpUyo1Sj1MaulgYoo1bl3RKdb4lem5x4iYpVBqVCQI6l6D//lQo4y/Ee5zW27gb8SD2LjJg/8kuvK+tbtQX+V1/I1ovgCsU9B+gCfM2hGhyWR6//33p02bdvLkyX379pmnANi8eXNxcfHTTz/t5ua2bt06u7UTAAAcgQOeI7R2RPiXv/yFWx4YGLhy5Up57QEAAKhUmGsUAACkUeiVVch/cZQRIQAAwB9zwHOE6AgdCZnRIC7R5OZiqLQIFUWhci7F9znlJUV6buUSXmXGWAlv46W8BA2jQzTcdQoNev5GyHUKZaxHKDTHmqz1CLlhGVc3saUEubkYD0/iOSzjf27wwzLEE2vizV3HGDOZJHwo0SEaTpnmw5Fqxpb1CLV+CXBtHwAAODWMCAEAQBpbpkzTekiIESEAADg1dIQAAODUcGgUAACksSUsY5+WWA8dYRXFjSVSYUVyTV3u2rlEOpQKdhYTQdDie5xysnKRQGq0jMipUgv2clOj3NV6Gb0Ab7mBeG65C/PaNTVKLswrsGAvtQAvlRrlT1NHpUN5Tzgjnltq56RTo9xiMdRUasRTS5Rq/vHsmBQmfI7QYZZhAgAAqJYwIgQAAHlsuKBe68E3OkIAAJCnCkyiLQodIQAASGPDdYQ4RwgAAKAljAg1Ri/xyvkDmQ4lkpDcNXXJuUOJwOf9u0Q5NzXKK2SMFd/nl3MbQ841SqVGeeXUnKLUc2WUkRoVXJhXLDWqI1KjrrzUqKsb/7miFublz9dKJmwFFjE2EhPhiqZDhZ5DqtxV4TxXOv5TwpAmdR7oCAEAQBqbDo3aqS3WwqFRAABwahgRAgCANDbMLIPLJwAAoHrRumMThY5Qc9QEVJxCoanUGGOlvJVsuQvqMjrnwg3FMMbu3y3jbISaYo2401JeeSkxCZyUsIyTTLHGXa2XMeZGZYi4a+pSz5XIrGnkc0WgBhMK72nhFjJ6cWPuc0ttRFFEXk74L0URvhwCl08AAABoCSNCAACQx4bVJ7QeY6MjBAAAeWyYa1RrODQKAABODSNCAACQRmEKtR5kRf+kKXSElYSeSo1fzg0rUlNeUXFK7gRmJeTauQIL8FL1qQV4qY3zp1gjUqPcdChVrhdNjRJJSG5S176pUSLx6EKVcxfmJVKjdGiWu6YuMZUasaYu+fB56NAsvz73aRF6TqiN0KlRgU9nzc9yVR22rFCv9bOHjhAAAKRysAEhzhECAIBzw4gQAACksWXSbTs1xWroCAEAQBpHnGsUh0YBAMCpYURYaQSmZ2REYM9QRkzCyZtTlFGpUWquUSrwSaypy50+lEqHUneK1OjvSUmNllOpUeJhcsupnVNs9lCRuUNZBQ9fZEpVspw3Xyv1xOqIAKvCX8hX60FN1WHDBfVaP3noCAEAQB5bzhHiOkIAAKguHHBAiHOEAADg3DAiBAAAaWy4fELzqWXQEcrHjRIITaXGiEyHvoyYYo1Il5QUc8pLePmUCsq5eRZGrakruBHunVIPhyrX8zJEwgvzVrOwjEEsLCO2/jCB+4DIdYYFZ0dzdeeFZdz5wRU3d/6u4sar7+rG3yWoF4KK+fAra37UTxOOdmwUh0YBAMCpYUQIAADS2DLptn1aYj10hAAAII2iiC3cwarAAWR0hAAAII8DXj+Bc4QAAODUMCKsJGRqlEjxGfSccm48ktEL83KnXivlRUkZY2Vk4JNY9ZeX4aRmR6PKuUFQcoo1YiY5bpiWSo1SixsbDRLilNSrzD3yQ69MS5TzpgdjjLm6cu6V2q9M/Ecvhmo5N06pc+F/2ybToa7EfsgLfLp78CuX1eCnSd157yBulJQx5sJ7YhljOl51zQ/uVR0OOCBERwgAABI54HWEODQKAABODSNCAACQSvNjnYLERoR6vT4xMXHMmDHmkrNnz/bu3TsgIKBfv34ZGRmymwcAAI5EvY5Q9EdbYiPC+fPnX7lyxbxMl8lkevzxx8eNG7dx48aFCxeOGjXq8OHDdmhkVUXNP8WLTNDrDhJhGRlTrHHLqcULycSNSLnQkoFUC6lQDNUSobAMOcWajLCMnCnWqDnJyvnl3F2ImrpPylKCQrOmkaEYYslAN95UaowxN95eUVbCz7mQs/F58nYVYpdwK+e3hP9eJj/KqfUYqfoOz4a5RoXPKcomMCI8derUli1bXn31VXNJSkpKbm7urFmzfHx83nrrrfPnzx89etQOjQQAALAXaztCg8HwzDPPfPrpp25ububCs2fPtm7d2sXFhTHm4eERFRV15swZuzQTAAAcgiL+ozVrO8KFCxfGxcV16dLFsjA3N7d27drmX729vW/fvk1twWiUcQUTAABopLycf8zZki0nCLXuC63qCC9fvvzZZ59NmjQpIyPj5s2bpaWlGRkZJpPJx8fn3r175mp37tzx9fUl70mHSzUAAByYevyv+rEqLHP79m1/f//JkyczxvLy8m7cuDFy5Mj9+/c3bdr07NmzJpNJURSDwXDhwoXw8HA7NxgAAKoudZgn+C8as6ojbN++fWpqqnr722+//fDDD9Vfe/furSjKV199NXny5C+++KJ+/fpxcXF2bGwVQ4ZGuYXUArxEarScF3qkplgjy7nBTiJ6SpaLBEFFU6Pclgu3UMrCvEKpUeLVFEKt70plMo1Gqtxea+oqxEEcFxmpUT2RGtV78Mcc3FdfaL+iyrnvNUa/N125Tzixii9F849+O7LDHGsmkyk1NbWgoCAuLq5WrVrcOhkZGRcuXGjevHlYWJhleVZWFmMsKCjIXPLbb7+Zj/F6eXk1bNhQ+HBlrVq1zFt0dXVdt27dhx9+6OXltWLFijVr1mieggUAgOrEYDAMHDhw4sSJH3/8cWRk5Llz535f57PPPouLi1uyZEn79u2//PJLtXDdunV+fn6NGzceOnSoZeV27doNHjx45MiRI0eOXLRoEbNhZplBgwYNGjTI/GtcXJz56KjopgAAoNoRv46wwiHhjz/+ePHixePHj3t6es6cOfOdd95Zu3atZYXCwsI33ngjOTm5bdu2Bw4cGDhw4JgxYzw9PTt16vTvf//7wIEDn3766QPb/P777yMiIsy/ygmwoBcEAABmh9Topk2bhg0b5unpyRgbN27c5s2bH7gGYceOHaGhoW3btmWMdenSxdvbe+/evYyxRo0aUbGVy5cvp6WllZSUqL8iyQkAAPLYcB1hhR3h1atXQ0ND1duhoaGlpaW3bt2yrHDt2jXL84JhYWHXrl2rYIPu7u4zZswYMWJEUFDQxo0bGSbdBgAAzV2+fHnZsmWWJa1atVLTl6WlpeaJXNzd3RljxcXFljVLS0tdXf+/L3N3dzcP9bjOnz+vXgG/du3aCRMm9OrVCx2hfNy8HhUzFJprlEpC6kXKDXqxoB25wi031yrYQu5GhO5R1kaoF0IoNSoy1ai01KhQPlRoglMdkYTUu/KfQxfec+sq+ELQ+yEvYCxlf6PWaqZ2CV6xTmRN5urNlrlGmVJQUHDkyBHLQh8fH/WGv7+/eaqWW7duKYoSEBBgWdPf3z83N9f8661btx6o8ADzPDCjR49+8cUXT506hY4QAACkse3qiTZt2ixdupT71y5duuzevXvmzJmMsX379rVt29bDw8OyQlxc3HPPPVdYWFinTp3c3NwzZ8506tTJmvu9cePGnTt3GjZsiI4QAACqrokTJy5YsGD27NnNmzefOXOmesEDY6xr166jRo166aWXmjVr1rdv31GjRk2aNGnZsmXDhg1TzylmZmauX7/+2LFj169f/+CDD8LDw4cPH75///41a9Z06NChuLh48eLF/fv3j4qKQlgGAADkkR2W8fPz+/XXX8vKypKTk5cvXz5q1Ci1/KmnnurcubN6e82aNY888si2bdsGDBiwcuVKy3+PjY194YUXzL9GRESEhobu378/LS1txowZCMsAAIBkig0X1P1R9fDw8IULFz5QqM76qfL09JwxY8YDFRo3bqweULXk7+8/a9asBwrRET4ManVWkYV5iUU5uDOBkXNByUiRUHOPUVECfgtFKlPlQvdIbkR0ijUiGWHiLnsrOIEZl474sDASERUXk8CHC/VBRCVxDNzZ0Qz8ykIvhJT9ionuzCLl1D1S703+e1l4l6i2KRobVpzX/LnAoVEAAHBqGBECAIA8dph0297QEQIAgDw2XEeo9eWWODQKAABODSNCAACQxoawDA6NOjKB0Ci5UKqRP58Uf3onchVfco4oXkaO3AjRQqKcW5/aCHmn3EymyD1S5dT8WFRLqDsVWvaWSg5yPxeoECg1WRdFUQReCB2xdRcXkf1NZBeSsl9Rd8rdw5lgNpicSo14b3JffTI0KjTtHmgEHSEAAEijyF6PsBLgHCEAADg1jAgBAEAaWy6o1/pAMTpCAACQxpZlmLTuCXFoFAAAnBpGhJWECpWRaVJeWFFo8ViqvvBGRMqFWkK20J4t4c4dKrpxSQvzErWpjVABY14QVM5zKPKqkRuR0RLqTuW8I6hXUzAbDP9hw8wyWkNHCAAA8uA6QgAAcGYKU0Qvh9C6H8Q5QgAAcG4YEQIAgDxYfaJ6onIuQvXJE+8CJ+TJyvz5pBixgKhoSwTKtdmISJ6FWlNXbOOiMQru+5x41aip14gHxHRCcZ6q/2qKbZxbl3xHiG2E/gOvkKjLL2a8efHUPxDljkMRvy5Q8weNQ6MAAODUMCIEAABpFEX8AnmtL6hHRwgAAPLgHCEAADgzW+YatU9LrIdzhAAA4NQwIpSPHxqlKov8gUyxyWmK9Q2RtnFuufAapw9bt6I7FYkAC26c+BosvHEZtfl3KuPVlLNfEeWi7wjiYYplgIWeKieESbcBAAAcDDpCAABwajg0CgAAMmHSbQAAcF62nCO0U1Osho5QPv7ic1RlkT9QlcU2TtYW2QhVX3Dj3HKxe5RQt6I75RcL5ly4GxfeJaj6MmpL2VWEdgnhV1lkVxF8RxAvssjGNf8or0Ic8DpCnCMEAACnhhEhAABI44iXT6AjBAAAaWyYWUZzODQKAABODR0hAAA4NRwatQIVHhNaWpNMpgkk1sjKxPcZsY3IKNdmIzreRniFjDEdsRFqOVz+c0us+0ruKvwtiz1MquVCD98BXk2xjXPrVvCOENkI/QdeIVGXX6x9TtJ+HPEcIUaEAADg1DAiBAAAmRxtXV50hAAAIJEDXlBvbUe4c+fOzZs35+Tk+Pn5Pfnkk926dVPL79y58/77758/f75Vq1YzZszw8vKyW1MBAKCqs+EcoeasPUeYkZERExMzadKkyMjIfv367d+/Xy0fMWLEhQsXpkyZcuTIkQkTJtitnQAAAHZh7Yhw6tSp6o3ExMTDhw/v2rWre/fup06dOnDgwM2bN2vWrNmpU6eAgICMjIwmTZrYrbUOTDSZpuOl/nQuApWp+sIbESkXagnZQnu2xEi1RGjVXylzjYq89IwxReSFk/Mcirxqdm0Jdady3hHCCVtuMfyHDRfUa/6UiqVGy8vLT548eejQoR49ejDGDh8+3K5du5o1azLGfHx8WrRokZqaapdmAgCAQ1Bs+tGUQEe4dOlSb2/v1q1bDx48OD4+njGWk5NTr149cwU/P5UejZsAAB9YSURBVL+cnBzq341G6torAABwANX1Y1ygI3z22Wfv3r2bkZGxb9++jz/+mDHm5eVVWlpqrlBcXFxBWEanwzWLAAAOzJqPcYUx5T/HRwV+KqHxFRDunBo3bjxu3LikpCTGWHBw8OXLl81/+u2330JCQiQ2DgAAHIvCHK0btD4sc/XqVbWTKy0t3blzZ6tWrRhjCQkJkydPPnjwYOfOnZOSkkpKSnr27GnHxlY1QhMwUckIF/5GuGf1XYhogIsr/wuNzoVTTm+EaCFRzq1PbYS8U254gdpIOZF/MXLKuYWMMRcq6EJQynn3KJKsoZBTpok8V4zaVUReNUY850KvGrVxKfsVdafcPZzR7whiI2LvTaF52rT/jK98tl1HKOM9ZTNrO8JevXrVqFHD19f33LlzMTExb775JmOsVq1an3zyycCBA1u3bn3ixInPP//cw8PDnq0FAACQzNqOMD09PT09vaCgICQkJDg42Fw+efLkgQMHXrp0KTIy0tfX1z6NBAAAB2HDBfWK4hgjQhcXl6ioKO6fGjRo0KBBA3lNAgAAR2XLdYT2aYn1kOQEAACnhkm3AQBAGpvWI7RTW6yFjvBhCKzFSa4TS1yXw429ubgRlYlyV145t5C6R8aYKxXA47ZQpDJV7urKP11gLCfKeSFOarYzikKss2zkzb2mI2KjYks1y5hjjBEvqOgLwX2VpbyaUvYrJrozi5SToWvivcl/4RAbdWToCAEAQBobzhFq/m0B5wgBAMCpYUQIAADSKEz8HKHWQ0J0hAAAII8DrlCPQ6MAAODUMCKUj3tUgEgC0lFAbpySyMK5iZS7uvGnUHRz55e7uvFm2yQaQ7WknCo3cMqF0qGMMZNIQpQ6ZlNuEImkyphsVDQ1Sk8fKrCrCJVL2t/EWkLvh5xyoZZQd0rlV8k5SIVCo87HERfmRUcIAADy2DDFmtbHRtERAgCAPDhHCAAAIN29e/dycnIqqKDX67OysgwGgzVbKy8vz8rK0uv16q/oCAEAQBrFpp+KzZo1Kzg4uH379l26dLl169bvK2zfvj04OLhHjx6hoaF79uxRC5OSkh555BFvb++OHTtaVv7Xv/4VFhbWo0ePwMDAzZs3MxwafRhCo3kqGUGuccqNohA5ArLcg5cvcKdCCkQ5byOMMTe98feF5QZOIWOsnMi/lPOjKNy6oqEYwQnMeFOpMXL+NqIlInOsibZQSlhGaBciX3qRcuGNiOyf9EYEysn52Ij3Jve9jAnWzBR63yb/pcLqKSkpq1atOnfunL+//5gxY955553PP//cskJZWdmkSZOWL18+aNCgdevWTZo0KSMjw8XFxdvb+6WXXrpw4cJ3331nrmwymSZPnjx37tynnnpq9+7dw4cP79u3L0aEAABQdX3zzTcjR4709/dnjL300kvffPPNA99Ed+3a5enpOWjQIMbYiBEjysrKUlJSGGPt27cfOnRoQECAZeXDhw/fuHHjySefZIzFx8cHBQX99NNP6AgBAEAe2cdGMzIyIiMj1duRkZF37tzJy8uzrJCZmWmuoNPpmjZtmpmZSW0tMzOzSZMmbm5u5g1mZmbi0CgAAEhjy6TbjN28eTMpKcmyJCIiIiwsjDF29+7dmjVrqoVeXl6MsTt37vj6+pprFhYWenp6mn/18vK6c+cOdUfcyugIAQBAGpvWI1QuXLjwwQcfWBY+9thj06ZNY4w1aNCgoKBALczPz2eMNWzY0LKmZQW1zgMVKq7cvn17dIQAAKCxrl27btq0ifunmJiYQ4cOqbcPHToUHh6ujgvNWrVqdfz4cb1e7+bmVlxcfPr06ZiYGOqOWrZsefHixYKCAm9vb6PRmJqa+tprr6EjfAjklx5eqIw4Gys0xRoVqHOvwc/Iccvda/BfdH0pfyo1fRk/xGmQkRrlz5pGhkP5f+B+/aSylwYqpusiITVKhUm534+Fc61EalRsNj4iTsnfVWSkRqmNCJfz9luhPZ8R7yDRKdbEFuZ10tio+L/QnnrqqdatW69du7ZZs2Zvvvnm888/r5ZPmDAhISFhzJgxnTp1ioyMnD59+tSpUxctWtShQ4eWLVsyxm7cuJGSknLo0KG8vLwNGzYEBAR069YtPDw8Pj7+5ZdfnjVr1tdff+3n59erVy+EZQAAQCL1LKHYTwWba9y48aZNm7788supU6eOHj365ZdfVsuDgoK8vb3V25s2bcrNzR03blxpaem6devUwhs3bmzYsOHGjRtt27bdsGHDL7/8opb/4x//cHd3Hzdu3G+//bZ161ZFUTAiBACAKi0+Pj4+Pv6Bwvnz55tvBwUFrV69+oEKMTEx69ev//3W/Pz8li9fblmCjhAAAKRRqsBqEqLQEQIAgDwOOOk2OsJKQn1FIsMybpxyMukgki/w8KRCMfydwUDlX3jlxnJ+ZWoBP6FZ0yiKjvOIqCfWxcAvF1uP0K5TrJFxHhlTrIlEVKjIiYdIOVnZk2iJJ38/5Nankl9kzIf3DuK+11gFYRmtP7WrONsun7BTY6yEsAwAADg1jAgBAEAaW1aot09LrIeOEAAA5HHAc4Q4NAoAAE4NI0IAAJDGhrCM5iNCdITyCe0D9CqsEqZYq1HGKTdQ6VBiKjWhWdP4U6bRc49xCa9Yy0v36V0lPBzGmJFXXvVTo26CqVFuOR34pIKdvJRyTYHKjLEaIuU1qOipyBRr3CeQ0fsbMWEet64zcsAjozg0CgAAzg0jQgAAkErzIZ4gdIQAACCPA15Qj44QAACkseU6Qq1HkDhHCAAATg0jwkpDLfLJr61z4U0gSa6qyk8rctfO5RYyOk7JjU0yavpQwalD+QE8kXQoVe4i+DCr0FyjZJC4qqdGPb3cOIVEatSzJqcyY8zTi0iN8rZDtYSag5T7DuK+1xi9HzreGbBKZkNsVGvoCAEAQBpFET/np/WxURwaBQAAp4YRIQAASGPDwrxaDwjREQIAgDwKUxTBk4Ra94PoCCuL6Fce7jxbrsRcUNQypOUGzusrNJcYY8zID5cIrqlLrUvMyyNQ4QVqKiz+yrTOHZaRsjAvGZYh8y+ccm6ChjHmWYsKy/DLa/DKybAM8TC57yDRBXg1H75UdQ44xxrOEQIAgFPDiBAAAOQRv6Bec+gIAQBAGhuWYVIURfAKZMms7Qjz8vJ27dp15cqV0NDQwYMH16hRQy03mUw//vjj2bNnW7du3a9fP7u1EwAAwC6sPUfYrl27NWvW5ObmfvHFF23atCkoKFDLX3zxxTlz5ty/f3/atGmzZ8+2WzsBAMARKDb9aMraEeHhw4f9/PwYY+Xl5S1btvz+++8nT56clZX15ZdfXrx4MSgoaPz48TExMa+99pqvr689G1z9UIk1zqECKt5GRQQ9eBNN8adGq6Bcxpq6QrOmUelQ6mHqeeXUTHIGaiY5GalR6rkilnIVS43qiNQoNwkpJzVKLHtLTrHGmzWNSofWpFKjVDl3ijVqKjXi4XPfQfRhPK0/nh2TLZNuC0/OKJm1I0K1F2SMubi4uLu7u7i4MMb27t3bqlWroKAgxljTpk0bNWq0f/9+OzUUAACqPuU/XaHAj+bfOYQvn1i/fv3t27cHDx7MGMvOzg4ICDD/yd/fPzs7m/pHscvOAACgiqmuH+NiqdHk5OQXXnjhhx9+8Pb2ZmrUx+J5MZlMmq+vCAAAGnO0C+oFOsIDBw6MHDly7dq1cXFxaklAQMD169fNFXJycgIDA6l/Rx8JAODQrPkYr84L8x49evTxxx9fuXJl7969zYW9e/c+ffr0tWvXGGPp6elXrlzp0aOHXZoJAABgH9aOCAcMGODl5bVq1apVq1YxxoYOHfrEE08EBgY+++yzffv2HTp06Pr161977TUfHx97trYaEkqs6fghPuZq4n+hMXlwE4/8F1304D+35TpynWGReTJdDdzKbu5EapSXeNSXlXMrO8lco27EGs7Uc8hNjboT6dAa5MK8dpxrlBtVdfMgAsZE8JhIjYrNNQoVs+WCejs1xWrWdoRffPFFefn/f7JERUWpNxYtWrRjx460tLSlS5daDhYBAMAZOeCk29Z2hI8//jj1p4SEhISEBEntAQAAR2bDXKNad4RYfQIAAJwaJt0GAABpFPFrBEQX8pUOHWEVRZ2/55aSIRoZI35yeVL+mroCgQ7GmCsvu1FKBT1K+CEafRkn/2LQ88My1NRrQusV2zcsQz2HIhPsuboJhmV4c5VRU6xRYZkavFnQuAvqMmLKNFbBWru8VYKFplJjxHOOUAzg0CgAADg1jAgBAEAa29YjtFNjrISOEAAApLFhZhmtTxGiIwQAAIm0X0xCGM4RAgCAU8OI0JGQBxyI7zMuclKjAmvtCiUbGWOuvICouwc/HVpWwg8x6ks5AVE9tTAvUW4sp8q5qVFuXdGFefmVqcSjzkVgvWI30YV5eZlMbpSUMVZDZMFeKgVKralLzprGe0QuxHOiEDu+1qeinIJN5wjt1BZroSMEAAB5qvHqEwAAANUSOkIAAHBqODQKAADSKMyGZZhwHSE8NHI30nHCG1SChtqI0Exg5FRqVHbDnZNzKeMFNxhjZSX8WdO4Sw9SYZlyqpxYp5C/HiG/rtAMa2Sgg16PkP8PLkJhGWL6On5YRiRZw4hwDbURMj8ltJQgFYrR+oPVmTnidYQ4NAoAAE4NI0IAAJCnGi/MCwAA8Iccca5RHBoFAACnhhEhAABIY0NYRusjo+gIqwdq7VzuH6iwIrERVyKZxw03UtODubryy7khRndeCpQxpvfkhzW5qVFyAV57LswrOMWanIV5ualROqYrsGCvUGVGTJgnlAKtoBxr6joSRztHiEOjAADg1DAiBAAAmRzuOk50hAAAII0t5wi17jfREQIAgDS2XD6h9QgS5wgBAMCpYUTodOjvalSKj6jNS41yC1lFC/ZyMpxUWNFATAfKDYKSlak5RanUKHeuUWpSUZHYKLkwL/EckmlSXiyTympyI6ZUfVc3sVlPuUsHUw+HXjuXSo3y60OVg5llAADAmWHSbQAAAMlWr17drFkzf3//qVOnlpaW/r5CWlpaz549fX194+Pj09PTzeXvvfdeWFhYSEjI3LlzzZf/xsfHt/+vN954g2FECAAAEklfjzAtLe3555/funVr8+bNhw0bNn/+/Llz51pWMJlMw4YNmzhx4tatWxcuXDhq1Khjx44xxjZv3rx48eKkpCR3d/dHH300IiJizJgxjLETJ06sWLEiJCSEMebj48MwIgQAgKps5cqVw4YN69atm5+f35///OcVK1Y8UCElJSUvL+/111+vVavW7NmzMzMzDx8+zBhbsWLF888/HxkZ2ahRo5dfftnyH6Ojo9u1a9euXbsmTZowjAjBTDREoyicYAg5bZhOYJ4tF1d+5MStnL8R7uxodPiFCstwi/mzqdl1ijV66jX+xrnPreg8bdz65CxoZP5FaBY0hGKqJ+nXEZ49ezYhIUG93bp16+zs7MLCwjp16lhWaNWqlYuLC2PM3d29RYsW586d69Chw9mzZ1944QXzP86fP9/8L8OHD1cUpVOnTu+8805AQAA6QgAA0FhZWVl+fr5liZeXl7u7O2MsNzfX3O3VrVuXMXb79m3LjjAvL6927drmX+vWrXv79m213PIf1ULG2KJFi2JjY0tKSt57772EhIQjR46gIwQAAHnEL6hnCtu1a1d4eLhl2dNPP71gwQLGWL169e7evasWFhYWMsZ8fX0ta9arV+/evXvmX+/cuaNWeOAfzf81btw49cbq1at9fX1PnjyJjhAAAOSx4TpCxh577LFNmzZx/xQREZGWlqbeTktLa9CggTouNGvatOmZM2dMJpOiKAaD4fz58xEREWp5Wlqaelg1LS2tadOmD2zZ3d3dzc2ttLQUYRkAAKi6Jk6cuGHDhrS0tPv373/wwQcTJ05Uy999992dO3cyxh555BEPD4+lS5eWl5d/+umnDRs27Ny5M2Ns0qRJixcvzsnJuX379qeffjpp0iTG2KVLlw4cOFBaWlpYWDhjxow6deq0adMGHSEAAEijhmVEfyrQtm3befPm9enTJyAgwNvb+6233lLLT58+nZOTwxhzcXFZv379kiVL6tSp880336xdu1Y9Njtq1Khhw4ZFR0dHRkb27dtX7UHv3LkzZcqUunXrhoaGnj59euvWrTVr1lSo/JtcxcXFz/9lwcCJUyrhvkBD9N4kELOkNmLizXZGlRN16Y0IBUGpyvz7FMtHUh8KZMqUu0IytREZgU+hIChSoNXJ480C/rDOP/dl3LpTIrTZ0//enX92D3VotBLgHCEAAEjlaN9+cGgUAACcGkaEAAAgjU2hUY2hIwQAAGlsWZhX6zPJODQKAABODSNCkEl0YkmxjRDlJu48mYJr54qlp6lIKlGdSGqK3GNFT4tI4FN46w9ZF5wPFuYFAACnpggf6tT80Cg6QgAAkMaW1Sfs0xLriZ0jzMnJUa/kt3T9+vV9+/bduHFDXqsAAAAqibUd4bJlywICAgIDA59++ukHylu1avX+++9HR0d/++23dmghAAA4FEX8R1PWHhrt1q3bnj17fvrpp+TkZHPh3bt3p0+fnpSU1LFjxz179qgTu3l4eNinqVANiZ5KsNuW6aDLw9YlCb/3pQSRJGwDoCLq7KFat0KMtSPCFi1aREVFPbAW9o4dO0JCQjp27MgY6927d82aNfft2ye9iQAAAPbzUGGZq1evNmrUyPxrWFjYlStXHrZFAADgsGwJy2g9gHyojrCkpMTd3d38q4eHR3FxMVXZaDQ+zH0BAIC2jEbjA8cFOarAOT9RD9UR+vv75+bmmn+9fft2QAC5SMcfP30AAFCFWfMx7nRTrHXs2PHIkSNFRUWMsby8vLNnz3bo0EFSwwAAACqDtSPCtLS0n376af/+/ZcuXfrggw9iYmL69+8fHR3do0ePsWPHTp48efHixYMGDbI8ZQggmei0YXbbiqJNbBTAAdhwjlDz94K1I8KysrL8/PyWLVsOHjw4Pz9fHQUyxtavXx8bG7t69eoePXqsWrXKbu0EAABH4GgXETLrR4SxsbGxsbG/L69du/bbb78ttUkAAACVB3ONAgCARI4XlkFHCAAA0jjgKkzoCAFsoPkbF6DKcsCeENf2AQCAU8OIEAAApLHhgnrNoSMEAABpHHGuURwaBQAAp4YRIQAASKMowpdDaH4oFR0hAABIpfWhTlE4NAoAAE4NI0IAAJDJ0QaE6AgBAEAeR1yPEB0hAADIg5llAAAAHAtGhAAAII0jXlCPjhAAAKRxxCnWcGgUAACcGkaEAAAgjwOGZdARAgCANIr4OT/Nj6SiIwQAAGkUpiiCQzzR+tLhHCEAADg1jAgBAEAeG84Rag0dIQAAyCN+HaHmHScOjQIAgFPDiBAAAKTBpNsAAODccI4QAACcmcJsGBHaqS3WwjlCAABwahgRAgCANA44wxo6QgAAkMgBe0IcGgUAAKeGESEAAEijKDZcDoHLJwAAoNrACvUAAODUcI4QAADAsWBECAAA0thyQb3WQ0J0hAAAII2C1ScAAAAcCzpCAABwajg0CgAA0ti0DJOd2mItdIQAACCNDecIte4HcWgUAACcG0aEAAAgjwNeUI+OEAAA5BE/R6j5SUIcGgUAAKeGESEAAEijiA/wtB4QoiMEAAB5bLl8QuuThJXUEbq6up5L3pG8/h+Vc3eVLDc318fHR6er5seZ8TCrk9u3b/v6+oqvG+dg8DDlivrpp6ioqIrrxDfyE91sk+LYX3UltjZKAsVkMlXOPeXl5RUUFFTOfQEAgHTBwcHu7u5at0K+yusIAQAAqqBqfvwHAACgYugIAQDAqaEjBAAAp4aOEAAAnBquIxR248aN1NTUrKys+Pj48PBwc/lvv/329ddfFxUVjRgxokOHDhq2UIqjR4/u3Lnz1q1bUVFRY8eO9fT0VMsLCwuXL1+elZXVq1evQYMGadvIh5eUlHTgwIGCgoLQ0NDx48f7+vqq5QUFBcuXL79+/XqfPn0GDBigbSMlWr9+vYeHx+DBg9VfS0tLv/zyy4sXL8bGxo4dO9bRLxrZsmVLTk6OertevXrDhw9Xb+fl5a1YsSInJ6dfv359+/bVroHS3Lx5c9WqVdnZ2Y0bN54wYULdunWZxU4bHx//2GOPad1GR+LY+70munfvPn/+/JkzZ6amppoLc3JyOnToUFBQ0KBBg0cffTQlJUXDFj68goKCxMTEW7duhYaGrl69unv37qWlpYwxo9HYu3fvAwcOhIeHT5s27W9/+5vWLX1Y69atMxqNTZo0+eWXX9q0aZOXl8cYMxgMPXr0SE1NbdKkyZ/+9KclS5Zo3Uw5tmzZ8swzzyxYsMBcMnr06O+++y4iIuLjjz+ePn26hm2TYsGCBTt37szIyMjIyLh27ZpaWFZW1q1bt5MnTzZu3Hjy5Mlff/21pm2U4MKFCzExMWlpaY0aNUpPT1cfqcFg6NmzZ2pqanh4+PPPP7948WKtm+lQTCCovLzcZDK1bt167dq15sK5c+cOGTJEvb1gwYIBAwZo0zhJysvLS0tL1dvFxcV169ZNSUkxmUzbtm1r1KiRXq83mUy7d+8ODAxUb1cDRqOxcePGmzZtMplMmzZtioyMVF9o9SGrtx1aQUFBdHT03Llzu3TpopacPn3ay8ursLDQZDJdunTJ09MzNzdX0zY+rK5du27ZsuWBwjVr1rRs2dJoNJpMpk2bNkVERKi3HVe/fv1mz579QKH60NQddfv27WFhYQaDQYvWOSSMCIVxDx+lpKSYD7k8+uijycnJldsoyXQ6nfmyWaPRWFZWVrt2bcZYcnLyI4884urqyhjr2bNnbm5uenq6lg2VJz09vaCgoHnz5oyx5OTk+Ph49YWOj4+/cuXK5cuXNW7fQ3vllVdeeeWVwMBAc0lKSkqnTp3UV7ZJkyZBQUGHDh3SroFybN++/aOPPtq2bZvpv1dIp6Sk9OnTR511pW/fvhcuXMjOzta0jQ9Fr9fv2rVr8ODBK1euXLJkiXng+8BOe+3atWqw01YadIRyXL9+vX79+urtBg0aFBUVFRYWatskWV5//fUePXq0adOGMZaTk2N+mC4uLr6+vtevX9e0dRLMmDEjKCgoJiZm4cKFakdo+Wq6u7v7+Pg4+sPcvXt3Zmbm5MmTLQstX03GWIMGDRy6h2CMtWjRwsPD4+bNmy+++GJiYqLRaGT/+2rWrFmzVq1aDv1qXr161Wg0Pvfcc5cvXz516lTr1q3PnDnD/vfVdHNzqwY7bWVCWEYOV1dXg8Gg3lZvuLm5adoiORYtWrRr1y7zKU9XV9fy8nLzX/V6fTWYb+ntt99+9dVX9+/fP3Xq1FatWnXo0MHNza06PcyioqIXX3zx+++/f2Auyur3ai5btky9MWvWrMjIyJ07d/br18/yvckYMxgMDv0wdTqdyWR67rnn1K81er3+o48++vLLL6vfq1mZMCKUIygoyPxtOisrq169euaYpeP6+9///tlnn+3Zs8ff318tCQoKysrKUm8XFxfn5+dbHmpzUF5eXv7+/iNGjOjXr9/mzZvZ/z7Mu3fv3r1716EfZnJyclZW1rhx49q3bz9v3ryTJ0+2b9/eaDRaPkzGWFZWlkM/TEs+Pj4tWrTIzMxk//vevH37dklJiUM/zICAAJ1O16JFC/XX6Ojo3377jf1upy0sLHToh1nJ0BHKkZiYuHHjRvVQzIYNGxITE7Vu0cNasWLFxx9/vGvXruDgYHNhYmLirl271MnTN23a1KxZM8sLSByOwWDQ6/Xq7bKyspMnT4aGhjLGEhMTf/7557t37zLGvvvuu9jY2KCgIC0b+nC6du26Z8+epUuXLl26dPz48U2aNFm6dKlOp+vXr9+xY8fUj9F//etfpaWlXbp00bqxttPr9eaR39WrV48fPx4dHc0YS0xM3LZt2/379xlj3333XVxcnJ+f8PIIVYeHh0f//v0PHjyo/nrw4EG1U0xMTNyxY4e6027cuLFNmzaW71z4A1qndRzP888/365dO09PzyZNmrRr1y41NdVkMt27d69t27Y9evQYOXJkw4YNz58/r3UzH0pWVpaiKKGhoe3+a+vWreqfxo4d26JFi4kTJ/r5+f3444/atvMhXblypWHDhkOGDBk7dmxYWFifPn2Ki4vVPw0fPrxly5YTJkzw8/PbsWOHtu2UaPny5ebUqMlkmj17dqNGjSZPnuzv77906VING/bwLl26FBAQ8Pjjj48cOdLHx+e5555Ty41GY2JiYps2bZ588klfX9+9e/dq2kwJjh492qBBgyeffHLAgAFNmzbNzs5Wy0eMGGHeaX/++WdtG+lYsPqEsAsXLlgGYSIjI9XcXWlp6Z49e+7du9enTx8fHx/tGihBWVnZqVOnLEsaNWqkXmxuMpn279+flZXVpUuXsLAwjRoozZUrV44fP15SUhIREREbG2suN5lMycnJOTk5Xbt2DQkJ0bCFct2+fTs3N7dZs2bmktTU1AsXLrRp0+YP15mr+s6ePXv27Fmj0RgTExMZGWkuNxqN+/btu3nzZvfu3R16cG+Wm5u7Z88eb2/vbt26mc/CmEymlJSU69evV7OdthKgIwQAAKeGc4QAAODU0BECAIBTQ0cIAABODR0hAAA4NXSEAADg1NARAgCAU0NHCAAATg0dIQAAODV0hAAA4NTQEQIAgFNDRwgAAE7t/wDKS2+FAfyCnAAAAABJRU5ErkJggg==",
"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
}