{
"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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter Function value Gradient norm \n",
" 0 3.029727e+01 7.277087e+00\n",
" * time: 0.004927873611450195\n",
" 1 2.714988e+01 6.090159e+00\n",
" * time: 0.014632940292358398\n",
" 2 2.710517e+01 1.353371e+01\n",
" * time: 0.03787398338317871\n",
" 3 1.732697e+01 5.975572e+00\n",
" * time: 0.06557393074035645\n",
" 4 1.242066e+01 2.145855e+00\n",
" * time: 0.09319901466369629\n",
" 5 1.177662e+01 3.740385e+00\n",
" * time: 0.11169099807739258\n",
" 6 1.079531e+01 2.428641e+00\n",
" * time: 0.1308889389038086\n",
" 7 1.027271e+01 2.621279e+00\n",
" * time: 0.15029287338256836\n",
" 8 9.464550e+00 1.770236e+00\n",
" * time: 0.16956090927124023\n",
" 9 8.948249e+00 1.327452e+00\n",
" * time: 0.18882989883422852\n",
" 10 8.806555e+00 1.037333e+00\n",
" * time: 0.2080700397491455\n",
" 11 8.680631e+00 7.621556e-01\n",
" * time: 0.22749090194702148\n",
" 12 8.617569e+00 4.089333e-01\n",
" * time: 0.24246001243591309\n",
" 13 8.595580e+00 3.253824e-01\n",
" * time: 0.25714993476867676\n",
" 14 8.579926e+00 2.478894e-01\n",
" * time: 0.27170896530151367\n",
" 15 8.560862e+00 1.780827e-01\n",
" * time: 0.28620100021362305\n",
" 16 8.558200e+00 1.705515e-01\n",
" * time: 0.3007829189300537\n",
" 17 8.552687e+00 1.503738e-01\n",
" * time: 0.31540393829345703\n",
" 18 8.547576e+00 1.989331e-01\n",
" * time: 0.32993197441101074\n",
" 19 8.543197e+00 1.458953e-01\n",
" * time: 0.34454798698425293\n",
" 20 8.542894e+00 1.006022e-01\n",
" * time: 0.3593180179595947\n",
" 21 8.538131e+00 1.465741e-01\n",
" * time: 0.3745739459991455\n",
" 22 8.533205e+00 1.881328e-01\n",
" * time: 0.38982605934143066\n",
" 23 8.530288e+00 8.305813e-02\n",
" * time: 0.4052770137786865\n",
" 24 8.524636e+00 9.840524e-02\n",
" * time: 0.4207789897918701\n",
" 25 8.522453e+00 1.079799e-01\n",
" * time: 0.43636488914489746\n",
" 26 8.520491e+00 6.014887e-02\n",
" * time: 0.45166993141174316\n",
" 27 8.518853e+00 6.773988e-02\n",
" * time: 0.46700000762939453\n",
" 28 8.518007e+00 6.731495e-02\n",
" * time: 0.5736539363861084\n",
" 29 8.517302e+00 7.673625e-02\n",
" * time: 0.5877430438995361\n",
" 30 8.516241e+00 6.200554e-02\n",
" * time: 0.6017050743103027\n",
" 31 8.515526e+00 4.584172e-02\n",
" * time: 0.6155829429626465\n",
" 32 8.514833e+00 6.136310e-02\n",
" * time: 0.6294019222259521\n",
" 33 8.514602e+00 2.985688e-02\n",
" * time: 0.6479389667510986\n",
" 34 8.514205e+00 3.781997e-02\n",
" * time: 0.6617748737335205\n",
" 35 8.513718e+00 2.953043e-02\n",
" * time: 0.6756250858306885\n",
" 36 8.513113e+00 3.301718e-02\n",
" * time: 0.6895098686218262\n",
" 37 8.512782e+00 2.778672e-02\n",
" * time: 0.7033388614654541\n",
" 38 8.512547e+00 2.329645e-02\n",
" * time: 0.7172040939331055\n",
" 39 8.512357e+00 1.945132e-02\n",
" * time: 0.7309520244598389\n",
" 40 8.512216e+00 1.881472e-02\n",
" * time: 0.7444770336151123\n",
" 41 8.512169e+00 1.793595e-02\n",
" * time: 0.7583260536193848\n",
" 42 8.512130e+00 1.458037e-02\n",
" * time: 0.7719740867614746\n",
" 43 8.512057e+00 1.078752e-02\n",
" * time: 0.7857198715209961\n",
" 44 8.512022e+00 1.296468e-02\n",
" * time: 0.7993419170379639\n",
" 45 8.511978e+00 1.117038e-02\n",
" * time: 0.813060998916626\n",
" 46 8.511954e+00 6.353316e-03\n",
" * time: 0.8264319896697998\n",
" 47 8.511929e+00 8.156709e-03\n",
" * time: 0.8399059772491455\n",
" 48 8.511903e+00 5.590825e-03\n",
" * time: 0.8535819053649902\n",
" 49 8.511890e+00 4.846164e-03\n",
" * time: 0.8672900199890137\n",
" 50 8.511880e+00 4.933957e-03\n",
" * time: 0.8811249732971191\n",
" 51 8.511870e+00 3.279049e-03\n",
" * time: 0.8949649333953857\n",
" 52 8.511865e+00 4.160331e-03\n",
" * time: 0.9086959362030029\n",
" 53 8.511861e+00 4.479635e-03\n",
" * time: 0.9227349758148193\n",
" 54 8.511856e+00 3.315607e-03\n",
" * time: 0.9366018772125244\n",
" 55 8.511854e+00 2.701586e-03\n",
" * time: 0.9504499435424805\n",
" 56 8.511853e+00 3.002228e-03\n",
" * time: 0.9644579887390137\n",
" 57 8.511851e+00 2.254416e-03\n",
" * time: 0.9780700206756592\n",
" 58 8.511849e+00 1.556914e-03\n",
" * time: 0.9918179512023926\n",
" 59 8.511848e+00 1.387774e-03\n",
" * time: 1.005734920501709\n",
" 60 8.511847e+00 1.052732e-03\n",
" * time: 1.0196490287780762\n",
" 61 8.511847e+00 9.997167e-04\n",
" * time: 1.0338470935821533\n",
" 62 8.511846e+00 7.443073e-04\n",
" * time: 1.048198938369751\n",
" 63 8.511846e+00 8.770664e-04\n",
" * time: 1.0626740455627441\n",
" 64 8.511846e+00 8.408076e-04\n",
" * time: 1.0770630836486816\n",
" 65 8.511846e+00 6.707615e-04\n",
" * time: 1.0914089679718018\n",
" 66 8.511846e+00 6.600849e-04\n",
" * time: 1.105733871459961\n",
" 67 8.511846e+00 5.658636e-04\n",
" * time: 1.1200649738311768\n",
" 68 8.511846e+00 4.944989e-04\n",
" * time: 1.1344590187072754\n",
" 69 8.511845e+00 3.950944e-04\n",
" * time: 1.1489908695220947\n",
" 70 8.511845e+00 3.166125e-04\n",
" * time: 1.1636199951171875\n",
" 71 8.511845e+00 2.550969e-04\n",
" * time: 1.1782269477844238\n",
" 72 8.511845e+00 2.210601e-04\n",
" * time: 1.1927909851074219\n",
" 73 8.511845e+00 1.844223e-04\n",
" * time: 1.2072288990020752\n",
" 74 8.511845e+00 1.614021e-04\n",
" * time: 1.221898078918457\n",
" 75 8.511845e+00 1.975962e-04\n",
" * time: 1.2363760471343994\n",
" 76 8.511845e+00 1.624466e-04\n",
" * time: 1.2508509159088135\n",
" 77 8.511845e+00 1.355871e-04\n",
" * time: 1.2658488750457764\n",
" 78 8.511845e+00 1.218384e-04\n",
" * time: 1.2809700965881348\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVxTx/o/8DkJi0hREapEFhEViyiiuO8VK264F6vWjdrWLnb1trb6vV2u3az12lWttr1W27rUYjet4oqttYq7FEUFNxYXFKGyJSS/P05vfrk6DzJxwiHk837xRxiGySQ5yWTOPOcZxWKxMAAAAFel07oDAAAAWsJACAAALg0DIQAAuDQMhAAA4NIwEAIAgEvDQAgAAC4NAyEAALg0DIQAAODSMBACAIBLc9O6AwAAAJXJyMj46KOPCgoKRo4cOXLkyFsr3LhxY+HChWlpae3atXvqqae8vLwYY4cOHVq/fv3p06cbN248bdq0e+65R638+uuvFxcXq7ejoqLGjx9fTQOh2Wye8+rrL86eXT13V81MJpObW+3/SoGHWZvgYdYm1fYw63u637bOiRNXS0pMQs2eOXPq4sWjjz76KPevly9f7tGjx2OPPda1a9cnn3yytLR03LhxN9UZN26c2WyeOnXqkiVLDh48uGrVKsbYCy+80KFDh0GDBh06dCgmJmbv3r2RkZGMsQ8++GDq1Kn+/v6MMW9vb8aYUj25RktKSp54fd7QKY9Uw30BAIAjjGpluG2d6Ojlhw9fFmz42IgR15OSkrh/mzdv3o4dOzZs2MAYW7ly5YIFCw4cOGBbIT09PSYm5uLFiz4+PlevXm3SpMnx48dDQ0PNZrNO9/fyX3x8fMeOHV955RXGmL+//++//96yZUtrC1gjBACAmmvPnj19+/ZVb/fp0+fQoUOlpaW2Ff7444/27dv7+Pgwxho2bBgZGbl3717GmHUUZIxdvHixUaNG1l8/+OCDWbNmJSUlqVNBDIQAACCNotjzU4m8vDw/Pz/1tr+/v8ViycvLoyqodXJzc20rLFmy5MqVK5MmTVJ/HT58eFhYmI+Pz/PPPz958mSGYBkAANDcjh07OnbsaFsyZsyYWbNmMcbq1KlTXl6uFpaVlTHG6tata1vTy8vLaDRafy0tLbWtsH79+ldeeWXLli3qciBj7LPPPlNvTJo0qVmzZnPmzMFACAAA0iiKolQ+xbuFxaJER0fPnz/ftjAwMFC9ERQUdO7cOfX2uXPn6tSpo8a52NY8e/as9dfz588HBQWptzds2DB9+vSff/65TZs2t95vcHBww4YNc3JycGoUAACkUoR/GjRoEPO/AgIC1MZGjx69bt069YKHFStWjBgxQl38S05OzsjIYIwNGDDgwoUL+/fvZ4z9+uuv169fV9cUk5OTp06d+v3338fExFi7VlhYaJ1fbty4saioKDIyEjNCAACouYYOHfrFF1906NAhKCjoxIkTW7duVcvnzJkzbty48PDwevXqzZs3b+DAgV26dPnjjz8WLFigXkc4ffr0kpKS8ePHq/UTEhLeeuut33///cEHH2zbtm1JSUl6evrixYvvvvtuDIQAACDP7YJfROn1+vXr1x86dOj69eudO3dWBznG2Pr1663LftOnTx86dOjx48cjIyMNhr+v8di+fbvJ9P+vaFTDSuPi4g4ePHj69GlPT8+IiIj69eszBMsAAIBEdqwRMsvt60dHR99UYh3wVEFBQdalQVVISAi3qVtrYo0QAABcGgZCAABwaTg1CgAA0tz2AvlbVeHMqGNhIAQAAGnsWCOUG1xjB5waBQAAl4YZIQAAyKNeI+9UMBACAIA0ivipTs1PjWIgBAAAaey5jlDrkRBrhAAA4NIwIwQAAHnsWCPEqVEAAKhNnG0cFBwIDx06dPLkybvvvrtr16516tRRCw8fPpyRkdG6devIyEgH9BAAAMCBqrpGWFFRMXHixGHDhn377bevvvrqypUr1fK33npryJAhP//8c//+/T/88EOH9RMAAJyAovwdLyNA6zlhVWeEixYtOnr0aFpamrqThcViYYzl5+fPnTs3NTU1IiIiNTW1f//+U6dOveuuuxzYXwAAqMmccI2wqjPClStXPv300/n5+Xv27Llx44YaHZucnNyyZcuIiAjGWMeOHf38/Hbu3OnAzgIAQM1mx4RQ+HIL2ao6Izx16tSqVasWL17s4+Nz/Pjxn376KTo6Ojs7Ozg42FonODj4woULVAvqJBIAAJyUxWLRfNByhKrOCIuLi318fPbs2bNly5bJkyfPnDmTMWY0GvV6vbWOm5ub0Wjk/rvFYsFACADg1Mxm823rqLtPCP1orqoDYZMmTWJjY9XvAvfdd9/hw4cZYwaD4fLly9Y6ly5dumnLYCtFUXQ6XLwPAODEbGc+JMWuH01VdXDq27dvVlaWejsrK6tJkyaMsZ49ex44cODatWuMsdzc3IyMjG7dujmoowAAAI5Q1TXC5557rk+fPg0aNKhXr94bb7zx73//mzHWvHnz4cOHjxgxYsKECV988cXEiRPVARIAAFyTer5T8F80VtUZYevWrXft2lVcXJyTk5OUlPTAAw+o5StWrBg3blxaWtrDDz+8ePFih/UTAACcgfgaoebLhAKZZe655565c+feVOju7j59+nSpXQIAAKdVA9b8RCGABQAAXBqSbgMAgDTqBfWC/+MkF9QDAADclh1rflqPgzg1CgAArg0zQgAAkMcJg2UwEAIAgDR2JNHWPH8pBkIAAJDGCSeEWCMEAADXhhkhgDgpO6k43ddmgKpwwo15MRACAIA8dqwROqgnVYZTowAA4NIwIwQAAGkU8QvktQ4axUAIAADyKIr45RBaj4QYCMF5ECEqgpErElqREysj2gr5WSHwISKjDYBKOWGwDNYIAQDApWFGCAAA0jjhhBADIQAAyGNPijWth0KcGgUAAJeGGSEAAMjjhOdGMRCClixCkZNEbW4p2TLViFhPhIqJt7msCHNeuXg4ukjoqdYfW1CT2XEdoeYDIU6NAgCAS8OMEAAA5HHCXKMYCAEAQB6sEQIAgCtTFOfLNYo1QgAAcGmYEYJMdOylQKwm1YjFTDTCKyfq0o0IRZOKxK8yMmpUqJhMZKzoOOW8MrIyY0zhfSWWEnqq+Zd9qGYKU0QvkNf8gnoMhAAAII8da4Raw6lRAABwaZgRAgCANHYEy2g+g8RACAAA0ihM/DpCrVeSMRDCbQjFv4jGuZi5cS4VRGWivIJXTjZiNhONc4v5QTRikTVEwAgZ/EKU6/T8xnU6zgKHTs9vRE+Uc+tTjehEIm5oYrFC4DSwRggAAOBcMCMEAAB5FOFTnZqfBsBACAAA0jjhmVGcGgUAANeGGSEAAMiDpNvgvMiAT5EISTrgkx+rWWHilJuM/EZMvMqMsQojrxGqMlFOR5kKBMcKhY1S6yJUTCYZCOrGOa/jxitkjOnd+eXc+m7uAvfIGNPpBeJXFUVKQjqocRR7tmHC5RMAAFBb2LFDvebfcjAQAgBAjVZYWLh69erCwsLBgwdHRETcWsFisfz0009//vlnVFTUoEGD1MKioqLk5OTMzMzGjRuPGDHCx8fHWn/Lli379+8PDw8fMWKEogheAQsAAFAZxa4fWnFxcZcuXTZs2HDp0qXu3bvv3Lnz1jpPP/30Sy+9VFxc/Nxzz82aNUst7Nev37Jlyy5duvT1119HRETk5OSo5a+//vr06dOLi4tff/31adOmMcwIAQBAIjvWCCs/N/r111/Xr1//u+++UxTFYDDMnTu3T58+thVycnKWLl166tSpwMDAyZMnt2nT5h//+Iefn9/GjRv9/f3VOr169fryyy9nzZpVWFg4f/783377rW3btjNmzAgJCZk9ezZmhAAAUHMlJycPGTJEHVyHDh26bds2k8lkW2H79u2RkZGBgYGMsbCwsLCwsJSUFMaYdRRkjHl5een1esbYnj17fH1927Ztq1bo3Lnzli1bMCN0OaLRoUKZP6lYTVM5P5WnsZxT3yhSmapv4oWSMiLElBEJSxnxMB2aa1Q0TSg3ENSNiA519+BnLHX34NQXqswYc+PVp+JXyUSmREpV7PrrLOzYfaLy6tnZ2f3791dvGwwGs9mcl5cXFBRkrZCTk2MwGKy/BgQEWM+CqrZv375///4vvvji1soGgyEnJwcDIQAAaOzQoUOPPvqobUmfPn3Gjx/PGNPpdNavnuqNm75EKopi+93UYrHYVjh06NC4ceNWrFihThm5lTEQAgCAPHZcR6iwBg0axMTE2Ba2aNFCvWEwGPLy8tTbubm5er2+cePGtjUNBkNubq7117y8vCZNmqi3jx07Nnjw4I8//njw4MHcyrm5ubGxsRIGwvLy8uzs7KCgIHd39ztvDQAAXE1oaOgjjzzC/VNcXNzixYvnzJmj0+l+/PHH2NhYNzc3xtiZM2fq1avXsGHDfv36TZs27fz588HBwadOnTpz5kzv3r0ZYxkZGYMGDXr33XdHjx5tba1bt24FBQWHDx9u167dpUuX9u3bt3z58qoGy8ycOVOxUVxcrJb/8MMPgYGBgwYNCg4O3rJlyx09EwAA4OTUNUKhn8oXCR944IGSkpLhw4c///zzc+fOnTNnjlo+duzYL7/8kjFmMBgee+yxuLi42bNnDx48+Nlnn23YsCFjbMyYMRaL5fvvv09ISEhISFi8eDFjzMfHZ9asWSNGjJg9e3b//v0nTZoUGhoqMCOcPXv23LlzbUvKysoeeuih5cuXDx48eM2aNYmJiVlZWWpkDlQrKv6Fu3cuP1iEDIoh4194USflZfw4l/JSE1HOqU83QgXRcMqNosEyVOo1boo14jkUShpGXcRLpVijEptxg2XcBYNlPOpwyj08BSozxjzqcJ4AqhEqnMfNQgXXcEu5demUXQiuqR5Sn+e6dev+/vvv3333XUFBwR9//NGyZUu1fP78+eqyH2NswYIFcXFxaWlpixYtio2NVQvffvvtGzduWNtp1qyZeuPll1/u0aPH/v373377bfXqe7FToyUlJV5eXtZfN2/eXK9ePfXc65gxY5566qmUlJR7773XzocLAABwi7vuumvSpEk3Ffbq1cv217i4uLi4ONsS67rgrfr06WN7MaLAdYTz58+/++67AwIC3nvvPbUkKysrPDz874Z0uhYtWpw5c6bqDQIAQC2j2EXbPld1IJw+ffq1a9f++uuv9evXv/HGG+vWrWOMFRUV2U4Q69atW1hYyP13i8VSUcE/rwUAAE7BaDTeto4da4Raj4NVHghbtGihjnldu3adOHHiL7/8whhr1KjRtWvXrHWuXbt2U1SrlaIoWDsEAHBqVbk0wI5xUPPFW3tSrF27ds3b25sx1q5du4MHD5aXlzPGbty4kZaW1q5dO8kdBAAAcKSqBsu88sorvXr18vHx2bFjx9q1a3/99VfGWOfOncPDw5999tlHH3104cKF3bp1426QAY7GjQ5lRHBjBbFHLpWTzFjGLy/jBYKWlfCjQ6ny0hJe1KhIiCljzMiLMqWiRqmHSW0dTKRY49YVTbHGr0znHiPCLIWiRkUCQT3q8D8f6pTzG+E+txUmfiOeROMWT/6T6Mb71q6nvsrr+I1ovgGsS9B+giesqjNCi8Xy9ttvP/PMM0eOHNmxY4c1BcD69etLSkqmTZvm7u6+evVqh/UTAACcgROuEVZ1Rvj6669zy5s0afL555/L6w8AAEC1Qq5RAACQRqF3ViH/xVlmhAAAALfnhGuEGAidCRmjQWT84sbFUNEiVCgKFedSUswpL73Bv8yolFeZMVbKa7yMF0HD6CAa7j6FJiO/EXKfQhn7EQrlWJO1HyE3WMbNXWwrQW5cjKcX8RyW8z83+MEyxBNr4eWuY4xZLBI+lOggGk6Z5tORWsae/Qi1fgmwQz0AALg0zAgBAEAae1KmaT0lxIwQAABcGgZCAABwaTg1CgAA0tgTLOOYnlQdBsIaihuWSAUrknvqcvfOJaJDqcDOEiIQtOQvTjlZ+YZA1Gg5EadKbdjLjRrl7tbL6A14K0zEc8vdmNehUaPkxrwCG/ZSG/BSUaP8NHVUdCjvCWfEc0sdnHTUKLdYDJVKjXhqiVLNP56dk8KE1widZhsmAACAWgkzQgAAkMeOC+q1nnxjIAQAAHlqQBJtURgIAQBAGjuuI8QaIQAAgJYwI9QYvcUr5w9kdCgRCcndU5fMHUoEfBYXEeXcqFFeIWOspJhfzu0MmWuUihrllVM5RannyiwjalRwY16xqFEdETXqxosadXPnP1fUxrz8fK1khK3AJsZmIhGuaHSo0HNIlbspnOdKx39KGKJJXQcGQgAAkMauU6MO6ktV4dQoAAC4NMwIAQBAGjsyy+DyCQAAqF20HthEYSDUHJWAilMolEqNMVbG28mWu6Euo+NcuEExjLHionJOI1SKNeJOy3jlZUQSOCnBMi6SYo27Wy9jzJ2KIeLuqUs9VyJZ08jnikBNJhTe08ItZPTmxtznlmpEUUReTvgvRRG+HAKXTwAAAGgJM0IAAJDHjt0ntJ5jYyAEAAB57Mg1qjWcGgUAAJeGGSEAAEijMIXaD7Kyf9IUBsJqQqdS45dzgxWplFdUOCU3gVkpuXeuwAa8VH1qA16qcX6KNSJqlBsdSpUbXTtq1J2IGqUfPndPXSKVGrGnLvnweeiHz6/PDQTVE9Gh3OeEaoSOGhX4dNZ8lavmsGeHeq2fPQyEAAAglZNNCLFGCAAArg0zQgAAkMaepNsO6kqVYSAEAABpnDHXKE6NAgCAS8OMsNoIpGdkRMCeqZxIwsnLKcqoqFEq1ygV8EnsqctNH0pFh1J3WoOiRon8mdz8ro6NGqUiIUUiJCuoqFHiYXLLqYNTLHuoSO5QVsnDF0mpSpbzIm+pJ1ZHBLAq/I18tZ7U1Bx2XFCv9ZOHgRAAAOSxZ40Q1xECAEBt4YQTQqwRAgCAa8OMEAAApLHj8gnNU8tgIJSPG0oglEqNETEdxnIixRoRXVJawikv5cWnVFLOjWdh1J66go1w75R6OFS5kRdDJLwxby0LljGJBcuIZZIjcB8QmTFOMDuamwcvWMaDH7ji7sE/VNx59d3c+YcE9UJQYT78ypqf9dOEs50bxalRAABwaZgRAgCANPYk3XZMT6oOAyEAAEijKGIbd7AacAIZAyEAAMjjhNdPYI0QAABcGmaE1YSMGiWi+ExGTjk3PJLRG/NyU6+V8UJJGWPlZMAnsesvL4aTyo5GlXMDQemoUX4PucG0VNQotbmxWcbGvNSrzD3zI7oxr47YmNfNjXOv1HFFZk3jlhKonnPDKXV6/rdtMjrUjTgOeQGfHp78yuV1+NGkHrx3EDeUlDGm5z2xjDEdr7rmJ/dqDiecEGIgBAAAiZzwOkKcGgUAAJeGGSEAAEil+blOQWIzQqPRGB8fP378eGtJenp6v379DAbDwIEDMzMzZXcPAACciXodoeiPtsRmhG+++ea5c+es23RZLJZRo0Y9+OCD69atmz9//tixY/ft2+eATtZUVIABL2SC3neQCJaRkWKNCEUh4guoiBuRcqEtA2X1UChYhkyxJiNYRk6KNSonWQW/nHsIUan7hLKmUZ9OQlnTyKAYYstAd14qNcaYO++oKC/lx7mQ2fi8eIcKcUi4V/B7wn8vkx/l1H6MVH2nZ0euUeE1RdkEZoRHjx794YcfnnvuOWtJSkpKfn7+rFmzfH19//nPf544ceLAgQMO6CQAAICjVHUgNJlMDz/88Icffuju7m4tTE9Pb9eunV6vZ4x5enpGRET8+eefDukmAAA4BUX8R2tVHQjnz5/frVu37t272xbm5+f7+PhYf23QoMGVK1e4/26xWEwm/vkuAABwCkaj8bZ17Fkg1HosrNJAeObMmY8++mjq1KmZmZmXLl0qKyvLzMy0WCy+vr5//fWXtdr169f9/Py4LSiK4uaGCFUAACdme0awNqnS4HTlypWAgIDExETG2NWrVy9evJiQkLBr164WLVqkp6dbLBZFUUwm08mTJ5s3b+7gDgMAQM2lTvME/0VjVRoIO3bsmJqaqt7++uuv3333XfXXfv36KYryxRdfJCYmfvLJJ3fffXe3bt0c2Nkahgwa5RZSG/ASUaMVvKBHKsUaWc4N7CRCT8lykUBQ0ahRbs+FeyhlY16hqFHi1RRC7e9KxWSazVS5lD11OY0rOv6rRu6pKxI1aiSiRo2e/EBQ7qsvdFxR5dz3GqPfm27cJ5x4Tiiaf/Q7kANyrJnN5tTU1IKCgm7dutmux9nKzMzMyMi45557QkNDbcsvXLigKEpgYKC15OzZsxUVfx8J3t7ejRs3Fs4sc9ddd1lbdHNzW7169bvvvuvt7b1s2bJvvvlG8yhYAACoTUwmU3x8fGJi4sKFC1u1anX8+PFb63z00UfdunX79NNPO3fu/Nlnn6mFq1at8vPzCwsLGzlypG3lmJiY4cOHJyQkJCQkLFy4kNmRWWbYsGHDhg2z/tqtWzfr2VHRpgAAoNYRv46w0inhjz/+eOrUqUOHDnl5ec2aNevVV19dtWqVbYXCwsKXXnopJSWlffv2u3fvHjp06Pjx4728vLp27bp3797du3d/+OGHN7X53XfftWzZ0vqrnFyjGAUBAIA5IGo0KSlp9OjRXl5ejLEJEyasX7/ebP6fs9mbNm0KCQlp3749Y6x79+6+vr7bt29njIWGhlJhK1lZWceOHSstLVV/RdJtAACQx47rCCsdCM+fPx8SEqLeDgkJKSsru3Tpkm2FCxcuNG3a1PprSEjI+fPnK2nQw8Nj1qxZCQkJgYGB69atY0i6DQAAmjtz5synn35qW9K2bVs1+rKsrMx62YaHhwdjzDqTU5WVldlenufh4VFWVlbJfZ04cUKNuFm1atXkyZP79u2LgVA+brweFWYolGuUioQ0ipSbjGKBduQOt9y4VsEechsRukdZjVAvhFDUqEiqUWlRo0LxoUIJTnVEJKTRjf8c6nnPrZvgC0Efh7wAYynHG7VXM3VI8Ip1Insy12725BplSkFBwf79+20LGzZsqN4wGAzWVC2XL19WFMVgMNjWNBgM+fn51l8vX758U4WbWONOH3jggRkzZhw9ehQDIQAASGPf1RPR0dFLlizh/rVbt25bt2598cUXGWM7duzo0KGDp6enbYWuXbtOnz69sLCwXr16+fn5f/75Z5cuXapyvxcvXiwoKGjcuDEGQgAAqLmmTJkyb968l19+uVWrVrNmzXr//ffV8h49eowdO/app55q1arVwIEDx44dO3ny5KVLl44ePVpdU8zKylqzZs3Bgwdzc3Pfeeed5s2bjxkzZteuXd98803Hjh1LSkoWLVo0ZMiQiIgIBMsAAIA8soNl/P39f//9d6PRmJKSsmzZsoSEBLV82rRpXbt2VW9//fXX/fr127Rp09ChQz///HPbf2/fvv2TTz5p/TU8PLxp06a//fZbenr6rFmzECwDAACSKXZcUHe76s2bN3/33XdvKpw6dar1tpeX1z/+8Y+bKjRr1kw9oWqrcePGtxZiILwT1O6sIhvzmgUygZG5oGREkVC5x6hQAn4PRSpT5UL3SDYimmKNiIywcLe9FUxgxqUjPizMVAIzi8CHC/VBREXimLjZ0Uz8ykIvhJTjiokezCLl1D1S703+e1n4kKi1UTR27Div+XOBU6MAAODSMCMEAAB5HJB029EwEAIAgDx2XEeo9eWWODUKAAAuDTNCAACQxo5gGZwadWYCQaPkRqlmfj4pfnonchdfMkcUL0aObIToIVHOrU81InSnQvdINkLdI5U3S2hjXuLVpCIHuZ8LVBAolayLoigCL4SOaF2vFzneRF5NKccVdafcI5wJxgaTqdSI9yb31SeDRoXS7oFGMBACAIA0iuz9CKsB1ggBAMClYUYIAADS2HNBvdYnijEQAgCANPZsw6T1SIhTowAA4NIwI6wmVFAZGU3KC1YU2jyWqi/ciEi5UE+EG5HRE27uUNHGJW3MS9SmGqECjHmBoHKeQzKcsrpfTepO5bwjqFdTMDYY/mZHZhmtYSAEAAB5cB0hAAC4MoUpopdDaD0OYo0QAABcG2aEAAAgD3afqJ2oOBeh+uTCu8CCPFmZn09KsBEZ5do0IhLPQu2pK9a4aBgF931OvWrEhwLxgJhOKJyn5r+aYo1z61byjhBphP4Dr5Coyy9mvLx46h+IcuehiF8XqPmDxqlRAABwaZgRAgCANIoifoG81hfUYyAEAAB5sEYIAACuzJ5co47pSdVhjRAAAFwaZoTy8YNGqcoifyCj2IQaF9xBVGzHUdHtSflxrYKN3Gndyu5UJARYsHHia7Bw4zJqSzlUhA4J4VdZ5FARfEeIxQALPVUuCEm3AQAAnAwGQgAAcGk4NQoAADIh6TYAALgue9YIHdSVKsNAKB9/8zmqssgfqMpijZO1RRqh6gs2zi0Xu0cJdSu7U36xYJwLt3HhQ4KqL6O2lENF6JAQfpVFDhXBdwTxIos0rvlHeQ3ihNcRYo0QAABcGmaEAAAgjTNePoGBEAAApLEjs4zmcGoUAABcGgZCAABwaTg1WgVU8JjQ1ppkZJpAxBpZmfg+I9aIjHJtGtHxGuEVMsZ0RCPUdrj855bY95U8VPgtiz1MqudCD98JXk2xxrl1K3lHiDRC/4FXSNTlF2sfJ+k4zrhGiBkhAAC4NMwIAQBAJmfblxcDIQAASOSEF9RXdSDcvHnz+vXr8/Ly/P39J02a1LNnT7X8+vXrb7/99okTJ9q2bfvCCy94e3s7rKsAAFDT2bFGqLmqrhFmZmZGRUVNnTo1PDx84MCBu3btUsvvv//+kydPPvLII/v37588ebLD+gkAAOAQVZ0RTp8+Xb0RHx+/b9++5OTkXr16HT16dPfu3ZcuXapbt26XLl0MBkNmZmZYWJjDeuvERCPTdLyoP51eoDJVX7gRkXKhnpA9dGRPzFRPhHb9lZJrVOSlZ4wpIi+cnOdQ5FVzaE+oO5XzjhCOsOUWw9/suKBe86dULGq0oqLiyJEje/fu7d27N2Ns3759MTExdevWZYz5+vq2bt06NTXVId0EAACnoNj1oymBgXDJkiUNGjRo1+33iT8AAB/KSURBVK7d8OHDY2NjGWN5eXkNGza0VvD398/Ly+P+r8ViqaiouMO+AgCAhkwmk9ZdcAiBgfDRRx8tKirKzMzcsWPHggULGGPe3t5lZWXWCiUlJVSwjKIoOh2uWQQAcGJubrdfTVMYU/4+PyrwUw2dr4Tw4NSsWbMHH3xwy5YtjLGgoKAzZ85Y/3T27Nng4GDqH50ujggAAEQpzNmGwaoHy5w/f14d5MrKyjZv3ty2bVvGWFxcXGJi4p49e7p27bply5bS0tI+ffo4sLM1jVACJioyQs9vhLuqrydCA/Ru/C80Oj2nnG6E6CFRzq1PNULeKTd4gWqkgoh/MXPKuYWMMT0V6EJQeGf0zSKRNRQyZZrIc8WoQ0XkVWPEcy70qlGNSzmuqDvlHuGMfkcQjYi9N4XytGn/GV/97LuOUMZ7ym5VHQj79u1bp04dPz+/48ePR0VF/d///R9j7K677vr3v/89dOjQdu3aHT58+OOPP/b09HRkbwEAACSr6kCYkZGRkZFRUFAQHBwcFBRkLU9MTBw6dOjp06fDw8P9/Pwc00kAAHASdlxQryjOMSPU6/URERHcPzVq1KhRo0byugQAAM7KnusIHdOTqkMkJwAAuDQk3QYAAGns2o/QQX2pKgyEd0JgL05yn1ji8kpu2JvenahMlLvxyrmF1D0yxtyoADxuD0UqU+VubvzlAnMFUc4L4qSynVEUYp9lMy/3mo4IGxXbqllGjjFGvKCiLwT3VZbyako5rpjowSxSTgZdE+9N/guHsFFnhoEQAACksWONUPNvC1gjBAAAl4YZIQAASKMw8TVCraeEGAgBAEAeJ9yhHqdGAQDApWFGKB/3rAARCUiHAnLDKYlYOHeRcjd3fgpFdw9+uZs7f/8sbmeonlRQ5SZOuVB0KGPMIhIhSp2zqTCJhKTKSDYqGjVKpw8VOFSEyiUdb2I9oY9DTrlQT6g7peJXyRykQkGjrscZN+bFQAgAAPLYkWJN63OjGAgBAEAerBECAABIV1RUlJeXV0kFo9F44cIFk8lUldZMJtOFCxfKy8vVXzEQAgCANIpdP5V78cUXg4ODO3Xq1L1798uXL99aYePGjUFBQX379g0ODt62bZtauGXLlnvvvbd+/fqdO3e2rfzbb7+Fhob27ds3MDBw/fr1DKdG74TQbJ6KjCD3OOWGohBxBGS5Jy++wIMKUiDKeY0wxtyN5lsLK0ycQsZYBRH/wi2nQlEEg2IEE5jxUqkxMn8b0RORHGuiPZQSLCP06gsdV1S5UGXhHpKNCJST+diI9yb3vYwEa1YKfWyT/1Jp9ZSUlC+//PL48eMBAQETJkx49dVXP/74Y9sK5eXlU6ZMWbp06bBhw1avXj116tTMzEy9Xt+gQYOnnnrq5MmT3377rbWyxWJJTEx87bXXHnrooa1bt44ZM2bAgAGYEQIAQM311VdfJSQkBAQEMMZmzJjx1Vdf3fRNNDk5uW7dusOGDWOM3X///eXl5SkpKYyxjh07jhw50mAw2Fbet2/fxYsXJ02axBiLjY0NDAz86aefMBACAIA8ss+NZmZmhoeHq7fDw8OvX79+9epV2wpZWVnWCjqdrkWLFllZWVRrWVlZYWFh7u7u1gazsrJwahQAAKSxJ+k2Y5cuXdqyZYttScuWLZs2bcoYKyoqqlu3rlro7e3NGLt+/bqfn5+1ZmFhoZeXl/VXb2/vwsJC6o64lTEQAgCANHbtR6icPHnynXfesS0cMmTIM888wxhr1KhRQUGBWnjt2jXGWOPGjW1r2lZQ6zRq1Ii6r1srd+zYEQMhAABorEePHklJSdw/RUVF7d27V729d+/e5s2bq/NC2wqHDh0qLy/38PAoKSk5duxYVFQUdUdt2rQ5depUQUFBgwYNzGZzamrqzJkzMRDeAfJLDy+ojFiNFUqxRgXUedThx8hxyz3q8F90Yxk/lZqxnB8IapIRNcqPyeRWpXG/flKxlyYqTFcvIWqUCiblfj8WjmslokbFsvFRMZm8o4I6roQCQT2IysLlIj0ke857+KIp1sQ25nXRsFHxf6E99NBD7dq1++abb1q1ajVnzpwnnnhCLZ88eXJcXNz48eM7d+4cHh4+c+bM6dOnv//++506dWrTpg1j7OLFiykpKXv37r169eratWsNBkPPnj2bN28eGxv71FNPvfjii8uXL/f39+/Tpw+CZQAAQCJ1lVDsp5LmmjVrlpSU9MUXXzz++OPjx49/+umn1fKgoKAGDRqot5OSkgoKCiZNmmQ0GlevXq0WXrx4ce3atRcvXuzQocPatWt//fVXtXzFihV16tSZPHnyhQsXfv75Z0VRMCMEAIAaLTY2NjY29qbCN954w3o7MDDwyy+/vKlCVFTUmjVrbm3Nz8/v008/tS3BQAgAANIoNWA3CVEYCAEAQB4nTLqNgbCaUF+RyGAZd045lThKKL7A04sKiuEfDCYq/oVXbq7gV5aTNY0q13EeEfXE6k38cjKch5sEzqEp1shwHikp1oTiqviVPUXKycpeRE+8+Mchtz4V+UW9I7gPn/teY5UEy2j9qV3D2Xf5hIM6U0UIlgEAAJeGGSEAAEhjzw71julJ1WEgBAAAeZxwjRCnRgEAwKVhRggAANLYESyj+YwQA6F8QscAvQurhBRrdco55SYqOpRIpSaUNY2bkIzRuce4hHes5UX3Gd0kPBzmtFGj7lTUqEh2NDrgkwrs5EUp1xWozBirI1Jehwo9FUmxxn0CGX28EQnzuHVdkROeGcWpUQAAcG2YEQIAgFSaT/EEYSAEAAB5nPCCegyEAAAgjT3XEWo9g8QaIQAAuDTMCKsNtcknv7ZOz0sgSSaK5EcrcvfO5RYyOpySGzbJqPShgpvq8gPwRKJDqXK94MOsMEnYmFdO1CgZSFzTo0a9vN05hUTUqFddTmXGmJc3ETXKa4fqCZWDlPsO4r7XGH0cOt8KWDWzI2xUaxgIAQBAGkURX/PT+twoTo0CAIBLw4wQAACksWNjXq0nhBgIAQBAHoUpiuAiodbjIAbC6iL6lYebZ8uNyAVFbUNaYeK8vkK5xBhjZn5widieutSRzg0MoYIXqFRY/J1pXTtYhtyYl9rDWShYhox/4ZRzI2gYY153UcEy/PI6vHIyWIZ4mNx3kOgGvJpPX2o6J8yxhjVCAABwaZgRAgCAPOIX1GsOAyEAAEhjxzZMiqIIXoEsWVUHwqtXryYnJ587dy4kJGT48OF16tRRyy0Wy48//pient6uXbuBAwc6rJ8AAAAOUdU1wpiYmG+++SY/P/+TTz6Jjo4uKChQy2fMmDF79uzi4uJnnnnm5Zdfdlg/AQDAGSh2/WiqqjPCffv2+fv7M8YqKiratGnz3XffJSYmZmdnf/bZZ6dOnQoMDJw4cWJUVNTzzz/v5+fnyA7XPlTEGudUARXeRoUIevISTfFTo1VSLmNPXaGsaVR0KPUwjbxyKpOcicokJyNqlHquiK1cxaJGdUTUKDcSUk7UKLHtLZlijZc1jYoOrUtFjVLl3BRrVCo14uFz30H0aTytP56dkz1Jt4WTM0pW1RmhOgoyxvR6vYeHh16vZ4xt3769bdu2gYGBjLEWLVqEhobu2rXLQR0FAICaT/l7KBT40fw7h/DlE2vWrLly5crw4cMZYzk5OQaDwfqngICAnJwc7n9ZLBYzdT0aAAA4g9r6MS4WNbpz584nn3zy+++/b9CgAVNDfWxOBlksFvLMmNOF0wIAgH2c7YJ6gYFw9+7dCQkJq1at6tatm1piMBhyc3OtFfLy8po0aUL9u06Hi/cBAJxYVT7Ga/PGvAcOHBg1atTnn3/er18/a2G/fv2OHTt24cIFxlhGRsa5c+d69+7tkG4CAAA4RlVnhIMHD/b29l6+fPny5csZYyNHjhw3blyTJk0effTRAQMGjBw5cs2aNc8//7yvr68je1sLCUWs6fhBfMzNwv9CY/HkRjzyX3Sh6FBG9Jz6vkjFu/LzZLqZuJXdPYioUV7Eo7G8glvZRXKNuhN7OFPPITdq1IOIDq1DbszrwFyj3FBVd08iwJgIPCaiRsVyjULl7Lmg3kFdqbKqDoSffPJJRcX//2SJiIhQbyxcuHDTpk1paWlLliyxnSwCAIArcsKk21UdCEeNGkX9KS4uLi4uTlJ/AADAmdmRa1TrgRABLAAA4NKQdBsAAKRRxK+XE93IVzoMhDUUtX7PLSWDaGTM+MntSfl76goEdDDG3HixG2VUoEcpP4jGWM6JfzEZ+cEyVOo1of2KHRssQz2HIgn23NwFg2V4ucqoFGtUsEwdXhY07oa6jEiZxirZa5e3S7BQKjVGPOcIigGcGgUAAJeGGSEAAEhj336EDupMFWEgBAAAaezILKP1EiEGQgAAkEj7zSSEYY0QAABcGmaEzoQ84UB8n9HLiRoV2GtXKLKRMebGCxD18ORHh5aX8oMYjWWcAFEjtTEvUW6uoMq5UaPcunSaOm46OuLVpCIedXqB/YrdRTfm5cVkckNJGWN1RDbspaJAqT11yaxpvEekJ54ThTjwtV6Kcgl2rRE6qC9VhYEQAADkqcW7TwAAANRKGAgBAMCl4dQoAABIozA7tmHCdYRwx8jDSMcJ3qAiaKhGhDKBkanUqNgND06cSzkvcIMxVl7Kz5rG3XqQCpapoMqJfQr5+xHy6wplWCMDOuj9CPn/oBcKliHS1/GDZUQiaxgRXEM1QsZPCW0lSAXFaP3B6sqc8TpCnBoFAACXhhkhAADIU4s35gUAALgtZ8w1ilOjAADg0jAjBAAAaewIltH6zCgGwtqB2juXn9qLX5nK+OVGROZxgxup9GBubvxybhCjBy8KlDFm9OIHa3KjRskNeB25MS+VYo37uSBrY15u1CgdpiuwYa9QZUYkzBOKAq2kHHvqOhNnWyPEqVEAAHBpmBECAIBMTncdJwZCAACQxp41Qq3HTQyEAAAgjT2XT2g9g8QaIQAAuDTMCF0O/V2NiuIjavOiRrmFrLINezkxnFSwoolIB8oNBCUrUzlFqahRbq5RKqmoSNgo9UJQuUbJaFJeWCYVq8kNMaXqu7mLZT3lbh1MR4Fyi+mNoLU+ewZVhcwyAADgypB0GwAAQLIVK1aEh4c3btx4+vTpZWVlt1ZIS0vr3bt3w4YN+/Xrl5GRYS1/6623QkJCgoODX331Vevlv7GxsR3/66WXXmKYEQIAgETS9yNMS0t78sknN2zYcM8994waNerNN9987bXXbCtYLJbRo0dPnTp148aN8+fPHzt27MGDBxlj69evX7Ro0datW93d3e+7777w8PDx48czxg4fPrxs2bLg4GDGmK+vL8OMEAAAarLPP/989OjRPXr08PPzmzNnzrJly26qkJKScvXq1ZkzZ3p7e7/00ktZWVn79u1jjC1btuyJJ55o2bJlaGjo008/bfuPkZGRMTExMTExYWFhDDNCsBINolEUTmAImTZMJ5BnS+/GDzlxr+A3ws2ORge/UMEy3GJ+NjWHplijU6/xG+c+t6J52rj1ySxoRDgPP35K+Lii6oNzkH4dYXp6elxcnHo7Ojo6JyensLCwXr16thXatm2r1+sZYx4eHq1btz5+/HinTp3S09OffPJJ6z+++eab1n8ZPXq0Tqfr0qXLq6++ajAYMBACAIDGysvLr127Zlvi7e3t4eHBGMvPz7cOe+qNK1eu2A6EV69e9fHxsf5av379K1euqOW2/6gWMsYWLlzYvn370tLSt956Ky4ubv/+/RgIAQBAHvEL6pnCkpOTmzdvbls2bdq0efPmMcYaNmxYVFSkFhYWFjLG/Pz8bGs2bNjwr7/+sv56/fp1tcJN/2j9rwcffFC9sXLlSj8/vyNHjmAgBAAAeey4jpCxIUOGJCUlcf/UsmXLtLQ09XZaWlqjRo3q169vW6FFixZ//vmnxWJRFMVkMp04caJly5ZqeVpamnpaNS0trUWLFje17OHh4e7uXlZWhmAZAACouaZMmbJ27dq0tLTi4uJ33nlnypQpavm//vWvzZs3M8buvfdeT0/PJUuWVFRUfPjhh40bN+7atStjbOrUqYsWLcrLy7ty5cqHH344depUxtjp06d3795dWlp6/fr1F154oV69etHR0RgIAQBAGjVYRvSnEh06dJg7d27//v0NBoOvr+8///lPtfzYsWN5eXmMMb1ev3bt2sWLF9erV++rr75atWqVem527Nixo0ePjoyMDA8PHzBggDqCXr9+/ZFHHmnQoEHTpk2PHTv2888/161bV6Hi3+QqKSl54vV5Q6c8Ug33BRqijyaBMEuqEQsv2xlVTtSlGxEKBKUq8+9TLD5SOPcYd4dkqhEZAZ9CgaCIAq1NRrUy3LbOih2Zl6+XCjV77I+t19K3UadGqwHWCAEAQCpn+/aDU6MAAODSMCMEAABp7Aoa1RgGQgAAkMaejXm1XknGqVEAAHBpmBGCTKKJJcUaIcot3DyZgnvnikVPUyGpRHUiUlPkHit7WkQCPoVbv8O64HqwMS8AALg0RfhUp+anRjEQAgCANPbsPuGYnlSd2BphXl6eeiW/rdzc3B07dly8eFFerwAAAKpJVQfCTz/91GAwNGnSZNq0aTeVt23b9u23346MjPz6668d0EMAAHAqiviPpqp6arRnz57btm376aefdu7caS0sKiqaOXPmli1bOnfuvG3bNjWxm6enp2O6CrWQ6FKCw1qmA13utC5J+L0vJRBJQhsAlVGzh2rdCzFVnRG2bt06IiLipr2wN23aFBwc3LlzZ8ZYv3796tatu2PHDuldBAAAcJw7CpY5f/58aGio9demTZueO3eOqlw92b0BAEBD9gTLaD2BvKOBsLS01MPDw/qrp6dnSUkJt6bFYjGbzXdyXwAAoC2TyeTmdrtRowas+Ym6o4EwICAgPz/f+uuVK1cMBv4mHYqi6PX6O7kvAADQ1u1HQRdMsda5c+f9+/ffuHGDMXb16tX09PROnTpJ6hgAAEB1qOqMMC0t7aefftq1a9fp06ffeeedqKioQYMGRUZG9u7de8KECYmJiYsWLRo2bJjtkiGAZKJpwxzWiqJN2CiAE7BjjVDz90JVZ4Tl5eXXrl1r06bN8OHDr127ps4CGWNr1qxp3779ypUre/fuvXz5cof1EwAAnIGzXUTIqj4jbN++ffv27W8t9/HxeeWVV6R2CQAAoPog1ygAAEjkfMEyGAgBAEAaJ9yFCQMhgB00f+MC1FhOOBJih3oAAHBpmBECAIA0dlxQrzkMhAAAII0z5hrFqVEAAHBpmBECAIA0iiJ8OYTmp1IxEAIAgFRan+oUhVOjAADg0jAjBAAAmZxtQoiBEAAA5HHG/QgxEAIAgDzILAMAAOBcMCMEAABpnPGCegyEAAAgjTOmWMOpUQAAcGmYEQIAgDxOGCyDgRAAAKRRxNf8ND+TioEQAACkUZiiCE7xROtLhzVCAABwaZgRAgCAPHasEWoNAyEAAMgjfh2h5gMnTo0CAIBLw4wQAACkQdJtAABwbVgjBAAAV6YwO2aEDupLVWGNEAAAXBpmhAAAII0TZljDQAgAABI54UiIU6MAAODSMCMEAABpFMWOyyFw+QQAANQa2KEeAABcGtYIAQAAnAtmhAAAII09F9RrPSXEQAgAANIo2H0CAADAuWAgBAAAl4ZTowAAII1d2zA5qC9VhYEQAACksWONUOtxEKdGAQDAtWFGCAAA8jjhBfUYCAEAQB7xNULNFwlxahQAAFwaZoQAACCNIj7B03pCiIEQAADksefyCa0XCatpINTpdCvmz9255svqubtqdu7cucDAQL1er3VHHOvcuXNBQUE6XS0/ne4iD/Ps2bPBwcGu8DBDQkLEt8dzMtX2MA+OH/+vf/2r8jqxof6izXplBpzx9LS3UxIoFouleu7p7NmzFRUV1XNf1aysrMxT01exeuBh1iZ4mLVJtT1Mg8Hg5eUlvVmLxVJeXq7hK1V9AyEAAEANVMtPjAAAAFQOAyEAALg0DIQAAODSMBACAIBLw3WEwi5evJiampqdnR0bG9u8eXNr+dmzZ//zn//cuHHj/vvv79Spk4Y9lOLAgQObN2++fPlyRETEhAkTrKFihYWFS5cuzc7O7tu377Bhw7Tt5J3bsmXL7t27CwoKQkJCJk6c6Ofnp5YXFBQsXbo0Nze3f//+gwcP1raTEq1Zs8bT03P48OHqr2VlZZ999tmpU6fat28/YcIEZ7+a4ocffsjLy1NvN2zYcMyYMertq1evLlu2LC8vb+DAgQMGDNCug9JcunRp+fLlOTk5zZo1mzx5cv369ZnNQRsbGztkyBCt++hMnPu410SvXr3efPPNF198MTU11VqYl5fXqVOngoKCRo0a3XfffSkpKRr28M4VFBTEx8dfvnw5JCRk5cqVvXr1KisrY4yZzeZ+/frt3r27efPmzzzzzPvvv691T+/U6tWrzWZzWFjYr7/+Gh0dffXqVcaYyWTq3bt3ampqWFjYY489tnjxYq27KccPP/zw8MMPz5s3z1rywAMPfPvtty1btlywYMHMmTM17JsU8+bN27x5c2ZmZmZm5oULF9TC8vLynj17HjlypFmzZomJif/5z3807aMEJ0+ejIqKSktLCw0NzcjIUB+pyWTq06dPampq8+bNn3jiiUWLFmndTadiAUEVFRUWi6Vdu3arVq2yFr722msjRoxQb8+bN2/w4MHadE6SioqKsrIy9XZJSUn9+vVTUlIsFsuGDRtCQ0ONRqPFYtm6dWuTJk3U27WA2Wxu1qxZUlKSxWJJSkoKDw9XX2j1Iau3nVpBQUFkZORrr73WvXt3teTYsWPe3t6FhYUWi+X06dNeXl75+fma9vFO9ejR44cffrip8JtvvmnTpo3ZbLZYLElJSS1btlRvO6+BAwe+/PLLNxWqD009UDdu3Ni0aVOTyaRF75wSZoTCuKePUlJSrKdc7rvvvp07d1ZvpyTT6XQeHh7qbbPZXF5e7uPjwxjbuXPnvffe6+bmxhjr06dPfn5+RkaGlh2VJyMjo6Cg4J577mGM7dy5MzY2Vn2hY2Njz507d+bMGY37d8eeffbZZ599tkmTJtaSlJSULl26qK9sWFhYYGDg3r17teugHBs3bnzvvfc2bNhg+e8V0ikpKf3791ezrgwYMODkyZM5OTma9vGOGI3G5OTk4cOHf/7554sXL7ZOfG86aC9cuFALDtpqg4FQjtzc3Lvvvlu93ahRoxs3bhQWFmrbJVn+8Y9/9O7dOzo6mjGWl5dnfZh6vd7Pzy83N1fT3knwwgsvBAYGRkVFzZ8/Xx0IbV9NDw8PX19fZ3+YW7duzcrKSkxMtC20fTUZY40aNXLqEYIx1rp1a09Pz0uXLs2YMSM+Pt5sNrP/fTXr1q171113OfWref78ebPZ/Pjjj585c+bo0aPt2rX7888/2f++mu7u7rXgoK1OCJaRw83NzWQyqbfVG+7u7pr2SI6FCxcmJydblzzd3Nxs8+QZjUbrxNF5vfLKK88999yuXbumT5/etm3bTp06ubu716aHeePGjRkzZnz33Xc35aKsfa/mp59+qt6YNWtWeHj45s2bBw4caPveZIyZTCanfpg6nc5isTz++OPq1xqj0fjee+999tlnte/VrE6YEcoRGBho/TadnZ3dsGFDR2Tkq2YffPDBRx99tG3btoCAALUkMDAwOztbvV1SUnLt2jXbU21OytvbOyAg4P777x84cOD69evZ/z7MoqKioqIip36YO3fuzM7OfvDBBzt27Dh37twjR4507NjRbDbbPkzGWHZ2tlM/TFu+vr6tW7fOyspi//vevHLlSmlpqVM/TIPBoNPpWrdurf4aGRl59uxZdstBW1hY6NQPs5phIJQjPj5+3bp16qmYtWvXxsfHa92jO7Vs2bIFCxYkJycHBQVZC+Pj45OTkwsKChhjSUlJrVq1sr2AxOmYTCaj0ajeLi8vP3LkSEhICGMsPj7+l19+KSoqYox9++237du3DwwM1LKjd6ZHjx7btm1bsmTJkiVLJk6cGBYWtmTJEp1ON3DgwIMHD6ofo7/99ltZWVn37t217qz9jEajdeZ3/vz5Q4cORUZGMsbi4+M3bNhQXFzMGPv222+7devm7y+8PULN4enpOWjQoD179qi/7tmzRx0U4+PjN23apB6069ati46Otn3nwm1oHa3jfJ544omYmBgvL6+wsLCYmJjU1FSLxfLXX3916NChd+/eCQkJjRs3PnHihNbdvCPZ2dmKooSEhMT8188//6z+acKECa1bt54yZYq/v/+PP/6obT/v0Llz5xo3bjxixIgJEyY0bdq0f//+JSUl6p/GjBnTpk2byZMn+/v7b9q0Sdt+SrR06VJr1KjFYnn55ZdDQ0MTExMDAgKWLFmiYcfu3OnTpw0Gw6hRoxISEnx9fR9//HG13Gw2x8fHR0dHT5o0yc/Pb/v27Zp2U4IDBw40atRo0qRJgwcPbtGiRU5Ojlp+//33Ww/aX375RdtOOhfsPiHs5MmTtoEw4eHhatxdWVnZtm3b/vrrr/79+/v6+mrXQQnKy8uPHj1qWxIaGqpebG6xWHbt2pWdnd29e/emTZtq1EFpzp07d+jQodLS0pYtW7Zv395abrFYdu7cmZeX16NHj+DgYA17KNeVK1fy8/NbtWplLUlNTT158mR0dHRERISGHZMiPT09PT3dbDZHRUWFh4dby81m844dOy5dutSrVy+nntxb5efnb9u2rUGDBj179rSuwlgslpSUlNzc3Fp20FYDDIQAAODSsEYIAAAuDQMhAAC4NAyEAADg0jAQAgCAS8NACAAALg0DIQAAuDQMhAAA4NIwEAIAgEvDQAgAAC4NAyEAALg0DIQAAODS/h/uNqfIGU5fbwAAAABJRU5ErkJggg==",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"execution_count": 1
}
],
"cell_type": "code",
"source": [
"using DFTK\n",
"using StaticArrays\n",
"using Plots\n",
"\n",
"# Unit cell. Having one of the lattice vectors as zero means a 2D system\n",
"a = 15\n",
"lattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];\n",
"\n",
"# Confining scalar potential, and magnetic vector potential\n",
"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...);\n",
"\n",
"\n",
"# Parameters\n",
"Ecut = 20 # Increase this for production\n",
"η = 500\n",
"C = η/2\n",
"α = 2\n",
"n_electrons = 1; # Increase this for fun\n",
"\n",
"# Collect all the terms, build and run the model\n",
"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": 1
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.5"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.5",
"language": "julia"
}
},
"nbformat": 4
}