{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Anyonic models"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We solve the almost-bosonic anyon model of https://arxiv.org/pdf/1901.10739.pdf"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter Function value Gradient norm \n",
" 0 8.410347e+01 1.624468e+01\n",
" * time: 0.003039836883544922\n",
" 1 6.368490e+01 1.253134e+01\n",
" * time: 0.008939981460571289\n",
" 2 5.741571e+01 1.676497e+01\n",
" * time: 0.022531986236572266\n",
" 3 4.275783e+01 1.139997e+01\n",
" * time: 0.0408778190612793\n",
" 4 3.487506e+01 1.017819e+01\n",
" * time: 0.05848383903503418\n",
" 5 1.490009e+01 3.157736e+00\n",
" * time: 0.14609384536743164\n",
" 6 1.138340e+01 3.520194e+00\n",
" * time: 0.15999794006347656\n",
" 7 9.900684e+00 2.191030e+00\n",
" * time: 0.17305684089660645\n",
" 8 9.126685e+00 3.320207e+00\n",
" * time: 0.18350481986999512\n",
" 9 8.811054e+00 3.575682e+00\n",
" * time: 0.19388580322265625\n",
" 10 8.504275e+00 3.066914e+00\n",
" * time: 0.20434999465942383\n",
" 11 8.062185e+00 1.735209e+00\n",
" * time: 0.21474194526672363\n",
" 12 7.559192e+00 1.312268e+00\n",
" * time: 0.22731399536132812\n",
" 13 7.172199e+00 1.502411e+00\n",
" * time: 0.23790788650512695\n",
" 14 6.895554e+00 1.094797e+00\n",
" * time: 0.24822092056274414\n",
" 15 6.719928e+00 8.923590e-01\n",
" * time: 0.2929229736328125\n",
" 16 6.607344e+00 1.118035e+00\n",
" * time: 0.3041698932647705\n",
" 17 6.502045e+00 7.087406e-01\n",
" * time: 0.3154029846191406\n",
" 18 6.349434e+00 7.500807e-01\n",
" * time: 0.32590389251708984\n",
" 19 6.119833e+00 7.781240e-01\n",
" * time: 0.3365139961242676\n",
" 20 5.941530e+00 9.128185e-01\n",
" * time: 0.3469219207763672\n",
" 21 5.830842e+00 6.289483e-01\n",
" * time: 0.35698390007019043\n",
" 22 5.767541e+00 5.504176e-01\n",
" * time: 0.36728882789611816\n",
" 23 5.757317e+00 1.224940e+00\n",
" * time: 0.3751859664916992\n",
" 24 5.724518e+00 7.870377e-01\n",
" * time: 0.3831508159637451\n",
" 25 5.714197e+00 6.780425e-01\n",
" * time: 0.3911128044128418\n",
" 26 5.697801e+00 7.782283e-01\n",
" * time: 0.42572784423828125\n",
" 27 5.677513e+00 5.120291e-01\n",
" * time: 0.43662190437316895\n",
" 28 5.661575e+00 4.481369e-01\n",
" * time: 0.4474318027496338\n",
" 29 5.650541e+00 6.530572e-01\n",
" * time: 0.45535778999328613\n",
" 30 5.624435e+00 5.653876e-01\n",
" * time: 0.46305394172668457\n",
" 31 5.608646e+00 3.829223e-01\n",
" * time: 0.4734818935394287\n",
" 32 5.600703e+00 4.283360e-01\n",
" * time: 0.4841279983520508\n",
" 33 5.592613e+00 3.799936e-01\n",
" * time: 0.4950439929962158\n",
" 34 5.590145e+00 3.047189e-01\n",
" * time: 0.5028908252716064\n",
" 35 5.581184e+00 2.546569e-01\n",
" * time: 0.5114209651947021\n",
" 36 5.575281e+00 2.352586e-01\n",
" * time: 0.5197818279266357\n",
" 37 5.571052e+00 2.531070e-01\n",
" * time: 0.5276718139648438\n",
" 38 5.568984e+00 2.161152e-01\n",
" * time: 0.5356478691101074\n",
" 39 5.567733e+00 1.375050e-01\n",
" * time: 0.5613539218902588\n",
" 40 5.566603e+00 1.622589e-01\n",
" * time: 0.569861888885498\n",
" 41 5.565739e+00 1.315509e-01\n",
" * time: 0.5783689022064209\n",
" 42 5.564717e+00 1.945912e-01\n",
" * time: 0.5863549709320068\n",
" 43 5.564227e+00 1.554160e-01\n",
" * time: 0.5944509506225586\n",
" 44 5.563492e+00 1.274773e-01\n",
" * time: 0.6049208641052246\n",
" 45 5.562688e+00 9.623832e-02\n",
" * time: 0.6153128147125244\n",
" 46 5.562032e+00 1.202397e-01\n",
" * time: 0.6230208873748779\n",
" 47 5.561544e+00 7.628104e-02\n",
" * time: 0.6310157775878906\n",
" 48 5.561288e+00 8.709403e-02\n",
" * time: 0.6388468742370605\n",
" 49 5.561095e+00 5.988173e-02\n",
" * time: 0.6469218730926514\n",
" 50 5.560959e+00 4.685654e-02\n",
" * time: 0.6572158336639404\n",
" 51 5.560826e+00 2.624886e-02\n",
" * time: 0.6827578544616699\n",
" 52 5.560742e+00 4.890421e-02\n",
" * time: 0.6913089752197266\n",
" 53 5.560731e+00 5.017957e-02\n",
" * time: 0.6994597911834717\n",
" 54 5.560645e+00 3.419509e-02\n",
" * time: 0.7072999477386475\n",
" 55 5.560596e+00 4.020835e-02\n",
" * time: 0.7176787853240967\n",
" 56 5.560555e+00 3.075647e-02\n",
" * time: 0.7279739379882812\n",
" 57 5.560543e+00 4.017573e-02\n",
" * time: 0.7355008125305176\n",
" 58 5.560533e+00 2.090276e-02\n",
" * time: 0.7430768013000488\n",
" 59 5.560515e+00 1.533790e-02\n",
" * time: 0.7506668567657471\n",
" 60 5.560502e+00 1.263486e-02\n",
" * time: 0.7585279941558838\n",
" 61 5.560487e+00 1.065857e-02\n",
" * time: 0.7664659023284912\n",
" 62 5.560478e+00 9.467035e-03\n",
" * time: 0.7745850086212158\n",
" 63 5.560473e+00 9.109492e-03\n",
" * time: 0.7848029136657715\n",
" 64 5.560471e+00 5.783410e-03\n",
" * time: 0.8103058338165283\n",
" 65 5.560469e+00 7.729054e-03\n",
" * time: 0.8187828063964844\n",
" 66 5.560468e+00 7.873965e-03\n",
" * time: 0.8270039558410645\n",
" 67 5.560468e+00 8.328648e-03\n",
" * time: 0.8348469734191895\n",
" 68 5.560467e+00 4.205446e-03\n",
" * time: 0.845383882522583\n",
" 69 5.560466e+00 6.067369e-03\n",
" * time: 0.8558318614959717\n",
" 70 5.560465e+00 4.957374e-03\n",
" * time: 0.863839864730835\n",
" 71 5.560464e+00 3.706725e-03\n",
" * time: 0.8739898204803467\n",
" 72 5.560464e+00 2.196989e-03\n",
" * time: 0.8816378116607666\n",
" 73 5.560463e+00 1.953301e-03\n",
" * time: 0.891960859298706\n",
" 74 5.560463e+00 2.417199e-03\n",
" * time: 0.899634838104248\n",
" 75 5.560463e+00 1.450785e-03\n",
" * time: 0.90982985496521\n",
" 76 5.560463e+00 1.426738e-03\n",
" * time: 0.9344789981842041\n",
" 77 5.560463e+00 1.438616e-03\n",
" * time: 0.9455118179321289\n",
" 78 5.560463e+00 9.532110e-04\n",
" * time: 0.9561269283294678\n",
" 79 5.560463e+00 1.237838e-03\n",
" * time: 0.9640858173370361\n",
" 80 5.560463e+00 7.899421e-04\n",
" * time: 0.9721689224243164\n",
" 81 5.560463e+00 7.256114e-04\n",
" * time: 0.9830179214477539\n",
" 82 5.560463e+00 5.471578e-04\n",
" * time: 0.9907839298248291\n",
" 83 5.560463e+00 2.145598e-04\n",
" * time: 1.0010159015655518\n",
" 84 5.560463e+00 4.118338e-04\n",
" * time: 1.008711814880371\n",
" 85 5.560463e+00 2.278516e-04\n",
" * time: 1.016535997390747\n",
" 86 5.560463e+00 1.909337e-04\n",
" * time: 1.026641845703125\n",
" 87 5.560463e+00 2.276338e-04\n",
" * time: 1.036787986755371\n",
"e(1,1) / (2π)= 1.7391793874123116\n"
]
},
{
"output_type": "display_data",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3daXgUVb748dNbNpCEkIRMQiSyyKYssogCipCgSNAIiDI6sukMeB0fnDtcl6uouOJ4x51xGZYRFBlU8CKuqICibKIsRmJkJwuLyBIgS6f7/6K8/e86lZxONRUqTX8/T150VZ2uOql0cnLOr37nOPx+vwAAIFo57a4AAAB2oiEEAEQ1GkIAQFSjIQQARDUaQgBAVKMhBABENRpCAEBUoyEEAEQ1GkIAQFSjIQQARDX3mbmMz+fbuHlr+06dzszlIp3P53M6+R/FHG6aWdyxMET5TUuM9YQsU1h4+NQpr6nTHj9+NCPD2bZt23DrdbrOUENYWVk5c+HivHGpZ+ZyAADLjejwu5Blbrhh6aZNB02eeGt+/tHFixeHV6vTF73/2gAAIM5YjxAAEA0cDuFw2F0Jk+gRAgCiGj1CAIBlHA6Hw2SX0O+3uQtJQwgAsBRDowAARBB6hAAA60TgwzI0hAAAy4QRIxR2xwgZGgUARDUaQgBAVGNoFABgmTAS6u0eGaUhBABYJ4wYoe0P1zA0CgCIavQIAQDWcUReQj0NIQDAMg7zQ522D43SEAIALBNOHqHdLSExQgBAVKNHCACwThgxwlDlDx8+PH/+/CNHjgwfPrxHjx7GAj6fb9GiRQUFBV27dh0xYoTWJX3vvff2798fKJOWlpafn1/r+ekRAgCs5DD/pVBeXt6nT5+1a9cKIQYNGrR8+XJjmUmTJj311FNNmjSZPn36XXfdpe0sLCz89v88/vjj77zzTl2XoEcIAGi85s+fn56ePn/+fIfD0aJFi8ceeywnJye4wN69e19//fVdu3alp6ffeOONHTp0uO+++9LS0v7rv/5LK1BRUZGRkTFhwoS6LkGPEABgGYfjt+dlTFD2CT///POhQ4dqo51Dhw5dtWpVdXV1cIEVK1Z07do1PT1dCHHuuee2a9fuq6++Ci7w7rvvJiUlXX755XVdgoYQAGCdMAZGlWOjJSUlaWlp2uuWLVv6fL6ysrLgAmVlZS1btgxspqenl5SUBBeYPXv2xIkTnc462zuGRgEAltE6hObe4nCsX79+9OjRwTuvuuoqbTDT5XL5fD5tp/bC5XIFl3Q6nYECQoiamprgArt27Vq5cuWcOXMUFaAhBADYLCMj4/rrrw/e07Fjx8Ch0tJS7XVJSYnb7Q50EI0FhBClpaUZGRmBzdmzZ+fm5mZlZSmuTkMIALBMGKtPCCEyMzOlhjBg6NChzz777AMPPOByuZYsWTJkyBC32y2E2LZtW3Jyclpa2uDBgydMmLBz587zzjvvxx9/3Lt3byAc6PP5/vWvf/39739XX50YIQDAOlbHCEePHu1wOK688srJkyfPmDHjgQce0PaPHTv2zTffFEKkpaXdddddubm5d91119ChQ++5556kpCStzCeffHLixIm8vDx1lekRAgAar7i4uNWrV7///vvHjh277777AoOcL774ovakqBDi8ccfHzp06I8//jh69OhLLrkk8N527dotX748NjZWfQkaQgCAZRwiRDpEbW8JIS4ubtSoUdLO3r17B28OGDBgwIABUpl27drVpwI0hAAA65iPEdo95zYNIQDAQhG4HiEPywAAoho9QgCAZcJIqLd9bJSGEABgmTDyCO1uBxkaBQBEN3qEAADrRODDMjSEAADL/Lawksm3NFBl6omGEABgmQjsEBIjBABEN3qEAADrhNEltLsLSUMIALBOGDHCBqpJvTE0CgCIavQIAQCWcZhPkLf7oVEaQgCAdRwO8+kQdreENIQAAOtE4MMyxAgBAFGNHiEAwDIR2CGkIQQAWCecKdbsbgoZGgUARDV6hAAA60Tg2CgNIQDAMmHkEdreEDI0CgCIavQIAQDWicC5RmkIAQDWIUYIAIhmDkfkzTVKjBAAENXoEQIALOMQDrMJ8hGTUH/ixIn777+/f//+AwcO/Nvf/lZTU6PtLy4uvuWWW3r16jVx4sQDBw40WD0BAJHAYf7LbvVtCO+5555Vq1a99NJLM2bMmD179nPPPaftv+6665o1azZ37lwhxJgxYxqolgAANJD6Do1u2LDhtttu69atmxBi9OjRGzZsEEKsXbu2qKho9erVHo/n+eefT01NLSgo6Ny5cwPWFwDQiIXxsIztncL69gjz8/MXLFiwa9euH3/88b333svPzxdCbN68uUePHh6PRwjRpEmTLl26bNq0qQErCwBo3By/Tbttjr11rm9D+Kc//amioqJnz54XX3xxu3btrr32WiHEgQMHkpKSAmWaN2+uCBP6fL7TrCsAwEaBp0NUzuIY4ahRo/r27Xvw4MFDhw653e477rhDCNGsWbOTJ08GypSXlycmJtZ5JSepGgAQwVwul91VaBD1apz8fv/KlSt///vfO53OmJiYG2+8ccWKFUKI7OzsoqIirYzP59uxY0d2dnaDVRUA0OiFMzZqc5Xr1RA6HI4OHTp88MEHQgifz/fBBx906tRJCDFkyJDjx48vW7ZMCLFo0aLY2NgBAwY0aHUBAI1ZBI6M1vup0X/+859/+MMf5s6dW1lZmZaW9tZbbwkhYmNj58yZM378+KSkpPLy8vnz55+tHWcAwNmqvg1h3759i4qKDh486PF4gh+QGTZsWHFxcVlZ2e9+9zu3m3lqACC6NcCk2/v3758zZ87x48eHDx/et29fYwGv1ztv3rzCwsIuXbrcdNNNgUdSvF7vwoULt2zZ0rx585EjR7Zr167W85t7gCU1NTW4FdR4PJ6srCxaQQBAWMkTqpbw2LFjffr0KSoqSkxMHDZs2Icffmgsc+utt77yyiuZmZnPPfec9iynEKKqqio3N3fmzJlJSUm//vrr119/XdclaL0AAJYJY4V6dfl58+ZlZ2fPmjVLCJGYmPj4448PHTo0uMDu3bsXLly4Z8+e1NTUESNGtGvXbtq0aenp6TNnziwvL//mm29C9tNIaQAANF4rVqy48sortddDhgxZvXp1VVVVcIFVq1Z17do1NTVVCJGZmXn++eevXr1aCLFs2bKxY8d+9NFHL7zwwvfff6+4BA0hAMA6YTw2quwRlpaWpqWlaa9btmzp9/vLysrqKqCVKSkpEULs3LnzpZdeevvtt/fv3z9kyJB//etfdV2CoVEAgGXCSQx0ONavXz969Ojgfbm5ubfddpsQwu12B2a08Xq9QghtXs8At9sdPHNZdXW1VsDhcFx44YXamhBdu3b961//Onbs2FqvT0MIALBZRkbG9ddfH7ynS5cugUNaD08IUVxc7Ha7g/t/WoHi4uLAZklJSUZGhhAiMzPzggsu0HZeeOGFxcXFXq+31nghDSEAwDJhrD7hECIzM1NqCAPy8vJmzJjxwAMPuN3ud955Z+jQoVrC+qZNm1JTUzMyMnJyciZMmFBUVNS+ffstW7aUlJQMHDhQCHHttdd+/vnn2knWrVvXrl27up6aoSEEADReI0eOfPHFF6+44oq2bdu+//77H3/8sbb/j3/845gxY6ZMmZKSknLPPffk5OQMHTp02bJl06ZNa9asmRBi4sSJc+fOvfrqq7OyshYvXkyMEABwRpiPEaqLx8bGrlixYvny5b/++uuTTz6Znp6u7X/ttde0J0WFENOmTbv66qsLCgr+9Kc/9ejRQ9vZrFmzdevWffLJJ6dOnXrwwQe18dJa0RACABq1mJiYq6++WtrZtWvX4M1evXr16tVLKhMbGzt8+PCQ56chBABYJhJXqKchBABYyu6GzSwS6gEAUY0eIQDAMmEk1Nu+Mi8NIQDAMuHkEdo9lEpDCACwjEOEWFap9jfZihghACCq0SMEAFgnjBXq7UZDCACwDjFCoLHzK7YaNflvhd1/O4CzBg0hAMAyDvPpEPQIAQBnEWKEgA3MjXYqj5/Wm80J/bdCX0K+tF91AsZRYZdIzCMkfQIAENXoEQIALBPGFGu2dwnpEQIAoho9QkQAf4jQnCpIaOqtoY+Hqkr9+Q3/BcvbfhOHpYNyAFEZUBT2/0cO2ImGEABgmXAelmmYmtQfDSEAwDIOEXnLMBEjBABENXqEaARCR/XMRAGlo8rD8lHDtf1mQo4hgnoSh3wuefEaaUsdUlRezCFVzVDYrw8ikoaI8IWRUG/3B4yGEABgHfMxQtvREAIALBNGHiExQgAA7ESPEDYIkYynDusZtv1+VdjP3FGfsmKG8ubmHpVjfsb/gvWpgVJ5pz6M55AKS0elS4eKV0qhU6m8FGGMtIEvQI2GEABgmbCGRhuoLvXF0CgAIKrRI0TDCDGCaDKHQT286dNt+3yqwvJR5aYQQt6hrIlayOFKp3osVL/pVI6Uykf1//HW8g+4vEs1SEuuBRTCmFnG9s8MDSEAwFJ2N2xm0RACACzjcJhOhyB9AgAAO9EjhEVCTHvmVxxVx/yMe3zqzRozR5WFjXvkmOJpxAilmJ8wBPbkTZdyU11YeWZhDCIaCuiOKnMt5JnjIm2UDKcrjNUn7P6Q0BACAKwTxlyjdmNoFAAQ1egRAgAs4xAOeR2V+rzJVjSECJepoKCZic1CBupqanRvqPHqj3rNHNWfSgoZCiFqapQxQp8q9innDSrzAoUQLimS59Ztuly68RuXW9p0mDjqkoeCnH7VpjxyJAUFTYUMhe1/9NCwwlmh3u6PBA0hAMBSEdYhpCEEADRuxcXFr7322pEjR/Lz8wcOHGgsUFVVNWvWrIKCgq5du44fP97tdgsh1q5du2LFikCZSZMmJSYm1np+HpYBAFjGEQblCY8cOdKnT5+DBw+2bdt21KhR//u//2ssM27cuAULFlx44YWzZ8+eNGmStnPFihULFiz49f/4fHUuLkOPEPWjni1UmAsKmsrtE4bAnlcf9vNW649W1+g2q3RHq+XCUkBR/lVRBxGl6KZfmqJTSEFBXWEp1U+EigK6PbpNj37THaMq7Pa4dPV0y3dYupZfSkMUUhqi9G4TIUNBouHZzvK5Rv/1r3917NjxpZdeEkLEx8c/+eST11xzTXCB7du3L168uLi4ODk5+ZprrsnOzp4+fXpGRoYQok+fPk8++WTI69MjBAA0Xl9++WVOTo72OicnZ82aNVVVVcEFVq9e3b179+TkZCFEenp6x44dv/76a+1QYWHho48+OmfOnCNHjiguQUMIALCOI6yvupWWlqampmqv09LS/H5/aWlpXQW0MlqBtLS0jh07er3et956q1OnTrt3767rEgyNol5qmUrMurFQKcPBaxiflAYwqyt1g5/VVbrNKulopXRUPzSqf680cCpCJWMY0ifUQ6NShoNhaNStHvzUDW/GxOoLx0pHdZueGP391x8VQkhjpW6/iX+RTY2UCmN+Rf2vhIgQxsK8wrF69erc3Nzgnddee+0dd9whhPB4PF6vV9tZXV0thIiNjQ0uGRMTEyiglYmJiRFCjB8/fvz48drOkSNHzpgxY+bMmbVWgIYQAGCZMGZYcwjRvn37u+++O3hnmzZttBeZmZnFxcXa6+LiYo/HE9z/E0JkZGTs27cvsLlv377MzEzpEpdccslnn31WVwUYGgUA2CwtLS1HL9AQXnPNNYsXL9bigv/+97/z8vJcLpcQYt26dXv27BFCDBkyZPv27QUFBUKIjRs3HjhwQEuxOHnypHYGr9f7/vvvX3jhhXVdnR4hAMAy9UiIqOU9ioMjRoz4xz/+0b9//+zs7BUrVnz66afa/j//+c9jxoyZMmVK8+bNp02blpubm5OT88knnzzyyCNNmzYVQlx22WXNmzdv0aLFhg0bmjdvfu+999Z1CRpC1M6wwFCI7dMKClZLMT85UCeF/aoqdJuVFV7d5ilpU3qv7qgcUDRcWp1fIc8Gp0wMcMoxQnk8Rk6Q0GdESGG/mDjdL29svO6oVz4aYrUpv08fNZR//HJMUUEOGRrTJ6T8Cr8+ZEjM8CxgafqEx+P57LPPVq5ceezYsZdffll7OlQI8frrrwde33333Xl5eQUFBffee2/Hjh21nR999NG3335bXl7+l7/8pXfv3ormmYYQANCouVyuQYMGSTs7dOgQvNmlS5cuXboE70lJSbnyyivrc34aQgCAZcKZdLthalJ/NIQAAMs4HMJ0+oTdLSENIf5PLamCQQeNU6xJ0TFpKSUzQUEpt08K4wlD2K9C2jyh2zx1Un/0ZLX+VKr4opR0KAxRwwaNEUpBQSk1MFYO++lqHlfpCd70JqhXm5J/8eWESL/6L4MqZCivRRfqb5z8R5AJ2CJdePkTtiJ9AgAQ1egRAgAsE4EdQhpCAICFrM4jPANoCPEb9bJKtSzD5DOROChNHyoF3qSgoBQCFEKcOqGL850q1xU4WV6lP6ovfEIKGaqSDqtqiRHq9kjBTp8+9ibdM+m32ymvsiT/8ns8UqagPkYYr/ttjUvQbUrRTW+1LmRYE2qhK2WAWDgcbmlbv+XUb+pO5jT8iXPIu1TLNtneV0A0MBEj/P777/Pz8zt27HjFFVesXr1a27lr166RI0d26tTpxhtvLCkpaZhKAgAih3VLT5wZ9W0It23blpubO3DgwCVLlkybNi0hIUHbn5+f37Zt2w8++CAlJeWGG25osHoCACKAlkdo9ste9R0anT59+k033TRlyhQhRGACm6+//nrv3r2PP/642+3+29/+lpKSsmXLFsXEpmhU1JOoqWdQE6GGRuU15fUTlUkjkNL4pDQQKoQ4eVy358Qx3VjoCf3Rk8d1R0/qR0qlXAtpGFaavE0YMj0MQ6NS1oF+GSb977e0JL1xaNStHBqN0w+NVlboh0ar9OkTZuppJNVMGsyUNqWRTrmw4W+cXz67coF7JmCLNGHMNWo6pmi1+vYI169f36FDh/Hjxw8dOvS5556rqakRQmzdurV79+5ut1sIER8f36VLly1btjRgZQEAsFp9e4R79+596qmnXnzxxcTExEmTJv36668PPfTQwYMHExMTA2WaN29+4MCBus7g88nTGQMAIkhNTY22BJJK4wj7mVLfHmFiYuKtt946bNiw/v37P/zww2+88YYQIikp6cSJE4Eyx48fb968eZ1XMqxjDQCIIKFbwbBihLY3nPXtEbZt2zbQ+UtKSiovLxdCnHfeeYWFhdpOn8/3888/B5ZSROSRo4L6LUNQSQ4K6hMJpBhVtXIdJSlQJ2VHCENQsPyoavOEFCPURxDVM64ZY4TVyhihYWYyVYzQIU+xZi59wrCAlD4oKC8XpaqnkTqcKSV+SHPFSZsOl64mUnKFMAYRQ6RuKI8CVqhvL23s2LFvvvnmyZMna2pqZs+era2IkZOTU1lZ+c477wgh5s2b17Rp0379+jVgZQEAjVsEPjRa7x7hxIkT169f37p1a4/H061bt7lz5wohYmJi5s+ff8stt9x5550ej+fNN99k/BMAoloEzrFW34bQ7Xb/85//fOGFFxwOR1xcXGD/4MGD9+7de/To0aSkJNsfgQUAwCxzU6zFx8cbdzqdTsUzMmhETiNx0BhkMiQOqhZakiYqk5dVOqGaMk0YMgUNMcJKXeFjqjxCacY1Q+DNmEcoLWCkyp4MMcWaUxWHE0K4PboBlZhKfYywUvfbaqiYfrI3ZVDQMMmZIcdR3lRNDidtyvFFp1wTp3TTzKQVyqX537sxMp9HaPcPkrlGAQCWCSfoZ/c/NDSEAADrRGCMkGdbAABRjR5hFAmx9I4ycdBnjBHWSDFCM5OL6iNzp05KeYTGuUZVmYLy1KPH9FOPlqvyCKXEQWl9KGFYQMoQitMVlu6RYcZO3aaUnCeEcLt1e7xVqlsqTx+qnLgpZAqjW79HmvXUHeNUbXqkCKKuKi5DKNSnjxqaSiskRNj4hTPXqN0/SRpCAIBlInBklKFRAEB0o0cYvQz5EtKmKjFAhEyf0A/rVcvrLumnWDupynAQhuFNadY0+aicLyHNqabKl/AahkalqeMM6xnpCvv1N1Ua8DEsWG+4ln6wtKZGt6lO1ZCoh2GlwUwhhCdGNxbqidXdtJhY/dxvsfqbpn+v26PPMPEYhuRDfNJUo5+2dx0QWgR2CWkIAQCWcYSxviANIQDgrBFGHqHd7SAxQgBAdKNHeFYLkTChOnqa6RPVcvqEflOfPmEIGcoxQmkONnVM0bDQklexKeVLSN+FEKJGzlJQB+pUCRTSv8mOGkNegUs5f5t61jQ5ACmlQ+gnb4uVV5WLjdPdlth43V+GqngpsKo7GhOnu2keOeFE/m9b+r6c8hxr6iQUZlxr9IgRAgCiWhh5hHYv2MDQKAAgqtEjBABYhkm30bioQ4LqOdXkYFjIGGGNaj4wr34ZJvWMa1IYTwhRIcf59GeTllKqUF2rWjlvmTFGqM7e86uz+aS0QmmBIcMt9ft0BUJlCuq+L2mNJ7dHumO6oGBFnHyH46Q7LK9OpV8BSv/TlO+hMvNShPpo+dUzrun/YqrDiUA90RACACzjiMD1CIkRAgCiGj1CAIBlwkmot3tQm4YwmpiIZ4UIGQpD8Mynn2tUnVZYrZyJVAryGfeYiwIqMwVr5KQ3w7epnnM1RIxQX1bKeTP+8hsurijvrdJtV7ukxE3lHQt1h6uVt7RaGViVPgnGDNQQuZhmPqV2j6ihFuEsw2R3S0hDCABo1Hbt2vWPf/zj+PHj11577ZVXXmksUFFRMXPmzMLCwi5dukyaNCkmJib46MKFCw8dOvQf//EfdZ2fGCEAwDoO819Khw8f7tu3b01NTY8ePW655ZZ33nnHWObmm2/+4IMP+vfv/84779x6663Bh7755ps77rjj8ccfV1yCHmH0MoxIqcakjNN7ycswyasyqUbPvOqR0mp54M6r3yMVkFM1vKqn+aUp0+QBXuOcaXI+hYm1kCQO+b0hlm7Xj4zKNXfWqL5r+Q7Ld0y+w1IBU0km8h32hRgaVU4VJ38O/dI0daq3onGwOo9w7ty53bp1e/rpp4UQHo9nxowZI0eODC5QVFS0bNmysrKyxMTEq666Kisr67HHHsvKyhJCVFZW3nHHHdOmTXvyyScVl6BHCACwjOO3x2VMfamsXr160KBB2utBgwZt2LChsrIyuMDXX3/do0ePxMREIURqamqnTp3WrFmjHXr44YdHjhzZqVMndZ1pCAEAjVdpaWlKSor2OjU11e/3l5WVBRcoKysLFBBCpKWllZaWCiG+//77jz/+eOrUqSEvwdAoAMA6Ya0+sXr16tzc3OB9w4YNmzJlihAiNja2qqpK26m9iIuLCy4ZGxtbXf3/15yprKyMi4vzer233nrriy++6PF4Ql6fhjCKGMJffsVhefKwWoJnqomy5AiiMrlCvVnb26XZ3fSXVk7xpX52v7aYX/hBQflE+vc65NnD5KihOoPFZ+a7Vt/A2gqof0CqmoScnE/IGSnKlZWkwg5Cho2dw3xeoEOI9u3b33333cE727dvr71o1arVvn37tNd79+6NiYlJTU0NLpmZmRkoIITYt29fZmbmTz/9tHnz5ptvvlkIcerUqQMHDrRt23bFihVa7FBCQwgAsFlaWlpOTk6th/Lz8++///4HHnggLi5uwYIF11xzjdPpFEKsWrWqVatWbdq0yc3NHT9+/ObNm7t27bp27dpffvnliiuucLvd27Zt086wevXqv/71r59++ml6enqtl6AhBABYxuEwnyCvLH/ttde++uqrffr0yc7OXr9+/fLly7X9U6dOHTNmzJQpU5KSkh599NEhQ4ZcdtllK1eufPLJJxMSEoQQbdq00Uru2LHD7XYHNo1oCAEA1rF6hXq32/3hhx+uW7fu8OHD8+fPb9asmbZ/4cKFgddTpkzJy8vbtm3bjBkzzjvvPOkM/fr1+/rrr1WXMFlfRI1QwbNQyzapjsrZe8rol3FPiHV8pEurp44LGQo9jaCgmvHMDrmAKu/QcP91m/LqUaHvsD7EqPwBmfpZh/zkhFotDBEmnLlGQxVwOp19+/aVdmZnZwdvtmvXrl27drW+PT4+vnXr1qrzh6oAAABnM3qEAADLROKk2/QIAQBRjR4hfnOakRr1zKXyRKbKKFFtUSVV5pkyH9I4h2pjDkmpFxlS3VL18lDqG2gsH+IHZKYmZjXmHw/OVjSEAAArWTvp9hlAQwgAsEw4McIGqkq90RDiN6f5WXTImw7FprG0bstQFen3yrBpoibq8Ua7qX8Iqlsqf5PyLVXdQGP5ED8gMzUxy/a/iThdVucRngE8LAMAiGr0CAEAlonE9AkaQgCAZcKYWcZ2NISoQ+i4nX7TqQ8UOVVHnepNlzxi73RJJ1e9Xb608huR/hX1G79N65Zhks9s/GMRIrCnerv0XRvuiXSH5WtL91z9AzL1sw75ybE9PgQQIwQARDV6hAAAy0RijJAeIQAgqtEjjCJyAp0cZdKFv+T/0YyZZ+oolH7T5ZY2ncpN+VqGt6siXnK4y6VbJchRo6+2T/qupSvLix85HOGHDI23UF1AHYUN9V2r73/IO2ziB6SO4EqbQsjfmOGTpiwsnwuNkaXr8p4JNIQAAOtEYEI9DSEAwDJhxAhtR4wQABDV6BFGr1BzcurCX8ZAjxwFlIOCqiCT26Pb9MToNz0u6Vpu/R6pgHQ2t/5aNfrgmc+lD/L59N+YIebnc0pBQSlkFX6Q0PhPs1O/S77DUtjPpfqu5Tss3bEY+Q579Hukn4h0NvmH65I2Vd+FqO2zFMzamUtx5oWRUG97B5KGEABgnQiMETI0CgCIavQIo4l6ASIzz+4L4yxcyufvPeqx0FjduFxMnDxwJ+2RNyt1m94qXb5ETY1+aNSnzn/wCeUOw8rvcnEF9ZRpopacE+XwcozqlsYob6l01LhH+okYxq5VFZPzOgyfHPUMbWaWokJj5KglzhLyLUy6DQA4WzjCyCNsmJrUHw0hAMA64cUIbV0hmxghACCq0SM8m6ljgqaCgrXECNVP80spDTGqGFVsnO5zGBsvfyzj9HuqKry6zUrdUW+1LqwnB3jk6S0AAB1ZSURBVAV9qv88jUM6NTWqt/vldZtUcdeQYVc5QcJ9GlFA+ZbqjsYZ7rB0z6WfiHQt6acpJ1coJ3sTtazTpAydhsi1QOMTRkK9w1wikuVoCAEAlgknj7BhalJ/DI0CAKIaPUIAgGXCWo+wgepSXzSEZzXDwkvStqK0tN5QLRNlqQNaHimgpd+UI1i6mF9cgvyxrKzQxwgrPcGbclCwxsRKSQ5nje5UVYZvs8bEyf36WyrPFqY/tzF4JsdZlUFBKaoXl6C7J/FNVEfjmhiisAmqmGKMHDIMP61QGD5LIVenUh21+w8ozg40hAAAy4QRI7T9HxpihACAqEaPEABgGYcwHyO0u0tIQxi91KEZQ3aXHGozTIyp25RnwpTjW/q0Nn3Mr7pSF7cTQlRXqYOCuk2/YbrQYHKgTv9dVLvkN3u9Jq6lnL1VnlTT6ZLHY9zqxEF9pqAU9ktoqtw8R3VUCBHfRB9E1J9c+nl55LRCKUYYcq5R/aYyj9DuMTOYd9avPuH1en//+99PmjQpsOenn3666qqrzjvvvOHDh+/evdvq6gEA0LDMNYQzZszYsmXL2rVrtU2/33/dddddfPHFa9as6dix4+jRoxughgCAiKE9LGP2y14mhka3bdu2aNGiv/zlL88//7y256uvvjpw4MC0adNcLtcjjzySmpr6/fffd+/evWGqitMVIplCKqwcQhSh0ifccvqEfqUk/aP/0sJJ3mp54M7r1Q06qqc9kxgGJPVDuB5ptjZ5VFadm6Fe1El9D415BeqbJuUwSAkS0mhnk2Yxus1zdJsJ+k0hRHxT/cnlbArljGse3Wbo9AlTc6oxUhpxwphize4fbH0bQp/Pd9ttt7300kvB458//PBD9+7dXS6XECIuLq5z587angapKQCg8WuAGOFPP/300ksvHTly5LrrrsvPzzcWOHHixHPPPVdQUNC1a9c777wzLi5OCLFu3bp333133759ycnJo0aNuuyyy+o6f32HRv/+97937969X79+wTt/+eWXc845J7CZlJR08ODBus7g8ymfYQAANG62/Bk/ePBgv379mjVrNmTIkNtvv33BggXGMmPGjFm9enV+fv7y5cvHjx+v7dyyZUtycnJeXl5GRsawYcPef//9ui5Rrx7hnj17nn322WXLlu3YsePAgQOVlZU7duw477zzmjdvXl5eHih27Nix5OTkuk7idJKzCAARrD5/xi3vEM6dO7d3796PPPKIEMLn8z399NNjxowJLrBt27bly5fv37//nHPOGTRoUEZGxu7du1u3bj1x4sRAmX379r3//vt5eXm1XqJeDWFJSUlsbOyIESOEEOXl5UeOHMnNzd26dWubNm22bdvm9/sdDkdNTU1RUVHbtm3rc0LYQxkklNdBkQI5hs9/iPQJfdzIE6M7eWycPuanT1GQQoDCOLFZiKCgFIpTBS8rT+nqWVVhyNyo1u2p8apXZdJPsSbfQ9UdE0J4PKqVleSckwRVgoQUFGyaKIUM5ShsQlNdAWkONmlVJo96GabQ6RPKZZik0gQJI43D8DMN/RZl8W+++WbgwIHa64EDB44dO7aiokIb/NSsWbOmR48e2vBkcnJyly5d1q5d27p160CBffv2ffXVV//5n/9Z1yXq1Uvr27fv9v/zzDPPdO7cefv27fHx8Tk5OT6f74033hBCvPbaa8nJyZdeeml9TggAQH2UlZW1aNFCe52SkuL3+8vKyuoqoJUpLS3VXi9YsMDhcGRlZV144YU333xzXZcwPVzZtGnTzMxM7bXb7X7rrbemTZuWlJT07LPPvvnmm+YfFgIAnEUc4XytWLGil96MGTO088XFxVVVVWmvKysrhRAJCQnBF4yPj6+urg5sVlRUBAqMGTPG7/f/9NNPP//884MPPlhXlU3PLHPNNddcc801gc3+/fvv2LFD6qgCAKJTeHmB3bt3f/rpp4P3BHpcrVq12rNnj/Z6z549cXFxKSkpUslAASHE3r17W7VqFVygffv2kydPfu6556ZPn17r1a2ZYo1W8GxgJq1QGGOE+lCc363Pt9Nnnknpd74at+KoMMTe5LrJNZGilarkvIo4fR6hIUboVcYIDasyqWKEhsxLYx6hKkYYJ627pMwjlDIFpaCglGUoDHmE0rWkmshTrLlVn4RaYoRmEgcRccJaj9CRlJTUs2fPWo+OHDly6tSpDzzwQEJCwvz58/Pz87Vndj799NPWrVuff/75Q4YMmTBhwrffftuzZ8+vvvrq6NGjWkzxwIEDaWlpQgifz/fZZ5+df/75dVWAuUYBAI1XXl7e7NmzL7roonPPPbegoOCzzz7T9t9///1jxow5//zzmzVrNmPGjKFDh/bt2/ebb775n//5n/j4eCHE1Vdf7ff709LSCgsLzznnnKVLl9Z1CRpCAIB1rM6fcLlcS5Ys+f77748dO9a7d+9A/G/JkiVNmjTRXk+ePDkvL6+wsPDll1/OyMjQdq5Zs2bLli2//PJLRkZGp06dFP1UGkIAgIUcZpdVClne4XD06NFD2vm73/0ueDMrKysrKyt4j9vtNr6rVjSE0Uv698hvJq1QGOfw1BeXpx7VH/X7XPpNfaSt9voG10UVezMsZqTPxtMHBeNOhYgRqvMI5VWZ1HONulT5diJ0HqE+jJegmno0Xl6VKUZ/VP7Fl5ZhCjW5qPSNqIKCjlpihCYSB3kOHWcADSEAwDKOCPz3hYYQAGCdCFyYl4YQvzHMvybvkMtLM5lJpfUDrW6/Uzqs3wrxOZSqYsjcUM7uFqu7tDToV1mhGxqtNizDVK1fIkqaDU7O9FAuUW+YlE6ezkJakl7KUpDmOVPPuCatoyTlWkjZEcI4FirPqaYf6FYutCRtGodGDfkS6pFSRJjw0icaqDL1xETYAICoRo8QAGCZMGaWsX0YgIYQAGAdYoSIYMqYYOh/8eRUAfUbXIpjDof8sVSvrCSnasSo5lSritfnS1TqriVFBIUQ3uozFyOU0hI8ym8kRhkyNAQUVVOmGU8uBwX1YVd5oSUpKCitm2QMEao/Gnb/TUQUoiEEAFgmjIdlbP/vh4YQAGCZCBwZ5alRAEB0o0eI2qknYBOh5mBzmvsXSx+yqmU6N+nkqvWMDOsu6VIDpaCgt0p3tLpajhFKQUFpijXD5HD6ZZik2cKU1RaGqKHHI8U+9THCWFXSoSHmp1pHSRhumpQpqA4KynOqhc4LZBK1s12k/UxpCAEA1onAhHoaQgCAZcLJI7S7B0mMEAAQ1egRol6M/7H55VCQKofOVMjQ4ZBLS9NVmooReqXQWpwu5qdOExRC1OgXWvLVSDFCXeFQMUJdYadhEk6XHJmTkvnUEURVYSkR0B0qPKmePjREUFA5laiIvPgRTArjsVG70RACACzjcJiP+dk9NsrQKAAgqtEjRP0YJ8qSJhc7jZFSw9haiCWfHC7diKQ0cOdy6466PbqzeaR0CGnk0ytfusanzJeQ0yf01ZY2Qy3d7pKyRPQDmC7lxHLS+HCIoy7D4LMyQcUR4uclRN3btfyrH2njZjAljIV57e4Q0hACAKzjEA6HyX927G4HaQgBABaKwDnWiBECAKIaPUKES1ptx1TIUJ7ATbdl/O/MsLKP/ll/p+79Ln3I0KePEfpqnPpN/VGfHCP0qYOCcnEV+bswxAjVgTpTKQ0hjhourY5fyjWX36wMCtr9zz7ONPMJ9bajIQQAWCaMZZgcxgfkziyGRgEAUY0eIQDAOhH4sAwNISxyOiFD6UyGCJZfGV1z6uN2PqcqjBci5meMEUo7/JbFCI2BFKfy21SvRSWdzJAIGGLlI3ksSw7KysXrLmv/HzXYK5xJtw05uGcYDSEAwDJh5BHa/t8TMUIAQFSjRwgAsFSEdQhpCNFAlJlnfr86ZGWIF0ghRjkAqS+rj4f59XE8p7602bxA6WzmIhty4C3Eb78679CQW2niaC0zx6qvbeIgoh0L8wIAEGHoEQIALBNOQn0DVaXeaAhhA/nXRJlrIUKlW4QYOPWrLuaXUhYM1/EbKqeoiHoZJlkt45OnkcMg52YoLxVqbSQyIhA+8ggBAFEtjLlG7W4IiRECAKIaPUIAgGUc9Xgi2vAWm7uENIRoBEI822+I88lb6giiajvkpPd+6f2mJlVTUwf9ailv5qDZc9s9NgUobN269fnnnz969Gh+fv6YMWOMBY4dO/bUU08VFhZ26dJl6tSpTZo0EUKsXbt28eLFO3fuTEtLu/XWW7t161bX+RkaBQA0Xvv377/sssuys7NHjx599913v/7668YyN9xwQ0FBwfjx4zds2DB27Fht5yOPPBIXFzd69Ojk5ORLLrlk06ZNdV2CHiEAwDLhrUeoODp79uz+/fvfd999QoiqqqoZM2bccsstwQV++OGHVatWHThwoEmTJpdeeml6evqOHTvatGmzdOlS7cwjR478/vvvlyxZUlenkB4hAMAy2swypr7UI/fr1q0bMGCA9nrAgAGbN2+uqKiQCvTo0UMbDk1KSrrgggs2bNgggtpXv9+/b9++9PT0ui5BjxARIMT/l9KsaSEijn7VYePkbvKlLAumhT7RacQMSQSEbcLII1QqKytr0aKF9jolJcXv95eVlWVnZwcK7N+/P1BACNGiRYvS0tLgM7z44ovHjx+X+pHBaAgBADZbsWJFr169gveMGDFCGw5NSEgIdAFPnTolhNA6fwHx8fGVlZWBzYqKiuACCxcufOKJJ7744ov4+Pi6rk5DCACwTFgxQtG9e/enn346eOe5556rvcjKytq9e7f2evfu3fHx8cH9P6mAVqZVq1ba68WLF0+ZMuXjjz/u0KGDogLECAEA1jEfI3Q4RFJSUk+91NRU7XyjRo16++23y8vLhRBz584dMWKE0+kUQixbtqygoEAIMWTIkJKSkjVr1gghvvjiixMnTlxxxRVCiI8++mjSpElLly7t2rWrusr0CBH5QqUh6kgBReWpQl3qjCLsh+h09dVXz5s3r2vXrpmZmfv27Vu+fLm2f/r06WPGjOncuXPTpk2feeaZvLy8Hj16fPfddy+88EJsbKwQ4vbbb6+srLzhhhu08qNHj37iiSdqvQQNIQCg8XI6nQsXLiwsLDx8+HDPnj1jYmK0/R9++GFcXJz2esKECXl5eUVFRR06dEhJSdF2fvfddz6fL3AerXWsFQ0hAMAyDhHGMkyhyxuDfMnJycGbaWlpaWlpwXsSExPrWQEaQkQZU+OoAEwKY4V6238PeVgGABDV6BECAKzDwrwAgGhm+VyjZwBDowCAqEaPEABgmTAelrF7ZJSGEABgrUiLETI0CgCIavQIAQBWqk+CfKNCQwgAsEw4MUK7200aQgCAZcJJn7C7B0mMEAAQ1egRAgCsE4Ezy9S3R/jEE0/06tUrLS2tZ8+eCxcuDOzfunXr5Zdf3qJFi5ycnKKiooapJAAgMoSxKm/ENIS//PLLc889V1BQ8MADD4wfP37dunVCCJ/PN3LkyKuvvnr37t39+vULrH8IAECkqG9D+PTTT/fr1y8lJSU/P793795r164VQqxaterIkSNTp05t2rTpf//3f2/fvn3Dhg0NWVsAQKPm+O1xGTPs7hKafljm8OHDW7Zsueiii4QQP/74Y9euXZ1OpxAiJiamU6dO27Zts76OAAA0GHMPy3i93rFjxw4fPrxfv35CiMOHD59zzjmBo4mJiYcOHarrvT6fL+xaAgBsV1NT43K51GUiMY/QRI/Q5/ONGzeuurr61Vdf1fYkJyeXl5cHChw9erRFixZ1XslJqgYARLCQrWCEqm+P0O/3T548ec+ePR9++GFsbKy2s127dj/88IPf73c4HF6v96effmrfvn2DVRUA0OiZT6i3O0RY7x7h5MmTv/zyy5deemn//v07duw4fPiwEGLQoEFut/vVV1/1+XzPP/98enr6xRdf3JC1BQA0bg7zX3arb0P41VdfVVZW5ufn5+bm5ubmzp49WwjhcrkWLVo0c+bMpk2bvvnmm2+99ZbtCw0DAGBKfYdGt27dWuv+Pn36bNq0ybr6AAAiWCQ+LMMUawAAy4Qz6bbdLSENIQDAUnb38MwipQEAENXoEQIALNM4ngM1h4YQAGCZSIwRMjQKAIhq9AgBANaJwIV5aQgBANZxmB7qtH1olIYQAGCZcBLqG6Ym9UeMEAAQ1egRAgAsRYwQABC1HMLhsL1lM4mhUQBAVKNHCACwDKtPAACiWwPMsbZx48ZnnnnmyJEj11133YQJE4wFfv3118cee+yHH37o1q3bfffd16xZMyFEWVnZ559/vnHjxhYtWtx7772K8zM0CgCwjCMsihOWlpYOHjz4oosuuv322x999NFZs2YZy1x//fXFxcVTp04tKiq6+eabtZ0fffTR/PnzCwoKFi9erK4zPUIAQOM1a9asgQMH3nXXXUKIJ598cvr06RMnTgwusHnz5jVr1hw8eDA+Pr5Xr14tW7b8+eef27VrN27cuHHjxs2bN++FF15QX4IeIQDAMlqM0NSXeih1w4YN/fr1017369fvhx9+OHXqlFSgR48e8fHxQohmzZpdcMEF3377rak60xACAKzjMP+ltH///uTkZO11ixYthBBlZWV1FdDKSAVCYmgUAGCzTz/9tG3btsF7br755ocfflgI0aRJk0AXUHtxzjnnBJds0qRJRUVFYPPkyZNNmzY1dXUaQgCAhcJZj/DSSy99+eWXg3empKRoL84999xdu3Zpr3fu3JmQkKD1CwOysrICBYQQu3btysrKMlUBhkYBAJYJY2TUIUSTJk3a6GkpEEKI0aNHL1q06OjRo0KIWbNmXX/99VpD++67727ZskUIceWVVx44cODLL78UQnzyySeVlZUDBw40VWd6hAAA61i9HuGQIUMGDhzYuXPn9PT048ePf/rpp9r+GTNmjBkz5sILL0xISHjxxRfz8/M7deq0bdu2V155JSYmRgjxxRdfjBw5sqqqqrKyMjk5+corr1ywYEGtl6AhBAA0Xk6nc+7cubt37z569GiXLl1cLpe2//PPP/d4PNrrm266KS8v7+effz7//PMDEcT+/ftv3749cJ5AYSMaQgCAZUImyIendevW0p4mTZoEbyYmJvbs2TN4j8fjad68eX1OTkMIALBMJM41ysMyAICoRo8QAGAZh0OEkT7RQJWpJxpCAICl7B7qNIuhUQBAVKNHCACwUqR1CGkIAQDWCSN9ghghAOAsYvXMMmcAMUIAQFSjRwgAsEwkJtTTEAIALNNAU6w1KIZGAQBRjR4hAMA6EfiwDA0hAMAyDvMxP9tHUmkIAQCWcQiHw2QXz2x5yxEjBABENXqEAADrhBEjtBsNIQDAOubzCG1vOBkaBQBENXqEAADLMOk2ACC6ESMEAEQzhwijR9hAdakvYoQAgKhGjxAAYJkInGGNhhAAYKEIbAkZGgUARDV6hAAAyzgcYaRDkD4BADhrsEI9ACCqESMEACCy0CMEAFgmnIR6u7uENIQAAMs4WH0CAIDIQkMIAIhqDI0CACwT1jJMDVSX+qIhBABYJowYod3tIEOjAIDoRo8QAGCdCEyopyEEAFjHfIzQ9iAhQ6MAgKhGjxAAYBmH+Q6e3R1CGkIAgHXCSZ+wO0h4hhpCt9u9beXHK//9+pm5XKQ7dOhQixYtzK/pFb2qqqpOnTqVmJhod0UiybFjx2JjY2NjY+2uSCQ5ePBgamqq3bWwTaf33+/UqZO6zODsFLOnbXOqxzfOinArZQGH3+8/M1c6fPjwkSNHzsy1AACWa9WqVUxMjN21sN6ZawgBAGiEeGoUABDVaAgBAFGNhhAAENVoCAEAUc310EMP2V2HqOb3+zds2PD666+/9957xcXFnTt3drt/y2n59ddfX3jhhXfeecfv97dr187eejZOBw4ceOONNxISEgJPtG/fvv3ZZ5/96KOPUlNT09PT7a1eY+P1et944425c+euX78+NTVVu2l+v3/BggVz5swpKirq2rVr4OMHIUR1dfWbb745b9689evXZ2VlNW/eXNt/+PBh7XdTCMHv5lmAHqHNSkpKRo0adfjw4aysrNdee23w4MFer1cIUV1dPWDAgO+++65Nmza33XbbP//5T7tr2hhNnjx56tSpK1eu1Db37NnTp0+fqqqq5OTkyy+/fN26dfZWr1Hxer1XXXXVyy+/3KpVq+rq6jVr1mj777vvvieeeKJ9+/ZLly4dMWKEvZVsbP7whz/MnDnzggsuOH78eLdu3YqKioQQVVVV/fv337Rp03nnnTdx4sQ5c+bYXU2cNj9sVV1dXV1drb0+fvx4XFzct99+6/f7Fy1a1KlTp5qaGr/fv3Tp0jZt2mivEbBw4cIbb7zxsssu+8c//qHtufvuu2+66Sbt9YMPPjhq1Cj7atfozJw5s1u3boEPm+bo0aNNmzbdunWr3+8/efJkYmLixo0bbapgo1NTU+PxeDZt2qRtDho06Pnnn/f7/QsWLLjgggt8Pp/f71+yZEn79u2114hc9Aht5na7A4NRXq/X6/U2bdpUCLFq1arBgwc7nU4hRG5u7s6dO/fu3WtnRRuZX3755aGHHnr22WeDd65atWrIkCHa69zc3EBPEUKIDz744Oabb/7ggw+eeeaZtWvXajs3btzYrFmzLl26CCHi4+P79+/PTQtwOp0dOnT47rvvhBBHjx7dvn27dqNWrVqVk5OjTfyUm5tbVFRUUlJic11xemgIG5EpU6YMHz78/PPPF0KUlpYG4l6xsbGJiYmlpaW21q5xufPOO6dOndqyZcvgncE3LS0t7dChQ1VVVXbUrjHauXPnq6++unjx4uPHj+fn58+cOVMIUVZWFjxhWMuWLfmbHuzdd9998MEHO3To0Lp168mTJw8aNEjoP2YJCQlNmzbldzPSERhvLB566KGNGzeuWLFC23S73TU1NYGj1dXVZ+XMRuFZtmxZaWnpuHHjpP1ut1uLsAohvF6vy+VyuVxnunKNldPp7NSpkxbQ6tmz57hx426//Xbjx4ypRwNqamrGjh2bn58/adKknTt3/vGPf7z44osHDhwY/DETQni9Xn43Ix0NYaMwY8aMf//731988UVycrK2JzMzs7i4WHt95MiREydOZGRk2FfBxmXevHm7du3q3bu3EKKwsHDXrl0lJSXTp0/PzMwMdGiKi4vT09NpCANatWrVuXNn7XWXLl0OHTqkfahKS0t9Pp82CF9cXNy3b19bq9mIbN68+dtvv121apXb7e7YseONN944Z86cgQMHBn/MDh06VFFRwe9mpGNo1H7PPvvsrFmzli9fHjzQN3z48A8//LC8vFwI8fbbb/fu3ZtkgIDHH3980aJFr7zyyiuvvNK+ffubbrppwoQJQojhw4cvWrTI7/cLIRYtWjR8+HC7a9qI5OfnB0KD33zzTXZ2dpMmTXr37u3xeD7//HMhRElJyZo1a4YNG2ZrNRuRFi1aVFdX79u3T9vcvn17SkqKEGL48OHLli07ceKEEOLtt9++5JJLtP2IXEy6bbOff/65ffv22dnZLVq00Pb87W9/u+KKK4QQI0aM2L59e/fu3ZctW7Zw4cLBgwfbWtNG6vLLLx8zZsykSZOEEEePHu3fv3/Lli2TkpLWrFnz1VdfZWdn213BxuLkyZODBg2Kj49v06bN0qVLZ82apf2jMHfu3LvvvnvYsGGrVq3Kz89/+umn7a5pI/LnP/95yZIleXl5O3bsKCgoWL169bnnnuv3+6+99to9e/Z069Zt2bJlixYt0n5hEbloCG1WUVHxww8/BO9p27ZtUlKSEMLn861cuXL//v39+/dv1aqVTRVs7AoLC5OTkwMPL1RUVCxfvryysjInJ4flCSXV1dVffPFFeXl53759g0fzCgsLN27c2KZNm4svvtjG6jVOW7duLSgoaN68ef/+/ePj47Wdgd/NAQMGZGZm2ltDnD4aQgBAVCNGCACIajSEAICoRkMIAIhqNIQAgKhGQwgAiGo0hACAqEZDCACIajSEAICoRkMIAIhqNIQAgKhGQwgAiGr/D+aCnL6t7ijSAAAAAElFTkSuQmCC",
"text/html": [
"\n",
"\n"
],
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {}
}
],
"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 = 14\n",
"lattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];\n",
"\n",
"# Confining scalar potential\n",
"pot(x, y, z) = ((x - a/2)^2 + (y - a/2)^2)\n",
"\n",
"# Parameters\n",
"Ecut = 50\n",
"n_electrons = 1\n",
"β = 5\n",
"\n",
"terms = [Kinetic(; scaling_factor=2),\n",
" ExternalFromReal(X -> pot(X...)),\n",
" Anyonic(1, β)\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-14) # Reduce tol for production\n",
"E = scfres.energies.total\n",
"s = 2\n",
"E11 = π/2 * (2(s+1)/s)^((s+2)/s) * (s/(s+2))^(2(s+1)/s) * E^((s+2)/s) / β\n",
"println(\"e(1,1) / (2π)= \", E11 / (2π))\n",
"display(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.1"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.1",
"language": "julia"
}
},
"nbformat": 4
}