{ "cells": [ { "cell_type": "markdown", "source": [ "# Gross-Pitaevskii equation with external magnetic field" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "We solve the 2D Gross-Pitaevskii equation with a magnetic field.\n", "This is similar to the\n", "previous example (Gross-Pitaevskii equation in one dimension),\n", "but with an extra term for the magnetic field.\n", "We reproduce here the results of https://arxiv.org/pdf/1611.02045.pdf Fig. 10" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using DFTK\n", "using StaticArrays\n", "using Plots" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "Unit cell. Having one of the lattice vectors as zero means a 2D system" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "a = 15\n", "lattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Confining scalar potential, and magnetic vector potential" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "pot(x, y, z) = ((x - a/2)^2 + (y - a/2)^2)/2\n", "ω = .6\n", "Apot(x, y, z) = ω * @SVector [y - a/2, -(x - a/2), 0]\n", "Apot(X) = Apot(X...);" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "Parameters" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "Ecut = 20 # Increase this for production\n", "η = 500\n", "C = η/2\n", "α = 2\n", "n_electrons = 1; # Increase this for fun" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "Collect all the terms, build and run the model" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter Function value Gradient norm \n", " 0 3.029827e+01 7.504715e+00\n", " * time: 0.0049610137939453125\n", " 1 2.603528e+01 4.489491e+00\n", " * time: 0.01473093032836914\n", " 2 1.857773e+01 4.669018e+00\n", " * time: 0.036540985107421875\n", " 3 1.584906e+01 4.664210e+00\n", " * time: 0.05862593650817871\n", " 4 1.223901e+01 1.681105e+00\n", " * time: 0.08082199096679688\n", " 5 1.048770e+01 7.141197e-01\n", " * time: 0.10336995124816895\n", " 6 1.006822e+01 9.057448e-01\n", " * time: 0.12119007110595703\n", " 7 9.642703e+00 5.654128e-01\n", " * time: 0.13910794258117676\n", " 8 9.481640e+00 5.290655e-01\n", " * time: 0.15692806243896484\n", " 9 9.446886e+00 1.276598e+00\n", " * time: 0.17073607444763184\n", " 10 9.396194e+00 9.363862e-01\n", " * time: 0.18424487113952637\n", " 11 9.331939e+00 8.767320e-01\n", " * time: 0.19732189178466797\n", " 12 9.320474e+00 6.245026e-01\n", " * time: 0.2104949951171875\n", " 13 9.248494e+00 5.587031e-01\n", " * time: 0.22507596015930176\n", " 14 9.202700e+00 3.513332e-01\n", " * time: 0.23935699462890625\n", " 15 9.164074e+00 3.670895e-01\n", " * time: 0.2534968852996826\n", " 16 9.149964e+00 1.772628e-01\n", " * time: 0.2683548927307129\n", " 17 9.141350e+00 2.215304e-01\n", " * time: 0.34978199005126953\n", " 18 9.140808e+00 3.628566e-01\n", " * time: 0.3627750873565674\n", " 19 9.138445e+00 1.909751e-01\n", " * time: 0.38048887252807617\n", " 20 9.138445e+00 1.909751e-01\n", " * time: 0.5199968814849854\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=1}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVyU1f4H8DPDDAiISqCyuIJCiIiEG67lQqihXk1z92oZ3sxbmW3en9cW29W8LZqVmuW9mWaa5pKaid5r/hLNHUXBJRFcWBSQZYaZ3x9Pv7mE3y/OGR94GObzfvnHeDw+nFngcJ7n83yPzmq1CgAAAFel13oAAAAAWsJECAAALg0TIQAAuDRMhAAA4NIwEQIAgEvDRAgAAC4NEyEAALg0TIQAAODSMBECAIBLw0QIAAAuzVAzX8ZisRw6erxtRETNfLkaZrFY9Pq6/ysFnmZdgqdZl9TY02zoYbxjn9Onc4uLzVKHLSi4ERSkDw0NdXRcd6uGJsLS0tLFX69/6M+Na+bLAQCA6oaHB96xzyOPbDpy5JrkgY8PG3Zj/fr13D8nJye/8847N27cGDZs2MyZM2+f9bOysv7nf/7n+PHj0dHRr7/+euPGjYUQmzZtWrNmTXp6ekBAwLRp0+Lj45XOjz76aEFBgfI4Li7umWeeqfu/KwEAgPO6cOFCYmLiqFGjFixYsGLFig8//PD2PsOHDzcajUuXLjWbzaNHj1Ya169f36dPn4ULFyYkJAwbNmzfvn1K+3fffdezZ8+RI0eOHDkyLi5O1NiKEAAAXIFOJ3Q6NQ+4bNmygQMHTpo0SQjx+uuvP//883/9618rdkhJSTl58mRycrK7u/uHH37o7+9/8uTJdu3aLV++XOnQrVu37du3b9u2rXv37krLwIED27ZtazsCVoQAAFB7HTlypGvXrsrjrl27njlzpqioqFKHjh07uru7CyG8vLyioqKOHDlSsUN5efnJkyfbtGlja3n22Wcffvjhd999t7i4WGBFCAAAKtLpdDrJJaHVWlX/K1euNGrUSHl8zz33KC0hISEVO/j6+tr+6uvre+XKlYpHmDt3rpeX19ixY5W/PvXUU9HR0cXFxYsWLdq6devOnTsxEQIAgKrkT41u3rxZmeRspk6d+vbbbwshGjRocOvWLaWxsLBQCNGwYcOKPSt2EEIUFRU1aNDA9td//OMfq1ev3rNnj8Hw+3w3Z84c5UF8fHxAQMCJEycwEQIAgMYGDBiwatWqii3e3t7Kg1atWqWnpyuP09PTGzRoUGnKbNmy5dmzZ5XHVqs1PT29VatWyl8//vjjRYsW7d69Oygo6PYv6uvrW79+/by8PFwjBAAA9eh+z8tI/XF3d/f9I+WanxBizJgxa9asuXbtmhDio48+GjNmjHLqdcWKFb/88osQIj4+vrCwcOvWrUKI9evXu7m59e7dW+nw5ptv7ty5s2XLlrbRZWVlZWVlKY8XL15ssViio6OxIgQAANU4cI1QVHmN8IEHHhg9enR4eHjDhg0bN268adMmpf2zzz575JFHunTp4uHhsWzZskmTJvn5+eXl5X3xxRfKWdBXX3314sWLtozMY4899umnn6anpz/00EM+Pj5lZWU+Pj5r165t2LAhJkIAAKjVFixYMGfOnKKiouDgYFvjf/7zH9vjxMTEzMzMS5cuNWvWzGj8vfzNuXPnbj9Uz549r1+/np2dXa9ePX9/f6UREyEAANR2jRo1smVHSUajsXXr1vYcymAwNGvW7A8tdzU0AACAChy4ob7KM6M1ARMhAACoxoFrhOpWonEAUqMAAODSsCIEAAD16By5oV5bmAgBAEA1OvlTnZqfGsVECAAAqnHkPkKtZ0JcIwQAAJeGFSEAAKjHgWuEODUKAAB1ibPNg5IT4eHDh8+cOdO4ceNu3brVq1dPaTxy5EhaWlq7du0iIyOrYYQAAADVyN5rhOXl5RMmTBgyZMg333zz8ssv2/bLePPNNwcPHrx58+b+/ft/8MEH1TZOAABwAjrd73kZCVqvCe1dES5ZsuTYsWMnTpzw8fERQlitViFETk7OvHnzUlJSIiIiUlJS+vfvP3ny5Pr161fjeAEAoDZzwmuE9q4IV61a9dRTT+Xk5Ozfv7+oqEhJx+7YsaNt27YRERFCiE6dOvn5+SUnJ1fjYAEAoHZzYEEofbuF2uxdEZ49e3b16tUff/yxj4/PqVOnvv/++44dO2ZmZjZv3tzWp3nz5pcuXeKOoCwiAQDASVmtVs0nrepg74rw1q1bPj4++/fv37lz56RJk2bNmiWEMJlMbm5utj4Gg8FkMlXLMAEAQGv2rGcc2J5ec/ZOhEFBQf369VN+FxgwYMCRI0eEEIGBgdeuXbP1uXr1amBgIHeEOvl7BACA69Dr7ZgydA790ZS9E+H9999v2+333LlzQUFBQoiePXseOnQoLy9PCJGVlZWWlhYXF1dNAwUAAKgO9l4jnDlzZp8+fRo1atSgQYPXX3/9vffeE0KEhoYOHTp02LBh48aNW7FixYQJE5QJEgAAXJNyvlPyv2jM3hVhu3bt9u7de+vWrcuXL69fv3706NFK+5dffjlmzJgTJ05MnTr1448/rrZxAgCAM5C/Rqj5dTOJyjL33nvvvHnzKjUajcZp06apOiQAAHBateCanyzsPgEAAC4NRbcBAEA1yg31kv/HSW6oBwAAuCMHrvlpPQ/i1CgAALg2rAgBAEA9ThiWwUQIAACqcaCItuZ1xzARAgCAapxwQYhrhAAA4NqwIoQ6iqmSX/c2A6N/+3a638mhznDCjXkxEQIAgHocuEZYTSOxG06NAgCAS8OKEAAAVKOTv0Fe69AoJkIAAFCPTid/O4TWMyEmQqh9pHMuxL9YZVMxMv1VSdxIf+ur87NC4iha/3QC5+SEYRlcIwQAAJeGFSEAAKjGCReEmAgBAEA9jpRY03oqxKlRAABwaVgRAgCAepzw3CgmQqgJfIZTJvDJpUmp/8AdhOxcZX+Z3hzq+5w9HcQ1s+3EP/Cd6XbJn0QyZd20/hkHNcyB+wg1/5Dg1CgAALg0rAgBAEA9TlhrFBMhAACoB9cIAQDAlel0zldrFNcIAQDApWFFCA7i0pcynYXVIhH4tFCdhRCWcuogXGfJdnIwXPSULUJKpkaZX4O53471eqY/1c51lmrXMb8nc+30CS6r3NMEZ6cTOtkb5DW/oR4TIQAAqMeBa4Raw6lRAABwaVgRAgCAahwIy2i+gsRECAAAqtEJ+fsItb5ijIkQ/p9k/oMubCYbUaFyLuVUoxDCYraQ7WT/cpnO3EgEM3LuacpkZfiwjBvd7sblXAxEu5sbfdXDjeoshHAzEP3dmJHouXaZxA0P4Ronh2uEAAAAzgUrQgAAUI9O+lSn5st9TIQAAKAaJzwzilOjAADg2rAiBAAA9aDoNtR+/La3TGySK3hGxSwt5XRW02xi2qlsZznXmWk3Ue1cZy5NyqZGyacpU4+No2dOx+glA58GKvBpMNIHMbgz7UY3iYNQX5EbIfd02OipjkvkSuw/DBrSObINE26fAACAusKBHeo1/4UG1wgBAMClYUUIAADqccJrhFgRAgCAanTy7nhudMmSJUFBQQ0bNhw7dmxRUdHtHQ4ePBgTE+Pt7d2pU6djx44pje+++25kZKSnp2erVq3mz59v63zq1Klu3bp5e3tHRUX9/PPPAhMhAADUZocPH549e/a2bdsuX76ck5Mzb968Sh0sFsvIkSMfe+yxmzdvjhkzZvTo0UoBSKvV+sUXX9y4cWPNmjVvvfXW+vXrlf7jx48fOHDgzZs3Z86cOWLECJPJpGO3GFVVcXHx9FffeejPj9fA1wIbqU1lufqZXHFOMpZpLisnO5eVMYHPUqK/iTkI2VkIUVZKpkbpzuVmFaqhsiVVZTbsJTfUFXyFTz4ISgQ+jUw61Eh15vpLdWZHwgzbjWtnnj75cslubgx3aXh44B37JD68IfV0rtRhi27+2u2+S7aJqpKnnnqquLj4k08+EUL89NNPY8aMyc7Ortjhxx9/HD9+fGZmpl6vN5vNTZs23bx5c7du3Sr2GTVqVHh4+GuvvXb06NHu3btfu3bN09NTCBEaGrpo0SKsCAEAoPZKS0tr37698jgqKurKlSv5+fmVOkRGRur1eiGEwWC49957T58+XbFDQUHBvn37lKkxLS0tNDRUmQWVA54+fRphGQAAUI8D9xHqRFFRUUZGRsXGJk2a1K9fXwiRm5vr4+OjNCoPcnNzGzVqZOuZl5en9FQ0aNAgN/e/S1KLxfLoo49269Zt8ODByv+9vbMKE2FZWVlmZmazZs2MRuPdHw0AAFzNvn37BgwYULFl7Nixr732mhDC39//5s2bSuONGzeEEI0bN67Y08/Pr6CgwPbX/Px8Wwer1fr4449fu3Zty5YtSou/v3+lzjExMfaeGp01a1bFkM+tW7eU9o0bNwYHBw8cOLB58+Y7d+60+1kDAEAdpIRApf4InRgwYED6HymzoBAiPDz86NGjyuOjR48GBgbaFoi2DsePH7dYLEIIk8l06tSp8PBwIYTVap0xY0ZqaurGjRtt50LDw8PPnj1rm8KOHTsWHh4usSL829/+VimuU1pa+uijj65cuXLQoEFr1qyZMmXKuXPn3NzoC+xQfdiIBvUPXC0xNhTD5V+o6EpZCd25tNhMt1P9y0okOgsmRMOVWGNDMVzVNKqd70w202+EfFiGia54EL/OunM5Fw8u/0LlXDzoF9ydOYi7B/H8LfXoHzJG7lPLPE038jDsb/II0WhK1dd58uTJvXr1mjp1alhY2Ouvv/7YY48p7S+99FLPnj0HDx7cu3fvhg0bLliw4Mknn5w/f36LFi06d+4shJg1a9Z333335ZdfpqWlCSH8/f1btmwZGRkZHR392muvzZkzZ8WKFeXl5fHx8XJhmeLi4op/3b59e4MGDQYNGiSEePjhh8vKyvbs2aPOUwcAABAiKirq/fffnzhxYkRERFhY2OzZs5X27Oxs5SSnXq//9ttvldOTu3btWrNmjdIhLS2tadOms2bNSkpKSkpK+vLLL5X2f/7znykpKc2bN1+1atV3331nMBgkVoTz589ftGhR/fr1n3vuuWeffVYIce7cubCwMOVf9Xp9mzZtzp8/r9JzBwAA5+NI0e079Z84ceLEiRMrNa5YscL2uH379nv37q3UYdOmTeTRQkJCduzYUbHF3hXhtGnT8vLyCgsLN2zY8Prrr69bt04IUVBQYDvxKoTw8vKyXdK8nXICFwAAnFR5OX22vCIHrhFqfsra3omwTZs2ypzXrVu3CRMmbNu2TQjRpEmTvLw8W5+8vLymTZuyX4nbdQYAAJyBPREQB+ZBzYuNOjI55eXleXt7CyGio6N//fXXsrIyIURRUdGJEyeio6NVHiAAAEB1svca4dy5c3v16uXj47N79+61a9f++9//FkJ06dIlLCzsmWeeSUpKWrRoUVxcXERERHWO1tVJpUMFExDldqY1MVXQSrkM5y2ivYRqrKqdSpOyEVOmnRw5t0Uwd4aefQ2pgKiV28WX28WYwl0X4XasZaumURlOE5cOLWMCn2RqlOlsZj4qZFKXSyNbLPQPHzZMSv3Wzm0RLJUm1fy8XF2j/QJPmr0rQqvV+tZbbz399NNHjx7dvXt3bGys0r5hw4bi4uLHHnvMaDR+/fXX1TZOAABwBk54jdDeFeGrr75KtgcFBS1fvly98QAAANQo1BoFAADV6Oy4HaLyf3GWFSEAAMCdOeE1QkyEtRQdGZAJdAgmF8Pt9ldSzLQXmcj2W4VEezHVKCRDNFxnLrZDZjekcitCNizDJG4kwzJ0uxsTAOE29jOSG0MyFebcmXaTO9Fu5DaAZOq3eVDhGi6cxZX64/Z6FFbqhxXTl9u7Ucc000dxtp/mtYQD1/w0XxHi3j4AAHBpWBECAIBqHCixpvmSECtCAABwaZgIAQDApeHUKAAAqMaRsEz1jMR+mAi1xsUMyT112XQo3U7WHitl0qFc4PNWQRnZXlRA9CejpIKPnhYXkalRujO36y+ZkORSoBy+fB3VyG3MK/NF9dzGvExq1MxUTSNrmHEfCS5NajAQry03EgNT7K2slPh5IrtDstzbxvwEY69RUd31TB1pnRVpUkfohPrbMFU3nBoFAACXhhUhAACox4Eb6rVeZGMiBAAA9dSCItqyMBECAIBqHLiPENcIAQAAtIQVocZkQqNseUaziakJSdWK5Cp5cqlRMh0qhCi6QaRJiwrpiGlxIVNrlEqTkrv1Cj41Wk7twSsZGmXfCTIIyiYemRqk5CUQN2YDXq5OpomJX5qp4rHcnrrcwckdbvWSVU/JMrYWbitkBrtfMdXMJW+5dh1VbJRdjjDLBGc77Qd3hokQAABU49Cp0Woai71wahQAAFwaVoQAAKAaByrLaH66GRMhAACoSuuJTRYmwhrCV/CS2A/WzOxxSu5MK4QopVInXAEzrjoaX2KNaL/FJGtuMSXWyLAMOWzBVIwTktvhsm8EcxAyoETGc6oYCXnJxM3AhWW4nAsTiaKiK0Zm71yuOhoZomG3CGYOYjYTP090zA9FNufCZIjIdr0bPRI989qSB9Hp6XeNGaAgf8xrfpWr9tDppG+HwO0TAAAAWsKKEAAA1OPA7hNar6cxEQIAgHocqDWqNZwaBQAAl4YVIQAAqEYndFxCqqr/pClMhDVGrigXHVZkUqNkKTXBxC/JvXBFFRvzMu3FVOCTbKzi4GRqlHs65A60ggl8culQfnNjJpFLFTbjOnPVxMiEJFvtzEgfhUuZktlOLjXKBT7Jg3AjdGe2CCZfW+l0qEx/NyY1ajAyiVyykpzkSLT+oV3bObJDvdYvKSZCAABQlZMtCHGNEAAAXBtWhAAAoBpHim5X01DshokQAABU44y1RnFqFAAAXBpWhOoj84qyJS7JXCJXbLOU2bGW3IOXDXYyNUhLmY18yW1yucCnVDv3NPkKn1QjlwKVrNdKp0aZkXDITGZ5Of1raLmZqZ/JhBvpg7MhWDrwSe4SzOVUTcyuv+SHnAt2StUUFVw4lqnLamRyrUZ34mXhSqpa3ZjKsUwz2ar5WkcDDtxQr/WrhIkQAADU48g1QtxHCAAAdYUTLghxjRAAAFwbVoQAAKAaB26f0Ly0DCbCGsKFZcqZqlxkpoONnDA72TIb83Kd5XIuZKSFjJZU1W6W2H+4XObg/LDpdm6EZK077t3UMWdYDNRByCMLIcrd6JHo2W1yieOwP1iYkZup/tyPs3LqXRNMTTKuqBu3RbCRSX6VeZC5Ku7dZKr0UVsHc4E17l2mOzPtmp/004aznRvFqVEAAHBpWBECAIBqHCm6XT0jsR8mQgAAUI1Ox55Ur+K/aAsTIQAAqMcJ75/ANUIAAHBpWBFWAyptxiXTLEwAj6z4ZeI24C2RSI2SpdGEZDpUMME8NsUnk+7jqp1xB5cKx3JPh3uDaFw5LaadzLvqmdpjZLUzIYQbs02usBLlxLhNZaWSkBzuIORGvlz9P/d6TLDTkyukRyVvmVJ37H7X1NDJRiGElWun3n6tlzS1iBMuCDERAgCAipzwPkKcGgUAAJeGFSEAAKhK83OdkuRWhCaTKTExcezYsbaW1NTUvn37BgYGJiQkZGRkqD08AABwJsp9hLJ/tCW3InzjjTcuXryo1/8+fVqt1uHDh48fP37dunXz589/5JFHDhw4UA2DrKXYLQapdvJSvxCinKmzZZIpG8blX8icApesKePa2YQOtR+hZBKH7M9nguh2codFLizDbdTHXdWgq6ZJRk50VHRFz5RS4zbwMzIfIXLg3NMxGJmDUCPkKsbpyYJszLvMfa5MZfQPH67yHPkl9fS2g2xWiEkQ0Qdhg0Vku9Y/ymsPB2qNSl9TVJvEivDYsWMbN26cOXOmrWXPnj05OTkvvviir6/v3//+99OnTx86dKgaBgkAAFBd7J0IzWbz1KlTP/jgA6PRaGtMTU2Njo52c3MTQnh4eERERJw8ebJahgkAAE5BJ/9Ha/ZOhPPnz4+Li+vevXvFxpycHB8fH9tfGzVqdP36de4IFu6WHwAAcAbl5fTliYocuUCo9Vxo10R4/vz5Dz/8cPLkyRkZGVevXi0tLc3IyLBarb6+voWFhbZuN27c8PPzY7+SHrdqAAA4MeX8X91jV1jm+vXrAQEBU6ZMEULk5uZeuXJl1KhRe/fubdOmTWpqqtVq1el0ZrP5zJkzoaGh1TxgAACovZRlnuR/0ZhdE2GnTp1SUlKUx//617/effdd5a99+/bV6XQrVqyYMmXK4sWLGzduHBcXV42DdRJkZSaughcXYqQLmEnGKcnaY2Sj4Dfm5fqTgU82HSqzTS73dLggKNnO5Ve5ull6prCZzkK1S6ZGyXAjF5Sz6plqfNxGvlSVPq4cHff0yagqV9SNi0CTOypzOx5zT9/Dk/6h5OltvL3Rox7dmdv1l36XJTcxRmj0Dpywxpr06cr69esHBwcrjw0Gw9dff/3uu+96e3t/9tlnX331leYpWAAAqHuysrJOnTpVRdaksLDwxIkTRUVFldpLS0sLCgoqtuTn5+f9P6W/9EQ4ZMiQ77//3vbXuLi41NTUwsLCo0ePxsTEyB4NAADqFp20KpeEVqs1KSmpY8eOI0eOjIqKunTp0u191q5d27Jly0mTJrVq1WrTpk1K4+bNm2NiYnx8fPr161exc5s2bVq3bh0aGhoaGjpjxgyhVq1RLAQBAEBUQ2p0586dW7ZsSU1NPXbsWPfu3V9++eVKHYqLi//yl7+sWbMmJSVlxYoVSUlJJpNJCNGiRYuFCxcuXLjw9mMeOHAgNzc3Nzd3+fLlAkW3AQBATQ7cR1jlRLh69epRo0bdc889QoikpKTVq1dXutS9fft2X19fZdk3ePBgvV6fnJwshIiKinrggQcaNmx4+zFLSkpu3rxp+ysmQgAAqL0uXLhgux8hNDS0qKio0g3rFTvodLqQkJDz589Xfcz+/fsHBwdHRETs2bNHYPeJu8Lu50k08pE/bhNamVqjMu1SAU7Bp/7IkbPb23LN1LG5+qtyI5HcgpbbylXoiOOQlTkFX+LSzUDVGmU6c7j9Y8vLqV1/y7kRMqFZpl0K+SWNHvSdZ43865HtAS3qM/09iYMz6VAu70puesxdoGJfEbmPlstdNnKk1qjQZWZmrl27tmJj+/btIyIihBCFhYX16v3+afH09FRaGjdubOtZVFTk4eFh+6uXl1fFG9xvd+DAgdatW1sslgULFowYMSI9PR0TIQAAqMaxuycuX75caSLMz89XJsKAgIC8vDylMTc3VwjRtGnTij2bNm1q66D0CQgIqOLLtW7dWgih1+ufe+65N99888iRI5gIAQBAY507d16zZg35TzExMT///LPy+Oeff46IiPDy8qrU4amnniopKalXr15hYeGxY8fsvIWhoKDg1q1bDRo0wEQIAADqUfuG+kcffTQyMnLJkiXh4eEvvfSSbQek4cOHJyYmTp48OSYm5r777ps2bVpSUtL7779///33h4eHCyEuXbq0ZcuW/fv3X7t27ZNPPmnRokVCQsKBAwd+/PHH2NjYkpKSBQsWdOrUqX379gjLAACAanTyNxJWPRE2a9Zs+/btP/300zvvvDNz5sykpCSlvUuXLi1atFAer1u3zsfHZ86cOQEBAV999ZXSePPmzYMHDxqNxvj4+IMHD549e1YI0bhx499++23+/PmffvppfHz81q1b3dzcsCK0g1z+Q52NecupYAibFpFJkZAVubjOgo/5kNgUiUGi5BUXIuEuwJMHcbPSX5GN8zDIy/5k5kII4cY8TRL3kbByLzix/bAQzDPikjhGJkVCb8wr+Xu9nqrTVr+BO9m5aXMfsr15m0Zku5cPUWKNq/9HZs2E5PemnnuDrMTrwgVo2Jew7mZoHNhx/o7du3btevuJ0xdffNH22N/f/4MPPqjUoV27dkuXLq3U2KpVq48++qhSI1aEAADg0rAiBAAA9Thh0W1MhAAAoB4H7iPUukgnTo0CAIBLw4oQAABU40BYBqdGnYB07SmZjXnZ0mtU3SyuM5s/JJtlKsMJvlaZVJ0p7htD6huGi6QaqKwmVzOMDc0yr62OOmnCjYR7TSzUu8nmDJmDc6+V2UQcyM3AxCbpumZ08pY7YcWWr6PauafDJWy5HZLJDzlXF5Db25kM0xqZbx+plDIX06USpkJo/5Mf/gATIQAAqEYnHKk1Wk2DsROuEQIAgEvDihAAAFTjyA31Wp8pxkQIAACqcWQbJq1nQpwaBQAAl4YV4V2QqUHKVbhkA59U6E8mrMe2S6dDubwr1S67/zBZPdXMdOZGQpa45FKjVmYHXm7kUlsHWy3c9r6EhvfQCc7Gwd5kOxmOFULkXy++vbG4iC7CyQc+JTpb6UgmvedzYX4p2Tn7YgHZzsUvPb2JH1bcR4KLqrpTuwRbvOgfg1J5V6sbMxKytep/cXYOVJbRGiZCAABQD+4jBAAAV6YTd9pXifgvGsM1QgAAcGlYEQIAgHqw+wQIwaRlpHMuVJ02yTwL2c5W9mLO63ORARIXOTEzG6WS8QoT1Sj4iAq5TS67RXA5UwqLeYfIkmxWqqqZ4KNP5Oa0oVH3kJ27DWhOtrt70N+tR/dn396Y9us1sjP5ggvujWPeee4jRG6HezOPDstkZtwk24sL6Q2IfRp53N7oWZ/YrVcwyRohhBfVX6p0nxDC6C7xbeWCG/bq5O8L1PxJ49QoAAC4NKwIAQBANTqd/A3yWt9Qj4kQAADUg2uEAADgyhypNVo9I7EfrhECAIBLw4pQfVIb+Ur9KsRUB2PDimQ798uam0FyA1U9tdks8+zZEVJhRa5OGRuapV5F9iqFZBKSLA7Hl1ij2w1GorJX02b1yc4RsU3Idg9P+rv1elbR7Y2/peWTncuY+C75BlnpMm3sG0G+tFxOtfAGnSbVM7+c0yNkRsIdhAyCGj24lDL9gltkEt3cB06ut1NB0W0AAAAng4kQAABcGk6NAgCAmlB0GwAAXJcj1wiraSh2w0SoPqZSGdOZOTkt90liq6YRjXqqIJngQzEGI7cLIBVRYbModDP5NLkd6Wbm2FUAACAASURBVKwWiYNzuRU2pSBVIotNRtDMZUQcIyf7Ftn5zNEcst3oTr9B2b8V3t5YUswGXZhmKgDCZYK4zTWpPSDZ/Sy5anxUTTshRDm5dSXVWFU7GX1ivqLcR0jyc6X9z/7q44T3EeIaIQAAuDSsCAEAQDXOePsEJkIAAFCNA5VlNIdTowAA4NIwEQIAgEvDqdFqQJ0WYAt+MQlJsp0rHCUVBGXToVxEkNqGVAjh7kGUDSvzoA/uZmTaqbpuBiNXlYpLkxKNUmXnRBVlw6ivyb2GnDIqNZp+IpfsXFRA70xL7j8shMi9QqRPi26W0QdhRm4l3ky2YJ5UFUE2YiqzmzTXn+8sNxhQkTNeI8SKEAAAXBpWhAAAoCZn25cXEyEAAKjICW+ot3ci3L59+4YNG7Kzs/39/SdOnNizZ0+l/caNG2+99dbp06ejoqKef/55b2/vahsqAADUdg5cI9ScvdcIMzIyOnToMHny5LCwsISEhL179yrtI0eOPHPmzOOPP37w4MFJkyZV2zgBAACqhb0rwmnTpikPEhMTDxw4sGPHjl69eh07dmzfvn1Xr1718vLq2rVrYGBgRkZGSEhItY22luGCoGSFTyYdyrWTcUoue8lVoSSDnVyNR24knvWNZHuDe4hG7wbuzMGJepiC2/aWqf3IRQHN1D9IbQVcBalqqG5M3tXNjXiDCvLonWlLb9FlQrmAsaWcePrcr+TcQcqpg7AfWurpCCaSylXTZROpbKzX/mPwIW2Z703utdL8PF4t58AN9ZovIOVSo+Xl5UePHv3ll1969+4thDhw4EBsbKyXl5cQwtfXt127dikpKdUyTAAAcAo6h/5oSmIiXLp0aaNGjaKjo4cOHdqvXz8hRHZ29j33/Hdd4O/vn52dzf13C/f7PAAAOIO6+mNcYiJMSkoqKCjIyMjYvXv3woULhRDe3t6lpf89vVNcXFxFWEbP3Q0OAADOwJ4f4zohdL+fH5X4UwODr4L05NS6devx48fv3LlTCNGsWbPz58/b/unChQvNmzdXcXAAAOBcdMLZpkH7wzK//fabMsmVlpZu3749KipKCPHggw9OmTJl//793bp127lzZ0lJSZ8+fapxsBrh3iT+Qr1MvIIreEa1G92pQlhCuNej38dyKhjCxSi4Om2N/OuR7Q39iPbCfLqyV9pReuRlJUTtsbJSolHwu62S+MyFZJUt8lVhPhPc78sGKuVE5qEcwEVXSORHQgih0xHteg96hEamkB4ZzuI+tFzyi/sckt9B0hk06rXivyLZLLfztvY/42ueY/cRalr9zt6J8P77769Xr56fn9+pU6c6dOgwZ84cIUT9+vXfe++9hx56KDo6+siRIx999JGHh0d1jhYAAEBl9k6EaWlpaWlp+fn5zZs3b9asma19ypQpDz30UHp6elhYmJ+fX/UMEgAAnIQDN9TrdM6xInRzc4uIiCD/qUmTJk2aNFFvSAAA4KwcuY+wekZiPyQ5AQDApaHoNgAAqMah/QiraSz2wkRoB67QEnNSm3xTuXJN3G6rBqpqGpnKE0J4eNLtZEKSzKMKITw86Q9DYEsf+9tv5tJlw/KuFZPtWedv3t7olkuPkC3WRb2EXDiUDc0ybxB9HK48mEwklS1UxsUpuYpf1Bdl73uWuRJDZixFVellop1sFPyH2ci0k98RXOiaS+SSAVHue1Cq9Jr0j36p3lDNMBECAIBqHLhGqPnvBbhGCAAALg0rQgAAUI1OyF8j1HpJiIkQAADU44Q71OPUKAAAuDSsCO8Cu/sn0cbVM+Rib2QwjwvUuZfR7eQmtNxZCy5ox22Te6vQdHtjUQFda7S0mN5slh4h8+sZWxNSprgrl+vk3gghJHb95ZBpUq4aKrd1MFvMknqm7BbBzMtCZjK5DZ+54DH5GhqNcqlRNk1KDUZqJIIt+spVPZWI70ptBSyE9mug6uOMG/NiIgQAAPU4UGJN698LMBECAIB6cI0QAADAuWBFCAAAqqmOBaHJZPrxxx9zc3P79evXtGlTss/BgwdTU1Pbt2/fsWNHW2NZWdmpU6fc3d3vvffeip2PHz9++PDhsLCwLl26CEyE1YEuscZFVLgSa9RVfTa8wO19Su1ky21vW8LkWbIvFpDt1y4X3d54I7eE7Hw1s5BsLyslBsO9JlzsotydrDHG7UBLP/1yejNgQZ41ocfBB44M1Mi5qnscLp9joL6ouwf9rc3tqUumTrh0ko5rl4mJkZ9wwddvI19Drhwdf3AqccN8W0nVaeOviml91q/G6RwoOFdl97KysgceeEAIERIS8te//nXHjh0xMTGV+sybN++TTz558MEHX3zxxZkzZ86cOVMI8fnnnyclJXl5ebVt2/aXX36xdV6yZMkrr7ySmJj497//ffTo0W+88QYmQgAAqL3WrVtXUFBw8OBBo9H46quvvvLKKxs2bKjYIScn58033zx06FB4ePiRI0d69eo1depUHx+fhISEK1eubNq06YMPPrB1Li4unjNnznfffdejR48LFy5EREQ8+eSTuEYIAADq0Tn0h7dp06Zhw4YZjUYhxMiRIzdv3lz+xxM4O3bsaNOmTXh4uBAiOjo6MDBw9+7dQoiAgIBGjRpVOtq+ffs8PDx69OghhGjZsmVMTMy2bdswEQIAgGqU+wil/lQtMzMzODhYedysWTOz2Xz16tVKHZo1a2b7a3BwcGZmpj1HUw6YmZmJU6MAAKAah/Yj1J04ceLFF1+s2NinT5+BAwcKIcxms5vb71eIlQcm0x+qeZSXl+v1/13UGQyGSh3u2BkTIQAAaMxoNPr6+lZsMRh+n54CAwOvXbumPL569apOpwsICKjYMyAgwNZBCHHlypWgoCDuC93euVevXpgIHSf1O49esmwYXa2KCdQZjXTk0UTViCpjflkzldIH4fbaLblF/M5VmE+XWCu8QbeTL4t7Pe5jKVPTjslkck9TMGFaHVVijXvXuE1o6dSoxD67QvDPiKxJVs+Lfg25EZIH5/YZ5uKu5DpAauvpKtrJwKc7E4KV2vXXyERM2dJrdC6c7OuSHLp/Iiws7IUXXiD/sU+fPt9+++1LL70khNi+fXv37t3d3d2FEIWFhUaj0cPDo1evXo8//nhOTo6fn19WVlZaWlr37t25L9WlS5esrKz09PTQ0NCCgoL9+/d/8MEHuEYIAAAqcuAqYVUz58SJE9PT05OSkhYuXDh79mxlRhRC9OvXb8mSJUKI1q1bjxo1asiQIYsXLx46dOikSZMCAwOFECdPnkxKSlqxYsX58+eTkpI+/PBDIYSfn9+0adNGjBjx0UcfDRkyJCEhISIiAitCAACovRo2bPjLL7+sWLHiypUrGzdutK32XnjhhbZt2yqPly9fvnLlylOnTk2fPn38+PFKY4MGDWJjY2NjY5W/tmjRQnmwYMGCr7766tdffx0zZszkyZMFbqgHAAAV6aphN4mAgADbQtBm+PDhtscGg+HRRx+t1KFZs2aPP/44MUKdbuzYsWPHjv3vf1dvqAAA4PKcsOg2JsK7wL55EtuV8WEZop0rHMXmDujd15j6WEwEwlJOp0jKSojUidlMd+ZyBx5UpoN9TaReK+YrljAHEUyFuXLqMNwb4eHJFDajUk7cfoTlzH6EfK0yKgAi85EQzOfTamVqjDEfFfKLsptoqtHO5aq4N4I8CFuhkKskR72EbIE1rX/E1zzHbp+opsHYCWEZAABwaVgRAgCAahzZob56RmI/TIQAAKAeJ7xGiFOjAADg0rAiBAAA1TgQltF8RYiJUH1M1IzpzATw3KjqaHz1Ka6gFBHsNJZyddpkD060e5i57WNpZBKSfO5VtJNJSC4eyWUv6bpZQpRToVnu4FyIkcwfculQHZO8Zfd2pkKz7NNkPm/k55D7acZ9JMi6blyxNy7YybZLHZx5I+jUKFdKjd1/mN56m+zsgpzwzChOjQIAgGvDihAAAFSl+RJPEiZCAABQjxPeUI+JEAAAVOPIfYRaryBxjRAAAFwaVoTVQCY2yu1xSuYM2bAiU57RbCLeX1MZHUp0L2WCdsxOtvWogwtuU1mZ8qFs4lGmnYvjcriDW6lXS7aSJ7nBramMfmHZsrTMM2IqysptNitVr5ULx3rVJ9o9vY1kZ7adCYKS/dn9hz2ZqDOVd+Wq77IJW4RGq+ZAbFRrmAgBAEA1Op38NT+tf4/AqVEAAHBpWBECAIBqHNiYV+sFISZCAABQj07odJIXCbWeBzER1hSZTXyFYLIbXHiB2/vUgyrWRSZohBBmEx2i4dotTIUwklTFLy7nIvXdQuZTqmiXqrMlu5UruV8xvycz/YJz2Q33elStO6pR8CMnC+xxnysuouJVn8izeNV3t7+zEMLLR6J/PS/6IPzTp4JFTOk+Nm+l+fqllnPCGmu4RggAAC4NK0IAAFCP/A31msNECAAAqnFgGyadTie3bY3a7J0Ic3Nzd+zYcfHixRYtWgwdOrRevXpKu9Vq3bRpU2pqanR0dEJCQrWNEwAAoFrYe40wNjb2q6++ysnJWbx4cceOHfPz85X2GTNm/O1vf7t169bTTz89e/bsahsnAAA4A51DfzRl74rwwIED/v7+Qojy8vL27dt/++23U6ZMyczMXLZs2dmzZ4ODgydMmNChQ4dnn33Wz8+vOgfstCRjozrqVxQ9E28zGunzChaqFBa3H6zFItdO4pOQMlWsmNfKyo2QekZSoVYhhGCqprmRcUome8mdECKTt1wKtJzZmJfrT2Y7uap7ctssM9lLTzbwKZEa9aY6Cz41KlVijXuaZIU5ndQGvAiN3okjRbfZ4ow1xN4VoTILCiHc3Nzc3d3d3NyEED/99FNUVFRwcLAQok2bNq1atdq7d281DRQAAGo/3e9TocQfzZeE0rdPrFmz5vr160OHDhVCXL58OTAw0PZPAQEBly9f5v6jlbuBCwAAnEFd/TEulxpNTk5+8sknv/vuu0aNGgkl6lPhdbFarZrvrwgAABpzthvqJSbCffv2jRo1avXq1XFxcUpLYGBgVlaWrUN2dnZQUBD33zFHAgA4NXt+jNfljXkPHTo0fPjw5cuX9+3b19bYt2/f48ePX7p0SQiRlpZ28eLF3r17V8swAQAAqoe9K8JBgwZ5e3uvXLly5cqVQog//elPY8aMCQoKSkpKio+P/9Of/rRmzZpnn33W19e3OkdbB0n9KsRmMpn9YI3UCX2rlX7T2ZP/TDs5ctk9dclyjuxAmNQo2coViuQymdzTJ3OGBiYEy72dZPlQ7jXh8q5cfzIhyaZGmXayfz1viZqiQghvKvDJdeaip2w7FRDlnqZUMVi2pCjdDHfgyA311TQUu9k7ES5evLi8/L8bakdERCgPFi1a9MMPP5w4cWLp0qUVF4sAAOCKnLDotr0T4fDhw7l/evDBBx988EGVxgMAAM7MgVqjWk+E2H0CAABcGopuAwCAanTy9wjIbuSrOkyEtRT5ydBzJdqZGIUwUpEByR1rOXRYhgkeSLVz30XctrfkbqumenQohgvLyFWSk8xXWCzECLmgBxvbkXn63MHZdk/iR4HUBrxcO1kaTfBJnHrUSIQQ7lS9QCNXGI9JM9ERKu6nttbn66DG4NQoAAC4NKwIAQBANY7tR1hNg7ETJkIAAFCNA5VlND8LjYkQAADUo/1mEtJwjRAAAFwaVoS1FZN6I1v13O8z1NtrFXRuUG4gkiXWyEJlQggDVRyOq4/F7RNbVko8T1MZnQ41m8rJ9nIzU7+NTJOyFdbof6CDoFw8lK3TpkqJNfqNcPcgXkMPT/og9byYKmhUEJSLnnpQKVAhhDvzLpMfFa64IBfrJd8gra9P1TUOXSOsprHYCxMhAACopw7vPgEAAFAnYSIEAACXhlOjAACgGp1wYBsm3EcIduM+XVYuWkMt+A3Me65jQjTcZ5qsVsUFOsikg2C2x+MyGmUl9NDLSon8i4lqFFWEaLjSa8zugCS2VpdMRkOVwBG37yBXk0wqcePBVEEj00xkaTTBx3aMZF1A5qPF7TrJfmi1vhblCpzxPkKcGgUAAJeGFSEAAKinDm/MCwAAcEfOWGsUp0YBAMClYUUIAACqcSAso/WZUUyEdQL/sSODdkxfdi9T5tBUO58a5Sp+UYFPpsgWF/iUTI1yJdYkUqNckJRLgZOvFZd45DbgZVOjVDu5W68QwsjEd8mydmzElNvdlzoINxJy2IIPzdJV05AOrZ2c7RohTo0CAIBLw4oQAADUpPkN8rIwEQIAgGocuUao9byJiRAAAFTjyO0TWq8gcY0QAABcGlaEdRn5axlXmFTvRmchdUxslPylT+9Gd3Yz0JlMMpdoZkKJZiY1aqL22uU6czVFuY156dQot6cug97EmEmN8q+hRJpUqjAp1852Zg5ObpPLhWDZvXO5lDL5udX6lBoQUFkGAABcmTMW3cZECAAAtVpWVtaKFSvy8/OHDh3ao0eP2zuYzebPP//81KlT7du3nzBhgpvb72eV0tLSVq1aZbFYxo4d265dO6Xxiy++KCkpUR6HhIT0798f1wgBAEA1ut/jMjKqXBLeuHGja9euFy5cCAgIGDJkyObNm2/vM2XKlBUrVoSEhCxevHj69OlK45kzZ7p06WKxWNzd3ePi4o4fP660z5w58+DBgxkZGRkZGVevXhVYEQIAQG32xRdfhIaGLl26VAhRv379N998c/DgwRU7nDt3bu3atZcuXfLz8xs2bFhISMjcuXMDAwPff//9Rx55ZN68eUKIGzduvPfee8uWLVP+y6xZs9q2bWs7AiZClyNVj00INlnMVLySS0aQmQ6DkcmzuNMRFXczEa4pZzbUlSqlJoSwWsiwDNmXL70mUR2sihpjTESF6i+7uy99EC6ew4xERx6E+7yhOlodpfp9hMnJyfHx8crjAQMGTJs2rayszN3d3dZh79690dHRfn5+QoigoKDw8PB9+/aNGDEiOTl57ty5Sp/4+PgnnnjC9l/++c9/NmrUqHPnzsqJVpwaBQCA2isrK6tx48bK46ZNm1qt1uzs7IodsrOzbR2UPpcvX1b+Y5MmTZTGJk2aZGVlKY/j4uJu3bqVkZExbNiwZ599VmBFCAAAapK/oV7oxIEDB0aNGlWxrX///o8//rgQwmg0ms1mpVF5YDQaK/Y0GAzl5f+9h8pkMikdDAZDxf9o+1+bNm1SHkyfPj0yMnLGjBmYCAEAQD0O3EcoRFBQ0MiRIyu2REVF2f5JWeEJITIzM41Go22dZ+uQmZlp++vly5eDg4OFEMHBwRX/Y1BQUKUvGh4e7ufnd+7cOZwaBQAAjQUHB4/8o3vvvVf5p8TExA0bNphMJiHE2rVrBw0apNwdcfDgwUuXLgkhBgwYcPbs2dOnTwshjhw5kpWVdf/99wshHnroobVr1yoH+eabb4YMGSKEKCkpsRXE+N///d/c3Nzw8HCsCAEAQDWqh2VGjBjx0Ucf9enTJyQk5Icffti+fbvS/sQTT4wZM+bpp5/28/N76aWXBgwYEB8fv3Xr1pdfftnHx0cIMX369Li4uEGDBrm7ux8+fHjfvn1CiF27dj399NOdOnUqKSnZuXPn22+/HRQUpJMtFuWY4uLi6a++89CfH6+BrwU1g/ngMK1czJJqJ4OaQggL00725zrzByGbmdQo3VcO983Pbdir52KW1GkdvoCZRKyX68yOnPoH2ZAy1GbDwwPv2OfL3RnXbpRIHfb4//6Yl7pr/fr1XAeTybRr167c3Nx+/frZzoueOHHCz88vICBA+euhQ4dSU1Pbt28fHR1t+49FRUU7duywWCwDBgxQZkez2XzkyJG0tDRPT8/Y2NjmzZsLhGUAAEBlav+WYzQaH3zwwUqNkZGRFf9633333XfffZX6eHt7Dxs2rGKLwWCIjY2NjY2t2IhrhAAA4NKwIgQAANU4FBrVGCZCAABQjSMb82pdTwinRgEAwKVhRQgO4jKMd38QK1eEkwkxkslnPqcql2tVJyFK4rKXbH+JDKdUsJPtL5n41PrXeqgdsDEvAAC4NJ30qU7NT41iIgQAANU4ckN99YzEfnLXCLOzsyuV/RZCZGVl7d69+8qVK+qNCgAAoIbYOxF+8skngYGBQUFBjz32WKX2qKiot956KzIy8l//+lc1jBAAAJyKTv6Ppuw9NdqzZ89du3Z9//33ycnJtsaCgoJZs2bt3LmzS5cuu3bteuSRR0aMGOHh4VE9QwUnJnWqRMd8W3AhGqY33dnKfc9J7rUrRZ1vc6mnr0ZlM60v3IBT0gkd9y1ca9m7ImzXrl1ERESlPbJ/+OGH5s2bd+nSRQjRt29fLy+v3bt3qz5EAACA6nNXYZnffvutVatWtr+2bNny4sWLdzsiAABwWqrvPlED7moiLCkpcXd3t/3Vw8OjuLiY62yxMOX9AQDAGVgslkrnBQm14JqfrLuaCAMCAnJycmx/vX79emAgu0nHnV8+AACoxez5Me5yJda6dOly8ODBoqIiIURubm5qamrnzp1VGhgAAEBNsHdFeOLEie+//37v3r3p6elvv/12hw4dBg4cGBkZ2bt373Hjxk2ZMmXJkiVDhgypeMkQwEGytcdUOREj/UVrNa1/wwbX5cA1Qs2/zexdEZaVleXl5bVv337o0KF5eXnKKlAIsWbNmpiYmFWrVvXu3XvlypXVNk4AAHAGznYTobB/RRgTExMTE3N7u4+Pz9y5c1UdEgAAQM1BrVEAAFCR84VlMBECAIBqnHAXJkyEUFfVrfALgNNwwpkQ9/YBAIBLw4oQAABU48AN9ZrDRAgAAKpxxlqjODUKAAAuDStCAABQjU4nfTuE5qdSMRECAICqtD7VKQunRgEAwKVhRQgAAGpytgUhJkIAAFCPM+5HiIkQAADUg8oyAAAAzgUrQgAAUI0z3lCPiRAAAFTjjCXWcGoUAABcGlaEAACgHicMy2AiBAAA1ejkr/lpfiYVEyEAAKhGJ3Q6ySWebH/V4RohAAC4NKwIAQBAPQ5cI9QaJkIAAFCP/H2Emk+cODUKAAAuDStCAABQDYpuAwCAa8M1QgAAcGU64cCKsJrGYi9cIwQAAJeGFSEAAKjGCSusYSIEAAAVOeFMiFOjAADg0rAiBAAA1eh0DtwOgdsnAACgzsAO9QAA4NJwjRAAAMC5YEUIAACqceSGeq2XhJgIAQBANTrsPgEAAOBcMBECAIBLw6lRAABQjUPbMFXTWOyFiRAAAFTjwDVCredBnBoFAADXhhUhAACoxwlvqMdECAAA6pG/Rqj5RUKcGgUAAJeGFSEAAKhGJ7/A03pBiIkQAADU48jtE1pfJKyhidBgMJxK/iF5zRc18+VqWE5Ojq+vr15fx88z42nWJdevX/fz85PfN87J4GmqK+L77yMiIqru06+Vv+xhQ4pjftaXODooFeisVmvNfKXc3Nz8/Pya+VoAAKC6Zs2aubu7az0K9dXcRAgAAFAL1fHzPwAAAFXDRAgAAC4NEyEAALg0TIQAAODScB+htCtXrqSkpGRmZvbr1y80NNTWfuHChc8//7yoqGjkyJGdO3fWcISqOHTo0Pbt269duxYRETFu3DhPT0+l/ebNm59++mlmZub9998/ZMgQbQd593bu3Llv3778/PwWLVpMmDDBz89Pac/Pz//000+zsrL69+8/aNAgbQepojVr1nh4eAwdOlT5a2lp6bJly86ePRsTEzNu3Dhnv2lk48aN2dnZyuN77rnn4YcfVh7n5uZ+9tln2dnZCQkJ8fHx2g1QNVevXl25cuXly5dbt249adKkhg0bigof2n79+g0ePFjrMToT5/7ca6JXr15vvPHGCy+8kJKSYmvMzs7u3Llzfn5+kyZNBgwYsGfPHg1HePfy8/MTExOvXbvWokWLVatW9erVq7S0VAhhsVj69u27b9++0NDQp59++h//+IfWI71bX3/9tcViCQkJ+fe//92xY8fc3FwhhNls7t27d0pKSkhIyF/+8pePP/5Y62GqY+PGjVOnTn3nnXdsLaNHj/7mm2/atm27cOHCWbNmaTg2Vbzzzjvbt2/PyMjIyMi4dOmS0lhWVtazZ8+jR4+2bt16ypQpn3/+uaZjVMGZM2c6dOhw4sSJVq1apaWlKc/UbDb36dMnJSUlNDR0+vTpS5Ys0XqYTsUKksrLy61Wa3R09OrVq22Nr7zyyrBhw5TH77zzzqBBg7QZnErKy8tLS0uVx8XFxQ0bNtyzZ4/Vat2yZUurVq1MJpPVav3xxx+DgoKUx3WAxWJp3br1+vXrrVbr+vXrw8LClDdaecrKY6eWn58fGRn5yiuvdO/eXWk5fvy4t7f3zZs3rVZrenq6p6dnTk6OpmO8Wz169Ni4cWOlxq+++qp9+/YWi8Vqta5fv75t27bKY+eVkJAwe/bsSo3KU1M+qFu3bm3ZsqXZbNZidE4JK0Jp5OmjPXv22E65DBgwIDk5uWYHpTK9Xm+7bdZisZSVlfn4+AghkpOTH3jgAYPBIITo06dPTk5OWlqalgNVT1paWn5+/r333iuESE5O7tevn/JG9+vX7+LFi+fPn9d4fHftmWeeeeaZZ4KCgmwte/bs6dq1q/LOhoSEBAcH//LLL9oNUB1bt25dsGDBli1brP9/h/SePXv69++vVF2Jj48/c+bM5cuXNR3jXTGZTDt27Bg6dOjy5cs//vhj28K30of20qVLdeBDW2MwEaojKyurcePGyuMmTZoUFRXdvHlT2yGp5bnnnuvdu3fHjh2FENnZ2ban6ebm5ufnl5WVpenoVPD8888HBwd36NBh/vz5ykRY8d10d3f39fV19qf5448/njt3bsqUKRUbK76bQogmTZo49QwhhGjXrp2Hh8fVq1dnzJiRmJhosVjEH99NLy+v+vXrO/W7+dtvv1kslieeeOL8+fPHjh2Ljo4+efKk+OO7aTQa68CHtiYhLKMOg8FgNpuVx8oDo9GoLlI8agAABCVJREFU6YjUsWjRoh07dtgueRoMhvLyctu/mkymOlBvae7cuTNnzty7d++0adOioqI6d+5sNBrr0tMsKiqaMWPGt99+W6kWZd17Nz/55BPlwYsvvhgWFrZ9+/aEhISK35tCCLPZ7NRPU6/XW63WJ554Qvm1xmQyLViwYNmyZXXv3axJWBGqIzg42PbbdGZm5j333GOLWTqv999//8MPP9y1a1dAQIDSEhwcnJmZqTwuLi7Oy8ureKrNSXl7ewcEBIwcOTIhIWHDhg3ij0+zoKCgoKDAqZ9mcnJyZmbm+PHjO3XqNG/evKNHj3bq1MlisVR8mkKIzMxMp36aFfn6+rZr1+7cuXPij9+b169fLykpceqnGRgYqNfr27Vrp/w1MjLywoUL4rYP7c2bN536adYwTITqSExMXLdunXIqZu3atYmJiVqP6G599tlnCxcu3LFjR7NmzWyNiYmJO3bsUIqnr1+/Pjw8vOINJE7HbDabTCblcVlZ2dGjR1u0aCGESExM3LZtW0FBgRDim2++iYmJCQ4O1nKgd6dHjx67du1aunTp0qVLJ0yYEBISsnTpUr1en5CQ8Ouvvyo/Rv/zn/+UlpZ2795d68E6zmQy2VZ+v/322+HDhyMjI4UQiYmJW7ZsuXXrlhDim2++iYuL8/eX3h6h9vDw8Bg4cOD+/fuVv+7fv1+ZFBMTE3/44QflQ7tu3bqOHTtW/M6FO9A6reN8pk+fHhsb6+npGRISEhsbm5KSYrVaCwsL77vvvt69e48aNapp06anT5/Weph3JTMzU6fTtWjRIvb/bd68WfmncePGtWvX7s9//rO/v/+mTZu0HeddunjxYtOmTYcNGzZu3LiWLVv279+/uLhY+aeHH364ffv2kyZN8vf3/+GHH7Qdp4o+/fRTW2rUarXOnj27VatWU6ZMCQgIWLp0qYYDu3vp6emBgYHDhw8fNWqUr6/vE088obRbLJbExMSOHTtOnDjRz8/vp59+0nSYKjh06FCTJk0mTpw4aNCgNm3aXL58WWkfOXKk7UO7bds2bQfpXLD7hLQzZ85UDMKEhYUpubvS0tJdu3YVFhb279/f19dXuwGqoKys7NixYxVbWrVqpdxsbrVa9+7dm5mZ2b1795YtW2o0QNVcvHjx8OHDJSUlbdu2jYmJsbVbrdbk5OTs7OwePXo0b95cwxGq6/r16zk5OeHh4baWlJSUM2fOdOzY8Y77zNV+qampqampFoulQ4cOYWFhtnaLxbJ79+6rV6/26tXLqRf3Njk5Obt27WrUqFHPnj1tV2GsVuuePXuysrLq2Ie2BmAiBAAAl4ZrhAAA4NIwEQIAgEvDRAgAAC4NEyEAALg0TIQAAODSMBECAIBLw0QIAAAuDRMhAAC4NEyEAADg0jARAgCAS8NECAAALu3/AA+ZuawKQBfhAAAAAElFTkSuQmCC", "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n" ], "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "execution_count": 5 } ], "cell_type": "code", "source": [ "terms = [Kinetic(),\n", " ExternalFromReal(X -> pot(X...)),\n", " LocalNonlinearity(ρ -> C * ρ^α),\n", " Magnetic(Apot),\n", "]\n", "model = Model(lattice; n_electrons, terms, spin_polarization=:spinless) # spinless electrons\n", "basis = PlaneWaveBasis(model; Ecut, kgrid=(1, 1, 1))\n", "scfres = direct_minimization(basis, tol=1e-5) # Reduce tol for production\n", "heatmap(scfres.ρ[:, :, 1, 1], c=:blues)" ], "metadata": {}, "execution_count": 5 } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.2" }, "kernelspec": { "name": "julia-1.7", "display_name": "Julia 1.7.2", "language": "julia" } }, "nbformat": 4 }