{
"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.099788e+01 6.116984e+00\n",
" * time: 0.004317045211791992\n",
" 1 2.691010e+01 3.922713e+00\n",
" * time: 0.012609004974365234\n",
" 2 2.254508e+01 6.512151e+00\n",
" * time: 0.03277111053466797\n",
" 3 1.335793e+01 1.979179e+00\n",
" * time: 0.05796694755554199\n",
" 4 1.205233e+01 1.203847e+00\n",
" * time: 0.07844901084899902\n",
" 5 1.114881e+01 8.827009e-01\n",
" * time: 0.09951591491699219\n",
" 6 1.076293e+01 1.207446e+00\n",
" * time: 0.11600708961486816\n",
" 7 1.042888e+01 5.836201e-01\n",
" * time: 0.13241910934448242\n",
" 8 1.013778e+01 5.526818e-01\n",
" * time: 0.14917898178100586\n",
" 9 9.955479e+00 7.352275e-01\n",
" * time: 0.16536903381347656\n",
" 10 9.818391e+00 6.283512e-01\n",
" * time: 0.18167710304260254\n",
" 11 9.666933e+00 4.394054e-01\n",
" * time: 0.19826793670654297\n",
" 12 9.540123e+00 6.213487e-01\n",
" * time: 0.21120500564575195\n",
" 13 9.449217e+00 4.928176e-01\n",
" * time: 0.22431206703186035\n",
" 14 9.338689e+00 5.947827e-01\n",
" * time: 0.23704004287719727\n",
" 15 9.209000e+00 5.497079e-01\n",
" * time: 0.2513151168823242\n",
" 16 9.061924e+00 4.757460e-01\n",
" * time: 0.2642250061035156\n",
" 17 8.907827e+00 5.779805e-01\n",
" * time: 0.3481109142303467\n",
" 18 8.815952e+00 4.871571e-01\n",
" * time: 0.3613889217376709\n",
" 19 8.784493e+00 5.918396e-01\n",
" * time: 0.37372899055480957\n",
" 20 8.749661e+00 5.343363e-01\n",
" * time: 0.386944055557251\n",
" 21 8.705717e+00 4.570461e-01\n",
" * time: 0.39934611320495605\n",
" 22 8.655241e+00 4.795930e-01\n",
" * time: 0.41168808937072754\n",
" 23 8.610672e+00 2.811325e-01\n",
" * time: 0.4238431453704834\n",
" 24 8.583581e+00 3.116545e-01\n",
" * time: 0.4360971450805664\n",
" 25 8.564811e+00 1.439107e-01\n",
" * time: 0.44854307174682617\n",
" 26 8.549999e+00 1.710718e-01\n",
" * time: 0.46181201934814453\n",
" 27 8.531525e+00 1.076546e-01\n",
" * time: 0.4742560386657715\n",
" 28 8.515105e+00 1.231673e-01\n",
" * time: 0.4868040084838867\n",
" 29 8.501319e+00 1.256690e-01\n",
" * time: 0.49924397468566895\n",
" 30 8.493459e+00 2.355565e-01\n",
" * time: 0.5119380950927734\n",
" 31 8.483944e+00 1.467024e-01\n",
" * time: 0.524104118347168\n",
" 32 8.473073e+00 1.680108e-01\n",
" * time: 0.5364689826965332\n",
" 33 8.459393e+00 2.017228e-01\n",
" * time: 0.5488741397857666\n",
" 34 8.447167e+00 1.013763e-01\n",
" * time: 0.5626459121704102\n",
" 35 8.441930e+00 1.229784e-01\n",
" * time: 0.5749590396881104\n",
" 36 8.431801e+00 1.112461e-01\n",
" * time: 0.5873849391937256\n",
" 37 8.411176e+00 1.735535e-01\n",
" * time: 0.6002750396728516\n",
" 38 8.385674e+00 1.209362e-01\n",
" * time: 0.6132221221923828\n",
" 39 8.359474e+00 1.099458e-01\n",
" * time: 0.6256780624389648\n",
" 40 8.337680e+00 1.329525e-01\n",
" * time: 0.6381571292877197\n",
" 41 8.314685e+00 1.160096e-01\n",
" * time: 0.6504909992218018\n",
" 42 8.287493e+00 9.857422e-02\n",
" * time: 0.6642551422119141\n",
" 43 8.257994e+00 2.990935e-01\n",
" * time: 0.6774659156799316\n",
" 44 8.235790e+00 1.849781e-01\n",
" * time: 0.6896250247955322\n",
" 45 8.220807e+00 2.764824e-01\n",
" * time: 0.7021269798278809\n",
" 46 8.210301e+00 3.131298e-01\n",
" * time: 0.714954137802124\n",
" 47 8.196624e+00 2.009719e-01\n",
" * time: 0.7275559902191162\n",
" 48 8.187040e+00 1.158681e-01\n",
" * time: 0.7400341033935547\n",
" 49 8.184173e+00 1.363978e-01\n",
" * time: 0.7528059482574463\n",
" 50 8.179223e+00 1.774656e-01\n",
" * time: 0.7675590515136719\n",
" 51 8.177658e+00 1.473323e-01\n",
" * time: 0.7814490795135498\n",
" 52 8.176886e+00 9.391334e-02\n",
" * time: 0.7966580390930176\n",
" 53 8.175235e+00 1.017094e-01\n",
" * time: 0.8118901252746582\n",
" 54 8.174395e+00 3.558007e-02\n",
" * time: 0.8943359851837158\n",
" 55 8.174395e+00 3.243617e-02\n",
" * time: 0.9147469997406006\n",
" 56 8.173836e+00 5.044451e-02\n",
" * time: 0.9272620677947998\n",
" 57 8.173836e+00 5.044451e-02\n",
" * time: 1.0395660400390625\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd1xUV/o/8DMDA2IFQaqogKKIYMHeCxpEEUvUGGPDlo3GGJON+UWjxjWJGqPJpukaSNGYKPYkbuwo6tqjKBZURBRBlC5tZpj5/XF35zvB8yBnvDAO83m/+GM4nLlzp8Dh3Pu5z1Ho9XoGAABgrZTm3gEAAABzwkAIAABWDQMhAABYNQyEAABg1TAQAgCAVcNACAAAVg0DIQAAWDUMhAAAYNUwEAIAgFXDQAgAAFbNtnoeRqfTLVyydP6CBdXzcNVMq9Xa2lbTK2lGeJo1CZ5mTVJtT7OBveqpfa5fzy4u1gptNiXl5oMHl2bOnGnqfj2ravqIlJaWZmgUB1MeVc/DAQCA7Ea29Hhqn7Fjf7148aHghi8PH55nxoEQh0YBAMCq1fyDBgAAUG0UCqZQmHsnBGFGCAAAVg0zQgAAkI1CoVAITgn1ejNPITEQAgCArHBoFAAAwIJgRggAAPKxwLAMBkIAAJCNCecImbnPEeLQKAAAWDUMhAAAYNVwaBQAAGRjwgX15j4yioEQAADkY8I5QrOHa3BoFAAArBpmhAAAIB+F5V1Qj4EQAABkoxA/1Gn2Q6MYCAEAQDamXEdo7pEQ5wgBAMCqYUYIAADyMeEcIQ6NAgBATWJp46DgQHjhwoUbN240atSoa9eutWrVkhovXryYlJTUunXrwMDAKthDAACAKlTZc4RlZWUTJkwYNmzY1q1blyxZsnHjRqn9448/HjJkyO+//x4aGvrFF19U2X4CAIAFUCj+m5cRYO45YWVnhN98882lS5cSExPr1avHGNPr9YyxrKysZcuWnT17NiAg4OzZs6GhoVOmTKlbt24V7i8AADzPLPAcYWVnhBs3bnzjjTeysrJOnjxZWFgopWP379/fokWLgIAAxljHjh2dnZ2PHDlShTsLAADPNxMmhE+93KKsrOzQoUM7d+7Mzs6m+iQmJsbGxl6/ft3QotfrL1y4sH379hMnTuh0OkN7Xl7e/v37d+/efe/ePamlsgPhzZs3f/nll7Fjxy5cuLBly5YXLlxgjKWlpXl7exv6eHt7G7b7JGkSCQAAFsosf8Y1Gk1oaOg777yzcePGVq1aJSQkPNln5cqVAwcO3LVrV58+fb788kupccCAAS+99NJPP/00bdq0Ll265OfnM8aOHz/epEmT5cuXx8TEBAYGfvPNN6zyh0aLiorq1av3xx9/KBSKBQsWvP322wcOHNBoNDY2NoY+tra2Go2Ge3e9Xo+BEADAoul0OuO/+VwmrD5RsZ07d2ZmZv755592dnaLFi364IMPtm3bZtwhJydn6dKlp06dCgwMPHfu3IABA6ZMmVKnTp3PP/88KCiIMabVajt27BgTEzN37lw/P78bN264uroyxvbs2fPiiy9Onz69sjNCT0/PAQMGSBPYgQMHXrx4kTHm4eHx8OFDQ5/MzEwPDw/u3RUKhVKJi/cBACzYU0dBxv53jlD0i7Zr164RI0bY2dkxxsaOHfvbb7+VlZUZd9i/f7+vr6902UJISEijRo0OHz7MGJNGQcaYra2tt7d3YWEhY8zd3V0aBRljPj4+arVarVZXdnDq27fv7du3pdu3b9/29PRkjPXs2fP8+fM5OTmMsfT09KSkpG7dulVygwAAAE917949wzk4b29vtVqdmZlJdWC8k3Tnz5+Pj48fM2ZMuS2vWLFi9OjRtWvXruyh0Xnz5vXp08fR0bF+/foffvjhmjVrGGN+fn6RkZHDhw8fP378d999N2HCBGmABAAA66RgwpdDKBhLSkpasWKFcWPHjh0HDBjAGDM+B6dSqRhjarXauOeTJ+mMO6Smpo4aNWrVqlUtWrQwvtenn376n//859ixY6zy5whbt24dHx+/cePGx48f79ixo2vXrlL7hg0boqOjExMTp0+fPmnSpEpuDQAAaibxc4QKBdNoNNLBRYOioiLphoeHx6NHj6TbmZmZCoWi3Dm4J0/SGaZk9+/fDw0NnTNnzrRp04zv8tVXX3399ddxcXGNGjViQpVlWrVqtWzZsnKNKpXq1VdfrfxGAACgJjNpPcLAwMDly5dzf9SzZ8/ffvvt3XffZYwdPHiwS5cu0vlCjUajVCptbGx69Ogxc+bMnJwcJyenBw8eXLt2TTpJl5mZOXDgwMmTJ7/55pvGG4yJiVm5cmVcXJzhgCpqjQIAwPNr0qRJK1eunDNnTqtWrRYvXrx+/XqpvWfPnuPGjZOCoJGRkSNGjJgwYUJMTMzLL7/s5eXFGBsxYsTjx4/z8/OlQbRTp06jRo2Kj4+fPn36qFGj1q1bJ23nrbfewkAIAACykS6oF7xPRf2dnJxOnz69fv36pKSk2NjYvn37Su1vvvlmy5Ytpds//vhjdHT0+fPnJ02aFBUVJTVOmDAhLy/PsJ3atWszxho1avTRRx8Zb9/GxkZRPZf3FRcXz1q6cujkGdXwWAAAUBVGtuRfIGesV99fLic+EtqsRp0wsP/DHTt2mLpfzwrX9gEAgFXDoVEAAJCPSWEZ88JACAAAsqlMEe0n71JFO1NJGAgBAEA2FjghxDlCAACwbpgRAlSIl6oWTVoL15sCsFwWuDAvBkIAAJCPCecIq2hPKg2HRgEAwKphRggAALJRPKVQDO8u5p4SYiAEAADZKBTil0OYeyTEQAiWg8ioiEVXBGsKylKBUJaHJP9U8H8g0lrRDwAEWWBYBucIAQDAqmFGCAAAsrHACSEGQgAAkI8pJdbMPRTi0CgAAFg1zAgBAEA+FnhsFAMhmBMR4SRaqfAltwoatWmRjTDG9EKRT6F4KPHLTx4mopp57dShKb3gtkU6m//PGTwPTLiO0OyfHBwaBQAAq4YZIQAAyMcCa41iIAQAAPngHCEAAFgzhcLyao3iHCEAAFg1zAhBTnQmUyDDSQU+dTpiIzqRztSe8DZC9ResV8pH/RdMnV9REP+1KpSc/kpyI/x2JW/j5P/pIrFRREytjYIpRC+QN/sF9RgIAQBAPiacIzQ3HBoFAACrhhkhAADIxoSwjNlnkBgIAQBANgomfh3h0/rfuXNn3bp1eXl5w4cPHzhw4JMdSktL165de+XKleDg4BkzZqhUKsZYUlLS7t27k5OT3d3dJ02a1LRpU6mzVquNiYk5d+6cv7//a6+95uDggIEQnkIo/0J1pqIr3HZdGb9zWRk/0KLTcvqXERsh94Tor+f1JxM33FayUplAboUxprTh/0DJy7/Y2PI3bkNtxIaXuOE1Uo/IqD0XXVDY3JMDeFZynyPMzs7u2rXruHHjgoODx48fv3bt2pEjR5brM2HChEePHk2aNCk6OvrMmTPff/89Y2z69OlBQUHBwcEJCQlBQUFnzpxp2bIlY2zOnDnnzp3729/+tnnz5kOHDv3+++8YCAEA4Pn1ww8/tGnTZvXq1YwxOzu7FStWlBsIb968+euvv6anpzs6Og4ePLhJkybLli1r3LjxgQMHpKkhYywlJWXLli3vv//+w4cPY2Jirl+/3rRp0zFjxri5uSUkJCAsAwAA8lH8d0lCERVt79ixYwMGDJBu9+/f/8yZM6WlpcYdjh8/3r59e0dHR8aYq6trQEDAyZMnGWOGUZAxlpeXJ3U4c+aMt7e3dJi0du3a3bp1O378OAZCAACQjUL8q2IZGRkuLi7SbVdXV71en5GRQXVgjDVq1Oj+/fvGHX788cfU1NRJkyZxO6enp+PQKAAAmNnx48fLpWCGDh36xhtvMMbs7Oy0Wq3UqFarGWO1atUy7mlnZ6fRaAzfqtVq4w579+79+9//vmfPnvr163M729vbYyAEAAD5mFR0u0WLFvPnzzdu8/Pzk254eXndvXtXun337l07OzvjKZ3U4d69e4Zv79275+XlJd0+fPjwxIkTd+zYERISYuiclpam1+ulA7L37t0bMmQIBkL4LyrwSdckE8hqlvGCnYyxMi0nCKrV8NOhVLvQRridWQUpU1676Kq/xNq5gplMKghqyznBYUt0tlXZ8NvteBtR8U+d2PIekTFmIxI9pSrGIU1q6Z5+0u/JuzCFq6traGgo96fDhw9ftGjRokWL7O3tf/nll4iICBsbG8bY8ePHPT09fXx8Bg0aFBUVdenSJSka+ujRo759+0odXnrppV9++aV79+6GrXXv3r2srCwuLq5fv343b968ePFieHg4BkIAAJCNCSvUV9w/MjJy3bp1Xbp08fHx+c9//nPgwAGpfd68eePGjZs7d66jo+PSpUsHDRrUt2/fw4cPf/zxx3Xq1GGMTZ48WafT/f3vf5f6Dx8+fOHChfb29itXrhw7duyAAQOOHj26YMECV1dXBfW/rbyKi4tnLV05dPKMangsMA1mhNx2zAh5G5dlRkjN/DAjfH6NbOnx1D5DRuy4ej1baLNFBX9265i2Y8cOqoNOpztx4kRubm7Pnj2l8CdjLDk5uUGDBs7OztK3SUlJV65cCQoKMhxTvXTpknROUeLi4mK4pv727dsXL15s2bJlQEAAQ2UZAACQUxUszKtUKnv27Fmu0dfX1/hbf39/f39/45agoCBqgz4+Pj4+PoZvMRACAIBsTDhHaPb5Pq4jBAAAq4YZYY3GXfaWqIhJrUxLFeHU8ip/lhFn5jRqqr2M01jKaWSMqYl27sa1Gn5n+txhtZ8jJCt2UucCBc4Rqnjn/BhjKjv+OUKVPafdjtdY0UZETjRyd5vRpVa5/7STC7rihKL5mLD6hNnfLswIAQDAqmFGCAAA8jHhOkJzTwllmBGq1erbt28bF60BAACwFJUdCN9++23jUuFFRUVS++7du728vAYPHuzt7W24zhEAAKyTdI5Q6MvsJwkFDo0uWLBg2bJlxi2lpaVTp0794YcfwsPDt2zZEhUVdfv2ban4DVQrkWvhiYVpmY640lxLtGtKOe1UnkVdouW2l5Zw+pcW8zuTYRleOxXPIS+oJ9q5C/lSwSIqiMTNdFA1xqgr0Kk1dblpFCqiwg3FMCIXY1eL39meWNDbvhanndqIyo7/WlF7ruT9165UyLIWMsjN0l5nsUOjxcXFxt/u27evfv364eHhjLEXX3xRrVYfPXpUzr0DAACoYgID4apVqxo1auTu7v7pp59KLbdv3zZcya9UKps3b56SkiL7LgIAgKUQX5VX/AJ8uVV2IHz11VdzcnIeP368c+fODz/8cNu2bYyxgoICBwcHQ5/atWvn5+dz767X68vK+Me1AADAIlQmFGnCOUJzj4OVHgibN28ujXldu3adMGHCH3/8wRhzdXXNyckx9MnJyXFzc+PeXaFQ4NwhAIBFU6lUT+1jwjho9pOKplw+kZOTIy1y0bZt2z///FMq711YWJiYmNi2bVuZdxAAAKAqVTY1unjx4l69etWrVy8uLi42NvbYsWOMsc6dO/v7+7/55pszZ8787LPPunXrJi1pAVVFcKUkHS/cWMYrjcYY0xIxSyqrWcoLgpYW8QOfJSLtVGq0VCQ1SpZSo9qJl0UsNSpSY406HMRdzIjRNcmqLjXKTYEyIjDMGNPW5lXd0/I3oiOip9RLaMudilB12qhfFe7Mw9zH5Woa80/whFV2RqjX65cvXz537tyEhIS4uDjDsvc7d+4sLi6eNm2aSqXavHlzle0nAABYAgs8R1jZGeHSpUu57Z6enjExMfLtDwAAQLVCrVEAAJCNgjGLqzWKgRAAAORjgecIMRBaEqqCF1U1jRsAoUIx3PALo/Mvxbz24sf8y4yKC/nt/LAMsSdqXj02Rq5HSAQ6iHYdEZYp465HSLziZFhGpMSacFiG125LLBloZ8/fiB0vF0MtDKkh1nrklqkTSiGxCl5DPWcPyb+3VIiGV5KNXNQQTGLKeoTmfgewHiEAAFg1zAgBAEA2ppRMM/eUEDNCAACwahgIAQDAquHQKAAAyMaUsEzV7EnlYSB8TnGjc1RlL3JNXV6ckspkUlXQioggaFGB+slGOjUqUmKN2EMyxMh7mlSykQrNctOhjEiT6ojOZGiUR6nk//qTC/OKlFhTCZZYs3fgLm7M//tArdVcpuU8f+qFpT7MZHE0LgV/D1XEX1YF9xgYcVzM3OetLJWCCZ8jtJhlmAAAAGokDIQAACAfhUlfFdq5c2fnzp1btWo1f/58rZZzxOjWrVvDhg3z8/MbOXJkamqq1Hj16tVFixZFRkbOnj3b0DMxMXHgX504cQIDIQAAyEfuittJSUkTJ05csmTJ77//fvjw4U8++eTJPiNHjmzTps3Ro0d9fHzGjh1ruOPjx48bNWp0+vRpQ08vL6/5/zN06NBjx461atUKAyEAAMhGYZIKNhgdHR0ZGRkeHu7n57dkyZK1a9eW63D8+PF79+4tXbrUy8vro48+unz58oULFxhjkZGRq1ev7tOnj3FnR0fH0P+5ffv2iBEjGjZsiIEQAACeX1euXDEs/BcSEpKamlpQUFCuQ3BwsK2tLWPM3t6+TZs2iYmJT92sWq3etGlTVFQUQ2rU7MjCirwfUOUZqRQfd01dqnYolQ4tzOekQxmRGqU2UkKlRnlr8KoFa42qubVG1URqlFyYl0qNylFrlPcPr5L4L1Rpw/8BWYNUKDVayk+NatScdm4cl1XwGvJSo8I1RanYKHdJXWImQU0wuP2J9C7xkEiTVrvMzMwGDRpItx0dHRljDx8+rFevnqHDw4cPDR2kPpmZmU/d7M6dO2vVqtWvXz+GsAwAAMjItCOjO3fuLNc0d+5caYNOTk6PHz+WbktzwYYNGxo/oqOjY2FhoeHb/Pz8ch24oqOjp06damNjwzAjBAAAsxs+fPiOHTu4P/L19b127Zp0+9q1aw0bNpTmheU66PV6hUKh0+lu3Ljh6+tb8cPdu3fv8OHD//rXv6RvMSMEAADZiEZGFU+7fGLixIlbtmxJSUkpKytbvXr1xIkTpfY1a9YcOXKEMRYaGqrX63/++WfG2HfffVe/fv0ePXowxjQaTU5OTlFRkVarzcnJMUwrGWMxMTH9+vVr2rSp9C1mhAAAICtZT6N27dp13rx57dq1UygU3bp1W7x4sdS+f//+unXr9unTx9bWdtOmTRMnTpwzZ46jo+PGjRuVSiVj7NixY6NGjZI6+/n5vfDCC9Jgqdfrd+3a9e677xoeAgOhuRGRAW4cg7vwKWNMUyqw1i53QV1GV0crKiBCNLx2aiNU/TZuWKaU18gqWJiXlwmiQzFiYRk9r50OgHCb+X8XhEusESEaW149MS1Rj01jJ5B/oUMx1OLGvGCRUN25CnIuvJdL9DXk9lcQSyRTKycTuR2qs9VRKIRLpj21///7f//vnXfeUavVDg4OhsY9e/YYbvfu3TslJeXx48d169Y1NPbr1y87O5v7cOfOnTNuwaFRAAB43tnY2BiPglzGo6AQzAgBAEA+Jqw+Ye75NAZCAACQTyVqhz5vcGgUAACsGmaEAAAgG+maCOE7mRUGwmpCl1Ljt3NDjFSKj1tKjRHV1IoLiXQotaYuWTWN006lQ6kH5S/MK5oa5b0sZaKpUV55MEYtzCuaGuWp0oV5uY2MMRXxNLkvi+j6w/x2kWpnjH5ZuKFZ6jUh23mvLfWIZJQRodEKmbJCvblfPgyEAAAgKwubEOIcIQAAWDfMCAEAQDZPXV+Qc5cq2pVKw0AIAACyMeEcodlHQhwaBQAAq4YZYXUhknZU/lDHS/dRC6WqSwUqeVIBTrK9iGoX2Di1MC9346XFRDpUZK1dKmErWj+T2y6aGuX+g0zGJuVIjars+J2pMK2ujLMwr57fl1o5l088HUo9fU47d1FiRodmue3UC6uw4T9RfrPZg4/PDxMuqDf3i4eBEAAA5GPKOUJcRwgAADWFBU4IcY4QAACsG2aEAAAgGxMunzD7GVYMhPLjJiaoGAUVu9DyMhpUWoRKl3ALmJFr5BLtVMEzbv02ssQakbjh9qcqxnEX4GWMaXkZIuq1ogJH5GKzvDeILqUmsDIvte6rUI0xxpiWlyLRajjhF1ZR5Tmhp8nH/WumJJ4mGYohniY3F6NS8Z+mHZEVUtlx+lPJGnJ1X+7zFMlJ1XyWdmwUh0YBAMCqYUYIAACyMaXodtXsSeVhIAQAANkoFPTCHfRdzAsDIQAAyMcCr5/AOUIAALBqmBFWBW4Aj1r3VSQ1SsQm1SUCgU962VtqI0QklbcdOnpK7DnvGWmIBXipNCl/I1TElKoxJlY1jVqalt/MX8qVOB6kIFKjOqLil66M8+8sWbqPV0qNMUZ0F8MNiNIV4/gfFXJJYV4Q1M6eSI2WUu2cT4XWXqxOm573RtBTGnNPdqqdBU4IMRACAICMLPA6QhwaBQAAq4YZIQAAyMrsxzoFic0INRpNRETEyy+/bGi5evVq//79PTw8wsLCkpOT5d49AACwJNJ1hKJfT5WdnZ2amlpBh5KSkuTk5NLSUhP2WWxG+NFHH6Wmpir/d05cr9ePHDnylVde2bZt26pVq8aOHXvmzBkTdsJSUYEJXruOWNqtrIwoecXLdJAl1oh0CbedjpxQix1SCR2BjZNLCfIKnqmJzkJ7otHwO3MXemSM6ah6Yvx3k0rWiIRoRBfq45VSY4zZqnjhLMH/cYkVE4nOVHE4Xi6GWu1PQ0RR1HYCH2Z74pNP5a00DrywDFWOjvjdtNFx9px6Tci/8JY2Z6o8E2qNPrX/3LlzN2zY4Ojo6Ojo+Pvvv7u7u5frsGvXrqlTp7q5uT169OjHH3984YUXGGN79+794IMPLly40KZNm9OnTxs6u7i4ZGVlSbfHjBmzefNmgd+WS5cu7d69e968eYaWo0ePZmVlvfvuu05OTosWLbp+/fr58+crv0EAAICKHT58ODY29tq1a7du3QoKClqyZEm5DqWlpdOnT9+wYUNiYuI333wzbdo0rVbLGHN1dV24cOGHH3745DaTkpL0er1er9+8eTOr/KFRrVY7ffr0L774QqVSGRqvXr3atm1bGxsbxpi9vX1AQMCVK1dMeqYAAFAjKMS/KrRp06YxY8Y0atSIMfbaa69t2rSp3KGXffv21a1bd/DgwYyxESNGlJWVHT16lDHWvn378PBwFxeXJ7ep0WiMD6JWdiBctWpVt27dunfvbtyYlZVVr149w7eOjo6PHj3i3l2v10tDNAAAWCiNhr+MjDFTThBWOBampKQ0b95cut2iRYuCggLDgU1DhxYtWvzv0RXNmze/fft2xTvZrVs3R0fHDh06nDx5klXyHGFKSsqXX365Z8+e5OTkzMzM0tLS5ORkHx8fJyenx48fG7rl5eU5Oztzt6BQKGxtkVAFALBgxkcE5ZWZmXngwAHjFl9fX19fX8ZYQUGBg4OD1Fi7dm2pxXieZ9xB6lNQUFDBYx09erR169ZlZWUfffTR8OHDb968WanB6dGjR+7u7lFRUYyx7OzsBw8ejBkzJj4+vnnz5levXtXr9QqFQqvV3rhxw8/Pr7LPGwAAahxpmid4F3bjxo0VK1YYN0ZERMyZM4cx5ubmlpubKzVmZ2dLLcY93dzccnJyDN/m5OSU61BO69atGWM2NjYLFy5cvXr1hQsXKjUQduzY8ezZs9LtTZs2ffLJJ9K3/fv3VygU3333XVRU1Ndff92oUaNu3bpVZoM1g9DirHqy5BVRYo2fGuXH2KhyYtx2blCTMaYlV7IVaKc2QhU2425E9GlqeQHRMiodSsV0iTeCmzKl8oRUmpSLSsqJpka5VdOo+Cr9oJxnpFXzO2ts+G+ELS8gquGthcsqeJepzxvv3SeDxCLLMpNV94iPkJ4X02V6/mtFNNfg0KiJNdZ69OixY8cO7g/btm176tQp6fapU6f8/f2leaFBu3bt5s2bp1ar7ezsCgsLL1261K5du8o8bFFRUUlJSd26dYUry9StW9fLy0u6bWtru3nz5k8++aROnTrffvvtzz//LFxZBwAAgDZt2rQ//vjj+++/P3ny5IIFC2bPni21jxs3bsOGDYyxkJCQwMDAOXPm/Pnnn7Nnz+7evXtAQABj7MGDB7GxsadPn87Ozo6NjT127Bhj7Ny5c5999tmxY8cOHDgwcuTI4ODgNm3aCJ+3GzZs2LBhwwzfduvWzXB0VJ4nDQAAFkz8OsIKp5BNmjT57bffPvnkk9zc3KlTpxoGwlatWrm6ukq3t2/f/v7778+aNSs4OPiXX36RGqWBkDHWoUOH2NjYDh069OzZ09HR8eLFi9u2batVq1bXrl3nzZtnK1eABaMgAAAwk1aof+qh1F69evXq1atc4+LFiw233d3d169fX65DcHDwli1byjX6+fl999135RqR5AQAAPlY4DpMWH0CAACsGmaE8uMWnKTyhFS4USg1SsXeuO10gFOsnZs+1RJPh1s6ldoZsjApsRFu4FPHW9mYMaYR2RPGWBlvO0LpUEacO1AoiJwqcVxJpyNyiSL7IlQ+lMqvaog0qY0t541TEbVDNcSaulSqmR8wJvPS1LvMrTUqGCTmvftkTJffXJOZUmvU3K8TBkIAAJCNBR4ZxaFRAACwbpgRAgCAfCxwSoiBEAAAZKMw4YI6DISWjFqdldsocOKdERkNocgJtRGhZA1jrIxY4ZZbZowqYKYloivcB6WeJrmIMS+hQz0iXadNoGoalSLhrkxL9VcQG6FQf1u41fuop0M9KHcPbZT8Uyc2tgLJLzqdRLwRWmKZXO67LLpx7ueNSngRnzfu77KeLrlIMPff/ipjwnWEZn8tcI4QAACsGmaEAAAgH5wjBAAAq2bCdYTmLtKJQ6MAAGDVMCMEAADZVEXR7aqGgfAZiMTE9PwAGrkwL7e8ExmbFFlsls5eCmQyGbHn1NMhnyY3HEuWtuI286OqouFYKr7LDYJSqVGhIzxUkJhcmNeGf/yGiqpyUYlcDTc1asvvbKvj7wn3NSQ/VxqBTz4jPrdkZyoIKrTMskiJNfKvAdVu7j/9YAwDIQAAyEYh93qE1QDnCAEAwKphRggAAGrBFuwAACAASURBVLIx5YJ6cx8oxkAIAACyMWUZJnOPhDg0CgAAVg0zwkogcl9kTIxbipCICHILRTIqgCeYyRQLdorsCdVOlk4V2UPqNRGq10qFFcmspkiZUGoVVupBuVEAopAn01PpUCX/UbnlQ6n0gY7YdX4RTrLYJreZj/yEk1lN4kF5zfTGBaKnVGfyc8jdE+oXnHgjiFWZa0Sa1ITKMuaGgRAAAOSD6wgBAMCaKZhC9HIIc4+DOEcIAADWDTNCAACQD1afAAqVLxAK0ZCdRTZOFnsTjKjw1x8WTNxQG+d3FilTR9Vjo4LaVGEzLioUY6viH2Kp52j/ZKOdPX8F2tJiLbddreYvNst9WRT8bdPJL270iUiRCKVlqK6i4SzuR4XK/pAbFwlnEc3UR5/obH0U4tcFmnscxKFRAACwbpgRAgCAbBQK8Qvkn9a/pKTk119/LSgoCA0NbdKkCbdPfHz81atX27Rp0717d0Njfn7+pUuXVCpV586dDY137949fvx4aWlpx44dAwMDGWaEAAAgJ4VJX7SSkpKePXuuXbv21KlT7dq1O3Xq1JN9FixYMHny5CtXrrz88ssffvih1Lh+/fpGjRqNGjVq9uzZhp5xcXHt2rXbvn17fHx87969P/74Y4YZIQAAyMiUWqMV/jQ2Nlan0+3bt8/GxsbHx+cf//jHb7/9ZtwhMzNz9erViYmJvr6+M2bM6Ny58+zZsxs0aDBq1KhXXnll69atX3zxhaFzUFDQnTt36tatyxgbP3784MGD33rrLcwIAQDg+bVnz55hw4bZ2NgwxkaMGLF3716t9i+BsoMHD7Zq1crX15cx1rp1a29v77i4OMZYw4YNHRwcym3N2dlZGgUZYy4uLjqdTqfTYUZYBaouPybHlvXCJeOENk79gCpAJUQsNMtFrqlLtHNXsqXqsXk0q8dtb9vd48nGBs61uJ1vX8nmtt+8zG8vfqzhtnMpqf/Vuc2Cb4/QNICMOhNxX35/0TA2N+8plA5FPvRpZC+6nZaW1q9fP+m2l5eXVqt98OCBl5eXcQfjbz09PdPS0irzuIsXL544cWKtWrUwEAIAgJklJia+++67xi3dunWLjIxkjOl0OsNIKd3Q/fXqKL1ebzyUKpVKHXX5lJEFCxakpKR89913DOcIAQDA7FQqlZOTk3FLgwYNpBuenp6ZmZnS7QcPHiiVSnd3d+Oe7u7uDx48MHybkZHh6elZ8cMtW7Zs9+7dhw4dkh4FAyEAAMjJhKLb/v7+8+fP5/6wf//+GzduXLBgAWPs3//+d69evVQqFWPs0aNHDg4OderU6du37/Tp0x88eODm5paamnrz5s2ePXtW8GirV6/esGFDXFxco0aNpBYMhAAAIBtTzhFW+NPx48evWrXqlVdeadmy5Zo1azZv3iy1DxkyZNy4cXPnzvX29p4wYcKQIUPGjBnz008/zZw509XVlTF25cqVzz///MaNGykpKTNnzgwKCpo9e/bBgwffeuutiIiIJUuWSNtZunQpBsIqUHX1guTYMlUYnvzoVuXi0fydIR9QYE/IUAjxA6FnaU+cXG/q78Rt7zPM58lGN++63M7H/32H2/4ovYjbzi3JRj5NahFEbgKE/EgItIu+mfTnkNoQPGfkrjVar16906dP//TTTzk5OQcPHmzfvr3UvmTJkmbNmkm3161bFxsbe+XKlUWLFo0cOVJqrF+/fkhISEhIiPStdCW+j4/PunXrjLdvb2+PgRAAAJ5rDRs2fP3118s1Dh482HBbqVSOHTu2XIfGjRvPmDGjXKN0rWG5RgyEAAAgG9kvn6gGGAgBAEA2JlSWMTtUlgEAAKuGgRAAAKwaDo1WE+EQI6/iF1Ufiyobxm2ncoNUjTEl0Z+7cWoPhbKaoolHbrtoOlTojaA2oivj17MoKeIEO4uI0mjqEv4CvNQitMS7LHiShrdxaiNirxX1uRL8qHA3Q+6hyIeZ/FwJ7aGlHQysOpZ4jhAzQgAAsGqYEQIAgJzkXpe3ymEgBAAA+ch9QX01qOxAuG/fvp07d2ZkZLi4uEycONFQyS0vL2/58uXXr18PCgp655136tSpU2W7CgAAzzsTzhGaXWXPESYnJwcHB0+ZMsXf3z8sLCw+Pl5qHz169I0bN2bMmHHu3LlJkyZV2X4CAABUicrOCF999VXpRkRExJkzZ/bv39+rV69Lly6dOHEiMzOzdu3aXbp08fDwSE5OlpYJrlGonCGxQCf3vyGhoB2jMpnEerBku9BGBAN4/Iig4B7yw7EirwnVTkUBKdS/sDa8Pddq+enQ1KQ8bvuh7beebKzbwI7bOS05n9tekFPKbefuOXe3WUUr1nL624i+EbxdITuLfphtOO8otXHq6fM3wmusYE+oX2WB1op+YPFMuKDe7BNIsb8WZWVlCQkJp0+f7t27N2PszJkzISEhtWvXZow5OTm1bt367NmzVbKbAABgERQmfZmVwEC4bt06R0fHtm3bRkZGDhgwgDGWkZHRsGFDQwcXF5eMjAzuffV6fVkZ/+ooAACwCFot54rYGkBgIJw5c2ZBQUFycnJcXNzq1asZY3Xq1Ckt/b/DNcXFxVRYRqFQKKkLswEAwBLY2j79bJqCMcV/j48KfFXDzldAeHDy8fF55ZVXDhw4wBhr3LhxSkqK4Ud37tzx9vam7mhxOSIAABClYJY2DFY+LHP37l1pkCstLd23b19QUBBj7IUXXoiKijp58mTXrl0PHDhQUlLSp0+fKtxZS0ZlN6gT8tyz/TaCZ/WVtjJshModEHsow0Z0ZVRn/p7b2nLay7T8VEgZUQWN+kfNVsVpLyvjbzwrk792btEJTjU1WxX/6Wg1/D3U6/gPyt0OlSLR8bfNrz1mQ+yhDe9zxYiPEPm54r1rjH5ZuA9K7YnQxm2ppyMUE6P+lpv9b3z1M+06QiJ7WD0qOxD27du3Vq1azs7O165dCw4Ofv/99xljdevWXbNmzdChQ9u2bXvx4sWvvvrK3t6+KvcWAABAZpUdCJOSkpKSknJzc729vRs3bmxoj4qKGjp06K1bt/z9/Z2dnatmJwEAwEKYcEG9QmEZM0IbG5uAgADuj1xdXV1dXeXbJQAAsFSmXEdYNXtSeUhyAgCAVUPRbQAAkI1J6xFW0b5UFgbCZ0Au8cprFKwFxY292fASjIyITVLt4mk9op235+SeqGz4G1FxUow6Ih7JLQPGGLO146VGiWCnXvBchFASkqhfxtTFnMuQNaViNcaoN4gKN3LR4VjOxlWCHwnunot+3oQ+h1RnFe8jQfUnP+EiaVL6T7m5/8ZDJWAgBAAA2ZhwjtDs/y3gHCEAAFg1zAgBAEA2CiZ+jtDcU0IMhAAAIJ8avEI9AACAWSQmJv7zn//Mz8+PjIx86aWXnuxQUFDwySefXLt2rU2bNm+//ba0OGBWVlZ8fPyFCxdcXFxmz55t6FxSUrJ69eoLFy40b978nXfecXR0xED4LATW7eSuWcoEI4JU9pIbm6Q2IhSoq6CdW4iS3AjxoCotp11HBT6JTCa3nSqqSaEelEt0EWNi6WBqeVv+g1L9dbznTx2bopKQdvacj5aK18gEY5ZU9FRlR3yYRT6HKuI3gto4t538hFNVefmpUXNPap4bsi/M++DBg969e7/55pstW7acN2+eRqOZMGFCuT5jx461s7ObOnXq2rVrJ06cuHXrVsbY9u3bf/rpJ1tb2/z8fOOBcMaMGWlpaW+88cZPP/0UGRl55MgRDIQAACAfE0qsVXhs9LvvvuvevfvChQsZYxqNZuXKleUGwsTExCNHjmRmZtapU6dXr17u7u63b9/28fGZPn369OnTN2zY8MUXXxg6379/f/PmzSkpKR4eHmFhYW5ubmfPnkVqFAAA5CP3CvWnTp3q3bu3dLt3794JCQklJSXGHU6fPt2+fXtpNVxHR8c2bdqcOXOG2tr58+ebNWvm4eHBGLOzs+vSpcupU6cwEAIAwPMrIyOjYcOG0m0XFxe9Xp+RkWHc4cGDB8ZLPjg7O6enp1NbK9fZxcUlPT0dh0YBAEA2poVG4+LiOnbsaNw4cuTI9957jzFWu3bt0tJSqVGaC0qTPwMHBwdDB6lPuQ5P7YyB0HTUm80tBEaXWJMh50JGA7gJCKIzNy7BGFOTuQNOHIVaDteOiKJwl8kVyq0wxn8nqJMU3BVoGWMaNT9dQ9V7EyOyh9QHS0+EhbjnY6jSfXa1iHffQeCjQudZBBI35J4Q/bk7o7InfiPIjQhk0LjrWjPqd5mKOPGbazKFeHRIoWDt2rVbtWqVcaO0FLx0486dO9LtlJQUBweHckv+NW7c2NCBMXbnzh3jtQLLady48d27d3U6nVKplDpHRkbi0CgAAJiZo6NjyF8ZVvcbNWrUtm3bHj9+zBj74YcfRo4cKY1he/bsuXr1KmPshRdeuH///qlTpxhjcXFxhYWFffv2pR6oe/fuKpVqz549jLGEhISrV68OGTIEM0IAAJCP3BfUDxky5McffwwODm7cuHFqaurBgwel9g8++GDcuHEBAQF169ZdvXr1kCFDOnTocP78+c8//7xWrVqMsWPHjk2aNOnx48d5eXl+fn79+/dfv369ra3tP//5z8mTJ3fo0OHPP/9cvnx5w4YNMRACAIBsTCm6XSGlUhkbG3v16tXc3NwOHTrY29tL7Xv27HFwcJBuT506dejQoUlJSa1atWrUqJHU2KVLl7Nnzxq2o1KppBujRo3q27fvlStXmjdvLsVHMRACAIBsTFqP8On9AwICyrWUO1Po5ubm5uZm3KJSqZycnLhbc3Z27tWrl+FbnCMEAACrhhnhMxDJiVFhRaESa1SgTiiAR3Yu5X8Y7NVl3Hathpca5aVAGWM6Hf9BuZlMPVkdjbO8LaNq2lEFzMiVkPmPWqYVyLVSReC4//BS/wQL1WljxGLIwllNXjuZySTSpPa8jXAbK2inPuTExokPrcjGxUuscRpFM8A1GYpuAwCAdVOILqtk9mWYcGgUAACsGmaEAAAgG8XTVpN4DmEgBAAA+eAcITDGuO8qWfGLym5wS6yROQL++6iuxcm5aEr54ReNmr9xrYa/cW5YhqyORtUp48ZLyMgJv12p5DwjLS/kwujaY9zICWOsjPeMqLprepF6bHTVPSrOI8uqkwK1x+xEQjGMMXsHzkfFvjb/82NHfGjteBsR3Tj1G8H9DaLCMtTvJpH1N/ff8udGFV0+UaVwjhAAAKwaZoQAACAbU1aor5o9qTwMhAAAIB8LPEeIQ6MAAGDVMCMEAADZmBCWMfuMEAOh/IQyZWJ1swRTfGpepk4jUjKNMVZGtfPW4CVrj5GxUR7it0hBlanjvYZKYqFdbsk0RgdBVdxYK1FLTWgNX7IIHJUmJdaJ5aZJhSKmjEiNUtnLWsSC3rV4n7daZAqU/6Gl+gttnLvOMCN+g6gXlnojuB9Pc8cenyMWeGQUh0YBAMC6YUYIAACyMvsUTxAGQgAAkI8FXlCPgRAAAGRjynWE5p5B4hwhAABYNcwIqwvxPw+1YC+34CRVQJIq21iLF/gk06FE4JNs58UsqewlhVixVjROyYtNEiVVtUSaVEcsKSxSDJVE5AzFnia9pDCnnSqdaksEj7kr1nKDmhW0O9RRPdlIlQOtVZvTmTHmQERSuRun9sTOnqg1yvsNohK21Bth/vnLc86E2Ki5YSAEAADZKBTi5/zM/b8FDo0CAIBVw4wQAABkY8LCvOaeEGIgBAAA+SiYQiF4ktDc4yAGwupCvtMiC/ZS9bHsiBJrZVrO+0vWGCPSIuQitESZMS7qPz7u4rSiNca4CQgqWKQRLb3GywqJPHXGGP9dpsMy/G0obYj8C+9loarx0Ws7c9rJnItIFTRuyIUxVrsu0V7PjtvuwOtPr/rLf5rc3yBqJWTyQ2v2P9vPOQussYZzhAAAYNUwIwQAAPmIX1BvdpgRAgCAbBQmqXib33//vZ+fn6ur67Rp00pKSp7skJCQ0L17d2dn5969e1+7ds3Q/o9//MPLy8vLy2vhwoXS2ZyLFy92/Kv4+PjKDoTZ2dmbN2/+5JNPNm/ebLwfer1+9+7dK1as+OOPPyq5KQAAgEq6dOnS3LlzN23adP369Vu3bn344YflOuh0uhdffHHkyJH37t0bNGjQ2LFjpfbt27dHR0fHx8efPHlyy5YtP/30E2PMz89v3f/MnDnz2rVrbdu2rexAGBIS8vPPP2dlZX399dft2rXLzc2V2l9//fUFCxYUFRXNnTv3vffek+mJAwCAZVKY9EWLiYkZNWpUly5dnJycFixYEB0dXa7D0aNH8/Ly5s2b5+DgMH/+/NTU1NOnTzPGoqOjX3vtNV9fX29v7zfeeEO6Y926dUP+59y5c6NHj65fv35lzxGeOXPGxcWFMVZWVtamTZvt27dHRUWlpaVFR0ffvHnTy8trwoQJwcHBb731lrOzswkvXc0nGBvlLkJL1YJS2fFTjDpeuo9MgfJTk3RCktcuXB2Nl4S0VWm5ndWl/KfPrZpGxSa1RDpULDVKvYbcVoJwapR693lJSNHUKLfEGrV2rj2RGnXgVU3jpj0ZY3Xqi6VJuaXXqKWDqcww9zeIG11mDPFQE5lSdLvC353r16+HhYVJt9u2bZuenp6Xl9egQQNDh2vXrgUFBSmVSsaYSqUKCAi4du1a586dr127NmfOHMMdly1bZrzZ4uLizZs37969m1X+HKE0CjLGbGxs7OzsbGxsGGOHDx8OCgry8vJijDVv3rxZs2bx8fGV3CAAANQ8iv8OhQJfjCnUanXOX6nVammDWVlZ9erVk27Xr19fajF+xOzs7Lp16xq+bdCggdQhOzvb+I7l7rV161ZXV9eePXsyE1KjW7ZsefToUWRkJGPs/v37Hh4ehh+5u7vfv3+fey+9Xq/TETMOAACwBDqdTkkdsng2+/fv9/PzM26ZPn36ihUrGGPOzs75+flSY15eHjOamEkaNmxYUFBg+DYvL0/q0LBhQ+M7ljtaGR0dPWXKFOnAjNhAeOTIkdmzZ+/atcvR0ZExplAojK+q1uv11NEes6+7CAAA1UT8gvohQ4bs2LGD+8MWLVokJiZKtxMTE93c3KR5oXGHK1euSIO0Vqu9fv16ixYtGGP+/v6JiYnSYdXExESpUZKcnHzixIlNmzZJ3wqM7SdOnBgzZswvv/zSrVs3qcXDwyM9Pd3QISMjw9PTk7p7Ff0fAQAA1aMyf8aFD4w+7ZzilClTYmNjExISCgsLly9fPmXKFKl9yZIle/fuZYz169evdu3a33zzjVar/fzzz93d3bt06SLd8euvv75//35mZuYXX3wRFRVl2GZ0dHRYWJhhwKrs4HT+/PmRI0fGxMT079/f0Ni/f//Lly/fu3ePMZaUlJSamtq7d+9KbhAAAOCp2rVrt3z58vDwcC8vr0aNGr3//vtS+/Xr1zMzMxljSqUyNjY2Ojrayclp8+bNmzdvlo5Bjh49+qWXXmrXrl3r1q3Dw8MnTpxo2Obly5dfe+01w7eKSlaMdHd3r1OnTkhIiPTtiBEjxo0bxxibO3fuvn37RowYsWXLlnHjxi1dupR79+Li4llLVw6dPMOUl8Eq8deDJcKKVOJRw1uDt7SYn8ksKeS3Fz1WE+2aJxuLiY2UFHI6M8aKizj9S3mNjDE1sdauRs1pJ9cfplKjxPrD3NecXH9YJDZKhRVFF+bl1s8kU6NEnJKbJqUymfTCvLxao3X5tUPr1CPSoURqlPugZGrUnkiN8lLK1BuB8zlPGtnS46l9vtx6JT2rSGizSRfitBnHqUOj1aCy5wi//vrrsrL/+3MTEBAg3fjss8/27t2bmJi4bt0648kiAABYIwssul3ZgXDkyJHUj1544YUXXnhBpv0BAABLZkKtUXMPhAiwAACAVcPqEwAAIBuF+PVyogv5yg4D4XOK/0ESXbKV26oXfNPJ+lO8NXWJPaHW1CWCHvyNqEv4YRluiIZbd43RJda4pdQYkYshI2ZUWEZsYV6xdYm5ZcO4ddcYY7ZEiIa7kq3QAryMWIOXCr84EBuh1trlLh1sS3xUqM8ht3IhQjGAQ6MAAGDVMCMEAADZVGZ9wSfvUkU7U0kYCAEAQDYmrD5h7lOEGAgBAEBGJlxHaG44RwgAAFYNM0JLQoWMlUqRtCIRtGMK/odBKNxILR3MTYcyohIYtXishiixxk2NaqjUqGDpNX5qVHBhXu4rKFpiTTA1Si3My38juLXKqAAnlRrlpkzpFCi/nbtEMCMCojZE2TmqNLTZY/rWwKRzhFW0L5WFgRAAAORjwgr15h4IcWgUAACsGgZCAACwajg0CgAAslEw8XOE5j53i4HQopCfFpHIAPEZpT66VDt3eTwyLEMkdOx4YRluMS1W0XqEnJwLd5FCZkJYhld6TXg9QpESa2QAhHptuWEZkXUHGfGa2xF5FnsHqk4bpz8VfqH2kMpVcaumCYdizH0uyhpY4nWEODQKAABWDTNCAACQTw1emBcAAOCpLLHWKA6NAgCAVcOMEAAAZGNCWMbcR0YxENYI9MeOt3augh9tVBCHB8RKrBGJRyoKaGfHyXaqS4kSa0QQlJsaFU2HCpVYq9LUKDeOyxizoZZfVnH621Il1qj4Li/bSUVMuYXxqI1Tbz0VMKYqzPHX1EU69PlkaecIcWgUAACsGmaEAAAgJ7NfIC8KAyEAAMjGlHOE5h43MRACAIBsTLl8wtwzSAyEAADwvLt161Z+fn5QUJCtLX/YysrKSklJ8fX1dXJyMm6/cuWKTqcLDAw0DM+5ubl6/X+DbXZ2dnXq1MFAWKNxw4p6qngo0UzWIOXWfhRYPJYRkUIqrKjVCJQPJVOjvNqhjKgpygQX5hUivDAvlSblvbaiKyRz2+nAp8ieELtNPX1qMmH2o2dQWXJXltHpdJMmTYqLi/Pw8MjOzj5w4ECzZs3K9dmwYcPcuXMDAwMTExPXrl07evRoxlhhYWF4ePiDBw9sbGzq16+/d+/e+vXrM8aaN2+u1WqVSiVjbPjw4TExMUiNAgCAbKRzhEJfFQ+Ee/fujY+PT0xMPH36dFhY2JIlS8p1KCwsnDNnzs6dO48ePfrzzz/Pnj1brVYzxtavX6/T6S5fvnzp0qX69et/+eWXhrucOXMmOzs7Ozs7JiaG4fIJAAB4nm3ZsmX06NHSZC4qKio2NtZwYFOyb98+V1fXXr16McYGDRpkb28fFxcn3XHSpEm2trZKpXLy5Mlbtmwx3CUnJycjI8PwLQZCAACQjeK/cRkRFU4JU1NTfX19pdu+vr5FRUUPHz6kOjDGfHx87ty58+QdU1NT/7uHCsWoUaPatm3btGnTAwcOMIRlAADA7NLS0mJjY41bWrVqFRQUxBgrKiqyt7eXGmvVqsUYKywsNO5ZVFRkZ2dn+NbBwaGoqIgxVlxcbHxHw70uXrzo6enJGPvqq6/GjBmTnJyMgdD6CBelEil5RSUdlPx0CTdJQWU0ysr47TotZ+N0KEaglBpjTMfrXu7IjNEP+M1iJdao1BIVruGVXqNX9xWo30YVeyP3hNdOhlwQiqmhTLuO8P79++UGwrCwMGkgdHd3z87OlhqzsrKkFuOexh2kPlIHNzc34zsa7iWNgoyxWbNmLVq0KCEhAQMhAACYWadOnYzP4Rnr2LHjsWPH3n77bcbY8ePH27Rp4+DgYNwhJCTk9ddfLy4udnBwyM/Pv3z5ckhIiOGOERER0h07depUbsu5ubmPHz92cnLCQAgAAPIRv6C+4tRoVFTUqlWrPv3001atWs2fP//999+X2sPDw0eMGDF9+vTg4OBu3bpNmTJl2rRpX3/99aBBg5o3b84Ye/3110NDQ9u0aaNSqdasWbNr1y7G2MmTJ3/77bfOnTsXFxd/9tlnPXr0CAwMRFgGAADkoxD/qpCHh8fhw4cvX768fv36JUuWREVFSe0DBw709/eXbsfGxjZt2vTzzz8PCAjYuHGj1NipU6ft27f//vvv27dv/+mnn6RYqZeXV0lJSUxMzI4dO8aOHfv7778rlUoFebZDVsXFxbOWrhw6eUY1PBbIi/yA8H5AdabOwHGvTKdO75VRp/dwjvAJOEcIVWFkS4+n9vnh0K2HeSVCm7186mDe9cM7duwwdb+eFQ6NAgCAbFB0G2ogoVV/qaMcxNSC6XlzCKUNf2JlqyNSo3YCa+dS1dH0/Ikif/JHH0YRmhISXakSeGQil/caksvbCkxDhaud8X5g9r9xUM1MKbpt7k8JBkIAAJCVpf33g7AMAABYNcwIAQBANiYsPmF2GAgBAEA2lniOEIdGAQDAqmFGCKYSCELSPxHZiJ66wI631DB1qZ/odbP87chy8S2ZGhX9b1pkI2JRVVzqB+LkXpi3GmAgBAAA+ShM+GcOl08AAEBNYcoF9VWzJ5Undo4wIyPDeFVfSXp6elxc3IMHD+TbKwAAgGpS2YHwX//6l4eHh6en57Rp08q1BwUFLV++PDAwcNOmTVWwhwAAYFHkrrtd1Sp7aLRnz56HDh367bffjhw5YmgsKCh4++23Dxw40Llz50OHDo0dO3bUqFGGFYEB/o9gdOPZN05uQjQsU+2/o0Il7Sr+gcCDmvsvEdQYCqZQmH1kE1TZGWHr1q0DAgKUf61IuHfvXm9v786dOzPG+vfvX7t27bi4ONl3EQAAoOo8U1jm7t27zZo1M3zbtGnT1NRUqnP1rPcEAABmZHWrT5SUlNjZ2Rm+tbe3Ly4u5vbU6/U67sJuAABgIbRara3t00aN5+Ccn6hnGgjd3d2zsrIM3z569MjDg79so0KhsLGxeZbHAgAA83r6KGiFJdY6d+587ty5wsJCxlh2dvbVq1c7deok044BAABUh8rOCBMTE3/77bf4+Phbt26tWLEiODh48ODBgYGBvXv3Hj9+fFRU1DfffDNs2DDjU4YA8hL9L1OomUqTcrtXZYU1Td/ncwAADl5JREFUGe8AYAYmnCM0+2e7sjNCtVqdk5PTpk2byMjInJwcaRbIGNuyZUv79u03btzYu3fvH374ocr2EwAALIGlXUTIKj8jbN++ffv27Z9sr1ev3uLFi2XdJQAAgOqDWqMAACAjywvLYCAEAADZWOAqTBgIASQiv4tm/70FeH5Z4EiIFeoBAMCqYUYIAACyMeGCerPDQAgAALKpilqjeXl5P//8c35+/pAhQwIDA5/soNfrd+7cef369cDAwIiICEP7/fv3Y2NjdTrdqFGjmjRpIjVevHjx2LFjpaWlnTt37tmzJ8OhUQAAeJ4VFhZ26dLlwIEDeXl5PXv25K5x9Prrry9evFij0cyfP/+dd96RGlNTU9u2bXvt2rWUlJT27dvfvHmTMbZ///6IiIhLly5lZmaOGTNG6owZIQAAyEahEL4couL+mzZtcnJyio2NVSgUrq6uy5Yt69u3r3GHtLS0b7/9Njk52dPT85VXXgkMDHznnXdcXFy++OKL8PDwb775hjGm0WjWrFnz1Vdfde3a9datWyqVijEWGRnZp0+fpUuXYkYIAACykrWyzP79+8PDw6XBcsiQIYcPH9ZoNMYdDh8+HBQU5OnpyRjz8fHx9fWNj4+X7jhkyBCpz9ChQ/ft28cYq1evnjQKMsbs7e2VEpmfPwAAgHzu37/v7u4u3fbw8NDpdBkZGVQHxpi7u3taWlq5dnd39/v37xvfS6fTvfvuu6+++qqdnR0OjQIAgJxMyIxeuHBh5syZxi3Sig6MMaVSaVjXXVrXttyifsYdGGN6vV7qoFAojO9oPPHT6/WzZs0qKSn56KOPGM4RAgCAjExbj9DR0TEkJMS4sWXLltINT09PwxQwPT3dxsbG1dXVuKenp2d6errh2/T0dGllXOP2jIwM6dipZN68eX/++ee+fftq167NMBACAICcTKos06xZsxkzZnB/GBYW9tVXXy1cuFCpVO7evTs0NFRaH/jWrVuOjo7Ozs79+/efOnVqampqkyZNkpKS7ty506dPH+mOu3fvfumllxhju3fvDgsLkzb43nvvxcXFHTx4sH79+lILBkIAAHh+jR07ds2aNUOHDvX39//xxx9//fVXqf3ll18eN27c3Llz3d3dZ82aNWjQoOHDh2/btu3tt992cnJijM2ePbtTp04TJ060s7P79ddfT548yRjbs2fPxx9/PGDAgFdffVXazmeffYaBEAAAZCP7BfUODg4nTpzYtWtXTk7OmTNn/Pz8pPY1a9YYjnauWrUqLCzsypUr69evN1xc4eXllZCQsGvXLp1O9+GHH7q5uTHGQkJC9u/fb7x9R0dHDIQAACCbqiixVqdOnZdffrlcY/fu3Y2/DQ0NDQ0NLdfHxcVl6tSpxi1ubm7SiGgMl08AAIBVw4wQAADkY4HLMGEgBAAA2SgqUUS7/F0wEAIAQI2hYAqF4BRPtL/scI4QAACsGmaEAAAgHxPOEZobBkIAAJCP+HWEZh84cWgUAACsGmaEAAAgG9OKblfRzlQSBkIAAJAPzhECAIA1UzATZoRVtC+VhXOEAABg1TAjBAAA2VhghTUMhAAAICMLHAlxaBQAAKwaZoQAACAbhcKEyyFw+QQAANQYcq9QXw0wEAIAgHxwjhAAAMCyYEYIAACyMeWCenNPCTEQAgCAbBRYfQIAAMCyYCAEAACrhkOjAAAgG5OWYaqifaksDIQAACAbE84RmnscxKFRAACwbpgRAgCAfCzwgnoMhAAAIB/xc4RmP0mIQ6MAAGDVMCMEAADZKMQneOaeEGIgBAAA+Zhy+YS5TxJW00CoVCo3rFp2ZMuP1fNw1Sw1NdXLy8vGxsbcO1K1UlNTGzdurFTW8MPpVvI079y54+3tbQ1Ps0mTJuLL41mYanuaf7788j/+8Y+K+wxo5iK6WYdk9xR7e1N3SgYKvV5fPY90586dsrKy6nmsalZaWmpv1nexeuBp1iR4mjVJtT1NDw8PBwcH2Ter1+vVarUZ36nqGwgBAACeQzX8wAgAAEDFMBACAIBVw0AIAABWDQMhAABYNVxHKOzBgwdnz55NS0sbMGCAn5+fof3OnTvff/99YWHh6NGjO3XqZMY9lMX58+f37dv38OHDgICA8ePHG6Ji+fn569evT0tL69u377Bhw8y7k8/uwIEDJ06cyM3NbdKkyYQJE5ydnaX23Nzc9evXp6enh4aGhoeHm3cnZbRlyxZ7e/vIyEjp29LS0ujo6Js3b7Zv3378+PGWfjXF7t27MzIypNsNGzZ88cUXpdvZ2dnffvttRkZGWFjYoEGDzLeDssnMzPzhhx/u37/v4+MzadKkBg0aMKMP7YABA4YMGWLufbQklv25N4tevXp99NFH8+fPP3v2rKExIyOjU6dOubm5rq6uAwcOPHr0qBn38Nnl5uZGREQ8fPiwSZMmGzdu7NWrV2lpKWNMp9P179//xIkTfn5+c+fO/fzzz829p89q8+bNOp3O19f32LFj7dq1y87OZoxptdrevXufPXvW19f3b3/729q1a829m/LYvXv39OnTV65caWh56aWXtm7d2qJFi9WrV7/99ttm3DdZrFy5ct++fcnJycnJyffu3ZMa1Wp1z549ExISfHx8oqKivv/+e7Puowxu3LgRHBycmJjYrFmzpKQk6Zlqtdo+ffqcPXvWz89v1qxZ33zzjbl306LoQVBZWZler2/btu0vv/xiaPzggw+GDx8u3V65cmV4eLh5dk4mZWVlpaWl0u3i4uIGDRocPXpUr9fv2bOnWbNmGo1Gr9cfPHjQ09NTul0D6HQ6Hx+fHTt26PX6HTt2+Pv7S2+09JSl2xYtNzc3MDDwgw8+6N69u9Ry+fLlOnXq5Ofn6/X6W7duOTg4ZGVlmXUfn1WPHj12795drvHnn39u06aNTqfT6/U7duxo0aKFdNtyhYWFvffee+UapacmfVD//e9/N23aVKvVmmPvLBJmhMK4h4+OHj1qOOQycODAI0eOVO9OyUypVNrZ2Um3dTqdWq2uV68eY+zIkSP9+vWztbVljPXp0ycrKyspKcmcOyqfpKSk3NzcVq1aMcaOHDkyYMAA6Y0eMGBAampqSkqKmffvmb355ptvvvmmp6enoeXo0aNdunSR3llfX18vL6/Tp0+bbwfl8e9///vTTz/ds2eP/n9XSB89ejQ0NFSqujJo0KAbN27cv3/frPv4TDQazf79+yMjI2NiYtauXWuY+Jb70N67d68GfGirDQZCeaSnpzdq1Ei67erqWlhYmJ+fb95dksvf//733r17t2vXjjGWkZFheJo2NjbOzs7p6elm3TsZvPPOO15eXsHBwatWrZIGQuN3087OzsnJydKf5sGDB2/fvh0VFWXcaPxuMsZcXV0teoRgjLVu3dre3j4zM/P111+PiIjQ6XTsr+9m7dq169ata9Hv5t27d3U63WuvvZaSknLp0qW2bdteuXKF/fXdVKlUNeBDW50QlpGHra2tVquVbks3VCqVWfdIHp999tn+/fsNpzxtbW2N6+RpNBrDxNFyLV68eN68efHx8a+++mpQUFCnTp1UKlVNepqFhYWvv/769u3by9WirHnv5r/+9S/pxrvvvuvv779v376wsDDj303GmFarteinqVQq9Xr9a6+9Jv1bo9FoPv300+jo6Jr3blYnzAjl4eXlZfhvOi0trWHDhlVRka+a/fOf//zyyy8PHTrk7u4utXh5eaWlpUm3i4uLc3JyjA+1Wag6deq4u7uPHj06LCxs586d7K9Ps6CgoKCgwKKf5pEjR9LS0l555ZWOHTsuW7YsISGhY8eOOp3O+GkyxtLS0iz6aRpzcnJq3br17du32V9/Nx89elRSUmLRT9PDw0OpVLZu3Vr6NjAw8M6dO+yJD21+fr5FP81qhoFQHhEREdu2bZMOxcTGxkZERJh7j57Vt99+u3r16v379zdu3NjQGBERsX///tzcXMbYjh07WrZsaXwBicXRarUajUa6rVarExISmjRpwhiLiIj4448/CgoKGGNbt25t3769l5eXOXf02fTo0ePQoUPr1q1bt27dhAkTfH19161bp1Qqw8LC/vzzT+nP6PHjx0tLS7t3727unTWdRqMxzPzu3r174cKFwMBAxlhERMSePXuKiooYY1u3bu3WrZuLi/DyCM8Pe3v7wYMHnzx5Uvr25MmT0qAYERGxd+9e6UO7bdu2du3aGf/mwlOYO61jeWbNmhUSEuLg4ODr6xsSEnL27Fm9Xv/48eMOHTr07t17zJgxbm5u169fN/duPpO0tDSFQtGkSZOQ//n999+lH40fP75169aTJ092cXH59ddfzbufzyg1NdXNzW348OHjx49v2rRpaGhocXGx9KMXX3yxTZs2kyZNcnFx2bt3r3n3U0br1683pEb1ev17773XrFmzqKgod3f3devWmXHHnt2tW7c8PDxGjhw5ZswYJyen1157TWrX6XQRERHt2rWbOHGis7Pz4cOHzbqbMjh//ryrq+vEiRPDw8ObN29+//59qX306NGGD+0ff/xh3p20LFh9QtiNGzeMgzD+/v5S7q60tPTQoUOPHz8ODQ11cnIy3w7KQK1WX7p0ybilWbNm0sXmer0+Pj4+LS2te/fuTZs2NdMOyiY1NfXChQslJSUtWrRo3769oV2v1x85ciQjI6NHjx7e3t5m3EN5PXr0KCsrq2XLloaWs2fP3rhxo127dgEBAWbcMVlcvXr16tWrOp0uODjY39/f0K7T6eLi4jIzM3v16mXRk3uDrKysQ4cOOTo69uzZ03AWRq/XHz16ND09vYZ9aKsBBkIAALBqOEcIAABWDQMhAABYNQyEAABg1TAQAgCAVcNACAAAVg0DIQAAWDUMhAAAYNUwEAIAgFXDQAgAAFYNAyEAAFg1DIQAAGDV/j+vuiM+zLcoQgAAAABJRU5ErkJggg==",
"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.2"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.2",
"language": "julia"
}
},
"nbformat": 4
}