{
"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.200289e+01 1.454467e+01\n",
" * time: 0.005261898040771484\n",
" 1 6.046576e+01 9.704580e+00\n",
" * time: 0.014676809310913086\n",
" 2 5.600278e+01 1.342225e+01\n",
" * time: 0.036894798278808594\n",
" 3 3.834285e+01 8.853759e+00\n",
" * time: 0.210892915725708\n",
" 4 2.860652e+01 7.242686e+00\n",
" * time: 0.24387001991271973\n",
" 5 2.110059e+01 5.436010e+00\n",
" * time: 0.2709629535675049\n",
" 6 2.077883e+01 6.471797e+00\n",
" * time: 0.29473090171813965\n",
" 7 9.501739e+00 2.286331e+00\n",
" * time: 0.3187689781188965\n",
" 8 6.752361e+00 2.244869e+00\n",
" * time: 0.3413078784942627\n",
" 9 6.242284e+00 1.328768e+00\n",
" * time: 0.3638739585876465\n",
" 10 6.080685e+00 1.415595e+00\n",
" * time: 0.382282018661499\n",
" 11 5.993464e+00 1.201395e+00\n",
" * time: 0.40092897415161133\n",
" 12 5.910314e+00 8.377229e-01\n",
" * time: 0.41904592514038086\n",
" 13 5.840276e+00 7.393215e-01\n",
" * time: 0.437175989151001\n",
" 14 5.797101e+00 5.593923e-01\n",
" * time: 0.45578598976135254\n",
" 15 5.773644e+00 6.088645e-01\n",
" * time: 0.5431559085845947\n",
" 16 5.752314e+00 5.112412e-01\n",
" * time: 0.5610578060150146\n",
" 17 5.737885e+00 1.144800e+00\n",
" * time: 0.5744550228118896\n",
" 18 5.705868e+00 7.519057e-01\n",
" * time: 0.588249921798706\n",
" 19 5.674153e+00 6.088696e-01\n",
" * time: 0.6020488739013672\n",
" 20 5.638859e+00 6.919618e-01\n",
" * time: 0.615670919418335\n",
" 21 5.618025e+00 3.521743e-01\n",
" * time: 0.6341109275817871\n",
" 22 5.611260e+00 5.395682e-01\n",
" * time: 0.6476008892059326\n",
" 23 5.598189e+00 4.058224e-01\n",
" * time: 0.665855884552002\n",
" 24 5.587904e+00 4.003910e-01\n",
" * time: 0.6794819831848145\n",
" 25 5.581449e+00 2.705804e-01\n",
" * time: 0.6933939456939697\n",
" 26 5.579960e+00 3.429019e-01\n",
" * time: 0.7083039283752441\n",
" 27 5.572904e+00 1.868303e-01\n",
" * time: 0.7227308750152588\n",
" 28 5.571181e+00 1.895262e-01\n",
" * time: 0.741567850112915\n",
" 29 5.567069e+00 2.424169e-01\n",
" * time: 0.7553949356079102\n",
" 30 5.566763e+00 2.805680e-01\n",
" * time: 0.7692439556121826\n",
" 31 5.566114e+00 1.559386e-01\n",
" * time: 0.7837178707122803\n",
" 32 5.564883e+00 1.838331e-01\n",
" * time: 0.798846960067749\n",
" 33 5.563777e+00 9.072058e-02\n",
" * time: 0.8781428337097168\n",
" 34 5.563245e+00 1.158905e-01\n",
" * time: 0.8939008712768555\n",
" 35 5.562361e+00 9.000053e-02\n",
" * time: 0.9076528549194336\n",
" 36 5.561655e+00 6.989616e-02\n",
" * time: 0.9212779998779297\n",
" 37 5.561441e+00 9.758331e-02\n",
" * time: 0.9347019195556641\n",
" 38 5.561137e+00 6.924282e-02\n",
" * time: 0.949638843536377\n",
" 39 5.561125e+00 8.751227e-02\n",
" * time: 0.9638168811798096\n",
" 40 5.560963e+00 3.948629e-02\n",
" * time: 0.9818358421325684\n",
" 41 5.560779e+00 5.584112e-02\n",
" * time: 0.9956657886505127\n",
" 42 5.560700e+00 6.867982e-02\n",
" * time: 1.008971929550171\n",
" 43 5.560636e+00 3.988749e-02\n",
" * time: 1.0221788883209229\n",
" 44 5.560602e+00 3.007347e-02\n",
" * time: 1.0353589057922363\n",
" 45 5.560597e+00 3.458839e-02\n",
" * time: 1.0494978427886963\n",
" 46 5.560570e+00 1.872054e-02\n",
" * time: 1.0687198638916016\n",
" 47 5.560548e+00 1.725754e-02\n",
" * time: 1.0878658294677734\n",
" 48 5.560518e+00 1.617112e-02\n",
" * time: 1.1019480228424072\n",
" 49 5.560514e+00 2.324838e-02\n",
" * time: 1.1154999732971191\n",
" 50 5.560497e+00 1.218610e-02\n",
" * time: 1.1290688514709473\n",
" 51 5.560486e+00 1.302545e-02\n",
" * time: 1.166562795639038\n",
" 52 5.560483e+00 1.638541e-02\n",
" * time: 1.18241286277771\n",
" 53 5.560479e+00 1.171836e-02\n",
" * time: 1.196631908416748\n",
" 54 5.560473e+00 8.094007e-03\n",
" * time: 1.215466022491455\n",
" 55 5.560470e+00 7.455271e-03\n",
" * time: 1.2339179515838623\n",
" 56 5.560467e+00 5.691004e-03\n",
" * time: 1.2475829124450684\n",
" 57 5.560466e+00 7.749086e-03\n",
" * time: 1.2617719173431396\n",
" 58 5.560466e+00 6.431027e-03\n",
" * time: 1.2750539779663086\n",
" 59 5.560465e+00 4.257302e-03\n",
" * time: 1.2938299179077148\n",
" 60 5.560465e+00 4.525134e-03\n",
" * time: 1.3072218894958496\n",
" 61 5.560464e+00 4.140827e-03\n",
" * time: 1.3203649520874023\n",
" 62 5.560464e+00 2.903045e-03\n",
" * time: 1.3344128131866455\n",
" 63 5.560464e+00 1.769375e-03\n",
" * time: 1.347625970840454\n",
" 64 5.560464e+00 1.896148e-03\n",
" * time: 1.365858793258667\n",
" 65 5.560463e+00 3.897075e-03\n",
" * time: 1.3795068264007568\n",
" 66 5.560463e+00 1.941512e-03\n",
" * time: 1.3986079692840576\n",
" 67 5.560463e+00 2.877500e-03\n",
" * time: 1.4125878810882568\n",
" 68 5.560463e+00 1.746326e-03\n",
" * time: 1.427060842514038\n",
" 69 5.560463e+00 1.291818e-03\n",
" * time: 1.4626379013061523\n",
" 70 5.560463e+00 1.029335e-03\n",
" * time: 1.4769179821014404\n",
" 71 5.560463e+00 9.204316e-04\n",
" * time: 1.4961438179016113\n",
" 72 5.560463e+00 6.169948e-04\n",
" * time: 1.516507863998413\n",
" 73 5.560463e+00 7.875155e-04\n",
" * time: 1.531303882598877\n",
" 74 5.560463e+00 1.278484e-03\n",
" * time: 1.545435905456543\n",
" 75 5.560463e+00 6.520756e-04\n",
" * time: 1.5588958263397217\n",
" 76 5.560463e+00 8.250427e-04\n",
" * time: 1.5722320079803467\n",
" 77 5.560463e+00 6.891213e-04\n",
" * time: 1.5858688354492188\n",
" 78 5.560463e+00 5.768669e-04\n",
" * time: 1.5992839336395264\n",
" 79 5.560463e+00 3.228600e-04\n",
" * time: 1.6176369190216064\n",
" 80 5.560463e+00 3.575937e-04\n",
" * time: 1.631910800933838\n",
" 81 5.560463e+00 5.129062e-04\n",
" * time: 1.6452419757843018\n",
" 82 5.560463e+00 4.093644e-04\n",
" * time: 1.6590509414672852\n",
" 83 5.560463e+00 9.684812e-04\n",
" * time: 1.6733458042144775\n",
" 84 5.560463e+00 6.799935e-04\n",
" * time: 1.6877129077911377\n",
" 85 5.560463e+00 4.643456e-04\n",
" * time: 1.7011868953704834\n",
" 86 5.560463e+00 2.734082e-04\n",
" * time: 1.7144629955291748\n",
" 87 5.560463e+00 6.052775e-04\n",
" * time: 1.7511019706726074\n",
" 88 5.560463e+00 5.340380e-04\n",
" * time: 1.765038013458252\n",
" 89 5.560463e+00 3.653612e-04\n",
" * time: 1.7833077907562256\n",
" 90 5.560463e+00 3.061180e-04\n",
" * time: 1.7968418598175049\n",
" 91 5.560463e+00 3.155435e-04\n",
" * time: 1.8107388019561768\n",
" 92 5.560463e+00 2.386411e-04\n",
" * time: 1.829455852508545\n",
" 93 5.560463e+00 1.996828e-04\n",
" * time: 1.8478968143463135\n",
" 94 5.560463e+00 4.678266e-04\n",
" * time: 1.8678460121154785\n",
" 95 5.560463e+00 4.255047e-04\n",
" * time: 1.8883638381958008\n",
" 96 5.560463e+00 3.006264e-04\n",
" * time: 1.9024338722229004\n",
" 97 5.560463e+00 2.568562e-04\n",
" * time: 1.9153308868408203\n",
" 98 5.560463e+00 2.517784e-04\n",
" * time: 1.9367499351501465\n",
" 99 5.560463e+00 1.650642e-04\n",
" * time: 1.9554708003997803\n",
" 100 5.560463e+00 1.588244e-04\n",
" * time: 1.9688389301300049\n",
" 101 5.560463e+00 7.881208e-05\n",
" * time: 1.986454963684082\n",
" 102 5.560463e+00 6.136900e-05\n",
" * time: 2.0039868354797363\n",
" 103 5.560463e+00 5.690330e-05\n",
" * time: 2.0254340171813965\n",
" 104 5.560463e+00 5.690336e-05\n",
" * time: 2.0908617973327637\n",
"e(1,1) / (2π)= 1.7391793984948878\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXwURd748eo5ciIJCQmRgEQIcrkcCwi7oCKQKJdGQDSru4DoI7jHw+4+PKyuootnVp/1QlZXERQVWVTwxaKuJyAol3gAEYzchIRD5CbJTGZ+f7Q7v+nqUJMeOnSG+bxf+WO6u6a70pmkUvXtb5UWDAYFAADxyuV0BQAAcBINIQAgrtEQAgDiGg0hACCu0RACAOIaDSEAIK7REAIA4hoNIQAgrtEQAgDiGg0hACCuec7OZQKBwF33Tp/65z+fncvFOr/f7/GcpR/NOaO2ttbtdjtdi1hSW1vrcrk0TXO6IrEkzn830xK9Ects2XLo1Cm/pdPu2PHdvn0bbrvttmjrdabO0k+0urq60qd9uOPg2bkcAMB2IzucH7HM9dcv/uqrAxZPvLGo6IiDDSFDowCAuBa/fXwAgO00TcTccDs9QgBAXKNHCACwjaZpVp/ACgYd7kLSEAIAbMXQKAAAMYQeIQDAPjH4sAwNIQDANlHECIXTMUKGRgEAcY2GEAAQ1xgaBQDYJoqEeqdHRmkIAQD2iSJG6PjDNQyNAgDiGj1CAIB9tNhLqKchBADYRrM+1On40CgNIQDANtHkETrdEhIjBADENXqEAAD7RBEjjFT+0KFDL7/88uHDh0eMGNGjRw9zgUAgsGDBgtLS0q5du44cOVLvkr711lv79u0LlcnOzi4qKqrz/PQIAQB20qx/KRw7dqx3796rV68WQgwcOPD99983l5k4ceJf//rX1NTU6dOn//73v9d3btmy5fP/ePDBB998883TXYIeIQCg8XrllVfOP//8l19+WdO0zMzMBx98sKCgILzA7t27X3rppR07duTk5Nxwww0dOnS48847s7Oz//d//1cvUFVV1bJly5tvvvl0l6BHCACwjab9+LyMBco+4UcffTRkyBB9tHPIkCHLly/3+XzhBZYuXdq1a9ecnBwhxAUXXJCfn79ixYrwAm+++WZ6evrll19+ukvQEAIA7BPFwKhybHTv3r3Z2dn66xYtWgQCgcrKyvAClZWVLVq0CG3m5OTs3bs3vMCsWbMmTJigeJaVoVEAgG30DqG1t2ja2rVrx4wZE77zqquu0gcz3W53IBDQd+ov3G53eEmXyxUqIISora0NL7Bjx47ly5fPmTNHUQEaQgCAw1q2bHndddeF7+nYsWPoUEVFhf567969Ho8n1EE0FxBCVFRUtGzZMrT5wgsvFBYWtm7dWnF1GkIAgG2iWH1CCJGbmys1hCFDhgx5/PHH7777brfbvWjRosLCQo/HI4T45ptvMjIyWrRoMWjQoAkTJmzbtq1t27bffPPN7t27Q+HAQCDw4osv/u1vf1NfnRghAMA+dscIx4wZ43K5rrzyyttuu62kpGTatGn6/nHjxs2bN08IkZ2d/fvf/76goGDy5MlDhgy544470tPT9TLvvffeiRMnhg8frq4yPUIAQOOVlJS0YsWKJUuWHD169K677goNcj799NOhZ2QeeOCBoUOHbtq0qbi4uE+fPqH35ufnf/jhh4mJiepL0BACAGyjiQjpEHW9JYKkpKRRo0ZJO3v16hW+2a9fv379+kll8vPz61MBGkIAgH2sxwidnnObhhAAYKMYXI+Qh2UAAHGNHiEAwDZRJNQ7PjZKQwgAsE0UeYROt4MMjQIA4hs9QgCAfWLwYRkaQgCAbX5cWMniWxqoMvVEQwgAsE0MdgiJEQIA4hs9QgCAfaLoEjrdhaQhBADYJ4oYYQPVpN4YGgUAxDV6hAAA22jWE+SdfmiUhhAAYB9Ns54O4XRLSEMIALBPDD4sQ4wQABDX6BECAGwTgx1CGkIAgH2imWLN6aaQoVEAQFyjRwgAsE8Mjo3SEAIAbBNFHqHjDSFDowCAuEaPEABgnxica5SGEABgH2KEAIB4pmmxN9coMUIAQFyjRwgAsI0mNKsJ8jGTUH/ixIm77rqrf//+AwYMeOSRR2pra/X95eXlv/rVr3r16jVhwoT9+/c3WD0BALFAs/7ltPo2hH/605+WL1/+9NNPl5SUvPDCC0888YS+/9prr23atOmcOXOEEMXFxQ1USwAAGkh9h0bXrVt36623duvWTQgxZsyYdevWCSFWr15dVla2cuVKr9f75JNPZmVllZaWdu7cuQHrCwBoxKJ4WMbxTmF9e4RFRUXz5s3bsWPHN99889ZbbxUVFQkhvv766x49eni9XiFEampqly5dvvrqqwasLACgcdN+nHbbGmfrXN+G8LbbbquqqurZs2efPn3y8/OvueYaIcT+/fvT09NDZZo1a3a6MGEwGAyFFQEAscjn80UudA7HCEePHt23b98DBw4cPHjQ4/H85je/EUI0bdr05MmToTLHjx9PS0ur8+2aprnd7jOvLgDAKfr437mnXg1hMBhctmzZL37xC5fLlZCQcMMNNyxdulQIkZeXV1ZWppcJBALbtm3Ly8trsKoCABq9aMZGHa5yvRpCTdM6dOjw9ttvCyECgcDbb7/dqVMnIURhYeGxY8eWLFkihFiwYEFiYuKll17aoNUFADRmMTgyWu+nRp9//vlf/vKXc+bMqa6uzs7Ofu2114QQiYmJs2fPHj9+fHp6+vHjx19++WXGPwEAsaW+DWHfvn3LysoOHDjg9XrDH5AZNmxYeXl5ZWXl+eef7/EwTw0AxLcGmHR73759s2fPPnbs2IgRI/r27Wsu4Pf7586du2XLli5dutx4440ulyu0f/78+Rs2bGjWrNmoUaPy8/PrPL+1uUazsrLCW0Gd1+tt3bo1rSAAIKrkCVVLePTo0UsuuaSsrCwtLW3YsGHvvPOOucyECROeffbZ3NzcJ554Qn+WUwhRU1NTUFAwc+bM9PT0H3744dNPPz3dJWi9AAC2iWKFenX5l156KS8vb9asWUKItLS0Bx98cMiQIeEFduzY8c9//nPXrl1ZWVkjR47Mz8+fNm1aTk7OzJkzjx8//tlnn0Xsp7H6BACg8Vq2bNmVV16pvy4sLFy5cmVNTU14gU8++aRr165ZWVlCiNzc3IsuumjlypVCiCVLlowdO/add9556qmnvvjiC8UlaAgBAPaJ4rFRZY+woqIiOztbf92iRYtgMFhZWXm6AnqZvXv3CiG2b98+Y8aMN998c9++fVdeeaU+J3adGBoFANgmmsRATVu7du2YMWPC9xUUFNx6661CCI/HE5qYzO/3C1Nev8fjCQQCoU2fz6cX0DStW7dus2fPFkJ069btj3/847hx4+q8Pg0hAMBhLVu2vO6668L3dOnSJXRI7+EJIcrLyz0eT3j/Ty9QXl4e2ty7d2/Lli2FELm5uaGTXHzxxeXl5X6/v854IQ0hAMA2Uaw+oQmRm5srNYQhw4cPLykpufvuuz0ezxtvvDF06FA9Yf3LL7/MysrKzc0dPHjwzTffXFZW1r59+w0bNlRUVAwYMEAIcc0113z00Uf6SdasWZOfn3+6p2ZoCAEAjdeoUaNmzJhxxRVXtGvX7l//+te///1vff9tt91WXFw8efLk5s2b/+lPfxo8ePCQIUOWLFkybdq0pk2bCiEmTJgwZ86coUOHtm7deuHChS+++OLpLkFDCACwj/UYobp4YmLi0qVLP/jggx9++OHhhx/OycnR9z/33HP6k6JCiGnTpg0dOrS0tPS2227r0aOHvrNp06Zr1qx57733Tp06dc899+jjpXWiIQQANGoJCQlDhw6Vdnbt2jV8s1evXr169ZLKJCYmjhgxIuL5aQgBALaJxRXqaQgBALZyumGzioR6AEBco0cIALBNFAn1jq/MS0MIALBNNHmETg+l0hACAGyjiQjLKtX9JkcRIwQAxDV6hAAA+0SxQr3TaAgBAPYhRgg0dkHFVqMm/61w+m8HcM6gIQQA2Eazng5BjxAAcA4hRgg4wNpop/L4Gb3Zmsh/K4wl5EsHVSdgHBVOicU8QtInAABxjR4hAMA2UUyx5niXkB4hACCu0SNEDAhGCM2pgoSW3hr5eKSq1F/Q9F+wvB20cFg6KAcQlQFF4fx/5ICTaAgBALaJ5mGZhqlJ/dEQAgBso4nYW4aJGCEAIK7RI0QjEDmqZyUKKB1VHpaPmq4dtBJyjBDUk2jyueTFa6QtdUhReTFNqpqpcNAYRCQNEdGLIqHe6Q8YDSEAwD7WY4SOoyEEANgmijxCYoQAADiJHiEcECEZTx3WM20Hg6qwn7WjAWXFTOWtzT0qx/zM/wUbUwOl8i5jGE+TCktHpUtHildKoVOpvBRhjLWBL0CNhhAAYJuohkYbqC71xdAoACCu0SNEw4gwgmgxh0E9vBkwbAeUo52BgHRUtSmEkHcoa6IWcbjSpR4LNW66lCOl8lHjf7x1/AMu71IN0pJrAYUoZpZx/DNDQwgAsJXTDZtVNIQAANtomuV0CNInAABwEj1C2CTCtGdBxVF1zM+8J6DerLVyVFnYvEeOKZ5BjFCK+QlTYE/edCs31YWVZxbmIKKpgOGoMtdCnjku1kbJcKaiWH3C6Q8JDSEAwD5RzDXqNIZGAQBxjR4hAMA2mtDkdVTq8yZH0RAiWpaCglamPZPidsIUqKutNeQG1vqNR/1WjhpPZb50ba0yRhhQxT7lvEFlXqAQwi1F8jyGTbfbMH7j9kibmoWjbnkoyBVUbcojR1JQ0FLIUDj+Rw8NK5oV6p3+SNAQAgBsFWMdQhpCAEDjtmfPnueff/7w4cNFRUUDBgwwF6ipqZk1a1ZpaWnXrl3Hjx/v8XiEEKtXr166dGmozMSJE9PS0uo8Pw/LAABso0VBecLDhw/36dPn4MGD+fn511133VtvvWUuM3bs2Ndee61r166zZ8+eOHGivnPp0qXz5s374T8CgdMuLkOPEPWjni1UWAsKWsrtE6bAnt8Y9vP7jEd9tYbNGsNRn1xYCijKvyrqIKI0kWlQmqJTSEFBQ2Ep1U9EigJ6vIZNr3HTk6Aq7PG6DfX0yHdYulZQSkMUUhqi9G4LIUNBouG5zva5Rl988cWOHTvOmDFDCJGcnFxSUnLNNdeEF9i6deuiRYvKy8szMjJGjBiRl5c3ffr0li1bCiEuueSShx9+OOL16RECABqvTz75ZPDgwfrrQYMGrVq1qqamJrzAypUru3fvnpGRIYTIycnp2LHjp59+qh/avHnz/fffP3v27MOHDysuQUMIALCPFtXX6VVUVGRlZemvs7Ozg8FgRUXF6QroZfQCLVq06Ny5c21t7WuvvdapU6edO3ee7hIMjaJe6phKzL6xUCnDwW8an5QGMH3VhsFPX41hs0Y6Wi0dNQ6NGt8rDZyKSMkYpvQJ9dColOFgGhr1qAc/DcObCYnGwonSUcOmN8F4/41HhRDSWKknaOFfZEsjpcKcX1H/KyEmRLEwr9BWrlxZUFAQvvPqq6/+7W9/K4Twer1+v1/f6fP5hBCJiYnhJRMSEkIF9DIJCQlCiHHjxo0bN07fOXr06JKSkpkzZ9ZZARpCAIBtophhTROiffv2U6dODd/Ztm1b/UVubm55ebn+ury83Ov1hvf/hBAtW7bcs2dPaHPPnj25ubnSJfr06fPxxx+frgIMjQIAHJadnT3YKNQQXnPNNQsXLtTjgvPnzx8+fLjb7RZCrFmzZteuXUKIwsLCrVu3lpaWCiHWr1+/f/9+PcXixIkT+hl8Pt+SJUsuvvji012dHiEAwDb1SIio4z2Kg9dee+3f//73/v375+XlLV269P3339f3//a3vy0uLp48eXKzZs2mTZtWUFAwePDg99577/7772/SpIkQ4vLLL2/WrFlmZua6deuaNWt2xx13nO4SNISom2mBoQjbZxQU9EkxPzlQJ4X9aqoMm9VVfsPmKWlTeq/hqBxQNF1anV8hL9ukTAxwyTFCeTxGTpAwZkRIYb+EJMMvb2Ky4ahfPhphtalgwBg1lH/8ckxRQQ4ZmtMnpPyKoDFkSMzwHGBr+oTX6/3ggw+WLVt29OjRZ555Rn86VAjx0ksvhV5PnTp1+PDhpaWld9xxR8eOHfWd//73vz///PNjx4794Q9/6N27t6J5piEEADRqbrd74MCB0s4OHTqEb3bp0qVLly7hezIzMwsLC+tzfhpCAIBtopl0u2FqUn80hAAA22iasJw+4XRLSEOI/6gjVTDsoHmKNSk6ZgyuWQoKSrl9UhhPmMJ+VdLmCcPmqZPGoyd9xlOp4otS0qEwRQ0bNEYoBQWl1MBEOexnqHlStTd805+iXm1K/sWXEyKD6r8MqpChvBZdpL9x8h9BJmCLddHlTziK9AkAQFyjRwgAsE0MdghpCAEANrI7j/AsoCHEj9TLKtWxDFPAQuKgNH2oFHiTgoJSCFAIceqEIc536rihwMnjNYqj0nurTqqSDmvqiBEa9kjBzoAx9ibdM+m32yWvsiT/8nu9UqagMUaYbPhtTUoxbErRTb/PEDKsjbTQlTJALDTNI20bt1zGTcPJXKY/cZq8S7Vsk+N9BcQDCzHCL7/8sqioqGPHjldcccXKlSv1nTt27Bg1alSnTp1uuOGGvXv3NkwlAQCxw76lJ86O+jaEmzdvLigoGDBgwKJFi6ZNm5aSkqLvLyoqateu3dtvv928efPrr7++weoJAIgBeh6h1S9n1XdodPr06TfeeOPkyZOFEKEJbD799NPdu3c/+OCDHo/nkUcead68+YYNG37yk580VGVhK/UkauoZ1ESkoVF5TXnjRGXSCKQ0PikNZgohTh4z7Dlx1DAWesJ49OQxw9GTx41DoydUw7DS5G3ClOlhGhqVsg6MyzAZf7+lJenNQ6Me5dBoknFotLrKODRaY0yfsFJPM6lm0mCmtCmNdMqFTX/jgvLZlQvcMwFbrIlirlHLMUW71bdHuHbt2g4dOowfP37IkCFPPPFEbW2tEGLjxo3du3f3eDxCiOTk5C5dumzYsKEBKwsAgN3q2yPcvXv3X//61xkzZqSlpU2cOPGHH3649957Dxw4kJaWFirTrFmz/fv31/n2YDCot50AgBjl8/m8Xm+EQo0j7GdJfXuEaWlpt9xyy7Bhw/r37/+Xv/zllVdeEUKkp6eHFnwSQhw7dqxZs2Z1vl3TNH0FKQBAjIrcCkYVI3S84axvj7Bdu3ahzl96evrx48eFEBdeeOGWLVv0nYFA4LvvvgstpYjYI0cFjVumoJIcFDQmEkgxKp9yHSUpUCflPwhTUPD4EdXmCSlGeExKn1DNuGaOEfqUMULTzGSqGKEmT7FmLX3CtICUMSgoLxelqqeZOpwpJX5Ic8VJm5rbUBMpuUKYg4gRUjeURwE71LdHOHbs2FdfffXkyZO1tbUvvPCCviLG4MGDq6ur33jjDSHE3LlzmzRp0q9fvwasLACgcYvBh0br3SOcMGHC2rVr27Rp4/V6u3XrNmfOHCFEQkLCyy+//Ktf/ep3v/ud1+t99dVXXfKinACAeBKDc6zVtyH0eDzPP//8U089pWlaUlJSaP+gQYN279595MiR9PR0xx+BBQDAKmtTrCUnJ5t3ulyu0z0jg8blDBIHzUEmU+KgaqElaaIyeVmlE6op04QpU9AUI6w2FD6qyiM8dUJKHIwwxZoUe5MS8gJyjNDwXnmKNZcqDieE8HgNAyoJ1cYYYbXht9VUMeNkb8qgoGmSM1OOo7ypmhxO2pTjiy65Ji7ppllJK5RL8793Y2Q9j9DpHyRzjQIAbBNN0M/pf2hoCAEA9onBGCHPtgAA4ho9wjgSYekdZeJgwBwjrJVihFYmFzVG5k6dlPIIzXONqjIF5alHjxqnHj2uyiOUEgel9aGEaQEpUyjOUFi6R6YZOw2bUnKeEMLjMezx16huqTx9qFxx6dKRpjk17pFmPfUkuFSbXimCaKiK2xQKDRijhpbSCgkRNn7RzDXq9E+ShhAAYJsYHBllaBQAEN/oEcYvU76EtKlKDBAR0yeMw3o+ed0l4xRrJ1UZDsI0vCnNmiYflfMlpDnVVPkSftPQqDR1nGk9I0PhoPGmSgM+pgXrTdcyDpbW1ho21akaEvUwrDSYKYTwJhjGQr2JhpuWkGic+y3ReNOM7/V4jRkmXtOQfIRPmmr00/GuAyKLwS4hDSEAwDZaFOsL0hACAM4ZUeQROt0OEiMEAMQ3eoTntAgJE6qjZ5g+4ZPTJ4ybxvQJU8hQjhFKc7CpY4qmhZb8ik0pX0L6LoQQtXKWgjpQp0qgkP5N1mpNeQVu5fxt6lnT5ACklA5hnLwtUV4cNDHJcFsSkw1/GWqSpcCq4WhCkuGmeeWEE/m/ben7cslzrKmTUJhxrdEjRggAiGtR5BE6vWADQ6MAgLhGjxAAYBsm3Ubjog4JqudUk4NhEWOEtar5wPzGZZjUM65JYTwhRJUc5zOeTVpKqUp1LZ9y3jJzjFCdvRdUZ/NJaYXSAkOmWxoMGApEyhQ0fF/SGk8er3THDEHBqiT5DidJd1hencq4ApTxpynfQ2XmpYj00QqqZ1wz/sVUhxOBeqIhBADYRovB9QiJEQIA4ho9QgCAbaJJqHd6UJuGMJ5YiGdFCBkKU/AsYJxrVJ1W6FPORCoF+cx7rEUBlZmCtXLSm+nbVM+5aiVTU855M//ymy6uKO+vMWz73FLipvKORbrDPuUt9SkDq9InwZyBGiEX08qn1OkRNdQhmmWYnG4JaQgBAI3a9u3bn3nmmSNHjhQVFV111VXmAlVVVU8//fSWLVu6dOkyadKkhISE8KPz588/ePDgr3/969OdnxghAMA+mvUvpUOHDv3sZz8LBAK9e/ceO3bsG2+8YS5z0003vfvuu5dffvnChQtvueWW8EOfffbZb37zmwcffFBxCXqE8cs0IqUakzJP7yUvwySvyqQaPfOrR0p98sCd37hHKiCnavhVT/NLU6bJA7zmOdPkfAoLayFJNPm9EZZuN46MyjV31aq+a/kOy3dMvsNSAUtJJvIdDkQYGlVOFSd/DoPSNHWqt6JxsDuPcM6cOd26dXvkkUeEEG63u6SkZNSoUeEFvv322yVLllRWVqalpRUWFrZu3fqBBx5o3bq1EKK6uvrXv/71tGnTHn74YcUl6BECAGyj/fi4jKUvlZUrVw4cOFB/PXDgwHXr1lVXV4cX+Oyzz3r06JGWliaEyMrK6tSp06pVq/RD99577+jRozt16qSuMw0hAKDxqqioaN68uf46KysrGAxWVlaGF6isrAwVEEJkZ2dXVFQIIb788sv33ntvypQpES/B0CgAwD5RrT6xcuXKgoKC8H3Dhg2bPHmyECIxMbGmpkbfqb9ISkoKL5mYmOjz/f81Z6qrq5OSkvx+/y233DJjxgyv1xvx+jSEccQU/goqDsuTh9URPFNNlCVHEJXJFaZN+Vqmt0uzuxkvrZziS/3sfl0xv+iDgvKJjO/V5NnD5KihOoMlYOW7Vt/AugpY+AFFuMPmkKCckaJcWUkqrBEybOw063mBmhDt27efOnVq+M727dvrL1q1arVnzx799e7duxMSErKyssJL5ubm7t69O7S5Z8+e3Nzcb7/99uuvv77pppuEEKdOndq/f3+7du2WLl2qxw4lNIQAAIdlZ2cPHjy4zkNFRUV33XXX3XffnZSUNG/evKuvvtrlcgkhli9f3qpVq7Zt2xYUFIwfP/7rr7/u2rXr6tWrv//++yuuuMLj8WzevFk/w8qVK//nf/7n/fffz8nJqfMSNIQAANtomvUEeWX5a6655rnnnrvkkkvy8vLWrl374Ycf6vunTJlSXFw8efLk9PT0Bx54oLCw8LLLLlu2bFlJSUlKSooQom3btnrJbdu2eTye0KYZDSEAwD52r1Dv8XjefvvtNWvWHDp06OWXX27atKm+f/78+fqTokKI//7v/x42bNjmzZtLSkouvPBC6Qz9+vX79NNPVZewWF/EjUjBs0jLNqmOytl7cvTLtBaSlSiUfGn11HERQ6FnEBRUM5/ZlGgolVdHEA2b8upRyviiMIcYlT8gSz/riJ+cSKuFIcZEM9dopAIul6tv377Szry8vPDN/Pz8/Pz8Ot+enJzcpk0b1fkjVQAAgHMZPUIAgG1icdJteoQAgLhGjxA/OsNIjXrmUnkiU2WUqK6okirzTJkPaZ5DNXZDUqpbKn+TylCoHBmNGLeT77CFmlgVuz8exC4aQgCAneyddPssoCEEANgmmhhhA1Wl3mgI8aMz/Cxq8qam2DSXNmyZqiL9Xpk2LdTEdLxRDcWpfwiqWyp/k/ItVd1Ac/kIPyArNbHK8b+JOFN25xGeBTwsAwCIa/QIAQC2icX0CRpCAIBtophZxnE0hDiNyHE746bLGChyqY661JtuecTe5ZZOrnq7fGnlNyL9Kxo0f5v2LcMkn9n8xyJCYE/1dum7Nt0T6Q7L15buufoHZOlnHfGT43h8CCBGCACIa/QIAQC2icUYIT1CAEBco0cYR+QEOjnKZAh/yf+jmTPP1FEo46bbI226lJvytUxvV0W85HCX27BKkFZrrHZA+q6lK8uLH2la9CFD8y1UF1BHYSN91+r7H/EOW/gBqSO40qYQ8jdm+qQpC8vnQmNk67q8ZwMNIQDAPjGYUE9DCACwTRQxQscRIwQAxDV6hPEr0pychvCXOdAjRwHloKAqyOTxGja9CcZNr1u6lse4Ryognc1jvFatMXgWcBuDfAHjN2aK+QVcUlBQCllFHyQ0/9PsMu6S77AU9nOrvmv5Dst3TL7D3gTjHU5QnU3+4bqlTdV3Ier6LIWzd+ZSnH1RJNQ73oGkIQQA2CcGY4QMjQIA4ho9wniiXoDIyrP7wjwLl/L5e696LDTRMC6XkCQP3El75M1qw6a/xpAvUVtrHBoNqPMfAkK5w7Tyu1xcQT1lmqgj50Q5vJyguqUJyltaxx02lpd+Iqaxa1XF5LwO0ydHPUOblaWo0BhpdcRZIr6FSbcBAOcKLYo8woapSf3REAIA7BNdjNDRFbKJEQIA4ho9wnOZOiZoKShYR4xQ/TS/lO4ePVcAAB2oSURBVNJgfDpfikglJhk+h4nJ8scyybinpspv2Kw2HPX7DGE9OSgYUP3naR7Sqa1VvT0or9ukirtGDLvKCRKeM4kCSrfUcDTJdIeley79RKRrST9NOblCOdmbqGOdJmXoNEKuBRqfKBLqNWuJSLajIQQA2CaaPMKGqUn9MTQKAIhr9AgBALaJaj3CBqpLfdEQntNMCy9J24rS0npDdUyUpQ5oeaWAlnFTjmAZYn5JKfLHsrrKGCOs9oZvykHBWgsrJWmuWsOpakzfZq2FkweNt1SeLcx4bnPwTI6zKoOCUlQvKcVwT5JTVUeTUk1R2BRVTDFBDhlGn1YoTJ+liKtTqY46/QcU5wYaQgCAbaKIETr+Dw0xQgBAXKNHCACwjSasxwid7hLSEMYvdWjGlN0lh9pME2MaNuWZMOX4ljGtzRjz81XLE376aoyRPDkoaNgMmqYLDScH6ozfhc8tv9nvt3At5eyt8qSaLrc8HuNRJw4aMwWlsF9KE+XmeaqjIlJMUfp5eeW0QilGGGmu0QifNONRgVhzzq8+4ff7f/GLX0ycODG059tvv73qqqsuvPDCESNG7Ny50+7qAQDQsKw1hCUlJRs2bFi9erW+GQwGr7322j59+qxatapjx45jxoxpgBoCAGKG/rCM1S9nWRga3bx584IFC/7whz88+eST+p4VK1bs379/2rRpbrf7vvvuy8rK+vLLL7t3794wVcWZipBMIRVWDiGKSOkTHjl9wrhSkvHRf2nhJL/PMBAqhPD7E8I31dOeSUwDksYhXK80W5vp0srcDPWiTuo51dymvAL1TZNyGKTBTGm0M7Wp4Y6lnmfYTDFuCiGSjW9PlrMplDOuGde7P9P0CeWm038wUQ9RTLHm9A+2vg1hIBC49dZbn3766fDxz02bNnXv3t3tdgshkpKSOnfurO9pkJoCABq/BogRfvvtt08//fThw4evvfbaoqIic4ETJ0488cQTpaWlXbt2/d3vfpeUlCSEWLNmzRtvvFFeXp6RkTF69OjLLrvsdOev79Do3/72t+7du/fr1y985/fff3/eeeeFNtPT0w8cOFDn24PBoN/vr/MQACAmOPJn/MCBA/369WvatGlhYeHtt98+b948c5ni4uKVK1cWFRV98MEH48eP13du2LAhMzNz+PDhLVu2HD58+L/+9a/TXaJePcJdu3Y9/vjjS5Ys2bZt2/79+6urq7dt23bhhRc2a9bs+PHjoWJHjx7NyMio8wyapnk8PKEKADGsPn/Gbe8Qzpkzp3fv3vfdd58QIhAIPProo8XFxeEFNm/e/MEHH+zbt++8884bOHBgy5Ytd+7c2aZNmwkTJoTK7Nmz51//+tfw4cPrvES9Gqe9e/cmJiaOHDlSCHH8+PHDhw8XFBRs3Lixbdu2mzdvDgaDmqbV1taWlZW1a9euPieEM5RBQnkdFOmhdtPYQYT0CWPcyJtgOHlikmGz1piiIIUAhXliswhBQePj+8r1oapPGepZUyXHCH3GgGWtX70qk3GKNfkequ6YEMLrVa2sJOecpKgSJKSgYJM0KWQop0+kNDEUkOZgk1Zl8qqXYYqYPqFehkkqTZAw1mimn2nktyiLf/bZZwMGDNBfDxgwYOzYsVVVVfrgp27VqlU9evTQhyczMjK6dOmyevXqNm3ahArs3r17xYoVf/zjH093iXoNjfbt23frfzz22GOdO3feunVrcnLy4MGDA4HAK6+8IoR47rnnMjIyfv7zn9fnhAAA1EdlZWVmZqb+unnz5sFgsLKy8nQF9DIVFRX663nz5mmadsEFF/zkJz+56aabTncJy1OsNWnSJDc3V3/t8Xhee+21adOmpaenP/7446+++qr1h4UAAOcQLZqvpUuX9jIqKSnRz5eUlFRTU6O/rq6uFkKkpKSEXzA5Odnn84U2q6qqQgWKi4uDwWBZWdl33313zz33nK7KluN2V1999dVXXx3a7N+//7Zt26SOKgAgPkWXF9i9e/dHH300fE+ox9WqVatdu3bpr3ft2pWUlNS8eXOpZKiAEGL37t2tWrUKL5Cfnz9p0qQnnnhi+vTpdV7dngdYaAXPBVbSCoU5RmgMxQU9xnw7Y+aZlH4XqPUojgpT7E2umzI/T52cV5VkzCM0xQj9yhihaVUmVYzQlHlpziNUxQiTpHWXlHmEUqagFBSUsgyFEMlNjCeXEgeTlFOseVSfhIhTrKkTBxFzolqPUEtPT+/Zs2edR0eNGjVlypS77747JSXl5ZdfLioqcrlcQoj333+/TZs2F110UWFh4c033/z555/37NlzxYoVR44c0WOK+/fvz87OFkIEAoEPP/zwoosuOl0FeJITANB4DR8+/IUXXvjpT3/aqlWrzZs3f/jhh/r+u+66q7i4+KKLLmratGlJSclVV13Vp0+f1atX/9///V9ycrIQYujQocFgMDs7e8uWLeedd97ixYtPdwkaQgCAfezOn3C73YsWLfryyy+PHj3au3fvUPxv0aJFqamp+utJkyYNHz58y5YtnTt3btmypb5z1apVGzduPHjwYMuWLTt16qTop9IQAgBspFldVilieU3TevToIe08//zzwzdbt27dunXr8D0ej6eeM53REMYv6d+joJW0QmGew9NYXJ561Hg0GHAbN42RtrrrG14XVezNtJiRMRvPGBRMOhUhRqjOI5RXZVLeM/XsrCJyHqG0UpJq6tFkeVWmBONR+Rc/OVVad0k9uaiUOKgKCmp1xAgtJA7yHDrOAhpCAIBttBj894WGEABgnxhcmJeGED8yzb8m75DLSzOZSaWNA62eoEs6bNyK8DmUqmLK3FDO7pZouLQ06FddZRga9ZmWYfIZl4iSZoOTMz2US9SbJqUzDY0al3qXshSkec7UM65J6yhJuRZSdoQwj4XKc6oZB7qVCy1Jm+ahUVO+hHqkFDEmuvSJBqpMPVmeWQYAgHMJPUIAgG2imFnG8WEAGkIAgH2IESKGKWOCkf/FU6YKmLgVxzRN/liqV1aSUzUSVHOq1SQb8yWqDdeSIoJCCL/v7MUIpbQEr/IbSVCGDE0BRdWUaeaTy0FBY9hVXmhJCgpK6yaZQ4Tqj4bTfxMRh2gIAQC2ieJhGcf/+6EhBADYJgZHRnlqFAAQ3+gRom7qCdhEpDnYXNb+xTKGrOqYzk06uWo9I9O6S4bUQCko6K8xHPX55BihFBSUplgzTQ5nXIZJmi1MWW1hihp6vVLs0xgjTFQlHZpifqp1lITppkmZguqgoDynWuS8QCZRO9fF2s+UhhAAYJ8YTKinIQQA2CaaPEKne5DECAEAcY0eIerF/B9bUA4FqXLoLIUMNU0uLU1XaSlG6JdCa0mGmJ86TdC8R0ocDBqLR4oRGgq7TJNwuuXInJTMp44gqgpLiYCeSOFJ9fShEYKCyqlERezFj2BRFI+NOo2GEABgG02zHvNzemyUoVEAQFyjR4j6MU+UJU0udgYjpaaxtQhLPmluw4ikaeV3w1GP13A2r5QOYVxiPuCXL10rj4Wq0yeM1ZY2Iy3d7payRIwDmG7lxHLS+HCEo27T4LMyQUWL8PMS4vTbdfyrH2vjZrAkioV5ne4Q0hACAOyjCU2z+M+O0+0gDSEAwEYxOMcaMUIAQFyjR4hoSavtWAoZyhO4GbbM/52ZVvYxPuvvMrzfbQwZBowxwkCty7hpPBqQY4SmfAnjplxcRf4uTDFCdaDOUkpDhKOmS6vjl3LN5Tcrg4JO/7OPs816Qr3jaAgBALaJYhkmzfyA3NnF0CgAIK7RIwQA2CcGH5ahIYRNziRkKJ3JFMEKKqNrLmPcLuBShfEixPzMMUJpR9C2GKE5kOJSfpvqtaikk5kSASOsfCSPZclBWbn46cs6/0cNzopm0m1TDu5ZRkMIALBNFHmEjv/3RIwQABDX6BECAGwVYx1CGkI0EGXmWTCoDlmZ4gVSiFEOQBrLGuNhQWMcz2UsbTUvUDqbtciGHHiL8Nuvzjs05VZaOFrHzLHqa1s4iHjHwrwAAMQYeoQAANtEk1DfQFWpNxpCOED+NVHmWohI6RYRBk6DqosFpZQF03WCpsopKqJehklWx/jkGeQwyLkZyktFWhuJjAhEjzxCAEBci2KuUacbQmKEAIC4Ro8QAGAbrR5PRJve4nCXkIYQjUCEZ/tNcT55Sx1BVG1HnPQ+KL3f0qRqauqgXx3lrRy0em6nx6YAhY0bNz755JNHjhwpKioqLi42Fzh69Ohf//rXLVu2dOnSZcqUKampqUKI1atXL1y4cPv27dnZ2bfccku3bt1Od36GRgEAjde+ffsuu+yyvLy8MWPGTJ069aWXXjKXuf7660tLS8ePH79u3bqxY8fqO++7776kpKQxY8ZkZGT87Gc/++qrr053CXqEAADbRLceoeLoCy+80L9//zvvvFMIUVNTU1JS8qtf/Sq8wKZNm5YvX75///7U1NSf//znOTk527Zta9u27eLFi/Uzjxo16ssvv1y4cOHpOoX0CAEAttFnlrH0pR65X7NmzaWXXqq/vvTSS7/++uuqqiqpQI8ePfTh0PT09IsvvnjdunUirH0NBoN79uw5//zzT3cJeoSIARH+v5RmTYsQcQyqDpsnd5MvZVswLfKJziBmSCIgHBNFHqFSZWVlZmam/rp58+bBYLCysjIvLy9UYN++faECQojMzMyKiorwM8yYMePYsWNSPzIcDSEAwGFLly7t1atX+J6RI0fqw6EpKSmhLuCpU6eEEHrnLyQ5Obm6ujq0WVVVFV5g/vz5Dz300Mcff5ycnHy6q9MQAgBsE1WMUHTv3v3RRx8N33nBBRfoL1q3br1z50799c6dO5OTk8P7f1IBvUyrVq301wsXLpw8efK///3vDh06KCpAjBAAYB/rMUJNE+np6T2NsrKy9PONHj369ddfP378uBBizpw5o0aNcrlcQoglS5aUlpYKIQoLCysqKlatWiWE+Pjjj0+cOHHFFVcIId59992JEycuXry4a9eu6irTI0Tsi5SGaCAFFJWninSps4qwH+LT0KFD586d27Vr19zc3D179nz44Yf6/unTpxcXF3fu3LlJkyaPPfbY8OHDe/To8cUXXzz11FOJiYlCiNtvv726uvr666/Xy48ZM+ahhx6q8xI0hACAxsvlcs2fP3/Lli2HDh3q2bNnQkKCvv+dd95JSkrSX48fP37YsGFlZWUdOnRo3ry5vvOLL74IBAKh8+itY51oCAEAttFEFMswRS5vDvJlZGSEb2ZnZ2dnZ4fvSUtLq2cFaAgRZyyNowKwKIoV6h3/PeRhGQBAXKNHCACwDwvzAgDime1zjZ4FDI0CAOIaPUIAgG2ieFjG6ZFRGkIAgL1iLUbI0CgAIK7RIwQA2Kk+CfKNCg0hAMA20cQInW43aQgBALaJJn3C6R4kMUIAQFyjRwgAsE8MzixT3x7hQw891KtXr+zs7J49e86fPz+0f+PGjZdffnlmZubgwYPLysoappIAgNgQxaq8MdMQfv/990888URpaendd989fvz4NWvWCCECgcCoUaOGDh26c+fOfv36hdY/BAAgVtS3IXz00Uf79evXvHnzoqKi3r17r169WgixfPnyw4cPT5kypUmTJn/+85+3bt26bt26hqwtAKBR0358XMYKp7uElh+WOXTo0IYNG376058KIb755puuXbu6XC4hREJCQqdOnTZv3mx/HQEAaDDWHpbx+/1jx44dMWJEv379hBCHDh0677zzQkfT0tIOHjxY5xuDwaDf7z+TigIAnOXz+bxer7pMLOYRWugRBgKBcePG+Xy+f/zjH/qejIyM48ePhwocOXIkMzOzzvdqmubx8IQqAMSwiK1gjKpv4xQMBidNmrRr16533nknMTFR35mfn79p06ZgMKhpmt/v//bbb9u3b99gVQUANHrWE+qdDhHWu0c4adKkTz755Omnn963b9+2bdsOHTokhBg4cKDH4/nHP/4RCASefPLJnJycPn36NGRtAQCNm2b9y2n1bQhXrFhRXV1dVFRUUFBQUFDwwgsvCCHcbveCBQtmzpzZpEmTV1999bXXXnN8oWEAACyp79Doxo0b69x/ySWXfPXVV/bVBwAQw2LxYRkeYAEA2CaaSbedbglpCAEAtnK6h2cVq08AAOIaPUIAgG0ax3Og1tAQAgBsE4sxQoZGAQBxjR4hAMA+MbgwLw0hAMA+muWhTseHRmkIAQC2iSahvmFqUn/ECAEAcY0eIQDAVsQIAQBxSxOa5njLZhFDowCAuEaPEABgG1afAADEtwaYY239+vWPPfbY4cOHr7322ptvvtlc4IcffnjggQc2bdrUrVu3O++8s2nTpkKIysrKjz76aP369ZmZmXfccYfi/AyNAgBso0VFccKKiopBgwb99Kc/vf322++///5Zs2aZy1x33XXl5eVTpkwpKyu76aab9J3vvvvuK6+8UlpaunDhQnWd6RECABqvWbNmDRgw4Pe//70Q4uGHH54+ffqECRPCC3z99derVq06cOBAcnJyr169WrRo8d133+Xn548bN27cuHFz58596qmn1JegRwgAsI0eI7T0pR5KXbduXb9+/fTX/fr127Rp06lTp6QCPXr0SE5OFkI0bdr04osv/vzzzy3VmYYQAGAfzfqX0r59+zIyMvTXmZmZQojKysrTFdDLSAUiYmgUAOCw999/v127duF7brrppr/85S9CiNTU1FAXUH9x3nnnhZdMTU2tqqoKbZ48ebJJkyaWrk5DCACwUTTrEf785z9/5plnwnc2b95cf3HBBRfs2LFDf719+/aUlBS9XxjSunXrUAEhxI4dO1q3bm2pAgyNAgBsE8XIqCZEampqWyM9BUIIMWbMmAULFhw5ckQIMWvWrOuuu05vaN98880NGzYIIa688sr9+/d/8sknQoj33nuvurp6wIABlupMjxAAYB+71yMsLCwcMGBA586dc3Jyjh079v777+v7S0pKiouLf/KTn6SkpMyYMaOoqKhTp06bN29+9tlnExIShBAff/zxqFGjampqqqurMzIyrrzyynnz5tV5CRpCAEDj5XK55syZs3PnziNHjnTp0sXtduv7P/roI6/Xq7++8cYbhw8f/t1331100UWhCGL//v23bt0aOk+osBkNIQDANhET5KPTpk0baU9qamr4ZlpaWs+ePcP3eL3eZs2a1efkNIQAANvE4lyjPCwDAIhr9AgBALbRNBFF+kQDVaaeaAgBALZyeqjTKoZGAQBxjR4hAMBOsdYhpCEEANgnivQJYoQAgHOI3TPLnAXECAEAcY0eIQDANrGYUE9DCACwTQNNsdagGBoFAMQ1eoQAAPvE4MMyNIQAANto1mN+jo+k0hACAGyjCU2z2MWzWt52xAgBAHGNHiEAwD5RxAidRkMIALCP9TxCxxtOhkYBAHGNHiEAwDZMug0AiG/ECAEA8UwTUfQIG6gu9UWMEAAQ1+gRAgBsE4MzrNEQAgBsFIMtIUOjAIC4Ro8QAGAbTYsiHYL0CQDAOYMV6gEAcY0YIQAAsYUeIQDANtEk1DvdJaQhBADYRmP1CQAAYgsNIQAgrjE0CgCwTVTLMDVQXeqLhhAAYJsoYoROt4MMjQIA4hs9QgCAfWIwoZ6GEABgH+sxQseDhAyNAgDiGj1CAIBtNOsdPKc7hDSEAAD7RJM+4XSQ8Cw1hC6Xa+6j9y/750tn53KxbufOnRdccIH1Nb3i16lTp44ePdqiRQunKxJLDhw4kJKSkpqa6nRFYsmOHTvy8vKcroVjvvjFL+677z51mUF5za2eNnlbzo7ExGgrZQMtGAyenSvt3Lmztrb27Fwr1lVXVyc6+rGIOcFg0OfzJSQkOF2RWFJTU+P1evl/y5I4/908//zzk5OTbT9tMBisqalx8MaevYYQAIBGiKdGAQBxjYYQABDXaAgBAHGNhhAAENfc9957r9N1iGvBYHDdunUvvfTSW2+9VV5e3rlzZ4/nx5yWH3744amnnnrjjTeCwWB+fr6z9Wyc9u/f/8orr6SkpGRlZel7tm7d+vjjj7/77rtZWVk5OTnOVq+x8fv9r7zyypw5c9auXZuVlaXftGAwOG/evNmzZ5eVlXXt2jX08YMQwufzvfrqq3Pnzl27dm3r1q2bNWum7z906JD+uymE4HfzHECP0GF79+4dPXr0oUOHWrdu/dxzzw0aNMjv9wshfD7fpZde+sUXX7Rt2/bWW299/vnnna5pYzRp0qQpU6YsW7ZM39y1a9cll1xSU1OTkZFx+eWXr1mzxtnqNSp+v/+qq6565plnWrVq5fP5Vq1ape+/8847H3roofbt2y9evHjkyJHOVrKx+eUvfzlz5syLL7742LFj3bp1KysrE0LU1NT079//q6++uvDCCydMmDB79mynq4kzFoSjfD6fz+fTXx87diwpKenzzz8PBoMLFizo1KlTbW1tMBhcvHhx27Zt9dcImT9//g033HDZZZf9/e9/1/dMnTr1xhtv1F/fc889o0ePdq52jc7MmTO7desW+rDpjhw50qRJk40bNwaDwZMnT6alpa1fv96hCjY6tbW1Xq/3q6++0jcHDhz45JNPBoPBefPmXXzxxYFAIBgMLlq0qH379vprxC56hA7zeDyhwSi/3+/3+5s0aSKEWL58+aBBg1wulxCioKBg+/btu3fvdrKijcz3339/7733Pv744+E7ly9fXlhYqL8uKCgI9RQhhHj77bdvuummt99++7HHHlu9erW+c/369U2bNu3SpYsQIjk5uX///ty0EJfL1aFDhy+++EIIceTIka1bt+o3avny5YMHD9YnIigoKCgrK9u7d6/DdcWZoSFsRCZPnjxixIiLLrpICFFRURGKeyUmJqalpVVUVDhau8bld7/73ZQpU6Q51cJvWnZ29sGDB2tqapyoXWO0ffv2f/zjHwsXLjx27FhRUdHMmTOFEJWVlaE7JoRo0aIFf9PDvfnmm/fcc0+HDh3atGkzadKkgQMHCuPHLCUlpUmTJvxuxjoC443Fvffeu379+qVLl+qbHo8nfEY65g8Lt2TJkoqKinHjxkn7PR6PHmEVQvj9frfb7Xa7z3blGiuXy9WpUyc9oNWzZ89x48bdfvvt5o9ZPM8fJqmtrR07dmxRUdHEiRO3b9/+X//1X3369BkwYED4x0wI4ff7+d2MdTSEjUJJSck///nPjz/+OCMjQ9+Tm5tbXl6uvz58+PCJEydatmzpXAUbl7lz5+7YsaN3795CiC1btuzYsWPv3r3Tp0/Pzc0NdWjKy8tzcnJoCENatWrVuXNn/XWXLl0OHjyof6gqKioCgYA+CF9eXt63b19Hq9mIfP31159//vny5cs9Hk/Hjh1vuOGG2bNnDxgwIPxjdvDgwaqqKn43Yx1Do857/PHHZ82a9cEHH4QP9I0YMeKdd945fvy4EOL111/v3bs3yQAhDz744IIFC5599tlnn322ffv2N95448033yyEGDFixIIFC4LBoBBiwYIFI0aMcLqmjUhRUVEoNPjZZ5/l5eWlpqb27t3b6/V+9NFHQoi9e/euWrVq2LBhjlazEcnMzPT5fHv27NE3t27d2rx5cyHEiBEjlixZcuLECSHE66+//rOf/Uzfj9jFpNsO++6779q3b5+Xl5eZmanveeSRR6644gohxMiRI7du3dq9e/clS5bMnz9/0KBBjta0kbr88suLi4snTpwohDhy5Ej//v1btGiRnp6+atWqFStWxPOKOZKTJ08OHDgwOTm5bdu2ixcvnjVrlv6Pwpw5c6ZOnTps2LDly5cXFRU9+uijTte0Efntb3+7aNGi4cOHb9u2rbS0dOXKlRdccEEwGLzmmmt27drVrVu3JUuWLFiwQP+FReyiIXRYVVXVpk2bwve0a9cuPT1dCBEIBJYtW7Zv377+/fu3atXKoQo2dlu2bMnIyAg9vFBVVfXBBx9UV1cPHjw4LS3N2bo1Nj6f7+OPPz5+/Hjfvn3DR/O2bNmyfv36tm3b9unTx8HqNU4bN24sLS1t1qxZ//79QysQhX43L7300tzcXGdriDNHQwgAiGvECAEAcY2GEAAQ12gIAQBxjYYQABDXaAgBAHGNhhAAENdoCAEAcY2GEAAQ12gIAQBxjYYQABDXaAgBAHHt/wGEur/X8/k+twAAAABJRU5ErkJggg==",
"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 = 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",
"# Collect all the terms, build and run the model\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",
"heatmap(scfres.ρ[:, :, 1, 1], c=:blues)"
],
"metadata": {},
"execution_count": 1
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.5"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.5",
"language": "julia"
}
},
"nbformat": 4
}