{
"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.510825e+01 1.433703e+01\n",
" * time: 0.0037479400634765625\n",
" 1 6.261635e+01 1.059977e+01\n",
" * time: 0.010985851287841797\n",
" 2 5.596146e+01 1.475653e+01\n",
" * time: 0.02825188636779785\n",
" 3 4.086197e+01 9.911974e+00\n",
" * time: 0.13123297691345215\n",
" 4 3.058236e+01 8.450051e+00\n",
" * time: 0.15758395195007324\n",
" 5 2.562763e+01 7.538135e+00\n",
" * time: 0.1791839599609375\n",
" 6 1.224479e+01 2.494525e+00\n",
" * time: 0.20121097564697266\n",
" 7 1.185749e+01 4.313602e+00\n",
" * time: 0.21592402458190918\n",
" 8 9.436821e+00 2.333889e+00\n",
" * time: 0.28729701042175293\n",
" 9 8.910663e+00 2.473528e+00\n",
" * time: 0.30191898345947266\n",
" 10 8.043010e+00 2.030384e+00\n",
" * time: 0.3199429512023926\n",
" 11 7.797338e+00 2.700208e+00\n",
" * time: 0.3357889652252197\n",
" 12 7.444471e+00 2.128942e+00\n",
" * time: 0.35129785537719727\n",
" 13 7.127486e+00 2.012159e+00\n",
" * time: 0.36583590507507324\n",
" 14 6.892565e+00 1.603986e+00\n",
" * time: 0.3803558349609375\n",
" 15 6.673610e+00 1.269919e+00\n",
" * time: 0.4379448890686035\n",
" 16 6.613841e+00 2.103153e+00\n",
" * time: 0.4540269374847412\n",
" 17 6.458632e+00 8.806043e-01\n",
" * time: 0.47229599952697754\n",
" 18 6.421592e+00 1.645461e+00\n",
" * time: 0.4867708683013916\n",
" 19 6.289115e+00 1.921790e+00\n",
" * time: 0.5013699531555176\n",
" 20 6.169125e+00 1.257745e+00\n",
" * time: 0.5160658359527588\n",
" 21 6.027653e+00 1.195836e+00\n",
" * time: 0.5306010246276855\n",
" 22 5.918568e+00 9.337476e-01\n",
" * time: 0.5638070106506348\n",
" 23 5.811785e+00 6.976539e-01\n",
" * time: 0.5778148174285889\n",
" 24 5.754954e+00 5.267787e-01\n",
" * time: 0.5922069549560547\n",
" 25 5.693027e+00 6.672500e-01\n",
" * time: 0.6029808521270752\n",
" 26 5.661399e+00 8.801142e-01\n",
" * time: 0.6139259338378906\n",
" 27 5.648746e+00 7.779108e-01\n",
" * time: 0.6248049736022949\n",
" 28 5.632913e+00 5.085560e-01\n",
" * time: 0.6399288177490234\n",
" 29 5.607723e+00 4.356511e-01\n",
" * time: 0.6522359848022461\n",
" 30 5.591711e+00 4.181805e-01\n",
" * time: 0.6629300117492676\n",
" 31 5.581046e+00 3.778050e-01\n",
" * time: 0.69038987159729\n",
" 32 5.577136e+00 2.760850e-01\n",
" * time: 0.7009599208831787\n",
" 33 5.573259e+00 2.752647e-01\n",
" * time: 0.7117528915405273\n",
" 34 5.570536e+00 2.239428e-01\n",
" * time: 0.7258179187774658\n",
" 35 5.568316e+00 1.454189e-01\n",
" * time: 0.7396500110626221\n",
" 36 5.565651e+00 1.651304e-01\n",
" * time: 0.7517080307006836\n",
" 37 5.564436e+00 2.092198e-01\n",
" * time: 0.7622559070587158\n",
" 38 5.564110e+00 2.471565e-01\n",
" * time: 0.7730698585510254\n",
" 39 5.563217e+00 1.175107e-01\n",
" * time: 0.7870948314666748\n",
" 40 5.562637e+00 1.292677e-01\n",
" * time: 0.8177399635314941\n",
" 41 5.561966e+00 1.095712e-01\n",
" * time: 0.8290579319000244\n",
" 42 5.561659e+00 1.091662e-01\n",
" * time: 0.8396730422973633\n",
" 43 5.561411e+00 8.586280e-02\n",
" * time: 0.8521609306335449\n",
" 44 5.561166e+00 5.440290e-02\n",
" * time: 0.8626489639282227\n",
" 45 5.561001e+00 3.317569e-02\n",
" * time: 0.8770749568939209\n",
" 46 5.560837e+00 4.469864e-02\n",
" * time: 0.8909668922424316\n",
" 47 5.560705e+00 2.592477e-02\n",
" * time: 0.9051389694213867\n",
" 48 5.560655e+00 3.082392e-02\n",
" * time: 0.932884931564331\n",
" 49 5.560612e+00 3.083745e-02\n",
" * time: 0.9481019973754883\n",
" 50 5.560581e+00 2.693332e-02\n",
" * time: 0.958873987197876\n",
" 51 5.560549e+00 2.192171e-02\n",
" * time: 0.9731018543243408\n",
" 52 5.560518e+00 1.826070e-02\n",
" * time: 0.9869489669799805\n",
" 53 5.560500e+00 1.779929e-02\n",
" * time: 1.001230001449585\n",
" 54 5.560486e+00 1.516361e-02\n",
" * time: 1.0153818130493164\n",
" 55 5.560485e+00 1.389700e-02\n",
" * time: 1.0261199474334717\n",
" 56 5.560479e+00 8.843074e-03\n",
" * time: 1.0542550086975098\n",
" 57 5.560475e+00 1.080395e-02\n",
" * time: 1.0683739185333252\n",
" 58 5.560471e+00 9.244356e-03\n",
" * time: 1.0824499130249023\n",
" 59 5.560467e+00 4.965850e-03\n",
" * time: 1.0965049266815186\n",
" 60 5.560466e+00 5.747469e-03\n",
" * time: 1.1073119640350342\n",
" 61 5.560465e+00 4.013168e-03\n",
" * time: 1.1219208240509033\n",
" 62 5.560465e+00 3.043593e-03\n",
" * time: 1.132699966430664\n",
" 63 5.560464e+00 2.298374e-03\n",
" * time: 1.146927833557129\n",
" 64 5.560464e+00 3.730200e-03\n",
" * time: 1.1752798557281494\n",
" 65 5.560464e+00 2.345714e-03\n",
" * time: 1.1857750415802002\n",
" 66 5.560463e+00 2.834053e-03\n",
" * time: 1.2000789642333984\n",
" 67 5.560463e+00 2.260962e-03\n",
" * time: 1.2109858989715576\n",
" 68 5.560463e+00 1.941850e-03\n",
" * time: 1.221524953842163\n",
" 69 5.560463e+00 1.110276e-03\n",
" * time: 1.2364559173583984\n",
" 70 5.560463e+00 1.174704e-03\n",
" * time: 1.2519810199737549\n",
" 71 5.560463e+00 9.657746e-04\n",
" * time: 1.2631080150604248\n",
" 72 5.560463e+00 7.808029e-04\n",
" * time: 1.2743608951568604\n",
" 73 5.560463e+00 6.031283e-04\n",
" * time: 1.3062939643859863\n",
" 74 5.560463e+00 9.867855e-04\n",
" * time: 1.3174388408660889\n",
" 75 5.560463e+00 6.201871e-04\n",
" * time: 1.328537940979004\n",
" 76 5.560463e+00 3.996932e-04\n",
" * time: 1.3397278785705566\n",
" 77 5.560463e+00 5.775584e-04\n",
" * time: 1.3507308959960938\n",
" 78 5.560463e+00 4.325642e-04\n",
" * time: 1.363562822341919\n",
" 79 5.560463e+00 4.241591e-04\n",
" * time: 1.3747169971466064\n",
" 80 5.560463e+00 3.225246e-04\n",
" * time: 1.3861379623413086\n",
" 81 5.560463e+00 3.093399e-04\n",
" * time: 1.4006319046020508\n",
" 82 5.560463e+00 1.713180e-04\n",
" * time: 1.4322278499603271\n",
" 83 5.560463e+00 2.974088e-04\n",
" * time: 1.442781925201416\n",
" 84 5.560463e+00 1.453810e-04\n",
" * time: 1.4582958221435547\n",
" 85 5.560463e+00 1.576106e-04\n",
" * time: 1.469101905822754\n",
" 86 5.560463e+00 1.187704e-04\n",
" * time: 1.4799468517303467\n",
" 87 5.560463e+00 1.043775e-04\n",
" * time: 1.4941918849945068\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -2.6967174168704244e-18\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -2.6967174168704244e-18\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
" 88 5.560463e+00 6.426286e-05\n",
" * time: 1.6433939933776855\n",
" 89 5.560463e+00 3.766925e-05\n",
" * time: 1.6932599544525146\n",
" 90 5.560463e+00 6.049655e-05\n",
" * time: 1.7052009105682373\n",
" 91 5.560463e+00 5.092156e-05\n",
" * time: 1.7168068885803223\n",
" 92 5.560463e+00 4.070035e-05\n",
" * time: 1.7282559871673584\n",
" 93 5.560463e+00 8.860389e-05\n",
" * time: 1.739408016204834\n",
" 94 5.560463e+00 5.879594e-05\n",
" * time: 1.7532639503479004\n",
" 95 5.560463e+00 6.064262e-05\n",
" * time: 1.7644739151000977\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -4.956352788505163e-19\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
" 96 5.560463e+00 5.787752e-05\n",
" * time: 1.7819039821624756\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -7.93016446160826e-18\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
" 97 5.560463e+00 4.820271e-05\n",
" * time: 1.7962238788604736\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -5.419679810653265e-18\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
" 98 5.560463e+00 3.891108e-05\n",
" * time: 1.8398759365081787\n",
" 99 5.560463e+00 8.089394e-05\n",
" * time: 1.8504128456115723\n",
" 100 5.560463e+00 5.232639e-05\n",
" * time: 1.8623018264770508\n",
" 101 5.560463e+00 4.646981e-05\n",
" * time: 1.8729689121246338\n",
" 102 5.560463e+00 5.396340e-05\n",
" * time: 1.8835728168487549\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -2.5861795497218356e-18\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -9.912705577010326e-19\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
" 103 5.560463e+00 5.396255e-05\n",
" * time: 1.9189469814300537\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -3.177095089133535e-18\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
" 104 5.560463e+00 5.369718e-05\n",
" * time: 1.9771080017089844\n",
" 105 5.560463e+00 5.070403e-05\n",
" * time: 1.9911308288574219\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -1.2390881971262907e-19\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -7.07107210761875e-19\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
" 106 5.560463e+00 5.070393e-05\n",
" * time: 2.0257468223571777\n",
"┌ Warning: Negative ρ detected\n",
"│ min_ρ = -9.46136862205366e-18\n",
"└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/densities.jl:5\n",
"e(1,1) / (2π)= 1.739179404220607\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXgUVbr48VO9ZGNJCElAwhIhrFFBkWUG1AgkIosEWTQjIyB6BZ1xuDOXy3VGkZFxQX1cEB03FkVERAUfBL0qsgjKJoosssgmhIRFZCdJd7p/f5S3f12nwul0U6HS9Pfz5I+uqtNVJ5VOTs556z1H8/v9AgCAWOWwuwIAANiJhhAAENNoCAEAMY2GEAAQ02gIAQAxjYYQABDTaAgBADGNhhAAENNoCAEAMY2GEAAQ01wX5zI+n++hiY+O/8c/Ls7lop3X63W5LtKP5pJRUVHhdDrtrkU0qaiocDgcmqbZXZFoEuO/m8nx7pBltm8/du6cN6zT7t3706FDm+69995I63WhLtJPtKysrMSjLdl79OJcDgBguVtbXxayzG23Ldy48UiYJ95cUHDCxoaQoVEAQEyL3T4+AMBymiaibridHiEAIKbRIwQAWEbTtHCfwPL7be5C0hACACzF0CgAAFGEHiEAwDpR+LAMDSEAwDIRxAiF3TFChkYBADGNhhAAENMYGgUAWCaChHq7R0ZpCAEA1okgRmj7wzUMjQIAYho9QgCAdbToS6inIQQAWEYLf6jT9qFRGkIAgGUiySO0uyUkRggAiGn0CAEA1okgRhiq/C+//PL2228fP368f//+11xzjbmAz+d77733tmzZ0r59+0GDBuld0o8++ujQoUOBMhkZGQUFBZWenx4hAMBKWvhfCqdOnercufPatWsdDkfPnj0///xzc5l777336aefrlOnzqRJk8aOHavv3L59+7f/5/HHH//www/Pdwl6hACAmmv27NmXXXbZ7NmzhRCpqamPPfZYXl5ecIH9+/fPmjVr3759DRo0uP3221u3bv2Pf/wjIyPjv//7v/UCpaWljRo1uuuuu853CXqEAADLaNpvz8uEQdknXLJkyc0336y/vvnmm7/66iuPxxNcYNmyZVdddVWDBg2EEE2bNs3Ozl65cmVwgQ8//DAlJeWGG2443yVoCAEA1olgYFQ5NlpcXJyRkaG/btCggc/nKykpCS5QUlKit4K6hg0bHjx4MLjAtGnTRo0apXiWlaFRAIBl9A5heG/RtHXr1g0dOjR4Z+/evfXBTKfT6fP59J36C6fTGVzS4XAECgghKioqggvs3bt3xYoVM2fOVFSAhhAAYLNGjRoNGTIkeE+bNm0Ch4qLi/XXBw8edLlcgQ5ioEBwF7C4uLhRo0aBzenTp+fn5zdp0kRxdRpCAIBlIlh9QgiRmZkpNYQBffr0ee655x5++GGn0zl//vz8/HyXyyWE+PHHH1NTUxs0aNCzZ8+77rpr9+7dzZs3//HHH/fv35+bm6u/1+fzvfnmm88++6z66sQIAQDWsTpGOGTIEIfDkZ+ff++99z711FMTJkzQ948YMWLOnDlCiIyMjL/+9a95eXljx469+eabH3zwweTkZL3MZ599dubMmX79+qmrTI8QAFBzJSQkrFy5ctGiRSdPnnzooYcCg5wvvfRS4BmZxx57rE+fPlu2bCksLOzSpUvgvdnZ2UuWLImPj1dfgoYQAGAZTYRIh6jsLSEkJCQMGjRI2nnttdcGb3br1q1bt25Smezs7KpUgIYQAGCd8GOEds+5TUMIALBQFK5HyMMyAICYRo8QAGCZCBLqbR8bpSEEAFgmgjxCu9tBhkYBALGNHiEAwDpR+LAMDSEAwDK/LawU5luqqTJVREMIALBMFHYIiRECAGIbPUIAgHUi6BLa3YWkIQQAWCeCGGE11aTKGBoFAMQ0eoQAAMto4SfI2/3QKA0hAMA6mhZ+OoTdLSENIQDAOlH4sAwxQgBATKNHCACwTBR2CGkIAQDWiWSKNbubQoZGAQAxjR4hAMA6UTg2SkMIALBMBHmEtjeEDI0CAGIaPUIAgHWicK5RGkIAgHWIEQIAYpmmRd9co8QIAQAxjR4hAMAymtDCTZCPmoT6M2fOPPTQQ927d8/NzX366acrKir0/UVFRXfeeee11147atSow4cPV1s9AQDRQAv/y25VbQj/53/+Z8WKFS+99NLkyZOnT5/+wgsv6PsHDhxYt27dmTNnCiEKCwurqZYAAFSTqg6Nrl+//p577mnfvr0QYujQoevXrxdCrFmzZufOnatWrXK73VOmTElPT9+6dWu7du2qsb4AgBosgodlbO8UVrVHWFBQMGfOnL179/74448fffRRQUGBEOKHH364+uqr3W63EKJWrVo5OTkbN26sxsoCAGo27bdpt8Njb52r2hDee++9paWlHTt27NKlS3Z29oABA4QQhw8fTklJCZSpV6/e+cKEfr8/EFYEAEQjj8cTutAlHCMcPHhw165djxw5cvToUZfL9ac//UkIUbdu3bNnzwbKnD59Ojk5udK3a5rmdDovvLoAALvo43+Xnio1hH6/f/ny5X/4wx8cDkdcXNztt9++bNkyIURWVtbOnTv1Mj6fb/fu3VlZWdVWVQBAjRfJ2KjNVa5SQ6hpWuvWrRcvXiyE8Pl8ixcvbtu2rRAiPz//1KlTixYtEkLMmzcvPj7+uuuuq9bqAgBqsigcGa3yU6NvvPHGH//4x5kzZ5aVlWVkZLz77rtCiPj4+BkzZowcOTIlJeX06dNvv/02458AgOhS1Yawa9euO3fuPHLkiNvtDn5Apm/fvkVFRSUlJZdddpnLxTw1ABDbqmHS7UOHDk2fPv3kyZO33HLL7373O3MBr9f71ltvbd++PScnZ9iwYQ6HI7D/3Xff3bRpU7169QYPHpydnV3p+cObazQ9PT24FdS53e4mTZrQCgIAIkqeULWEJ0+e7Ny5865du+rXr9+vXz89SCcZNWrU66+/3qRJkylTpujPcgohysvLe/Xq9corr6Smph4/fvzrr78+3yVovQAAlolghXp1+bfeeisrK+uNN94QQtSpU+eJJ57o06dPcIG9e/fOnTt3//796enpAwcOzM7OnjBhQsOGDV9++eWzZ89+/fXXIftprD4BAKi5li9fftNNN+mv8/PzV61aVV5eHlzgq6++6tChQ3p6uhAiMzOzVatWq1atEkIsWrTozjvv/OSTT6ZMmfLdd98pLkFDCACwTgSPjSp7hMXFxRkZGfrrBg0a+P3+kpISqYDeCgbKHDx4UAixZ8+eqVOnfvjhh4cPH77pppv0ObErxdAoAMAykSQGatq6deuGDh0avC8vL++ee+4RQrhcrsDEZF6vV5jy+oMLCCE8Ho9eQNO09u3bz5gxQwjRvn37v/3tbyNGjKj0+jSEAACbNWrUaMiQIcF7cnJyAoeKior010VFRS6XK9BBDBTQu4C6gwcPNmrUSAiRmZkZOMkVV1xRVFTk9XorjRfSEAIALBPB6hOaEJmZmVJDGNCvX7/JkydPmDDB5XJ98MEHffr00RPWv//++/T09MzMzF69et111107duxo1arVpk2biouLc3NzhRADBgxYsmSJfpK1a9dmZ2ef76kZGkIAQM01aNCgl156KTc3t0WLFosWLfrf//1fff+9995bWFg4duzYtLS0Bx98MC8vr3fv3osXL54wYULdunWFEKNGjXrzzTd79+7dtGnTBQsWvPnmm+e7BA0hAMA64ccI1cXj4+OXLl26ZMmSX3/9dfLkyQ0bNtT3v/HGG2lpafrrhx9+uE+fPlu3bh0zZkyHDh30nXXr1l2zZs0XX3xx5syZiRMn6uOllaIhBADUaHFxcTfffLO088orrwze7NixY8eOHaUy8fHxffv2DXl+GkIAgGWicYV6GkIAgKXsbtjCRUI9ACCm0SMEAFgmgoR621fmpSEEAFgmkjxCu4dSaQgBAJbRRIhllSp/k62IEQIAYho9QgCAdSJYod5uNIQAAOsQIwRqOr9iq0aT/1bY/bcDuGTQEAIALKOFnw5BjxAAcAkhRgjYILzRTuXxC3pzeEL/rTCWkC/tV52AcVTYJRrzCEmfAADENHqEAADLRDDFmu1dQnqEAICYRo8QUcAfIjSnChKG9dbQx0NVper8pv+C5W1/GIelg3IAURlQFPb/Rw7YiYYQAGCZSB6WqZ6aVB0NIQDAMpqIvmWYiBECAGIaPULUAKGjeuFEAaWjysPyUdO1/eGEHEME9SSafC558RppSx1SVF5Mk6pmKuw3BhFJQ0TkIkiot/sDRkMIALBO+DFC29EQAgAsE0EeITFCAADsRI8QNgiRjKcO65m2/X5V2C+8oz5lxUzlw5t7VI75mf8LNqYGSuUdxjCeJhWWjkqXDhWvlEKnUnkpwhhtA1+AGg0hAMAyEQ2NVlNdqoqhUQBATKNHiOoRYgQxzBwG9fCmz7Dt86kKy0eVm0IIeYeyJmohhysd6rFQ46ZDOVIqHzX+x1vJP+DyLtUgLbkWUIhgZhnbPzM0hAAAS9ndsIWLhhAAYBlNCzsdgvQJAADsRI8QFgkx7ZlfcVQd8zPv8ak3K8I5qixs3iPHFC8gRijF/IQpsCdvOpWb6sLKMwtzENFUwHBUmWshzxwXbaNkuFARrD5h94eEhhAAYJ0I5hq1G0OjAICYRo8QAGAZTWjyOipVeZOtaAgRqbCCguFMeybF7YQpUFdRYZgJrcJrPOoN56jxVOZLV1QoY4Q+VexTzhtU5gUKIZxSJM9l2HQ6DeM3Tpe0qYVx1CkPBTn8qk155EgKCoYVMhS2/9FD9YpkhXq7PxI0hAAAS0VZh5CGEABQsx04cOCNN944fvx4QUFBbm6uuUB5efm0adO2bt161VVXjRw50uVyCSHWrFmzbNmyQJnRo0cnJydXen4elgEAWEaLgPKEx48f79Kly9GjR7Ozs4cMGfLRRx+ZywwfPvzdd9+96qqrZsyYMXr0aH3nsmXL5syZ8+v/8fnOu7gMPUJUjXq2UBFeUDCs3D5hCux5jWE/r8d41FNh2Cw3HPXIhaWAovyrog4iSss2+aUpOoUUFDQUllL9RKgooMtt2HQbN11xqsIut9NQT5d8h6Vr+aU0RCGlIUrvDiNkKEg0vNRZPtfom2++2aZNm6lTpwohEhMTn3zyyQEDBgQX2LVr14IFC4qKilJTU/v375+VlfXoo482atRICNG5c+cnn3wy5PXpEQIAaq6vvvqqV69e+uuePXuuWbOmvLw8uMCqVas6dOiQmpoqhGjYsGGbNm2+/vpr/dC2bdsmTZo0ffr048ePKy5BQwgAsI4W0df5FRcXp6en668zMjL8fn9xcfH5Cuhl9AINGjRo166dz+d777332rRps2/fvvNdgqFRVEklU4lZNxYqZTh4TeOT0gCmp8ww+OkpN2yWS0fLpKPGoVHje6WBUxEqGcOUPqEeGpUyHExDoy714KdheDMu3lg4Xjpq2HTHGe+/8VRCCGms1OUP41/ksEZKhTm/oupXQlSIYGFeoa1atSovLy945y233PLnP/9ZCOF2u71er77T4/EIIeLj44NLxsXFBQroZeLi4oQQI0aMGDFihL5z8ODBkydPfvnllyutAA0hAMAyEcywpgnRsmXL8ePHB+9s3ry5/iIzM7OoqEh/XVRU5Ha7g/t/QohGjRodOHAgsHngwIHMzEzpEl26dFm6dOn5KsDQKADAZhkZGb2MAg3hgAED5s+fr8cF586d269fP6fTKYRYs2aNPtqZn5+/a9eurVu3CiE2bNhw+PBhPcXizJkz+hk8Hs+iRYuuuOKK812dHiEAwDJVSIio5D2KgwMHDvz3v//dvXv3rKysZcuWffHFF/r+Bx54oLCwcOzYsfXq1ZswYUJeXl6vXr0+++yzf/3rX7Vr1xZCXH/99fXq1UtLS1u/fn29evUefPDB812ChhCVMy0wFGL7goKCHinmJwfqpLBfealhs6zUa9g8J21K7zUclQOKpkur8yvkZZuUiQEOOUYoj8fICRLGjAgp7BeXYPjljU80HPXKR0OsNuX3O6Vt43E5pqgghwzN6RNSfoXfGDIkZngJsDR9wu12f/HFF8uXLz958uQrr7yiPx0qhHjrrbcCr8ePH9+vX7+tW7c++OCDbdq00Xd+9tln33777alTp/7617926tRJ0TzTEAIAajSn09mjRw9pZ+vWrYM3c3JycnJygvfUr18/Pz+/KuenIQQAWCaSSberpyZVR0MIALCMpomw0yfsbglpCPF/KkkVDDponmJNio4Zg2thBQWl3D4pjCdMYb9SafOMYfPcWePRsx7jqVTxRSnpUJiihtUaI5SCglJqYLwc9jPUPKHMHbzpTVKvNiX/4ssJkX71XwZVyFBeiy7U3zj5jyATsEW7yPInbEX6BAAgptEjBABYJgo7hDSEAAALWZ1HeBHQEOI36mWVKlmGyRdG4qA0fagUeJOCglIIUAhx7owhznfutKHA2dPlxqPGwmekkKEq6bC8khihYY8U7PQZY2/SPZN+ux3yKkvyL7/bLWUKGmOEiYbf1oQkw6YU3fR6DCHDilALXSkDxELTXNK2ccth3DSczGH6E6fJu1TLNtneV0AsCCNG+P333xcUFLRp0+bGG29ctWqVvnPv3r2DBg1q27bt7bfffvDgweqpJAAgeli39MTFUdWGcNu2bXl5ebm5uQsWLJgwYUJSUpK+v6CgoEWLFosXL05LS7vtttuqrZ4AgCig5xGG+2Wvqg6NPvroo3fcccfYsWOFEIEJbL7++uv9+/c//vjjLpfr6aefTktL27Rp05VXXlldlYWl1JOoqWdQE6GGRuU15Y0TlUkjkNL4pDQQKoQ4e8qw58xJw1joGePRs6cMR88aR0qlXAtpGFaavE2YMj1MQ6NS1oFxGSbj77e0JL15aNSlHBpNMA6NlpUah0bLjekT4dTTTKqZNJgpf1+aurB8cr98duUC90zAFm0imGs07Jii1araI1y3bl3r1q1Hjhx58803v/DCCxUVFUKIzZs3d+jQweVyCSESExNzcnI2bdpUjZUFAMBqVe0R7t+//6mnnpo6dWpycvLo0aN//fXXiRMnHjlyJDk5OVCmXr16hw8frvTtfr9fbzsBAFHK4/G43e4QhWpG2C8sVe0RJicn33333X379u3evfs///nP2bNnCyFSUlICCz4JIU6dOlWvXr1K365pmr6CFAAgSoVuBSOKEdrecFa1R9iiRYtA5y8lJeX06dNCiMsvv3z79u36Tp/P99NPPwWWUkT0kaOCxi1TUEkOChoTCaQYlUe5jpIUqJOyI4QpKHj6hGrzjBQjNEYQ1TOumWOEHmWM0DQzmSpGqMlTrIWXPmFaQMoYFJSXi1LV00wdzpQSP6S54qTCmtNQEym5QpiDiCFSN5RHAStUtUc4fPjwd9555+zZsxUVFdOnT9dXxOjVq1dZWdkHH3wghJg1a1bt2rW7detWjZUFANRsUfjQaJV7hKNGjVq3bl2zZs3cbnf79u1nzpwphIiLi3v77bfvvPPOBx54wO12v/POOw55UU4AQCyJwjnWqtoQulyuN95448UXX9Q0LSEhIbC/Z8+e+/fvP3HiREpKiu2PwAIAEK7wplhLTEw073Q4HOd7RgY1ywUkDpqDTKbEQdVCS9JEZfKySmdUU6YJU6agKUZYZih8UpVHKM24Zgq8mfMIpQWMVNmTIaZYU4bWhBAut2FAJa7MGCMsM/y2mipmnOxNGRQ0TXJmynGUN1WTw0mbcnzRIdfEId20cNIK5dL8710ThZ9HaPcPkrlGAQCWiSToZ/c/NDSEAADrRGGMkGdbAAAxjR5hDAmx9I4ycdBnjhFWSDHCcCYXNUbmzp2V8gjNc42qMgXlqUdPGqcePa3KI5QSB6X1oYRpASlTKM5QWLpHphk7DZtScp4QwuUy7PGWq26pPH2oXHHp0qGmOTXukWY9dcU5VJtuKYJoqIrTFAr1GaOGYaUVEiKs+SKZa9TunyQNIQDAMlE4MsrQKAAgttEjjF2mfAlpU5UYIEKmTxiH9TzyukvGKdbOqjIchGl4U5o1TT4q50tIc6qp8iW8pqFRaeo403pGhsJ+402VBnxMC9abrmUcLK2oMGyqUzUk6mFYaTBTCOGOM4yFuuMNNy0u3jj3W7zxphnf63IbM0zcpiH5EJ801ein7V0HhBaFXUIaQgCAZbQI1hekIQQAXDIiyCO0ux0kRggAiG30CC9pIRImVEcvMH3CI6dPGDeN6ROmkKEcI5TmYFPHFE0LLXkVm1K+hPRdCCEq5CwFdaBOlUAh/ZusVZjyCpzK+dvUs6bJAUgpHcI4eVu8vDhofILhtsQnGv4ylCdKgVXD0bgEw01zywkn8n/b0vflkOdYUyehMONajUeMEAAQ0yLII7R7wQaGRgEAMY0eIQDAMky6jZpFHRJUz6kmB8NCxggrVPOBeY3LMKlnXJPCeEKIUjnOZzybtJRSqepaHuW8ZeYYoTp7z6/O5pPSCqUFhky31O8zFAiVKWj4vqQ1nlxu6Y4ZgoKlCfIdTpDusLw6lXEFKONPU76HysxLEeqj5VfPuGb8i6kOJwJVREMIALCMFoXrERIjBADENHqEAADLRJJQb/egNg1hLAkjnhUiZChMwTOfca5RdVqhRzkTqRTkM+8JLwqozBSskJPeTN+mes7VcDI15Zw38y+/6eKK8t5yw7bHKSVuKu9YqDvsUd5SjzKwKn0SzBmoIXIxw/mU2j2ihkpEsgyT3S0hDSEAoEbbs2fPK6+8cuLEiYKCgt69e5sLlJaWvvTSS9u3b8/JyRkzZkxcXFzw0blz5x49evT+++8/3/mJEQIArKOF/6V07Nixrl27+ny+Tp06DR8+/P333zeXGTZs2KeffnrDDTfMnz//7rvvDj70zTff/OlPf3r88ccVl6BHGLtMI1KqMSnz9F7yMkzyqkyq0TOveqTUIw/ceY17pAJyqoZX9TS/NGWaPMBrnjNNzqcIYy0kiSa/N8TS7caRUbnmjgrVdy3fYfmOyXdYKhBWkol8h30hhkaVU8XJn0O/NE2d6q2oGazOI5w5c2aHDh2efvppIYTT6XzqqacGDx4cXGDHjh2LFi0qKSlJTk7Oz89v0qTJY4891qRJEyFEWVnZ/fffP2HChCeffFJxCXqEAADLaL89LhPWl8qqVat69Oihv+7Ro8f69evLysqCC3zzzTdXX311cnKyECI9Pb1t27arV6/WD02cOHHw4MFt27ZV15mGEABQcxUXF6elpemv09PT/X5/SUlJcIGSkpJAASFERkZGcXGxEOL777//7LPPxo0bF/ISDI0CAKwT0eoTq1atysvLC97Xt2/fsWPHCiHi4+PLy8v1nfqLhISE4JLx8fEez/9fc6asrCwhIcHr9d59991Tp051u90hr09DGENM4S+/4rA8eVglwTPVRFlyBFGZXKHerOzt0uxuxksrp/hSP7tfWcwv8qCgfCLjezV59jA5aqjOYPGF812rb2BlBdQ/IFVNQk7OJ+SMFOXKSlJhjZBhTaeFnxeoCdGyZcvx48cH72zZsqX+onHjxgcOHNBf79+/Py4uLj09PbhkZmbm/v37A5sHDhzIzMzcsWPHDz/8MGzYMCHEuXPnDh8+3KJFi2XLlumxQwkNIQDAZhkZGb169ar0UEFBwUMPPfTwww8nJCS88847AwYMcDgcQogVK1Y0bty4efPm+fn5I0eO3LhxY/v27desWfPLL7/ceOONLpdr27Zt+hlWrVr1X//1X59//nnDhg0rvQQNIQDAMpoWfoK8svyAAQNef/31zp07Z2VlrVu3bsmSJfr+cePGFRYWjh07Njk5+bHHHrvpppuuv/765cuXT548OSkpSQjRvHlzveTu3btdLldg04yGEABgHatXqHe5XIsXL167du2xY8fefvvtunXr6vvnzp2rPykqhPjLX/7St2/fbdu2TZ48+fLLL5fO0K1bt6+//lp1iTDri5gRKngWatkm1VE5e08Z/TLvCbGOj3Rp9dRxIUOhFxAUVDOf2ZRoKJVXRxANm/LqUaHvsDHEqPwBhfWzDvnJCbVaGKJMJHONhirgcDi6du0q7czKygrezM7Ozs7OrvTtiYmJzZo1U50/VAUAALiU0SMEAFgmGifdpkcIAIhp9AjxmwuM1KhnLpUnMlVGiSqLKqkyz5T5kOY5VKM3JKW6pfI3qQyFypHRkHE7+Q6HUZNwRe+PB9GLhhAAYCVrJ92+CGgIAQCWiSRGWE1VqTIaQvzmAj+LmrypKTbNpQ1bpqpIv1emzTBqYjoeRUNxqlsqf5PyLVXdQHP5ED+gcGoSLtv/JuJCWZ1HeBHwsAwAIKbRIwQAWCYa0ydoCAEAlolgZhnb0RDiPELH7YybDmOgyKE66lBvOuURe4dTOrnq7fKlld+I9K+o3/xtWrcMk3xm8x+LsEKh8v03bJruiXSH5WtL91z9AwrrZx3yk2N7fAggRggAiGn0CAEAlonGGCE9QgBATKNHGEPkBDo5ymQIf8n/o5kzz9RRKOOm0yVtOpSb8rVMb1dFvORwl9OwSpBWYay2T/qupSvLix9pWuQhQ/MtVBdQR2FDfdfq+x/yDofxA1JHcKVNIUKEQtUZkcQTo4Kl6/JeDDSEAADrRGFCPQ0hAMAyEcQIbUeMEAAQ0+gRxq5Qc3Iawl/mQI8cBZSDgqogk8tt2HTHGTfdTulaLuMeqYB0NpfxWhXG4JnPaQzy+YzfmCnm53NIQUEpZBV5kND8T7PDuEu+w1LYz6n6ruU7LN8x+Q6744x3OE51NvmH65Q2Vd+FqOyzFMzamUtx8UWQUG97B5KGEABgnSiMETI0CgCIafQIY4l6AaJwnt0X5lm4lM/fu9VjofGGcbm4BHngTtojb5YZNr3lhnyJigrj0KhPnf/gE8odppXf5eIK6inTRCU5J8rh5TjVLY1T3tJK7rCxvPQTMY1dqyom53WYPjnqGdpCdA7s7jogJK2SOEvItzDpNgDgUqFFkEdYPTWpOhpCAIB1IosR2rpCNjFCAEBMo0d4KVPHBMMKClYSI1Q/zS+lNBifzpciUvEJhs9hfKL8seCJBzwAAB2QSURBVEww7ikv9Ro2ywxHvR5DWE8OCvpU/3mah3QqKlRv98vrNqniriHDrnKChOtCooDSLTUcTTDdYemeSz8R6VrST1NOrlBO9iYqWadJGToNkWuBmieChHotvEQky9EQAgAsE0keYfXUpOoYGgUAxDR6hAAAy0S0HmE11aWqaAgvaaaFl6RtRWlpvaFKJspSB7TcUkDLuClHsAwxv4Qk+WNZVmqMEZa5gzfloGBFGCslaY4Kw6nKTd9mRRgn9xtvqTxbmPHc5uCZHGdVBgWlqF5CkuGeJNZSHU2oZYrCJqliinFyyDDytEJh+iyFXJ1KddTuP6C4NNAQAgAsE0GM0PZ/aIgRAgBiGj1CAIBlNBF+jNDuLiENYexSh2ZM2V1yqM00MaZhU54JU45vGdPajDE/T5k84aen3BjJk4OChk2/abrQYHKgzvhdeJzym73eMK6lnL1VnlTT4ZTHY1zqxEFjpqAU9kuqrdysozoqQsUUpZ+XW04rlGKEoeYaDfFJMx4ViDaX/OoTXq/3D3/4w+jRowN7duzY0bt378svv7x///779u2zunoAAFSv8BrCyZMnb9q0ac2aNfqm3+8fOHBgly5dVq9e3aZNm6FDh1ZDDQEAUUN/WCbcL3uFMTS6bdu2efPm/fWvf50yZYq+Z+XKlYcPH54wYYLT6Zw0aVJ6evr333/foUOH6qkqLlSIZAqpsHIIUYRKn3DJ6RPGlZKMj/5LCyd5PYaBUCGE1xsXvKme9kxiGpA0DuG6pdnaTJdW5maoF3VSz6nmNOUVqG+alMMgDWZKo5216hruWK06hs0k46YQItH49kQ5m0I545pxvfsLTZ9Qbtr9BxNVEMEUa3b/YKvaEPp8vnvuueell14KHv/csmVLhw4dnE6nECIhIaFdu3b6nmqpKQCg5quGGOGOHTumTp16/PjxgQMHDhw40FzgzJkzzz///JYtW9q3b/+Xv/wlISFBCLF27doPPvjgwIEDqampQ4YMuf766893/qoOjT777LMdOnTo1q1b8M5ffvmlTp06gc2UlJQjR45U+na/3+/1eis9BACICrb8GT9y5Ei3bt2Sk5Nvuumm+++/f86cOeYyhYWFX3/99a233rpkyZIRI0boOzdt2lS/fv1bbrmlcePG/fr1+/jjj893iSr1CH/++efnn39+0aJFu3fvPnz4cFlZ2e7duy+//PJ69eqdPn06UOzkyZOpqamVnkHTNJeLJ1QBIIpV5c+45R3CGTNmdOrUadKkSUIIn8/3zDPPFBYWBhfYtm3bF198cejQoTp16vTo0aNRo0b79u1r1qzZqFGjAmX279//8ccf9+vXr9JLVKlxOnjwYHx8/K233iqEOH369PHjx/Py8jZv3ty8efNt27b5/X5N0yoqKnbu3NmiRYuqnBD2UAYJ5XVQpIfaTWMHIdInjHEjd5zh5PEJhs0KY4qCFAIU5onNQgQFjY/vK9eHKjtnqGd5qRwj9BgDlhVe9apMxinW5HuoumNCCLdbtbKSnHOSpEqQkIKCtZOlkKGcPpFU21BAmoNNWpXJrV6GKWT6hHoZJqk0QcJoo5l+pqHfoiy+evXq3Nxc/XVubu7w4cNLS0v1wc9AgauvvlofnkxNTc3JyVmzZk2zZs0CBfbv379y5cq//e1v57tElYZGu3btuuv/PPfcc+3atdu1a1diYmKvXr18Pt/s2bOFEK+//npqaurvf//7qpwQAICqKCkpqV+/vv46LS3N7/eXlJScr4Bepri4WH89Z84cTdOaNm165ZVXDhs27HyXCHuKtdq1a2dmZuqvXS7Xu+++O2HChJSUlOeff/6dd94J/2EhAMAlRIvka9myZdcaTZ48WT9fQkJCeXm5/rqsrEwIkZSUFHzBxMREj8cT2CwtLQ0UKCws9Pv9O3fu/Omnnx555JHzVTnsuN0tt9xyyy23BDa7d+++e/duqaMKAIhNkeUFdujQ4ZlnngneE+hxNW7c+Oeff9Zf//zzzwkJCWlpaVLJQAEhxP79+xs3bhxcIDs7e8yYMS+88MKjjz5a6dWteYCFVvBSEE5aoTDHCI2hOL/LmG9nzDyT0u98FS7FUWGKvcl1U+bnqZPzShOMeYSmGKFXGSM0rcqkihGaMi/NeYSqGGGCtO6SMo9QyhSUgoJSlqEQIrG28eRS4mCCcoo1l+qTEHKKNXXiIKJOROsRaikpKR07dqz06KBBg8aNG/fwww8nJSW9/fbbBQUFDodDCPH55583a9asVatW+fn5d91117ffftuxY8eVK1eeOHFCjykePnw4IyNDCOHz+ZYsWdKqVavzVYAnOQEANVe/fv1mzJhxzTXXNG7cePv27UuWLNH3P/TQQ4WFha1atapbt+5TTz3Vu3fvLl26rFmz5tlnn01MTBRC9OnTx+fzpaen79ixo27dugsXLjzfJWgIAQDWsTp/wul0Lliw4Pvvvz9x4kTnzp31Rk4IsWDBglq1aumvR48e3a9fv23btuXk5Fx22WX6ztWrV2/evPmXX35p1KhRmzZtFP1UGkIAgIW0cJdVqkp585xlgQZP17hxYyk06HK5qjjTGQ1h7JL+PfKHk1YozHN4GovLU48aj/p9TuOmMdJWeX2D66KKvZkWMzJm4xmDggnnQsQI1XmE8qpMynumnp1VhM4jlFZKUk09miivyhRnPCr/4ifWktZdUk8uKiUOqoKCWiUxwjASB3kOHRcBDSEAwDJaFP77QkMIALBOFC7MS0OI35jmX5N3yOWlmcyk0saBVpffIR02boX4HEpVMWVuKGd3izdcWhr0Kys1DI16TMsweYxLREmzwcmZHsol6k2T0pmGRo1LvUtZCtI8Z+oZ16R1lKRcCyk7QpjHQuU51YwD3cqFlqRN89CoKV9CPVKKKBNZ+kQ1VaaKwp5ZBgCASwk9QgCAZSKYWcb2YQAaQgCAdYgRIoopY4Kh/8VTpgqYOBXHNE3+WKpXVpJTNeJUc6qVJxrzJcoM15IigkIIr+fixQiltAS38huJU4YMTQFF1ZRp5pPLQUFj2FVeaEkKCkrrJplDhOqPht1/ExGDaAgBAJaJ4GEZ2//7oSEEAFgmCkdGeWoUABDb6BGicuoJ2ESoOdgc4f2LZQxZVTKdm3Ry1XpGpnWXDKmBUlDQW2446vHIMUIpKChNsWaaHM64DJM0W5iy2sIUNXS7pdinMUYYr0o6NMX8VOsoCdNNkzIF1UFBeU610HmBTKJ2qYu2nykNIQDAOlGYUE9DCACwTCR5hHb3IIkRAgBiGj1CVIn5Pza/HApS5dCFFTLUNLm0NF1lWDFCrxRaSzDE/NRpguY9UuKg31g8VIzQUNhhmoTTKUfmpGQ+dQRRVVhKBHSFCk+qpw8NERRUTiUqoi9+hDBF8Nio3WgIAQCW0bTwY352j40yNAoAiGn0CFE15omypMnFLmCk1DS2FmLJJ81pGJE0rfxuOOpyG87mltIhjEvM+7zypSvksVB1+oSx2tJmqKXbnVKWiHEA06mcWE4aHw5x1GkafFYmqGghfl5CnH+7kn/1o23cDGGJYGFeuzuENIQAAOtoQtPC/GfH7naQhhAAYKEonGONGCEAIKbRI0SkpNV2wgoZyhO4GbbM/52ZVvYxPuvvMLzfaQwZ+owxQl+Fw7hpPOqTY4SmfAnjplxcRf4uTDFCdaAurJSGEEdNl1bHL+Way29WBgXt/mcfF1v4CfW2oyEEAFgmgmWYNPMDchcXQ6MAgJhGjxAAYJ0ofFiGhhAWuZCQoXQmUwTLr4yuOYxxO59DFcYLEfMzxwilHX7LYoTmQIpD+W2q16KSTmZKBAyx8pE8liUHZeXi5y9r/x812CuSSbdNObgXGQ0hAMAyEeQR2v7fEzFCAEBMo0cIALBUlHUIaQhRTZSZZ36/OmRlihdIIUY5AGksa4yH+Y1xPIexdLh5gdLZwotsyIG3EL/96rxDU25lGEcrmTlWfe0wDiLWsTAvAABRhh4hAMAykSTUV1NVqoyGEDaQf02UuRYiVLpFiIFTv+pifillwXQdv6lyioqol2GSVTI+eQE5DHJuhvJSodZGIiMCkSOPEAAQ0yKYa9TuhpAYIQAgptEjBABYRqvCE9Gmt9jcJaQhRA0Q4tl+U5xP3lJHEFXbISe990vvD2tSNTV10K+S8uEcDPfcdo9NAQqbN2+eMmXKiRMnCgoKCgsLzQVOnjz51FNPbd++PScnZ9y4cbVq1RJCrFmzZv78+Xv27MnIyLj77rvbt29/vvMzNAoAqLkOHTp0/fXXZ2VlDR06dPz48W+99Za5zG233bZ169aRI0euX79++PDh+s5JkyYlJCQMHTo0NTX1d7/73caNG893CXqEAADLRLYeoeLo9OnTu3fv/ve//10IUV5ePnny5DvvvDO4wJYtW1asWHH48OFatWr9/ve/b9iw4e7du5s3b75w4UL9zIMGDfr+++/nz59/vk4hPUIAgGX0mWXC+lKP3K9du/a6667TX1933XU//PBDaWmpVODqq6/Wh0NTUlKuuOKK9evXi6D21e/3Hzhw4LLLLjvfJegRIgqE+P9SmjUtRMTRrzpsntxNvpRlwbTQJ7qAmCGJgLBNBHmESiUlJfXr19dfp6Wl+f3+kpKSrKysQIFDhw4FCggh6tevX1xcHHyGqVOnnjp1SupHBqMhBADYbNmyZddee23wnltvvVUfDk1KSgp0Ac+dOyeE0Dt/AYmJiWVlZYHN0tLS4AJz58594oknli5dmpiYeL6r0xACACwTUYxQdOjQ4Zlnngne2bRpU/1FkyZN9u3bp7/et29fYmJicP9PKqCXady4sf56/vz5Y8eO/eyzz1q3bq2oADFCAIB1wo8RappISUnpaJSenq6fb/Dgwe+///7p06eFEDNnzhw0aJDD4RBCfPzxx1u3bhVC5OfnFxcXr169WgixdOnSM2fO3HjjjUKITz/9dPTo0QsXLrzyyivVVaZHiOgXKg3RQAooKk8V6lIXFWE/xKY+ffrMmjXrqquuyszMPHDgwJIlS/T9kyZNKiwsbNeuXe3atZ977rl+/fpdffXV33333YsvvhgfHy+EuO+++8rKym677Ta9/NChQ5944olKL0FDCACouRwOx9y5c7dv337s2LGOHTvGxcXp+z/55JOEhAT99ciRI/v27btz587WrVunpaXpO7/77jufzxc4j946VoqGEABgGU1EsAxT6PLmIF9qamrwZkZGRkZGRvCe5OTkKlaAhhAxJqxxVABhimCFett/D3lYBgAQ0+gRAgCsw8K8AIBYZvlcoxcBQ6MAgJhGjxAAYJkIHpaxe2SUhhAAYK1oixEyNAoAiGn0CAEAVqpKgnyNQkMIALBMJDFCu9tNGkIAgGUiSZ+wuwdJjBAAENPoEQIArBOFM8tUtUf4xBNPXHvttRkZGR07dpw7d25g/+bNm2+44Yb69ev36tVr586d1VNJAEB0iGBV3qhpCH/55ZcXXnhh69atDz/88MiRI9euXSuE8Pl8gwYN6tOnz759+7p16xZY/xAAgGhR1YbwmWee6datW1paWkFBQadOndasWSOEWLFixfHjx8eNG1e7du1//OMfu3btWr9+fXXWFgBQo2m/PS4TDru7hGE/LHPs2LFNmzZdc801Qogff/zxqquucjgcQoi4uLi2bdtu27bN+joCAFBtwntYxuv1Dh8+vH///t26dRNCHDt2rE6dOoGjycnJR48erfSNfr/f6/VeSEUBAPbyeDxut1tdJhrzCMPoEfp8vhEjRng8ntdee03fk5qaevr06UCBEydO1K9fv9L3aprmcvGEKgBEsZCtYJSqauPk9/vHjBnz888/f/LJJ/Hx8frO7OzsLVu2+P1+TdO8Xu+OHTtatmxZbVUFANR44SfU2x0irHKPcMyYMV999dVLL7106NCh3bt3Hzt2TAjRo0cPl8v12muv+Xy+KVOmNGzYsEuXLtVZWwBAzaaF/2W3qjaEK1euLCsrKygoyMvLy8vLmz59uhDC6XTOmzfv5Zdfrl279jvvvPPuu+/avtAwAABhqerQ6ObNmyvd37lz540bN1pXHwBAFIvGh2V4gAUAYJlIJt22uyWkIQQAWMruHl64WH0CABDT6BECACxTM54DDQ8NIQDAMtEYI2RoFAAQ0+gRAgCsE4UL89IQAgCso4U91Gn70CgNIQDAMpEk1FdPTaqOGCEAIKbRIwQAWIoYIQAgZmlC02xv2cLE0CgAIKbRIwQAWIbVJwAAsa0a5ljbsGHDc889d/z48YEDB951113mAr/++utjjz22ZcuW9u3b//3vf69bt64QoqSk5Msvv9ywYUP9+vUffPBBxfkZGgUAWEaLiOKExcXFPXv2vOaaa+67775//etf06ZNM5cZMmRIUVHRuHHjdu7cOWzYMH3np59+Onv27K1bt86fP19dZ3qEAICaa9q0abm5uf/5n/8phHjyyScfffTRUaNGBRf44YcfVq9efeTIkcTExGuvvbZBgwY//fRTdnb2iBEjRowYMWvWrBdffFF9CXqEAADL6DHCsL7UQ6nr16/v1q2b/rpbt25btmw5d+6cVODqq69OTEwUQtStW/eKK6749ttvw6ozDSEAwDpa+F9Khw4dSk1N1V/Xr19fCFFSUnK+AnoZqUBIDI0CAGz2+eeft2jRInjPsGHD/vnPfwohatWqFegC6i/q1KkTXLJWrVqlpaWBzbNnz9auXTusq9MQAgAsFMl6hL///e9feeWV4J1paWn6i6ZNm+7du1d/vWfPnqSkJL1fGNCkSZNAASHE3r17mzRpElYFGBoFAFgmgpFRTYhatWo1N9JTIIQQQ4cOnTdv3okTJ4QQ06ZNGzJkiN7Qfvjhh5s2bRJC3HTTTYcPH/7qq6+EEJ999llZWVlubm5YdaZHCACwjtXrEebn5+fm5rZr165hw4anTp36/PPP9f2TJ08uLCy88sork5KSpk6dWlBQ0LZt223btr366qtxcXFCiKVLlw4aNKi8vLysrCw1NfWmm26aM2dOpZegIQQA1FwOh2PmzJn79u07ceJETk6O0+nU93/55Zdut1t/fccdd/Tr1++nn35q1apVIILYvXv3Xbt2Bc4TKGxGQwgAsEzIBPnINGvWTNpTq1at4M3k5OSOHTsG73G73fXq1avKyWkIAQCWica5RnlYBgAQ0+gRAgAso2kigvSJaqpMFdEQAgAsZfdQZ7gYGgUAxDR6hAAAK0Vbh5CGEABgnQjSJ4gRAgAuIVbPLHMRECMEAMQ0eoQAAMtEY0I9DSEAwDLVNMVatWJoFAAQ0+gRAgCsE4UPy9AQAgAso4Uf87N9JJWGEABgGU1oWphdvHDLW44YIQAgptEjBABYJ4IYod1oCAEA1gk/j9D2hpOhUQBATKNHCACwDJNuAwBiGzFCAEAs00QEPcJqqktVESMEAMQ0eoQAAMtE4QxrNIQAAAtFYUvI0CgAIKbRIwQAWEbTIkiHIH0CAHDJYIV6AEBMI0YIAEB0oUcIALBMJAn1dncJaQgBAJbRWH0CAIDoQkMIAIhpDI0CACwT0TJM1VSXqqIhBABYJoIYod3tIEOjAIDYRo8QAGCdKEyopyEEAFgn/Bih7UFChkYBADGNHiEAwDJa+B08uzuENIQAAOtEkj5hd5DwIjWEDodj1jP/Wv7eWxfnctFu3759TZs2DX9Nr9h17ty5kydPNmjQwO6KRJMjR44kJSXVqlXL7opEk71792ZlZdldC9t894c/TJo0SV2mZ1ZauKdN3N1wb3x8pJWygOb3+y/Olfbt21dRUXFxrhXtysrK4m39WEQdv9/v8Xji4uLsrkg0KS8vd7vd/L8Vlhj/3bzssssSExMtP63f7y8vL7fxxl68hhAAgBqIp0YBADGNhhAAENNoCAEAMY2GEAAQ05wTJ060uw4xze/3r1+//q233vroo4+KioratWvncv2W0/Lrr7+++OKLH3zwgd/vz87OtreeNdPhw4dnz56dlJSUnp6u79m1a9fzzz//6aefpqenN2zY0N7q1TRer3f27NkzZ85ct25denq6ftP8fv+cOXNmzJixc+fOq666KvDxgxDC4/G88847s2bNWrduXZMmTerVq6fvP3bsmP67KYTgd/MSQI/QZgcPHhw8ePCxY8eaNGny+uuv9+zZ0+v1CiE8Hs9111333XffNW/e/J577nnjjTfsrmlNNGbMmHHjxi1fvlzf/Pnnnzt37lxeXp6amnrDDTesXbvW3urVKF6vt3fv3q+88krjxo09Hs/q1av1/X//+9+feOKJli1bLly48NZbb7W3kjXNH//4x5dffvmKK644depU+/btd+7cKYQoLy/v3r37xo0bL7/88lGjRs2YMcPuauKC+WErj8fj8Xj016dOnUpISPj222/9fv+8efPatm1bUVHh9/sXLlzYvHlz/TUC5s6de/vtt19//fX//ve/9T3jx4+/44479NePPPLI4MGD7atdjfPyyy+3b98+8GHTnThxonbt2ps3b/b7/WfPnk1OTt6wYYNNFaxxKioq3G73xo0b9c0ePXpMmTLF7/fPmTPniiuu8Pl8fr9/wYIFLVu21F8jetEjtJnL5QoMRnm9Xq/XW7t2bSHEihUrevbs6XA4hBB5eXl79uzZv3+/nRWtYX755ZeJEyc+//zzwTtXrFiRn5+vv87Lywv0FCGEWLx48bBhwxYvXvzcc8+tWbNG37lhw4a6devm5OQIIRITE7t3785NC3A4HK1bt/7uu++EECdOnNi1a5d+o1asWNGrVy99IoK8vLydO3cePHjQ5rriwtAQ1iBjx47t379/q1athBDFxcWBuFd8fHxycnJxcbGttatZHnjggXHjxklzqgXftIyMjKNHj5aXl9tRu5poz549r7322vz580+dOlVQUPDyyy8LIUpKSgJ3TAjRoEED/qYH+/DDDx955JHWrVs3a9ZszJgxPXr0EMaPWVJSUu3atfndjHYExmuKiRMnbtiwYdmyZfqmy+UKnpGO+cOCLVq0qLi4eMSIEdJ+l8ulR1iFEF6v1+l0Op3Oi125msrhcLRt21YPaHXs2HHEiBH33Xef+WMWy/OHSSoqKoYPH15QUDB69Og9e/b8x3/8R5cuXXJzc4M/ZkIIr9fL72a0oyGsESZPnvzee+8tXbo0NTVV35OZmVlUVKS/Pn78+JkzZxo1amRfBWuWWbNm7d27t1OnTkKI7du379279+DBg48++mhmZmagQ1NUVNSwYUMawoDGjRu3a9dOf52Tk3P06FH9Q1VcXOzz+fRB+KKioq5du9pazRrkhx9++Pbbb1esWOFyudq0aXP77bfPmDEjNzc3+GN29OjR0tJSfjejHUOj9nv++eenTZv2xRdfBA/09e/f/5NPPjl9+rQQ4v333+/UqRPJAAGPP/74vHnzXn311VdffbVly5Z33HHHXXfdJYTo37//vHnz/H6/EGLevHn9+/e3u6Y1SEFBQSA0+M0332RlZdWqVatTp05ut/vLL78UQhw8eHD16tV9+/a1tZo1SP369T0ez4EDB/TNXbt2paWlCSH69++/aNGiM2fOCCHef//93/3ud/p+RC8m3bbZTz/91LJly6ysrPr16+t7nn766RtvvFEIceutt+7atatDhw6LFi2aO3duz549ba1pDXXDDTcUFhaOHj1aCHHixInu3bs3aNAgJSVl9erVK1eujOUVcyRnz57t0aNHYmJi8+bNFy5cOG3aNP0fhZkzZ44fP75v374rVqwoKCh45pln7K5pDfLnP/95wYIF/fr1271799atW1etWtW0aVO/3z9gwICff/65ffv2ixYtmjdvnv4Li+hFQ2iz0tLSLVu2BO9p0aJFSkqKEMLn8y1fvvzQoUPdu3dv3LixTRWs6bZv356amhp4eKG0tPSLL74oKyvr1atXcnKyvXWraTwez9KlS0+fPt21a9fg0bzt27dv2LChefPmXbp0sbF6NdPmzZu3bt1ar1697t27B1YgCvxuXnfddZmZmfbWEBeOhhAAENOIEQIAYhoNIQAgptEQAgBiGg0hACCm0RACAGIaDSEAIKbREAIAYhoNIQAgptEQAgBiGg0hACCm0RACAGLa/wOTRMAf286kXQAAAABJRU5ErkJggg==",
"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.2"
},
"kernelspec": {
"name": "julia-1.8",
"display_name": "Julia 1.8.2",
"language": "julia"
}
},
"nbformat": 4
}