{
"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.170912e+01 6.360302e+00\n",
" * time: 0.0043621063232421875\n",
" 1 2.931172e+01 6.265971e+00\n",
" * time: 0.012782096862792969\n",
" 2 2.195605e+01 5.139085e+00\n",
" * time: 0.033094167709350586\n",
" 3 1.424589e+01 1.711165e+00\n",
" * time: 0.057192087173461914\n",
" 4 1.138878e+01 8.367846e-01\n",
" * time: 0.07734203338623047\n",
" 5 1.123285e+01 1.273805e+00\n",
" * time: 0.1677260398864746\n",
" 6 1.050670e+01 1.369371e+00\n",
" * time: 0.18364310264587402\n",
" 7 1.005320e+01 9.154337e-01\n",
" * time: 0.19940519332885742\n",
" 8 9.856926e+00 7.046263e-01\n",
" * time: 0.21506619453430176\n",
" 9 9.736397e+00 4.559933e-01\n",
" * time: 0.23055815696716309\n",
" 10 9.622823e+00 3.811662e-01\n",
" * time: 0.24585914611816406\n",
" 11 9.559846e+00 5.473489e-01\n",
" * time: 0.2577500343322754\n",
" 12 9.522914e+00 5.880825e-01\n",
" * time: 0.2696340084075928\n",
" 13 9.485757e+00 3.468251e-01\n",
" * time: 0.2852480411529541\n",
" 14 9.454671e+00 4.759135e-01\n",
" * time: 0.29700517654418945\n",
" 15 9.428250e+00 2.982504e-01\n",
" * time: 0.3087351322174072\n",
" 16 9.405865e+00 2.834709e-01\n",
" * time: 0.32050609588623047\n",
" 17 9.374254e+00 2.122735e-01\n",
" * time: 0.33239006996154785\n",
" 18 9.367758e+00 2.112867e-01\n",
" * time: 0.34413909912109375\n",
" 19 9.366211e+00 2.807264e-01\n",
" * time: 0.35590505599975586\n",
" 20 9.359579e+00 2.298018e-01\n",
" * time: 0.3674910068511963\n",
" 21 9.348747e+00 1.824376e-01\n",
" * time: 0.3790290355682373\n",
" 22 9.344617e+00 1.416999e-01\n",
" * time: 0.39046311378479004\n",
" 23 9.340467e+00 1.116812e-01\n",
" * time: 0.4022071361541748\n",
" 24 9.336310e+00 1.668980e-01\n",
" * time: 0.4138951301574707\n",
" 25 9.330638e+00 1.219833e-01\n",
" * time: 0.42561912536621094\n",
" 26 9.328471e+00 7.769024e-02\n",
" * time: 0.4372432231903076\n",
" 27 9.327201e+00 1.202868e-01\n",
" * time: 0.44901418685913086\n",
" 28 9.323275e+00 5.881263e-02\n",
" * time: 0.46082305908203125\n",
" 29 9.321338e+00 6.906074e-02\n",
" * time: 0.47248411178588867\n",
" 30 9.320866e+00 4.084347e-02\n",
" * time: 0.48805999755859375\n",
" 31 9.320345e+00 6.575464e-02\n",
" * time: 0.49985718727111816\n",
" 32 9.319015e+00 3.246313e-02\n",
" * time: 0.5115740299224854\n",
" 33 9.318330e+00 2.500694e-01\n",
" * time: 0.5197131633758545\n",
" 34 9.313040e+00 1.324316e-01\n",
" * time: 0.5319662094116211\n",
" 35 9.308644e+00 9.831239e-02\n",
" * time: 0.5438470840454102\n",
" 36 9.307580e+00 8.185243e-02\n",
" * time: 0.555927038192749\n",
" 37 9.307414e+00 4.341591e-02\n",
" * time: 0.5722291469573975\n",
" 38 9.307206e+00 5.370488e-02\n",
" * time: 0.5845181941986084\n",
" 39 9.306217e+00 4.127214e-02\n",
" * time: 0.5965940952301025\n",
" 40 9.303663e+00 9.985026e-02\n",
" * time: 0.6088690757751465\n",
" 41 9.301147e+00 2.859551e-01\n",
" * time: 0.6171572208404541\n",
" 42 9.295117e+00 1.973940e-01\n",
" * time: 0.6293151378631592\n",
" 43 9.292691e+00 1.684452e-01\n",
" * time: 0.6456341743469238\n",
" 44 9.289944e+00 1.457288e-01\n",
" * time: 0.6577341556549072\n",
" 45 9.282724e+00 8.454584e-02\n",
" * time: 0.6701111793518066\n",
" 46 9.282081e+00 6.871866e-02\n",
" * time: 0.6822412014007568\n",
" 47 9.281218e+00 5.156265e-02\n",
" * time: 0.6943860054016113\n",
" 48 9.277345e+00 9.731595e-02\n",
" * time: 0.7065331935882568\n",
" 49 9.272827e+00 7.534179e-02\n",
" * time: 0.7186992168426514\n",
" 50 9.271000e+00 6.589796e-02\n",
" * time: 0.7306911945343018\n",
" 51 9.267663e+00 7.901974e-02\n",
" * time: 0.7431690692901611\n",
" 52 9.260631e+00 2.819699e-01\n",
" * time: 0.7518370151519775\n",
" 53 9.251532e+00 2.216085e-01\n",
" * time: 0.7646300792694092\n",
" 54 9.247834e+00 2.295872e-01\n",
" * time: 0.7772471904754639\n",
" 55 9.241000e+00 1.250288e-01\n",
" * time: 0.7899470329284668\n",
" 56 9.236274e+00 1.444614e-01\n",
" * time: 0.8023531436920166\n",
" 57 9.231389e+00 1.084476e-01\n",
" * time: 0.8148891925811768\n",
" 58 9.227486e+00 1.118260e-01\n",
" * time: 0.8277740478515625\n",
" 59 9.222604e+00 9.252583e-02\n",
" * time: 0.8406310081481934\n",
" 60 9.215846e+00 1.135831e-01\n",
" * time: 0.8533680438995361\n",
" 61 9.202626e+00 1.671377e-01\n",
" * time: 0.9213471412658691\n",
" 62 9.195539e+00 1.501892e-01\n",
" * time: 0.9332292079925537\n",
" 63 9.192899e+00 1.098381e-01\n",
" * time: 0.94482421875\n",
" 64 9.189271e+00 1.353576e-01\n",
" * time: 0.9566349983215332\n",
" 65 9.179416e+00 1.323246e-01\n",
" * time: 0.9682102203369141\n",
" 66 9.169957e+00 1.149294e-01\n",
" * time: 0.9797391891479492\n",
" 67 9.167443e+00 1.102427e-01\n",
" * time: 0.9914782047271729\n",
" 68 9.158277e+00 1.766439e-01\n",
" * time: 1.0033190250396729\n",
" 69 9.148612e+00 1.337518e-01\n",
" * time: 1.0150411128997803\n",
" 70 9.135499e+00 1.433978e-01\n",
" * time: 1.0266201496124268\n",
" 71 9.124734e+00 2.055086e-01\n",
" * time: 1.03810715675354\n",
" 72 9.116504e+00 2.294350e-01\n",
" * time: 1.049926996231079\n",
" 73 9.111598e+00 2.020176e-01\n",
" * time: 1.0618531703948975\n",
" 74 9.087792e+00 1.637168e-01\n",
" * time: 1.0735700130462646\n",
" 75 9.078883e+00 1.523720e-01\n",
" * time: 1.0853102207183838\n",
" 76 9.074586e+00 1.406730e-01\n",
" * time: 1.096912145614624\n",
" 77 9.067875e+00 1.607302e-01\n",
" * time: 1.1086561679840088\n",
" 78 9.064213e+00 1.994306e-01\n",
" * time: 1.1203861236572266\n",
" 79 9.058364e+00 1.699841e-01\n",
" * time: 1.1319911479949951\n",
" 80 9.044933e+00 2.498340e-01\n",
" * time: 1.1437511444091797\n",
" 81 9.032408e+00 1.621983e-01\n",
" * time: 1.1551721096038818\n",
" 82 9.017918e+00 1.947078e-01\n",
" * time: 1.166896104812622\n",
" 83 9.003244e+00 1.474720e-01\n",
" * time: 1.1785449981689453\n",
" 84 8.991946e+00 1.343248e-01\n",
" * time: 1.1903181076049805\n",
" 85 8.979491e+00 1.508666e-01\n",
" * time: 1.2020840644836426\n",
" 86 8.968311e+00 1.334005e-01\n",
" * time: 1.2136991024017334\n",
" 87 8.940874e+00 1.688551e-01\n",
" * time: 1.2255311012268066\n",
" 88 8.926509e+00 1.947394e-01\n",
" * time: 1.237192153930664\n",
" 89 8.901758e+00 2.167859e-01\n",
" * time: 1.2491512298583984\n",
" 90 8.882959e+00 2.351602e-01\n",
" * time: 1.2612621784210205\n",
" 91 8.869064e+00 2.053316e-01\n",
" * time: 1.2735660076141357\n",
" 92 8.857656e+00 1.344689e-01\n",
" * time: 1.2855541706085205\n",
" 93 8.842519e+00 1.821028e-01\n",
" * time: 1.2976040840148926\n",
" 94 8.834062e+00 1.667926e-01\n",
" * time: 1.3098289966583252\n",
" 95 8.816064e+00 2.249383e-01\n",
" * time: 1.3218111991882324\n",
" 96 8.790293e+00 2.172484e-01\n",
" * time: 1.3339672088623047\n",
" 97 8.768550e+00 1.726694e-01\n",
" * time: 1.3462011814117432\n",
" 98 8.746495e+00 1.578501e-01\n",
" * time: 1.358199119567871\n",
" 99 8.714193e+00 2.055451e-01\n",
" * time: 1.370474100112915\n",
" 100 8.682707e+00 1.934821e-01\n",
" * time: 1.382620096206665\n",
" 101 8.662307e+00 6.809440e-01\n",
" * time: 1.390726089477539\n",
" 102 8.617488e+00 3.714069e-01\n",
" * time: 1.4028210639953613\n",
" 103 8.574112e+00 4.137236e-01\n",
" * time: 1.4148712158203125\n",
" 104 8.557076e+00 2.238488e-01\n",
" * time: 1.4268381595611572\n",
" 105 8.549236e+00 1.320331e-01\n",
" * time: 1.4388420581817627\n",
" 106 8.542633e+00 1.256028e-01\n",
" * time: 1.451146125793457\n",
" 107 8.536839e+00 1.454259e-01\n",
" * time: 1.463395118713379\n",
" 108 8.530770e+00 1.345589e-01\n",
" * time: 1.4754371643066406\n",
" 109 8.522545e+00 8.686581e-02\n",
" * time: 1.4878921508789062\n",
" 110 8.518893e+00 8.526019e-02\n",
" * time: 1.5006051063537598\n",
" 111 8.516610e+00 7.538590e-02\n",
" * time: 1.5131750106811523\n",
" 112 8.515162e+00 5.385389e-02\n",
" * time: 1.5258491039276123\n",
" 113 8.514090e+00 7.610562e-02\n",
" * time: 1.5383970737457275\n",
" 114 8.514016e+00 6.616039e-02\n",
" * time: 1.55122709274292\n",
" 115 8.513848e+00 4.204365e-02\n",
" * time: 1.5637531280517578\n",
" 116 8.513359e+00 3.889917e-02\n",
" * time: 1.576620101928711\n",
" 117 8.512972e+00 3.307892e-02\n",
" * time: 1.5891990661621094\n",
" 118 8.512772e+00 3.856756e-02\n",
" * time: 1.6510391235351562\n",
" 119 8.512628e+00 3.468493e-02\n",
" * time: 1.6629531383514404\n",
" 120 8.512485e+00 3.493368e-02\n",
" * time: 1.6748130321502686\n",
" 121 8.512295e+00 3.057295e-02\n",
" * time: 1.6866340637207031\n",
" 122 8.512248e+00 1.633541e-02\n",
" * time: 1.6984121799468994\n",
" 123 8.512187e+00 1.575410e-02\n",
" * time: 1.710076093673706\n",
" 124 8.512101e+00 1.708457e-02\n",
" * time: 1.7219312191009521\n",
" 125 8.512074e+00 1.346665e-02\n",
" * time: 1.7335050106048584\n",
" 126 8.512057e+00 7.816918e-03\n",
" * time: 1.7450881004333496\n",
" 127 8.512014e+00 9.067979e-03\n",
" * time: 1.7567050457000732\n",
" 128 8.511993e+00 8.302427e-03\n",
" * time: 1.7684612274169922\n",
" 129 8.511965e+00 8.005615e-03\n",
" * time: 1.7802441120147705\n",
" 130 8.511950e+00 6.628085e-03\n",
" * time: 1.7920451164245605\n",
" 131 8.511941e+00 9.270250e-03\n",
" * time: 1.8037691116333008\n",
" 132 8.511916e+00 7.158347e-03\n",
" * time: 1.8152611255645752\n",
" 133 8.511903e+00 6.406994e-03\n",
" * time: 1.8269431591033936\n",
" 134 8.511891e+00 4.684444e-03\n",
" * time: 1.8388011455535889\n",
" 135 8.511882e+00 4.126750e-03\n",
" * time: 1.8505051136016846\n",
" 136 8.511872e+00 6.334414e-03\n",
" * time: 1.8623671531677246\n",
" 137 8.511868e+00 3.683836e-03\n",
" * time: 1.873946189880371\n",
" 138 8.511865e+00 3.934017e-03\n",
" * time: 1.8856120109558105\n",
" 139 8.511859e+00 2.641955e-03\n",
" * time: 1.8974270820617676\n",
" 140 8.511856e+00 3.600353e-03\n",
" * time: 1.9091150760650635\n",
" 141 8.511853e+00 2.037216e-03\n",
" * time: 1.9207031726837158\n",
" 142 8.511853e+00 2.115360e-03\n",
" * time: 1.9324491024017334\n",
" 143 8.511852e+00 2.490188e-03\n",
" * time: 1.944092035293579\n",
" 144 8.511849e+00 1.629815e-03\n",
" * time: 1.9560530185699463\n",
" 145 8.511848e+00 1.600314e-03\n",
" * time: 1.9678380489349365\n",
" 146 8.511848e+00 1.436800e-03\n",
" * time: 1.9793200492858887\n",
" 147 8.511848e+00 1.375760e-03\n",
" * time: 1.9911072254180908\n",
" 148 8.511847e+00 1.410468e-03\n",
" * time: 2.0032410621643066\n",
" 149 8.511847e+00 1.050250e-03\n",
" * time: 2.0152721405029297\n",
" 150 8.511846e+00 9.749387e-04\n",
" * time: 2.027245044708252\n",
" 151 8.511846e+00 9.498488e-04\n",
" * time: 2.0394680500030518\n",
" 152 8.511846e+00 7.896550e-04\n",
" * time: 2.051664113998413\n",
" 153 8.511846e+00 5.867535e-04\n",
" * time: 2.063836097717285\n",
" 154 8.511846e+00 7.477354e-04\n",
" * time: 2.076181173324585\n",
" 155 8.511846e+00 6.299342e-04\n",
" * time: 2.0884170532226562\n",
" 156 8.511845e+00 3.499972e-04\n",
" * time: 2.100548028945923\n",
" 157 8.511845e+00 4.126935e-04\n",
" * time: 2.1123812198638916\n",
" 158 8.511845e+00 4.131686e-04\n",
" * time: 2.1246612071990967\n",
" 159 8.511845e+00 2.335790e-04\n",
" * time: 2.1367311477661133\n",
" 160 8.511845e+00 3.042421e-04\n",
" * time: 2.1487061977386475\n",
" 161 8.511845e+00 1.987897e-04\n",
" * time: 2.160861015319824\n",
" 162 8.511845e+00 1.966069e-04\n",
" * time: 2.1728310585021973\n",
" 163 8.511845e+00 1.707244e-04\n",
" * time: 2.18471622467041\n",
" 164 8.511845e+00 1.530225e-04\n",
" * time: 2.196885108947754\n",
" 165 8.511845e+00 1.407033e-04\n",
" * time: 2.2090752124786377\n",
" 166 8.511845e+00 1.228115e-04\n",
" * time: 2.221372127532959\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVxU1f8/8HNnWEREQVDZQRFcUNwX3BMTNXHJLfelRcsyK1PTT4vlx8qs7GPlmmUfyy2Dj6Wl4AL2MVPMFRdUcENBZRFFGGaY+f1x+85vPnbeyBnvcBnm9Xzwx3A43Dkzc2fOnHtf9xzJZDIxAAAAR6VRuwEAAABqQkcIAAAODR0hAAA4NHSEAADg0NARAgCAQ0NHCAAADg0dIQAAODR0hAAA4NDQEQIAgENzUrsBAAAA5UlPT//iiy8KCgqGDh06ZMiQv1coKir67LPPTp8+HRUVNWPGjBo1ajDGjh8/npCQcPHixQYNGjz99NNNmzaVK7/77rv379+Xb0dFRY0ZM6aSOkKj0fiPd96dM39+5dxdJTMYDE5O1f8rBR5mdYKHWZ1U2sOs4+r80DrnzuUVFxuENnvp0oWcnJNTp07l/vXWrVtdu3adNm1ax44dX3jhheLi4tGjRz9QZ/To0WVlZZMnT161atXRo0c3bNjAGHv99dfbtGkTGxt77Nix9u3b//HHH5GRkYyxf/3rX5MnT/bx8WGMubu7M8akyplrtLi4ePq7iwdOeq4S7gsAAGzhySZ+D63TuvW648dvCW741JAhd+Lj47l/++ijj/bu3btjxw7G2L///e+lS5ceOXLEssLZs2fbtm2bk5Pj4eGRl5fn7+9/7ty5kJAQo9Go0fx1+i8uLq5Dhw5vvfUWY8zHx+f3338PDw83bwHnCAEAoOr6/fffe/XqJd/u1avX0aNHS0pKLCscPHiwTZs2Hh4ejLG6detGRkb+8ccfjDFzL8gYy8nJqVevnvnXZcuWzZ07NyEhQR4KoiMEAADFSJI1P+XIzs729vaWb/v4+JhMpuzsbKqCXOfGjRuWFVauXHn79u0JEybIvw4aNCg0NNTDw+PVV1+dNGkSQ1gGAABUt2/fvvbt21uWjBgxYs6cOYyxGjVqlJaWyoU6nY4xVrNmTcuabm5uer3e/GtJSYllhYSEhLfffjspKUk+HcgYW7t2rXxjwoQJDRs2nD9/PjpCAABQjCRJUvlDvL8xmaTWrVsvWbLEsjAgIEC+ERgYeOXKFfn2lStXatSoIedcLGuaKzDGrl69GhgYKN/esWPH1KlTd+zY0aJFi7/fb1BQUN26da9fv45DowAAoChJ+MfT07Pd//L19ZU3NmzYsK1bt8oXPKxfv37IkCHyyb+kpKT09HTGWN++fa9evSonaH777bc7d+7I5xQTExMnT568bdu2du3amZtWWFhoHl/+8ssvd+/ejYyMxIgQAACqroEDB65du7Zt27bBwcGnT5/evXu3XD5//vzRo0dHRETUrl37ww8/7N+/f+fOnX///fePP/7Yzc2NMTZt2rTi4uIxY8bI9UeOHPn+++///vvv48aNa9mypU6nS0tLW758eb169dARAgCAch4WfhGl1WoTEhKOHTtWWFjYoUMH8/m/hIQE82m/559/fuDAgefOnVuxYoW/v79cuHfvXoPh/1/RKMdKY2Njjx49euHCBVdX1+bNm9epU4chLAMAAAqy4hwhMz2kviRJbdq0eaDQz+9/LmoMCgoKCgqyLAkODuZuLTAw0HwSUYZzhAAA4NDQEQIAgEPDoVEAAFDMQy+Q/7uHHRm1OXSEAACgGCvOESobrrECDo0CAIBDw4gQAACUI18jb1fQEQIAgGIk8UOdqh8aRUcIAACKseY6QrV7QpwjBAAAh4YRIQAAKMeKc4Q4NAoAANWJvfWDgh3hsWPHzp8/X69evc6dO9eoUUMuPH78eHp6evPmzSMjI23QQgAAABuq6DnCsrKy8ePHDxo06IcffnjnnXfWr18vl7///vtPPPHE9u3b+/Tps2zZMpu1EwAA7IAk/ZWXEaD2mLCiI8Lly5efPHkyLS1NXsnCZDIxxnJzcxcuXJiamtqsWbPU1NQ+ffpMnjy5Vq1aNmwvAABUZXZ4jrCiI8L169e//PLLubm5Bw8eLCoqktOxiYmJ4eHhzZo1Y4y1b9/e29s7OTnZho0FAICqzYoBofDlFkqr6IjwwoULGzduXLFihYeHx9mzZ3/++efWrVtnZWVZrv8UFBR07do1agvyIBIAAOyUyWRSvdOyhYqOCO/fv+/h4XHw4MGkpKSJEyfOmjWLMabX67VarbmOk5OTXq/n/rvJZEJHCABg14xG40PryKtPCP2orqIdob+/f0xMjPxd4PHHHz9+/DhjzM/P79atW+Y6N2/efGDJYDNJkjQaXLwPAGDHLEc+JMmqH1VVtHPq1atXZmamfDszM9Pf358x1q1btz///DM/P58xduPGjfT09OjoaBs1FAAAwBYqeo7w1Vdf7dmzp6enZ+3atf/5z39++umnjLGwsLDBgwcPGTJk7NixX3/99fjx4+UOEgAAHJN8vFPwX1RW0RFh8+bN9+/ff//+/evXr8fHxz/11FNy+b///e/Ro0enpaU9++yzK1assFk7AQDAHoifI1T9NKHAzDJNmzZduHDhA4XOzs7Tpk1TtEkAAGC3qsA5P1EIsAAAgEPDpNsAAKAY+YJ6wf+xkwvqAQAAHsqKc35q94M4NAoAAI4NI0IAAFCOHYZl0BECAIBirJhEW/X5S9ERAgCAYuxwQIhzhAAA4NgwIgQQp8hKKnb3tRmgIuxwYV50hAAAoBwrzhHaqCUVhkOjAADg0DAiBAAAxUjiF8irHRpFRwgAAMqRJPHLIdTuCdERgv0gIiqCyRUFtqJMVkZ0K+RnhcCHiBLbACiXHYZlcI4QAAAcGkaEAACgGDscEKIjBAAA5VgzxZraXSEOjQIAgEPDiBAAAJRjh8dG0RGCmkxCyUmiNreU3DK1EbGWCBUTb3OlEua8cvE4ukj0VO2PLajKrLiOUPWOEIdGAQDAoWFECAAAyrHDuUbREQIAgHJwjhAAAByZJNnfXKM4RwgAAA4NI0JQEp29FMhqUhsxGYmN8MqJuvRGhNKkIvlVRqZGhYrJiYwlDaecV0ZWZoxJvK/EikRPVf+yD5VMYpLoBfKqX1CPjhAAAJRjxTlCteHQKAAAODSMCAEAQDFWhGVUH0GiIwQAAMVITPw6QrXPJKMjhIcQyr+I5lyM3JxLGVGZKC/jlZMbMRqJjXOL+SEasWQNERghwy9EuUbL37hGwznBodHyN6Ilyrn1qY1oRBI3NLGsENgNnCMEAACwLxgRAgCAciThQ52qHwZARwgAAIqxwyOjODQKAACODSNCAABQDibdBvtFBj6phCQvfclNgTLGjGX8rGaZgVNu0PM3YuBVZoyV6XkboSoT5XTKVCAcKxQbpc6LUJlMMgjqxDmu48QrZIxpnfnl3PpOzgL3yBjTaAXyq5KkyIR0UOVI1izDhMsnAACgurBihXrVv+WgIwQAgCrt7t27mzZtKiwsHDBgQNOmTf9ewWQybd++PS0trVWrVv369TP/V1JS0sWLF319fYcMGVKrVi1z/d27d6empjZp0mTw4MGSJHgFLAAAQHkkq35o9+/f79Sp0/bt27Ozs6Ojo5OTk/9e5+WXX547d25RUdHMmTPfeOMNubB3796rVq3KyclZv359s2bNrl+/Lpe/9957zz33XFFR0YIFC5599lmGESEAACjIinOE5R8b/f777z08PH788UdJkvz8/BYuXNizZ0/LCtevX1+9evX58+cDAwMnTpzYsmXLWbNmeXt7//LLLz4+PnKd7t27f/vtt3Pnzi0sLPzoo49+++23qKiol156KSQkZN68eRgRAgBA1ZWYmDhw4EC5c42Li9uzZ4/BYLCssHfv3sjIyMDAQMZYWFhYw4YNU1JSGGPmXpAx5ubmptVqGWMHDx708vKKiopijNWrV69jx45JSUkYEToc0XSo0MyfVFbTUMqfylNfyqmvF6lM1TfwoqSMiJgyYsJSRjxMm841KjpNKDcI6kSkQ51d+DOWOrtw6gtVZow58epT+VVyIlNiSlWs+msvrFh9ovzqWVlZffr0kW/7+fkZjcbs7Gy525Ndv37dz8/P/Kufn5/5KKhs7969qampX3/9NVUZHSEAAKjs2LFjU6dOtSzp2bPnmDFjGGMajcb81VO+8cCXSEmSLL+bmkwmywrHjh0bPXr0+vXrAwICqMroCAEAQDlWXEcoMU9Pz3bt2lkWNm7cWL7h5+eXnZ0t375x44ZWq23QoIFlTcsKjLHs7Gx/f3/59qlTp/r37//5558PGDCAqhwTE6NAR1haWpqVlRUYGOjs7PzoWwMAAEcTGhr63HPPcf8UGxu7cuXKf/zjHxqN5qeffurdu7eTkxNj7NKlS7Vr165bt27v3r2feeaZq1evBgUFXbhw4dKlSz169GCMpaen9+vXb8mSJcOHDzdvLTo6Oj8///jx461atbp58+ahQ4e++eabioZlZs2aJVm4f/++XL5t27aAgID+/fsHBQUlJSU90jMBAAB2Tj5HKPRT/knCp556qqioaPDgwbNmzXrvvffefPNNuXzUqFHffvstY8zPz2/atGmxsbHz589/4oknZs6cWbduXcaY3P/95z//GTly5MiRI1esWMEY8/DwmDNnztChQ+fPn//444+PHz8+NDRUYEQ4f/78hQsXWpbodLqnn3563bp1AwYM2Lx585QpUzIzM+VkDlQqKv/CXTuXHxYhQzFk/oWXOinV8XMupSUGopxTn94IFaLhlOtFwzLU1GvcKdaI51Bo0jDqIl5qijVqYjNuWMZZMCzjUoNT7uIqUJkx5lKD8wRQG6HiPE4mKlzDLeXWpafsQrimcij6PNesWfPgwYNbt269c+fOoUOHwsPD5fKPPvrIHJn59NNPd+7ceerUqS+++MKcrPnggw+KiorM22nYsKF8Y/78+V27dj1y5MiiRYvkQ6Zih0aLi4vd3NzMv+7atat27dryhoYPHz5jxoyUlJTHHnvMyocLAADwN7Vq1Zo4ceIDhfLxT7PY2NjY2FjLEvN5wb/r1atXr169zL8KXEe4ZMmSevXq+fr6fvzxx3JJZmZmRETEXxvSaBo3bnzp0qWKbxAAAKoZySrqtrmiHeG0adPy8/Pv3buXkJDwz3/+c+vWrYyxu3fvWg4Qa9asWVhYyP13k8lUVsY/rgUAAHZBr9c/tI4V5wjV7gcr3BE2btxY7vM6d+48fvz4X3/9lTFWv379/Px8c538/PwHUq1mkiTh3CEAgF2ryKUBVvSDqp+8tWaKtfz8fHd3d8ZYq1atjh49WlpayhgrKiqSZ/5WuIEAAAC2VNGwzNtvv929e3cPD499+/Zt2bLlt99+Y4x17NgxIiLilVdemTp16tKlS6Ojo5s1a2bL1gIfNx3KiHBjGbFGLjUnmV7HL9fxgqC6Yn46lCovKealRkUipowxPS9lSqVGqYdJLR1MTLHGrSuYGqUW5iXnHiNilkKpUZEgqEsN/udDjVL+RrjPbZmBvxFXYuMmV/6z6MT71q6lvspr+BtRfQFYh6D+AE9YRUeEJpPpgw8+mDlz5okTJ/bt22eeAiAhIaG4uPiZZ55xdnbetGmTzdoJAAD2wA7PEVZ0RPjuu+9yy/39/deuXatcewAAACoV5hoFAADFSPTKKuS/2MuIEAAA4OHs8BwhOkJ7QmY0iBm/uLkYKi1CRVGonEvxfU55SRH/MqMSXmXGWAlv4zpegobRIRruOoUGPX8j5DqFSqxHKJSWUWo9Qm5YxslZbClBbi7G1Y14Dkv5nxv8sAzxxJp4c9cxxkwmBT6U6BANp0z14Ug1Y816hGq/BFihHgAAHBpGhAAAoBhrpkxTe0iIESEAADg0dIQAAODQcGgUAAAUY01YxjYtqTh0hFUUN5ZIhRXJNXW5a+cS6VAq2FlMBEGL73HKycpFAqnRUiKnSi3Yy02NclfrZfQCvGUG4rnlLsxr09QouTCvwIK91AK8VGqUP00dlQ7lPeGMeG6pnZNOjXKLxVBTqRFPLVGq+sezfZKY8DlCu1mGCQAAoFrCiBAAAJRjxQX1ag++0RECAIByqsAk2qLQEQIAgGKsuI4Q5wgBAADUhBGhysj4Ie8PZDqUSEJy19Ql5w4lAp/37xLl3NQor5AxVnyfX85tDDnXKJUa5ZVTc4pSz5VRidQo9Wpyv++KpkY1RGrUiZcadXLmP1fUwrz8+VrJhK3AIsZGYiJc0XSo0HNIlTtJnOdKw39KGNKkjgMdIQAAKMaqQ6M2aktF4dAoAAA4NIwIAQBAMVbMLIPLJwAAoHpRu2MThY5QddQEVJxCoanUGGM63kq23AV1GZ1z4YZiGGP375ZyNkJNsUbcqY5XriMmgVMkLOMgU6xxV+tljDlTGSLumrrUcyUyaxr5XBGowYTEe1q4hYxe3Jj73FIbkSSRlxP+jyQJXw6ByycAAADUhBEhAAAox4rVJ9QeY6MjBAAA5Vgx16jacGgUAAAcGkaEAACgGIlJ1HqQ5f2TqtARVhJ6KjV+OTesSE15RcUpuROYlZBr5woswEvVpxbgpTbOn2KNSI1y06FUud6xU6PORGqUfvjcNXWJqdSINXXJh89DP3x+fW4QVEukQ7nPCbUROjUq8Oms+lmuqsOaFerVfvbQEQIAgKLsbECIc4QAAODYMCIEAADFWDPpto2aUmHoCAEAQDH2ONcoDo0CAIBDw4iw0ghMz8iIwJ6hlJiEkzenKKNSo9Rco1Tgk1hTlzt9KJUOpe60CqVGifkzufO72jY1SiUhRRKSZVRqlHiY3HJq5xSbPVRk7lBWzsMXmVKVLOclb6knVkMEWCX+Qr5qD2qqDisuqFf7yUNHCAAAyrHmHCGuIwQAgOrCDgeEOEcIAACODSNCAABQjBWXT6g+tQw6QuVxowRCU6kxItOhLyWmWCPSJSXFnPISXj6lnHJunoVRa+oKboR7p9TDcZSFeRUJyxiIsAwVzrLZmrrkjHGCs6M5ufDCMi784IqzC39XcebVd3Lm7xLUC0HFfPiVVT/qpwp7OzaKQ6MAAODQMCIEAADFWDPptm1aUnHoCAEAQDGSJLZwB6sCB5DREQIAgHLs8PoJnCMEAACHhhFhJSFTo8SUVwY9p1xPTbFGxCm5U6/peFFSxlgpGfgkVv3lZTup2dGocm5AlEyNEjPJccO0VGqUWtzYqERqlHqVuUd+RBfm1RAL8zo5ce6V2q9M/EevzJq63DilRsv/tk2mQ52I/ZAX+HRx5VcurcFPk7rw3kHcKCljTMt7YhljGl511Q/uVR12OCBERwgAAAqyw+sIcWgUAAAcGkaEAACgKNWPdQoSGxHq9fq4uLgxY8aYS86cOdO7d28/P79+/fplZGQo3TwAALAn8nWEoj/qEhsRLlq06MqVK+Zlukwm05NPPjlu3LitW7cuWbJk1KhRhw8ftkEjqyoqX8CLTNDrDhJhGSWmWCOiKES+gErciJQLLRmoVAuFwjJ2MMUaNSdZGb+cuwtRU/cpspSg0KxpZCiGWDLQmTeVGmPMmbdXlJbwcy7UO0LvxttViF3CuYzfEv57mfwop9ZjpOrbPSvmGhU+p6g0gRHhyZMnt23b9uqrr5pLUlJScnNz586d6+Xl9dZbb507d+7PP/+0QSMBAABspaIdocFgePbZZ5ctW+bs7GwuPHPmTKtWrbRaLWPM1dW1WbNmp0+ftkkzAQDALkjiP2qraEe4ZMmS6OjoLl26WBbm5uZ6eHiYf/X09Lx9+zb3300mk8HAP94FAAB2Qa/XP7SONScI1e4LK9QRXrp06fPPP588eXJGRsbNmzd1Ol1GRobJZPLy8rp375652p07d7y9vblbkCTJyQkJVQAAO2Z5RLA6qVDndPv2bV9f3ylTpjDG8vLycnJyRo4cuX///saNG585c8ZkMkmSZDAYzp8/HxYWZuMGAwBA1SUP8wT/RWUV6gjbt2+fmpoq3/7+++8/+ugj+dfevXtLkvT1119PmTLlyy+/rFevXnR0tA0bW8WQoVFuIbUAL5EaLeOFHqkp1shybrCTiJ6S5SJBUNHUKLflwi3kbcS2qVHi1RRCre9KZTKNRqrcVmvqSsQBI60SqVE9kRrVu/KDoNxXX2i/osq57zVGvzeduE84sYovRfWPfhuywRxrJpMpNTW1oKAgOjq6Vq1a3DoZGRnnz59v2rRpSEiIZXlWVhZjLCAgwFxy+fLlsrK/9gR3d/cGDRoIzyxTq1Yt8xadnJw2bdr00Ucfubu7r1mzZsOGDaqnYAEAoDoxGAwDBw6cNGnSJ598EhERcfbs2b/X+fzzz6Ojo1esWNG+ffuvvvpKLty4caOPj0/Dhg2HDh1qWbldu3aDBw8eOXLkyJEjly5dyqyYWWbQoEGDBg0y/xodHW0+Oiq6KQAAqHbEryMsd0j4008/Xbhw4dixY25ubnPmzHnnnXc2btxoWaGwsPCNN95ITk5u27btgQMHBg4cOGbMGDc3t86dO//xxx8HDhxYtmzZA9v88ccfw8PDzb8qM9coekEAAGA2SI3Gx8cPGzbMzc2NMTZu3LiEhASj8X+OZu/cuTM4OLht27aMsS5dunh6eu7du5cxFhoaSsVWLl26dOrUqZKSEvlXTLoNAADKseI6wnI7wqtXrwYHB8u3g4ODdTrdrVu3LCtcu3bN8rxgSEjItWvXytmgi4vL7NmzR44cGRAQsHXrVoZJtwEAQHWXLl1atWqVZUnLli3l9KVOpzNftuHi4sIYKy4utqyp0+ksL89zcXExD/W4zp07J18Bv3HjxokTJ/bq1QsdofK4eT0qZig01yiVhNSLlBv0YkE7coVbbq5VsIXcjQjdo1IboV4IodSoyFSjiqVGhfKhQhOcaogkpN6J/xxqec+tk+ALQe+HvICxIvsbtVYztUvwijUiazJXb9bMNcqkgoKCI0eOWBZ6eXnJN3x9fc1Ttdy6dUuSJD8/P8uavr6+ubm55l9v3br1QIUHmOeBeeqpp1566aWTJ0+iIwQAAMVYd/VE69atV65cyf1rly5ddu/ePWfOHMbYvn372rZt6+rqalkhOjr6hRdeKCwsrF27dm5u7unTpzt16lSR+83Jyblz506DBg3QEQIAQNU1adKkxYsXz5s3r2nTpnPmzJEveGCMde3addSoUTNmzGjSpEnfvn1HjRo1efLkVatWDRs2TD6nmJmZuXnz5qNHj964cePDDz8MCwsbPnz4/v37N2zY0KFDh+Li4uXLl/fv379Zs2YIywAAgHKUDsv4+Pj8/vvvpaWlycnJq1evHjVqlFz+9NNPd+7cWb69YcOGxx57bMeOHQMGDFi7dq3lv7dp0+bFF180/xoeHh4cHLx///60tLTZs2cjLAMAAAqTrLig7mHVw8LClixZ8kChPOunzM3Nbfbs2Q9UaNiwoXxA1ZKvr+/cuXMfKERH+Cio1VlFFuY1CswERs4FpUSKhJp7jIoS8FsoUpkqF7pHciOiU6wRyQgTd9lbwQnMuDTEh4WRiKhoTQIfLtQHEZXEMXBnRzPwKwu9EIrsV0x0ZxYpp+6Rem/y38vCu0S1TdFYseK86s8FDo0CAIBDw4gQAACUY4NJt20NHSEAACjHiusI1b7cEodGAQDAoWFECAAAirEiLINDo/ZMIDRKLpRq5M8nxZ/eiVzFl5wjipeRIzdCtJAo59anNkLeKTeTKXKPVDk1PxbVEupOhZa9pZKD3M8FKgRKTdZFkSSBF0JDbF2rVeLVFNklhF9l7q7C28OZYDaYnEqNeG9yX30yNCo07R6oBB0hAAAoRlJ6PcJKgHOEAADg0DAiBAAAxVhzQb3aB4rREQIAgGKsWYZJ7Z4Qh0YBAMChYURYSahQGZkm5YUVhRaPpeoLb0SkXKglZAtt2RLu3KGiG1doYV6iNrURKmDMC4Kq82oKbUS0XJGdWWR/E80Gw1+smFlGbegIAQBAObiOEAAAHJnEJNHLIdTuB3GOEAAAHBtGhAAAoBysPlE9UTkXofrkiXeBE/JkZf58UoIbUaJcnY2I5FmoNXXFNi4ao+C+z6lXjfhQIB4Q0wjFear+qym2cW7dct4RIhuh/8ArJOryixlvXjz5D0S5/ZDErwtU/UHj0CgAADg0jAgBAEAxkiR+gbzaF9SjIwQAAOXgHCEAADgya+YatU1LKg7nCAEAwKFhRKg8fmiUqizyBzLFJrRxwRVExVYcFV2elJ9rFdzIo9Yt705FIsCCGye+BgtvXInaiuwqQruE8KsssqsIviPEMsBCT5UDwqTbAAAAdgYdIQAAODQcGgUAACVh0m0AAHBc1pwjtFFTKgwdofL4i89RlUX+QFUW2zhZW2QjVH3BjXPLxe5Rgbrl3Sm/WDDnwt248C5B1VeitiK7itAuIfwqi+wqgu8I4kUW2bjqH+VViB1eR4hzhAAA4NAwIgQAAMXY4+UT6AgBAEAxVswsozocGgUAAIeGjhAAABwaDo1WABUeE1pak0ymCSTWyMrE9xmxjShRrs5GNLyN8AoZYxpiI9RyuPznllj3ldxV+FsWe5hUy4Uevh28mmIb59Yt5x0hshH6D7xCoi6/WP2cpO3Y4zlCjAgBAMChYUQIAABKsrd1edERAgCAguzwgvqKdoS7du1KSEjIzs728fGZMGFCt27d5PI7d+588MEH586da9my5ezZs93d3W3WVAAAqOqsOEeouoqeI8zIyIiKipo8eXJERES/fv32798vl48YMeL8+fPPPffckSNHJk6caLN2AgAA2ERFR4TTpk2Tb8TFxR0+fDgxMbF79+4nT548cODAzZs3a9as2alTJz8/v4yMjEaNGtmstXZMNJmm4aX+NFqBylR94Y2IlAu1hGyhLVtipFoitOqvInONirz0jDFJ5IVT5jkUedVs2hLqTpV5RwgnbLnF8BcrLqhX/SkVS42WlZWdOHHi0KFDPXr0YIwdPny4Xbt2NWvWZIx5eXk1b948NTXVJs0EAAC7IFn1oyqBjnDlypWenp6tWrUaPHhwTEwMYyw7O7tu3brmCsyPiUYAAB+qSURBVD4+PtnZ2dz/NZlMZWVlj9hWAABQkcFgULsJNiHQEU6dOvXu3bsZGRn79u375JNPGGPu7u46nc5cobi4mArLSJKk0eCaRQAAO+bk9PCzaRJj0l/HRwV+KqHx5RDunBo2bDhu3LikpCTGWGBg4KVLl8x/unz5clBQEPWPdpcjAgAAURKzt26w4mGZq1evyp2cTqfbtWtXy5YtGWOxsbFTpkw5ePBg586dk5KSSkpKevbsacPGVjVCEzBRyQgtfyPcs/paIhqgdeJ/odFoOeX0RogWEuXc+tRGyDvlhheojZQR+Rcjp5xbyBjTUkEXgsQ7om8USdZQyCnTRJ4rRu0qIq8aI55zoVeN2rgi+xV1p9w9nNHvCGIjYu9NoXna1P+Mr3zWXUeoxHvKahXtCHv16lWjRg1vb++zZ89GRUW9+eabjLFatWp9+umnAwcObNWq1fHjx7/44gtXV1dbthYAAEBhFe0I09PT09PTCwoKgoKCAgMDzeVTpkwZOHDgxYsXIyIivL29bdNIAACwE1ZcUC9J9jEi1Gq1zZo14/6pfv369evXV65JAABgr6y5jtA2Lak4JDkBAMChYdJtAABQjFXrEdqoLRWFjvBRCKzFSa4TS1xeyY29aZ2JykS5E6+cW0jdI2PMiQrgcVsoUpkqd3Liny4wlhHlvBAnNdsZRSLWWTby5l7TELFRsaWalZhjjBEvqOgLwX2VFXk1FdmvmOjOLFJOhq6J9yb/hUNs1J6hIwQAAMVYcY5Q9W8LOEcIAAAODSNCAABQjMTEzxGqPSRERwgAAMqxwxXqcWgUAAAcGkaEyuMeFSCSgHQUkBunJLJwziLlTs78KRSdXfjlTs789bO4jaFaUkaVGzjlQulQxphJJCFKHbMpM4hEUpWYbFQ0NUpPHyqwqwiVK7S/ibWE3g855UItoe6Uyq+Sc5AKhUYdjz0uzIuOEAAAlGPFFGtqHxtFRwgAAMrBOUIAAADF3bt3Lzs7u5wKer0+KyvLYDBUZGtlZWVZWVl6vV7+FR0hAAAoRrLqp3xz584NDAxs3759ly5dbt269fcKv/zyS2BgYI8ePYKDg/fs2SMXJiUlPfbYY56enh07drSs/N///jckJKRHjx7+/v4JCQkMh0YfhdBonkpGkGuccqMoRI6ALHfl5QtcqJACUc7bCGPMWW/8e2GZgVPIGCsj8i9l/CgKt67YrGnkSshkREVo/jaiKSJzrIm2UJGwjNCrL7RfUeVClYVbSG5EoJycj414b3Lfy5hgzUyi923yX8qtnpKSsm7durNnz/r6+o4ZM+add9754osvLCuUlpZOnjx59erVgwYN2rRp0+TJkzMyMrRaraen54wZM86fP//DDz+YK5tMpilTpixYsODpp5/evXv38OHD+/btixEhAABUXd99993IkSN9fX0ZYzNmzPjuu+8e+CaamJjo5uY2aNAgxtiIESNKS0tTUlIYY+3btx86dKifn59l5cOHD+fk5EyYMIExFhMTExAQ8PPPP6MjBAAA5Sh9bDQjIyMiIkK+HRERcefOnby8PMsKmZmZ5goajaZx48aZmZnU1jIzMxs1auTs7GzeYGZmJg6NAgCAYqyZdJuxmzdvJiUlWZaEh4eHhIQwxu7evVuzZk250N3dnTF2584db29vc83CwkI3Nzfzr+7u7nfu3KHuiFsZHSEAACjGqvUIpfPnz3/44YeWhU888cTMmTMZY/Xr1y8oKJAL8/PzGWMNGjSwrGlZQa7zQIXyK7dv3x4dIQAAqKxr167x8fHcP0VFRR06dEi+fejQobCwMHlcaNayZctjx47p9XpnZ+fi4uJTp05FRUVRd9SiRYsLFy4UFBR4enoajcbU1NTXXnsNHeEjIL/08EJlxNlYoSnWqECdSw1+Ro5b7lKD/6Lrdfyp1PSl/BCnQYnUKH/WNOE1dTmFVPbSQMV0tQqkRqkwKbeFwrlWIjUqNhsflcnk7RXUfiUUBHUhKguXi7SQbDnv4YtOsSa2MK+DxkbF/4X29NNPt2rVasOGDU2aNHnzzTenT58ul0+cODE2NnbMmDGdOnWKiIh47bXXpk2b9tlnn3Xo0KFFixaMsZycnJSUlEOHDuXl5W3ZssXPz69bt25hYWExMTEzZsyYM2fOt99+6+Pj06tXL4RlAABAQfJZQrGfcjbXsGHD+Pj4tWvXPv/880899dTLL78slwcEBHh6esq34+Pj8/Lyxo8fr9PpNm3aJBfm5ORs2bIlJyenbdu2W7Zs+e233+Tyb7/91tXVdcKECZcvX96+fbskSRgRAgBAlRYTExMTE/NA4aJFi8y3AwIC1q9f/0CFqKiozZs3/31rPj4+q1evtixBRwgAAIqRqsBqEqLQEQIAgHLscNJtdISVhPqKRIZlnDnl1MRRQvkCVzcqFMPfGQxU/oVXbizjV6YW8BNaSpDCDSJRT6zWwC8n4zy8cttOsUbGeRSZYk0oV8Wv7CpSTlZ2I1rixt8PufWp5Bf1juA+fO57jZUTllH7U7uKs+7yCRs1poIQlgEAAIeGESEAACjGmhXqbdOSikNHCAAAyrHDc4Q4NAoAAA4NI0IAAFCMFWEZ1UeE6AiVJ7QP0KuwKjDFWo1STrmBSocSU6kJzZrGnzLNxmvqcuOUeicFHg6z29SoM5UaFZkdjQ58UsFOXkq5pkBlxlgNkfIaVPRUZIo17hPI6P2NmDCPW9cR2eGRURwaBQAAx4YRIQAAKEr1IZ4gdIQAAKAcO7ygHh0hAAAoxprrCNUeQeIcIQAAODSMCCsNtcgnv7ZGy5tAkpwokp9W5K6dyy1kdJySG5tk1PShSqypy1/4tJzYJK9cK/gwywwKLMyrTGqUDBJX9dSom7szp5BIjbrV5FRmjLm5E6lR3naollBzkHLfQdz3GqP3Q/s7A1bJrIiNqg0dIQAAKEaSxM/5qX1sFIdGAQDAoWFECAAAirFiYV61B4ToCAEAQDkSkyTBk4Rq94PoCCuL6Fce7jxbTsRcUNQypGUGzusrNJcYY8zID5cIrqlLrUvMyyNQ4QVqKiz+yrSOHZYhF+al1nAWCsuQ+RdOOTdBwxhzq0WFZfjlNXjlZFiGeJjcd5DoAryqD1+qOjucYw3nCAEAwKFhRAgAAMoRv6BedegIAQBAMVYswyRJkuAVyAqraEeYl5eXmJh45cqV4ODgwYMH16hRQy43mUw//fTTmTNnWrVq1a9fP5u1EwAAwCYqeo6wXbt2GzZsyM3N/fLLL1u3bl1QUCCXv/TSS/Pnz79///7MmTPnzZtns3YCAIA9kKz6UVVFR4SHDx/28fFhjJWVlbVo0eLHH3+cMmVKVlbWV199deHChYCAgPHjx0dFRb322mve3t62bHD1QyXWOIcKqHgbFRF05U00xZ8arZxyJdbUFZo1jUqHUg9TzyunZpIzUDPJKZEapZ4rYilXsdSohkiNcpOQyqRGiWVvySnWeLOmUenQmlRqlCrnTrFGTaVGPHzuO4g+jKf2x7N9smbSbeHJGRVW0RGh3AsyxrRarYuLi1arZYzt3bu3ZcuWAQEBjLHGjRuHhobu37/fRg0FAICqT/qrKxT4Uf07h/DlE5s3b759+/bgwYMZY9evX/fz8zP/ydfX9/r169z/MplMRup6NAAAsAfV9WNcLDWanJz84osv/uc///H09GRy1MfiYJDJZCKPjNldnBYAAKxjbxfUC3SEBw4cGDly5MaNG6Ojo+USPz+/GzdumCtkZ2f7+/tT/67R4OJ9AAA7VpGP8eq8MO+ff/755JNPrl27tnfv3ubC3r17nzp16tq1a4yx9PT0K1eu9OjRwybNBAAAsI2KjggHDBjg7u6+bt26devWMcaGDh06evRof3//qVOn9u3bd+jQoZs3b37ttde8vLxs2dpqSCixpuGH+JiTif+FxuTKTTzyX3ShdCgjWk59X6Tyrvx5Mp0M3MrOLkRqlJd41JeWcSs7yFyjzsQaztRzyE2NuhDp0Brkwrw2nGuUG1V1diUCxkTwmEiNis01CuWz5oJ6GzWlwiraEX755ZdlZf//k6VZs2byjaVLl+7cuTMtLW3lypWWg0UAAHBEdjjpdkU7wieffJL6U2xsbGxsrELtAQAAe2bFXKNqd4QIsAAAgEPDpNsAAKAYSfx6OdGFfBWHjrCKos7fc0vJEI0SI35yeVL+mroCgQ7GmBMvu6Gjgh4l/BCNvpSTfzHo+WEZauo1ofWKbRuWoZ5DkQn2nJwFwzK8ucqoKdaosEwN3ixo3AV1GTFlGitnrV3eKsFCU6kx4jlHKAZwaBQAABwaRoQAAKAY69YjtFFjKggdIQAAKMaKmWXUPkWIjhAAABSk/mISwnCOEAAAHBpGhPaEPOBAfJ/RKpMaFVhrVyjZyBhz4gVEXVz56dDSEn6IUa/jBET11MK8RLmxjCrnpka5dUUX5uVXphKPGq3AesXOogvz8jKZ3CgpY6yGyIK9VAqUWlOXnDWN94i0xHMiETu+2qeiHIJV5wht1JaKQkcIAADKqcarTwAAAFRL6AgBAMCh4dAoAAAoRmJWLMOE6wjhkZG7kYYT3qASNNRGhGYCI6dSo7IbLpycSykvuMEYKy3hz5rGXXqQCsuUUeXEOoX89Qj5dYVmWCMDHfR6hPx/0AqFZYjp6/hhGZFkDSPCNdRGyPyU0FKCVChG7Q9WR2aP1xHi0CgAADg0jAgBAEA51XhhXgAAgIeyx7lGcWgUAAAcGkaEAACgGCvCMmofGUVHWD1Qa+dy/0CFFYmNOBHJPG64kZoezMmJX84NMbrwUqCMMb0bP6zJTY2SC/DacmFewSnWlFmYl5sapWO6Agv2ClVmxIR5QinQcsqxpq49sbdzhDg0CgAADg0jQgAAUJLdXceJjhAAABRjzTlCtftNdIQAAKAYay6fUHsEiXOEAADg0DAidDj0dzUqxUfU5qVGuYWsvAV7ORlOKqxoIKYD5QZBycrUnKJUapQ71yg1qahIbJRcmJd4Dsk0KS+WSWU1uRFTqr6Ts9isp9ylg+kUKLeYXgha7aNnUFGYWQYAABwZJt0GAABQ2Pr165s0aeLr6ztt2jSdTvf3CmlpaT179vT29o6JiUlPTzeXv//++yEhIUFBQQsWLDBf/hsTE9P+/7zxxhsMI0IAAFCQ4usRpqWlTZ8+ffv27U2bNh02bNiiRYsWLFhgWcFkMg0bNmzSpEnbt29fsmTJqFGjjh49yhhLSEhYvnx5UlKSi4vL448/Hh4ePmbMGMbY8ePH16xZExQUxBjz8vJiGBECAEBVtnbt2mHDhnXr1s3Hx+cf//jHmjVrHqiQkpKSl5f3+uuv16pVa968eZmZmYcPH2aMrVmzZvr06REREaGhoS+//LLlP0ZGRrZr165du3aNGjViGBGCmWiIRpI4wRBy2jCNwDxbWid+5MS5jL8R7uxodPiFCstwi/mzqdl0ijV66jX+xrnPreg8bdz65CxoRJyHn58S3q+o+mAfFL+O8MyZM7GxsfLtVq1aXb9+vbCwsHbt2pYVWrZsqdVqGWMuLi7Nmzc/e/Zshw4dzpw58+KLL5r/cdGiReZ/GT58uCRJnTp1euedd/z8/NARAgCAykpLS/Pz8y1L3N3dXVxcGGO5ubnmbq9OnTqMsdu3b1t2hHl5eR4eHuZf69Spc/v2bbnc8h/lQsbY0qVL27RpU1JS8v7778fGxh45cgQdIQAAKEf8gnomscTExLCwMMuyZ555ZvHixYyxunXr3r17Vy4sLCxkjHl7e1vWrFu37r1798y/3rlzR67wwD+a/2vcuHHyjfXr13t7e584cQIdIQAAKMeK6wgZe+KJJ+Lj47l/Cg8PT0tLk2+npaXVr19fHheaNW7c+PTp0yaTSZIkg8Fw7ty58PBwuTwtLU0+rJqWlta4ceMHtuzi4uLs7KzT6RCWAQCAqmvSpElbtmxJS0u7f//+hx9+OGnSJLn8vffe27VrF2Pssccec3V1XblyZVlZ2bJlyxo0aNC5c2fG2OTJk5cvX56dnX379u1ly5ZNnjyZMXbx4sUDBw7odLrCwsLZs2fXrl27devW6AgBAEAxclhG9Kccbdu2XbhwYZ8+ffz8/Dw9Pd966y25/NSpU9nZ2YwxrVa7efPmFStW1K5d+7vvvtu4caN8bHbUqFHDhg2LjIyMiIjo27ev3IPeuXPnueeeq1OnTnBw8KlTp7Zv316zZk2Jyr8pq7i4ePq7iwdOeq4S7gtURO9NAjFLaiMm3mxnVDlRl96IUBCUqsy/T7F8pPDcY9wVkqmNKBH4FAqCIgVanTzZxO+hdf69L+PWnRKhzZ76Y3f+mT3UodFKgHOEAACgKHv79oNDowAA4NAwIgQAAMVYFRpVGTpCAABQjDUL86p9JhmHRgEAwKFhRAhKEp1YUmwjRLmJO0+m4Nq5YulpKpJKVCeSmiL3WN7TIhL4FN76I9YFx4OFeQEAwKFJwoc6VT80io4QAAAUY83qE7ZpScWJnSPMzs6Wr+S3dOPGjX379uXk5CjXKgAAgEpS0Y5w1apVfn5+/v7+zzzzzAPlLVu2/OCDDyIjI7///nsbtBAAAOyKJP6jqooeGu3WrduePXt+/vnn5ORkc+Hdu3dnzZqVlJTUsWPHPXv2yBO7ubq62qapUA2Jnkqw2ZbpoMuj1iUJv/cVCSIpsA2A8sizh6rdCjEVHRE2b968WbNmD6yFvXPnzqCgoI4dOzLGevfuXbNmzX379ineRAAAANt5pLDM1atXQ0NDzb+GhIRcuXKFqlw5s3sDAICKrAnLqD2AfKSOsKSkxMXFxfyrq6trcXExt6bJZDIajY9yXwAAoC6DweDk9LBeowqc8xP1SB2hr69vbm6u+dfbt2/7+fEX6ZAkSavVPsp9AQCAuh7eCzrgFGsdO3Y8cuRIUVERYywvL+/MmTMdOnRQqGEAAACVoaIjwrS0tJ9//nn//v0XL1788MMPo6Ki+vfvHxkZ2aNHj7Fjx06ZMmX58uWDBg2yPGUIoDDRacNsthVJndgogB2w4hyh6u+Fio4IS0tL8/PzW7RoMXjw4Pz8fHkUyBjbvHlzmzZt1q9f36NHj3Xr1tmsnQAAYA/s7SJCVvERYZs2bdq0afP3cg8Pj7ffflvRJgEAAFQezDUKAAAKsr+wDDpCAABQjB2uwoSOEMAKqr9xAaosO+wJsUI9AAA4NIwIAQBAMVZcUK86dIQAAKAYe5xrFIdGAQDAoWFECAAAipEk4cshVD+Uio4QAAAUpfahTlE4NAoAAA4NI0IAAFCSvQ0I0RECAIBy7HE9QnSEAACgHMwsAwAAYF8wIgQAAMXY4wX16AgBAEAx9jjFGg6NAgCAQ8OIEAAAlGOHYRl0hAAAoBhJ/Jyf6kdS0RECAIBiJCZJgkM80fqKwzlCAABwaBgRAgCAcqw4R6g2dIQAAKAc8esIVe84cWgUAAAcGkaEAACgGEy6DQAAjg3nCAEAwJFJzIoRoY3aUlE4RwgAAA4NI0IAAFCMHc6who4QAAAUZIc9IQ6NAgCAQ8OIEAAAFCNJVlwOgcsnAACg2sAK9QAA4NBwjhAAAMC+YEQIAACKseaCerWHhOgIAQBAMRJWnwAAALAv6AgBAMCh4dAoAAAoxqplmGzUlopCRwgAAIqx4hyh2v0gDo0CAIBjw4gQAACUY4cX1KMjBAAA5YifI1T9JCEOjQIAgEPDiBAAABQjiQ/w1B4QoiMEAADlWHP5hNonCSupI9RoNP9esjB587eVc3eV7MqVKwEBAVqtVu2G2NaVK1cCAwM1mmp+ON1BHubly5eDgoIc4WEGBweLL49nZyrtYR4dM+a9994rv05MqI/oZt0yfC+5ulrbKAVIJpOpcu7p8uXLZWVllXNflUyn07mq+ipWDjzM6gQPszqptIfp5+fn5uam+GZNJlNpaamKr1TldYQAAABVUDU/MAIAAFA+dIQAAODQ0BECAIBDQ0cIAAAODdcRCsvJyUlNTc3KyoqJiQkLCzOXX758+ZtvvikqKhoxYkSHDh1UbKEi/vzzz127dt26datZs2Zjx441R8UKCwtXr16dlZXVq1evQYMGqdvIR5eUlHTgwIGCgoLg4ODx48d7e3vL5QUFBatXr75x40afPn0GDBigbiMVtHnzZldX18GDB8u/6nS6r7766sKFC23atBk7dqy9X02xbdu27Oxs+XbdunWHDx8u387Ly1uzZk12dna/fv369u2rXgMVc/PmzXXr1l2/fr1hw4YTJ06sU6cOs9hpY2JinnjiCbXbaE/se79XRffu3RctWjRnzpzU1FRzYXZ2docOHQoKCurXr//444+npKSo2MJHV1BQEBcXd+vWreDg4PXr13fv3l2n0zHGjEZj7969Dxw4EBYWNnPmzM8++0ztlj6qTZs2GY3GRo0a/fbbb61bt87Ly2OMGQyGHj16pKamNmrU6Pnnn1+xYoXazVTGtm3bnn322cWLF5tLnnrqqR9++CE8PPyTTz6ZNWuWim1TxOLFi3ft2pWRkZGRkXHt2jW5sLS0tFu3bidOnGjYsOGUKVO++eYbVduogPPnz0dFRaWlpYWGhqanp8uP1GAw9OzZMzU1NSwsbPr06cuXL1e7mXbFBILKyspMJlOrVq02btxoLlywYMGQIUPk24sXLx4wYIA6jVNIWVmZTqeTbxcXF9epUyclJcVkMu3YsSM0NFSv15tMpt27d/v7+8u3qwGj0diwYcP4+HiTyRQfHx8RESG/0PJDlm/btYKCgsjIyAULFnTp0kUuOXXqlLu7e2FhoclkunjxopubW25urqptfFRdu3bdtm3bA4UbNmxo0aKF0Wg0mUzx8fHh4eHybfvVr1+/efPmPVAoPzR5R/3ll19CQkIMBoMarbNLGBEK4x4+SklJMR9yefzxx5OTkyu3UQrTaDQuLi7ybaPRWFpa6uHhwRhLTk5+7LHHnJycGGM9e/bMzc1NT09Xs6HKSU9PLygoaNq0KWMsOTk5JiZGfqFjYmKuXLly6dIlldv3yF555ZVXXnnF39/fXJKSktKpUyf5lW3UqFFAQMChQ4fUa6Ayfvnll48//njHjh2m/7tCOiUlpU+fPvKsK3379j1//vz169dVbeMj0ev1iYmJgwcPXrt27YoVK8wD3wd22mvXrlWDnbbSoCNUxo0bN+rVqyffrl+/flFRUWFhobpNUsrrr7/eo0eP1q1bM8ays7PND1Or1Xp7e9+4cUPV1ilg9uzZAQEBUVFRS5YskTtCy1fTxcXFy8vL3h/m7t27MzMzp0yZYllo+WoyxurXr2/XPQRjrHnz5q6urjdv3nzppZfi4uKMRiP731ezZs2atWrVsutX8+rVq0aj8YUXXrh06dLJkydbtWp1+vRp9r+vprOzczXYaSsTwjLKcHJyMhgM8m35hrOzs6otUsbSpUsTExPNpzydnJws58nT6/XmgaP9evvtt1999dX9+/dPmzatZcuWHTp0cHZ2rk4Ps6io6KWXXvrxxx8fmIuy+r2aq1atkm/MnTs3IiJi165d/fr1s3xvMsYMBoNdP0yNRmMymV544QX5a41er//444+/+uqr6vdqViaMCJUREBBg/jadlZVVt25dW8zIV8n+9a9/ff7553v27PH19ZVLAgICsrKy5NvFxcX5+fmWh9rslLu7u6+v74gRI/r165eQkMD+92HevXv37t27dv0wk5OTs7Kyxo0b1759+4ULF544caJ9+/ZGo9HyYTLGsrKy7PphWvLy8mrevHlmZib73/fm7du3S0pK7Pph+vn5aTSa5s2by79GRkZevnyZ/W2nLSwstOuHWcnQESojLi5u69at8qGYLVu2xMXFqd2iR7VmzZpPPvkkMTExMDDQXBgXF5eYmFhQUMAYi4+Pb9KkieUFJHbHYDDo9Xr5dmlp6YkTJ4KDgxljcXFxv/766927dxljP/zwQ5s2bQICAtRs6KPp2rXrnj17Vq5cuXLlyvHjxzdq1GjlypUajaZfv35Hjx6VP0b/+9//6nS6Ll26qN1Y6+n1evPI7+rVq8eOHYuMjGSMxcXF7dix4/79+4yxH374ITo62sdHeHmEqsPV1bV///4HDx6Ufz148KDcKcbFxe3cuVPeabdu3dq6dWvLdy48hNppHfszffr0du3aubm5NWrUqF27dqmpqSaT6d69e23btu3Ro8fIkSMbNGhw7tw5tZv5SLKysiRJCg4Obvd/tm/fLv9p7NixzZs3nzRpko+Pz08//aRuOx/RlStXGjRoMGTIkLFjx4aEhPTp06e4uFj+0/Dhw1u0aDFx4kQfH5+dO3eq204FrV692pwaNZlM8+bNCw0NnTJliq+v78qVK1Vs2KO7ePGin5/fk08+OXLkSC8vrxdeeEEuNxqNcXFxrVu3njBhgre39969e1VtpgL+/PPP+vXrT5gwYcCAAY0bN75+/bpcPmLECPNO++uvv6rbSPuC1SeEnT9/3jIIExERIefudDrdnj177t2716dPHy8vL/UaqIDS0tKTJ09aloSGhsoXm5tMpv3792dlZXXp0iUkJESlBirmypUrx44dKykpCQ8Pb9OmjbncZDIlJydnZ2d37do1KChIxRYq6/bt27m5uU2aNDGXpKamnj9/vnXr1s2aNVOxYYo4c+bMmTNnjEZjVFRURESEudxoNO7bt+/mzZvdu3e368G9WW5u7p49ezw9Pbt162Y+C2MymVJSUm7cuFHNdtpKgI4QAAAcGs4RAgCAQ0NHCAAADg0dIQAAODR0hAAA4NDQEQIAgENDRwgAAA4NHSEAADg0dIQAAODQ0BECAIBDQ0cIAAAODR0hAAA4tP8H62Kf4MGzGy4AAAAASUVORK5CYII=",
"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.3"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.3",
"language": "julia"
}
},
"nbformat": 4
}