{
"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.373262e+01 1.473173e+01\n",
" * time: 0.004369020462036133\n",
" 1 6.181286e+01 9.424847e+00\n",
" * time: 0.012845039367675781\n",
" 2 5.730653e+01 1.332171e+01\n",
" * time: 0.1429281234741211\n",
" 3 4.111002e+01 9.094423e+00\n",
" * time: 0.16821718215942383\n",
" 4 3.278998e+01 8.708089e+00\n",
" * time: 0.19269013404846191\n",
" 5 3.207671e+01 8.130955e+00\n",
" * time: 0.21332502365112305\n",
" 6 1.794137e+01 5.222238e+00\n",
" * time: 0.2336280345916748\n",
" 7 1.392088e+01 3.460109e+00\n",
" * time: 0.2509341239929199\n",
" 8 1.081712e+01 4.927425e+00\n",
" * time: 0.26448702812194824\n",
" 9 8.871162e+00 3.611963e+00\n",
" * time: 0.281447172164917\n",
" 10 7.899794e+00 2.446670e+00\n",
" * time: 0.29865503311157227\n",
" 11 7.183731e+00 2.423015e+00\n",
" * time: 0.3133699893951416\n",
" 12 6.723927e+00 3.022203e+00\n",
" * time: 0.38468408584594727\n",
" 13 6.426257e+00 1.836654e+00\n",
" * time: 0.3999931812286377\n",
" 14 6.245604e+00 1.149805e+00\n",
" * time: 0.4139900207519531\n",
" 15 6.183646e+00 1.576459e+00\n",
" * time: 0.42772912979125977\n",
" 16 6.044556e+00 1.806774e+00\n",
" * time: 0.4418652057647705\n",
" 17 5.921497e+00 9.543571e-01\n",
" * time: 0.45578813552856445\n",
" 18 5.886266e+00 1.399724e+00\n",
" * time: 0.46642017364501953\n",
" 19 5.804765e+00 1.120524e+00\n",
" * time: 0.47691917419433594\n",
" 20 5.769294e+00 1.292946e+00\n",
" * time: 0.4880549907684326\n",
" 21 5.729341e+00 6.018522e-01\n",
" * time: 0.5027279853820801\n",
" 22 5.684643e+00 6.545033e-01\n",
" * time: 0.516801118850708\n",
" 23 5.647523e+00 8.595996e-01\n",
" * time: 0.5271501541137695\n",
" 24 5.615312e+00 5.179389e-01\n",
" * time: 0.5410439968109131\n",
" 25 5.600388e+00 5.167113e-01\n",
" * time: 0.5513961315155029\n",
" 26 5.593794e+00 3.383487e-01\n",
" * time: 0.56581711769104\n",
" 27 5.587878e+00 4.156380e-01\n",
" * time: 0.6231989860534668\n",
" 28 5.580488e+00 2.350075e-01\n",
" * time: 0.6379520893096924\n",
" 29 5.574371e+00 2.759970e-01\n",
" * time: 0.6486921310424805\n",
" 30 5.571000e+00 1.704074e-01\n",
" * time: 0.6628050804138184\n",
" 31 5.567104e+00 1.623976e-01\n",
" * time: 0.6734540462493896\n",
" 32 5.566883e+00 2.346724e-01\n",
" * time: 0.6840620040893555\n",
" 33 5.565906e+00 1.551490e-01\n",
" * time: 0.6978662014007568\n",
" 34 5.564864e+00 1.586927e-01\n",
" * time: 0.7082531452178955\n",
" 35 5.563443e+00 1.820743e-01\n",
" * time: 0.7188570499420166\n",
" 36 5.562844e+00 1.079697e-01\n",
" * time: 0.7294712066650391\n",
" 37 5.562135e+00 8.784016e-02\n",
" * time: 0.739915132522583\n",
" 38 5.561503e+00 6.406086e-02\n",
" * time: 0.7503809928894043\n",
" 39 5.561314e+00 3.956035e-02\n",
" * time: 0.7640790939331055\n",
" 40 5.561110e+00 6.085595e-02\n",
" * time: 0.7777981758117676\n",
" 41 5.560936e+00 2.773294e-02\n",
" * time: 0.7918150424957275\n",
" 42 5.560786e+00 2.922346e-02\n",
" * time: 0.8212339878082275\n",
" 43 5.560705e+00 3.828851e-02\n",
" * time: 0.8326292037963867\n",
" 44 5.560678e+00 3.178743e-02\n",
" * time: 0.8432400226593018\n",
" 45 5.560640e+00 2.659538e-02\n",
" * time: 0.8572499752044678\n",
" 46 5.560596e+00 2.331114e-02\n",
" * time: 0.8711750507354736\n",
" 47 5.560566e+00 1.981264e-02\n",
" * time: 0.885310173034668\n",
" 48 5.560538e+00 2.399216e-02\n",
" * time: 0.8954579830169678\n",
" 49 5.560527e+00 2.180800e-02\n",
" * time: 0.905879020690918\n",
" 50 5.560514e+00 1.809218e-02\n",
" * time: 0.9161181449890137\n",
" 51 5.560502e+00 1.169501e-02\n",
" * time: 0.9300141334533691\n",
" 52 5.560490e+00 7.590705e-03\n",
" * time: 0.9437880516052246\n",
" 53 5.560485e+00 1.731625e-02\n",
" * time: 0.9543640613555908\n",
" 54 5.560480e+00 8.731067e-03\n",
" * time: 0.9646711349487305\n",
" 55 5.560474e+00 6.518268e-03\n",
" * time: 0.9750311374664307\n",
" 56 5.560471e+00 5.395783e-03\n",
" * time: 0.9884591102600098\n",
" 57 5.560469e+00 4.344324e-03\n",
" * time: 1.0215489864349365\n",
" 58 5.560467e+00 2.985214e-03\n",
" * time: 1.0368411540985107\n",
" 59 5.560467e+00 6.343350e-03\n",
" * time: 1.0477252006530762\n",
" 60 5.560465e+00 3.726513e-03\n",
" * time: 1.05885910987854\n",
" 61 5.560464e+00 2.442238e-03\n",
" * time: 1.0702450275421143\n",
" 62 5.560464e+00 1.714623e-03\n",
" * time: 1.0852010250091553\n",
" 63 5.560464e+00 1.888507e-03\n",
" * time: 1.0998461246490479\n",
" 64 5.560464e+00 4.058083e-03\n",
" * time: 1.11051607131958\n",
" 65 5.560463e+00 3.619347e-03\n",
" * time: 1.1248271465301514\n",
" 66 5.560463e+00 2.193256e-03\n",
" * time: 1.1353859901428223\n",
" 67 5.560463e+00 1.921064e-03\n",
" * time: 1.1458940505981445\n",
" 68 5.560463e+00 1.479551e-03\n",
" * time: 1.1598000526428223\n",
" 69 5.560463e+00 1.212007e-03\n",
" * time: 1.173539161682129\n",
" 70 5.560463e+00 1.008847e-03\n",
" * time: 1.187148094177246\n",
" 71 5.560463e+00 6.111069e-04\n",
" * time: 1.2009592056274414\n",
" 72 5.560463e+00 7.708128e-04\n",
" * time: 1.2297241687774658\n",
" 73 5.560463e+00 5.619155e-04\n",
" * time: 1.2411301136016846\n",
" 74 5.560463e+00 4.407222e-04\n",
" * time: 1.255547046661377\n",
" 75 5.560463e+00 4.871960e-04\n",
" * time: 1.2694401741027832\n",
" 76 5.560463e+00 3.619096e-04\n",
" * time: 1.2802400588989258\n",
" 77 5.560463e+00 4.568508e-04\n",
" * time: 1.2907171249389648\n",
" 78 5.560463e+00 3.385982e-04\n",
" * time: 1.3043551445007324\n",
" 79 5.560463e+00 3.215527e-04\n",
" * time: 1.3148210048675537\n",
" 80 5.560463e+00 4.307813e-04\n",
" * time: 1.325383186340332\n",
" 81 5.560463e+00 2.834495e-04\n",
" * time: 1.3393089771270752\n",
" 82 5.560463e+00 2.935447e-04\n",
" * time: 1.3495080471038818\n",
" 83 5.560463e+00 2.665229e-04\n",
" * time: 1.3598260879516602\n",
" 84 5.560463e+00 2.144549e-04\n",
" * time: 1.3698961734771729\n",
" 85 5.560463e+00 2.273653e-04\n",
" * time: 1.3802452087402344\n",
" 86 5.560463e+00 2.506765e-04\n",
" * time: 1.3904991149902344\n",
" 87 5.560463e+00 1.336142e-04\n",
" * time: 1.403954029083252\n",
" 88 5.560463e+00 1.801719e-04\n",
" * time: 1.4324381351470947\n",
" 89 5.560463e+00 1.550423e-04\n",
" * time: 1.4435451030731201\n",
" 90 5.560463e+00 1.283101e-04\n",
" * time: 1.4543280601501465\n",
" 91 5.560463e+00 1.988984e-04\n",
" * time: 1.4646780490875244\n",
" 92 5.560463e+00 1.219513e-04\n",
" * time: 1.4751310348510742\n",
" 93 5.560463e+00 1.133006e-04\n",
" * time: 1.4860990047454834\n",
" 94 5.560463e+00 1.775406e-04\n",
" * time: 1.4969711303710938\n",
" 95 5.560463e+00 1.140080e-04\n",
" * time: 1.5076391696929932\n",
" 96 5.560463e+00 3.138961e-04\n",
" * time: 1.5181941986083984\n",
" 97 5.560463e+00 2.873798e-04\n",
" * time: 1.5287971496582031\n",
" 98 5.560463e+00 2.310047e-04\n",
" * time: 1.5395431518554688\n",
" 99 5.560463e+00 1.805115e-04\n",
" * time: 1.5537981986999512\n",
" 100 5.560463e+00 1.717529e-04\n",
" * time: 1.5680911540985107\n",
" 101 5.560463e+00 1.470132e-04\n",
" * time: 1.581780195236206\n",
" 102 5.560463e+00 9.742409e-05\n",
" * time: 1.5956029891967773\n",
" 103 5.560463e+00 7.471977e-05\n",
" * time: 1.6096620559692383\n",
" 104 5.560463e+00 7.286124e-05\n",
" * time: 1.646496057510376\n",
" 105 5.560463e+00 7.286123e-05\n",
" * time: 1.6977920532226562\n",
" 106 5.560463e+00 7.286123e-05\n",
" * time: 1.7524051666259766\n",
"e(1,1) / (2π)= 1.7391794030083887\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "Plot{Plots.GRBackend() n=1}",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXwURd748eo5cnEk5AIJl9yHGhQEFNQIJCKHRDk0KyuXPsK6B8/uw7q6iqw36mtVPFZXORQRERV8IeijsByCEkEEOeSQSwgJh0g4k8xk5vdH+8xvujrUZIYOnWE+71f+mO6u6a50JqlUfftbpfn9fgEAQKxy2F0BAADsREMIAIhpNIQAgJhGQwgAiGk0hACAmEZDCACIaTSEAICYRkMIAIhpNIQAgJhGQwgAiGmuC3MZn8/30ORH7//73y/M5aKd1+t1uS7Qj+aiUVlZ6XQ67a5FNKmsrHQ4HJqm2V2RaBLjv5vJ8e6QZbZvP3b2rDes0+7d++OhQ5vuvffeSOt1vi7QT7S8vLzEoy3de/TCXA4AYLnb2l0Sssztty/cuPFImCfenJ9famNDyNAoACCmxW4fHwBgOU0TUTfcTo8QABDT6BECACyjaVq4T2D5/TZ3IWkIAQCWYmgUAIAoQo8QAGCdKHxYhoYQAGCZCGKEwu4YIUOjAICYRkMIAIhpDI0CACwTQUK93SOjNIQAAOtEECO0/eEahkYBADGNHiEAwDpa9CXU0xACACyjhT/UafvQKA0hAMAykeQR2t0SEiMEAMQ0eoQAAOtEECMMVf7nn39+5513jh8/PmjQoKuuuspcwOfzvf/++1u2bMnOzh4yZIjeJf34448PHToUKJOZmZmfn1/l+ekRAgCspIX/pXDy5Mlu3bp98803DoejT58+X3zxhbnMvffe++yzz9arV++xxx6bMGGCvnP79u3f/p8nn3zyo48+Otcl6BECAGqv2bNnX3LJJbNnzxZCpKamPvHEE7m5ucEF9u/fP2vWrH379jVs2PCOO+5o167d3//+98zMzL/+9a96gbKyssaNG48ZM+Zcl6BHCACwjKb9+rxMGJR9wqVLl958883665tvvvnLL7/0eDzBBZYvX37FFVc0bNhQCNGsWbPWrVuvWrUquMBHH32UkpJyww03nOsSNIQAAOtEMDCqHBstLi7OzMzUXzds2NDn85WUlAQXKCkp0VtBXaNGjQ4ePBhcYNq0aWPHjlU8y8rQKADAMnqHMLy3aNratWuHDx8evLNfv376YKbT6fT5fPpO/YXT6Qwu6XA4AgWEEJWVlcEF9u7du3LlypkzZyoqQEMIALBZ48aNhw0bFrynffv2gUPFxcX664MHD7pcrkAHMVAguAtYXFzcuHHjwOb06dPz8vKaNm2quDoNIQDAMhGsPiGEyMrKkhrCgP79+z///PMPP/yw0+mcP39+Xl6ey+USQmzdujUtLa1hw4Z9+vQZM2bM7t27W7Zs+cMPP+zfvz8nJ0d/r8/ne+utt/75z3+qr06MEABgHatjhMOGDXM4HHl5effee+8zzzwzadIkff/o0aPnzJkjhMjMzPzzn/+cm5s7YcKEm2+++YEHHkhOTtbLfP7556dPnx44cKC6yvQIAQC1V0JCwqpVqxYtWnTixImHHnooMMj5yiuvBJ6ReeKJJ/r3779ly5aCgoLu3bsH3tu6deulS5fGx8erL0FDCACwjCZCpENU9ZYQEhIShgwZIu3s2rVr8GbPnj179uwplWndunV1KkBDCACwTvgxQrvn3KYhBABYKArXI+RhGQBATKNHCACwTAQJ9baPjdIQAgAsE0Eeod3tIEOjAIDYRo8QAGCdKHxYhoYQAGCZXxdWCvMtNVSZaqIhBABYJgo7hMQIAQCxjR4hAMA6EXQJ7e5C0hACAKwTQYywhmpSbQyNAgBiGj1CAIBltPAT5O1+aJSGEABgHU0LPx3C7paQhhAAYJ0ofFiGGCEAIKbRIwQAWCYKO4Q0hAAA60QyxZrdTSFDowCAmEaPEABgnSgcG6UhBABYJoI8QtsbQoZGAQAxjR4hAMA6UTjXKA0hAMA6xAgBALFM06JvrlFihACAmEaPEABgGU1o4SbIR01C/enTpx966KFevXrl5OQ8++yzlZWV+v6ioqK77rqra9euY8eOPXz4cI3VEwAQDbTwv+xW3Ybwb3/728qVK1955ZUpU6ZMnz79xRdf1Pffeuut9evXnzlzphCioKCghmoJAEANqe7Q6Lp16+65557s7GwhxPDhw9etWyeEKCws3Llz5+rVq91u99SpUzMyMrZu3dqxY8carC8AoBaL4GEZ2zuF1e0R5ufnz5kzZ+/evT/88MPHH3+cn58vhPj++++vvPJKt9sthKhTp06nTp02btxYg5UFANRu2q/TbofH3jpXtyG89957y8rKunTp0r1799atWw8ePFgIcfjw4ZSUlECZBg0anCtM6Pf7A2FFAEA08ng8oQtdxDHCoUOH9ujR48iRI0ePHnW5XL///e+FEPXr1z9z5kygzKlTp5KTk6t8u6ZpTqfz/KsLALCLPv538alWQ+j3+1esWPGb3/zG4XDExcXdcccdy5cvF0K0aNFi586dehmfz7d79+4WLVrUWFUBALVeJGOjNle5Wg2hpmnt2rVbvHixEMLn8y1evLhDhw5CiLy8vJMnTy5atEgIMW/evPj4+Ouuu65GqwsAqM2icGS02k+Nvvnmm7/97W9nzpxZXl6emZn53nvvCSHi4+NnzJgxevTolJSUU6dOvfPOO4x/AgCiS3Ubwh49euzcufPIkSNutzv4AZkBAwYUFRWVlJRccsklLhfz1ABAbKuBSbcPHTo0ffr0EydO3HLLLddcc425gNfrffvtt7dv396pU6cRI0Y4HI7A/vfee2/Tpk0pKSnDhg1r3bp1lecPb67RjIyM4FZQ53a7mzZtSisIAIgoeULVEp44caJbt267du1KS0sbOHCgHqSTjB079o033mjatOnUqVPvu+8+fWdFRUXfvn1fe+211NTU0tLSr7766lyXoPUCAFgmghXq1eXffvvtSy+99M033xRC1KtX76mnnurfv39wgb17977//vv79+9PT0+/7bbbWrVq9cgjjzRq1OjVV189c+bMV199FbKfxuoTAIDaa8WKFXl5efrrvLy81atXV1RUBBf48ssvs7Oz09PThRCNGzdu27bt6tWrhRCLFi266667Fi9ePHXq1O+++05xCRpCAIB1InhsVNkjPHjwYGZmpv66YcOGfr+/pKQkuEBxcXFGRkZgs2HDhgcPHhRC7Nmz5+WXX54/f/7hw4dvuukmfU7sKjE0CgCwTCSJgZq2du3a4cOHB+/Lzc295557hBBut9vr9eo79RdSXr/L5Qqeuczj8egFNE3Lzs6eMWOGECI7O/svf/nLqFGjqrw+DSEAwGaNGzceNmxY8J7LLrsscEjv4QkhioqKXC5XoINoLiCEOHjwYFZWlhAiKyurU6dOgbMVFRV5vd4q44U0hAAAy0Sw+oQmRFZWltQQBgwcOHDKlCmTJk1yuVwffvhh//799YT1DRs2ZGRkZGVl5ebmjhkzZseOHW3btv3++++Li4tzcnKEEIMHD166dKl+ksLCwjZt2pzrqRkaQgBA7TVkyJBXXnklJyenVatWixYt+t///V99/7333ltQUDBhwoS0tLQHHnggNze3X79+ixcvnjRpUr169YQQY8eOfeutt/r169esWbMFCxa89dZb57oEDSEAwDrhxwjVxePj45ctW7Z06dJffvllypQpjRo10ve/+eab+pOiQoiHH364f//+W7duHT9+fOfOnfWd9evXLywsXLJkyZkzZyZPnty4ceNzXYKGEABQq8XFxd18883Szssvvzx4s0uXLl26dJHKxMfHDxgwIOT5aQgBAJaJxhXqaQgBAJayu2ELFwn1AICYRo8QAGCZCBLqbV+Zl4YQAGCZSPII7R5KpSEEAFhGEyGWVar6TbYiRggAiGn0CAEA1olghXq70RACAKxDjBCo7fyKrVpN/lth998O4KJBQwgAsIwWfjoEPUIAwEWEGCFgg/BGO5XHz+vN4Qn9t8JYQr60X3UCxlFhl2jMIyR9AgAQ0+gRAgAsE8EUa7Z3CekRAgBiGj1CRAF/iNCcKkgY1ltDHw9Vlerzm/4Llrf9YRyWDsoBRGVAUdj/HzlgJxpCAIBlInlYpmZqUn00hAAAy2gi+pZhIkYIAIhp9AhRC4SO6oUTBZSOKg/LR03X9ocTcgwR1JNo8rnkxWukLXVIUXkxTaqaqbDfGEQkDRGRiyCh3u4PGA0hAMA64ccIbUdDCACwTAR5hMQIAQCwEz1C2CBEMp46rGfa9vtVYb/wjvqUFTOVD2/uUTnmZ/4v2JgaKJV3GMN4mlRYOipdOlS8UgqdSuWlCGO0DXwBajSEAADLRDQ0WkN1qS6GRgEAMY0eIWpGiBHEMHMY1MObPsO2z6cqLB9Vbgoh5B3KmqiFHK50qMdCjZsO5UipfNT4H28V/4DLu1SDtORaQCGCmWVs/8zQEAIALGV3wxYuGkIAgGU0Lex0CNInAACwEz1CWCTEtGd+xVF1zM+8x6ferAznqLKweY8cUzyPGKEU8xOmwJ686VRuqgsrzyzMQURTAcNRZa6FPHNctI2S4XxFsPqE3R8SGkIAgHUimGvUbgyNAgBiGj1CAIBlNKHJ66hU5022oiFEpMIKCoYz7ZkUtxOmQF1lpWEmtEqv8ag3nKPGU5kvXVmpjBH6VLFPOW9QmRcohHBKkTyXYdPpNIzfOF3SphbGUac8FOTwqzblkSMpKBhWyFDY/kcPNSuSFert/kjQEAIALBVlHUIaQgBA7XbgwIE333zz+PHj+fn5OTk55gIVFRXTpk3bunXrFVdcMXr0aJfLJYQoLCxcvnx5oMy4ceOSk5OrPD8PywAALKNFQHnC48ePd+/e/ejRo61btx42bNjHH39sLjNy5Mj33nvviiuumDFjxr333qvvXL58+Zw5c375Pz7fOReXoUeI6lHPFirCCwqGldsnTIE9rzHs5/UYj3oqDZsVhqMeubAUUJR/VdRBRGnZJr80RaeQgoKGwlKqnwgVBXS5DZtu46YrTlXY5XYa6umS77B0Lb+UhiikNETp3WGEDAWJhhc7y+cafeutt9q3b//yyy8LIRITE59++unBgwcHF9i1a9eCBQuKiopSU1MHDRrUokWLRx99NCsrSwjRrVu3p59+OuT16RECAGqvL7/8sm/fvvrrPn36FBYWVlRUBBdYvXp1586dU1NThRCNGjVq3779119/rR/atm3bY489Nn369OPHjysuQUMIALCOFtHXuRUXF2dkZOivMzMz/X5/cXHxuQroZfQCDRs27Nixo8/ne//999u3b79v375zXYKhUVRLFVOJWTcWKmU4eE3jk9IApqfcMPjpqTBsVkhHy6WjxqFR43ulgVMRKhnDlD6hHhqVMhxMQ6Mu9eCnYXgzLt5YOF46ath0xxnvv/FUQghprNTlD+Nf5LBGSoU5v6L6V0JUiGBhXqGtXr06Nzc3eOctt9zyhz/8QQjhdru9Xq++0+PxCCHi4+ODS8bFxQUK6GXi4uKEEKNGjRo1apS+c+jQoVOmTHn11VerrAANIQDAMhHMsKYJ0aZNm/vvvz94Z8uWLfUXWVlZRUVF+uuioiK32x3c/xNCNG7c+MCBA4HNAwcO6AHCYN27d1+2bNm5KsDQKADAZpmZmX2NAg3h4MGD58+fr8cF586dO3DgQKfTKYQoLCzURzvz8vJ27dq1detWIcT69esPHz6sp1icPn1aP4PH41m0aNHll19+rqvTIwQAWKYaCRFVvEdx8NZbb/3Xv/7Vs2fPFi1arFixYsmSJfr+P/7xjwUFBRMmTGjQoMEjjzySm5vbu3fvJUuWPP7443Xr1hVCXH/99Q0aNEhPT1+3bl1qaurf/va3c12ChhBVMy0wFGL7vIKCHinmJwfqpLBfRZlhs7zMa9g8K21K7zUclQOKpkur8yvkZZuUiQEOOUYoj8fICRLGjAgp7BeXYPjljU80HPXKR0OsNuX3O6Vt43E5pqgghwzN6RNSfoXfGDIkZngRsDR9wu12L1myZOXKlaWlpa+//rr+dKgQYtasWQ0aNNBf//Wvfx04cOCWLVseeuihdu3a6Ts///zzb7/99uTJk3/5y1+6du2qaJ5pCAEAtZrT6bzxxhulnW3btg3e7NixY8eOHYP3pKWl5eXlVef8NIQAAMtEMul2zdSk+mgIAQCW0TQRdvqE3S0hDSH+TxWpgkEHzVOsSdExY3AtrKCglNsnhfGEKexXJm2eNmyePWM8esZjPJUqviglHQpT1LBGY4RSUFBKDYyXw36GmieUu4M3vUnq1abkX3w5IdKv/sugChnKa9GF+hsn/xFkArZoF1n+hK1InwAAxDR6hAAAy0Rhh5CGEABgIavzCC8AGkL8Sr2sUhXLMPnCSByUpg+VAm9SUFAKAQohzp42xPnOnjIUOHOqwnjUWPi0FDJUJR1WVBEjNOyRgp0+Y+xNumfSb7dDXmVJ/uV3u6VMQWOMMNHw25qQZNiUoptejyFkWBlqoStlgFhomkvaNm45jJuGkzlMf+I0eZdq2Sbb+wqIBWHECDds2JCfn9++ffsbb7xx9erV+s69e/cOGTKkQ4cOd9xxx8GDB2umkgCA6GHd0hMXRnUbwm3btuXm5ubk5CxYsGDSpElJSUn6/vz8/FatWi1evDg9Pf3222+vsXoCAKKAnkcY7pe9qjs0+uijj955550TJkwQQrRv317f+dVXX+3fv//JJ590uVzPPvtsenr6pk2bFBObolZRT6KmnkFNhBoaldeUN05UJo1ASuOT0kCoEOLMScOe0ycMY6GnjUfPnDQcPWMcKZVyLaRhWGnyNmHK9DANjUpZB8ZlmIy/39KS9OahUZdyaDTBODRaXmYcGq0wpk+EU08zqWbSYKb8fWnqwvLJ/fLZlQvcMwFbtIlgrtGwY4pWq26PcO3ate3atRs9evTNN9/84osvVlZWCiE2b97cuXNnl8slhEhMTOzUqdOmTZtqsLIAAFituj3C/fv3P/PMMy+//HJycvK4ceN++eWXyZMnHzlyJDk5OVCmQYMGhw8frvLtfr9fbzsBAFHK4/G43e4QhWpH2C8s1e0RJicn33333QMGDOjVq9c//vGP2bNnCyFSUlICCz4JIU6ePBmYC1yiaZq+ghQAIEqFbgUjihHa3nBWt0fYqlWrQOcvJSXl1KlTQohLL710+/bt+k6fz/fjjz8GllJE9JGjgsYtU1BJDgoaEwmkGJVHuY6SFKiTsiOEKSh4qlS1eVqKERojiOoZ18wxQo8yRmiamUwVI9TkKdbCS58wLSBlDArKy0Wp6mmmDmdKiR/SXHFSYc1pqImUXCHMQcQQqRvKo4AVqtsjHDly5LvvvnvmzJnKysrp06f37t1bCNG3b9/y8vIPP/xQCDFr1qy6dev27NmzBisLAKjdovCh0Wr3CMeOHbt27drmzZu73e7s7OyZM2cKIeLi4t5555277rrrj3/8o9vtfvfddx3yopwAgFgShXOsVbchdLlcb7755ksvvaRpWkJCQmB/nz599u/fX1pampKSYvsjsAAAhCu8KdYSExPNOx0Ox7mekUHtch6Jg+YgkylxULXQkjRRmbys0mnVlGnClCloihGWGwqfUOURSjOumQJv5jxCaQEjVfZkiCnWlKE1IYTLbRhQiSs3xgjLDb+tpooZJ3tTBgVNk5yZchzlTdXkcNKmHF90yDVxSDctnLRCuTT/e9dG4ecR2v2DZK5RAIBlIgn62f0PDQ0hAMA6URgj5NkWAEBMo0cYQ0IsvaNMHPSZY4SVUowwnMlFjZG5s2ekPELzXKOqTEF56tETxqlHT6nyCKXEQWl9KGFaQMoUijMUlu6RacZOw6aUnCeEcLkMe7wVqltqilYKhZApjC7jHmnWU1ecQ7XpliKIhqo4TaFQnzFqGFZaISHC2i+SuUbt/knSEAIALBOFI6MMjQIAYhs9wthlypeQNlWJASJk+oRxWM8jr7tknGLtjCrDQZiGN6VZ0+Sjcr6ENKeaKl/CaxoalaaOM61nZCjsN95UacDHtGC96VrGwdLKSsOmOlVD4lAOw0qDmUIId5xhLNQdb7hpcfHGud/ijTfN+F6X2zhm6zYNyYf4pKlGP23vOiC0KOwS0hACACyjRbC+IA0hAOCiEUEeod3tIDFCAEBso0d4UQuRMKE6ep7pEx45fcK4aUyfMIUM5RihNAebOqZoWmjJq9iU8iWk70IIUSkFBUME6lQJFNK/yVqlKa/AqZy/TT1rmhyAlNIhjJO3xcuLg8YnGG5LfKLhL0NFohRYNRyNSzDcNLeccCL/ty19Xw55jjV1EgozrtV6xAgBADEtgjxCuxdsYGgUABDT6BECACzDpNuoXdQhQfWcanIwLGSMsFI1H5jXuAyTesY1KYwnhCiT43zGs0lLKZWpruVRzltmjhGqs/f86mw+Ka1QWmDIdEv9PkMB9bk1h+H7ktZ4crmlO2YICpYlyHc4QbrD8upUxhWgjD9N+R4qMy9FqI+WXz3jmvEvpjqcCFQTDSEAwDJaFK5HSIwQABDT6BECACwTSUK93YPaNISxJIx4VoiQoTAFz3zGuUbVaYUe5UykUpDPvCe8KKAyU7BSTnozfZvqOVdDxAiNZaWcN/Mvv+niivLeCsO2xyklbirvWKg77FHeUo8ysCp9EswZqCFyMcP5lNo9ooYqRLIMk90tIQ0hAKBW27Nnz2uvvVZaWpqfn9+vXz9zgbKysldeeWX79u2dOnUaP358XFxc8NG5c+cePXr0vvvuO9f5iRECAKyjhf+ldOzYsR49evh8vquvvnrkyJEffPCBucyIESM+++yzG264Yf78+XfffXfwoa+//vr3v//9k08+qbgEPcLYZRqRUo1Jmaf3kpdhkldlUo2eedUjpR554M5r3CMVkFM1vKqn+Svlhd2Nm+Y50+R8ijDWQpJo8ntDLN1uHBmVa+6oVH3X8h2W75h8h6UCYSWZyHfYF2JoVDlVnPw59EvT1KneitrB6jzCmTNndu7c+dlnnxVCOJ3OZ555ZujQocEFduzYsWjRopKSkuTk5Ly8vKZNmz7xxBNNmzYVQpSXl993332TJk16+umnFZegRwgAsIz26+MyYX2prF69unfv3vrr3r17r1u3rry8PLjA119/feWVVyYnJwshMjIyOnTosGbNGv3Q5MmThw4d2qFDB3WdaQgBALVXcXFxenq6/jojI8Pv95eUlAQXKCkpCRQQQmRmZhYXFwshNmzY8Pnnn0+cODHkJRgaBQBYJ6LVJ1avXp2bmxu8b8CAARMmTBBCxMfHV1RU6Dv1FwkJCcEl4+PjPZ7/v+ZMeXl5QkKC1+u9++67X375ZbfbHfL6NIQxxBT+8isOy5OHVRE8U02UJUcQlckV6s2q3i7N7ma8tHKKL/Wz+1XF/CIPCsonMr5Xk2cPk6OG6gwWXzjftfoGVlVA/QNS1STk5HxCzkhRrqwkFdYIGdZ2Wvh5gZoQbdq0uf/++4N3tmnTRn/RpEmTAwcO6K/3798fFxeXkZERXDIrK2v//v2BzQMHDmRlZe3YsWPTpk0jRowQQpw9e/bw4cOtWrVavny5HjuU0BACAGyWmZnZt2/fKg/l5+c/9NBDDz/8cEJCwrvvvjt48GCHwyGEWLlyZVZWVqtWrfLy8saMGbNx48bs7OzCwsKff/75xhtvdLlcP/zwg36G1atX/8///M8XX3zRqFGjKi9BQwgAsIymhZ8gryw/ePDgN954o1u3bi1atFi7du3SpUv1/RMnTiwoKJgwYUJycvLjjz9+0003XX/99StWrJgyZUpSUpIQomXLlnrJ3bt3u1yuwKYZDSEAwDpWr1DvcrkWL178zTffHDt27J133qlfv76+f+7cufqTokKIP/3pTwMGDNi2bduUKVMuvfRS6Qw9e/b86quvVJcIs76IGaGCZ6GWbVIdlbP3lNEv854Q6/hIl1ZPHRcyFHoeQUE185lNiYZSeXUE0bAprx4V+g4bQ4zKH1BYP+uQn5xQq4UhykQy12ioAg6Ho0ePHtLOFi1aBG+2bt26devWVb49MTGxefPmqvOHqgAAABczeoQAAMtE46Tb9AgBADGNHiF+dZ6RGvXMpfJEpsooUVVRJVXmmTIf0jyHavSGpFS3VL08lPoGmsuH+AGFU5NwRe+PB9GLhhAAYCVrJ92+AGgIAQCWiSRGWENVqTYaQvzqPD+LmrypKTbNpQ1bpqpIv1emzTBqYjoeRUNxqlsqf5PyLVXdQHP5ED+gcGoSLtv/JuJ8WZ1HeAHwsAwAIKbRIwQAWCYa0ydoCAEAlolgZhnb0RDiHELH7YybDmOgyKE66lBvOuURe4dTOrnq7fKlld+I9K+o3/xtWrcMk3xm8x+LsEKh8v03bJruiXSH5WtL91z9AwrrZx3yk2N7fAggRggAiGn0CAEAlonGGCE9QgBATKNHGEPkBDo5ymQIf8n/o5kzz9RRKOOm0yVtOpSb8rVMb1dFvORwl9OwSpBWaay2T/qupSvLix9pWuQhQ/MtVBdQR2FDfdfq+x/yDofxA1JHcKVNIUKEQtUZkcQTo4Kl6/JeCDSEAADrRGFCPQ0hAMAyEcQIbUeMEAAQ0+gRxq5Qc3Iawl/mQI8cBZSDgqogk8tt2HTHGTfdTulaLuMeqYB0NpfxWpXG4JnPaQzy+YzfmCnm53NIQUEpZBV5kND8T7PDuEu+w1LYz6n6ruU7LN8x+Q6744x3OE51NvmH65Q2Vd+FqOqzFMzamUtx4UWQUG97B5KGEABgnSiMETI0CgCIafQIY4l6AaJwnt0X5lm4lM/fu9VjofGGcbm4BHngTtojb5YbNr0VhnyJykrj0KhPnf/gE8odppXf5eIK6inTRBU5J8rh5TjVLY1T3tIq7rCxvPQTMY1dqyom53WYPjnqGdpCdA7s7jogJK2KOEvItzDpNgDgYqFFkEdYMzWpPhpCAIB1IosR2rpCNjFCAEBMo0d4MVPHBMMKClYRI1Q/zS+lNBifzpciUr8/iAYAAB2WSURBVPEJhs9hfKL8sUww7qko8xo2yw1HvR5DWE8OCvpU/3mah3QqK1Vv98vrNqniriHDrnKChOt8ooDSLTUcTTDdYemeSz8R6VrST1NOrlBO9iaqWKdJGToNkWuB2ieChHotvEQky9EQAgAsE0keYc3UpPoYGgUAxDR6hAAAy0S0HmEN1aW6aAgvaqaFl6RtRWlpvaEqJspSB7TcUkDLuClHsAwxv4Qk+WNZXmaMEZa7gzfloGBlGCslaY5Kw6kqTN9mZRgn9xtvqTxbmPHc5uCZHGdVBgWlqF5CkuGeJNZRHU2oY4rCJqliinFyyDDytEJh+iyFXJ1KddTuP6C4ONAQAgAsE0GM0PZ/aIgRAgBiGj1CAIBlNBF+jNDuLiENYexSh2ZM2V1yqM00MaZhU54JU45vGdPajDE/T7k84aenwhjJk4OChk2/abrQYHKgzvhdeJzym73eMK6lnL1VnlTT4ZTHY1zqxEFjpqAU9kuqq9yspzoqhEisYwwiGk8u/bzcclqhFCMMOdeocVOZR2j3mBnCd9GvPuH1en/zm9+MGzcusGfHjh39+vW79NJLBw0atG/fPqurBwBAzQqvIZwyZcqmTZsKCwv1Tb/ff+utt3bv3n3NmjXt27cfPnx4DdQQABA19Idlwv2yVxhDo9u2bZs3b96f//znqVOn6ntWrVp1+PDhSZMmOZ3Oxx57LCMjY8OGDZ07d66ZquJ8hUimkAorhxBFqPQJl5w+YVwpyfjov7RwktdjGAgVQni9ccGb6mnPJKYBSeMQrluarc10aWVuhnpRJ/Wcak5TXoH6pkk5DFKChDTaWae+4Y7VqWfYTDJuCiES6xpPLmdTKGdcM653Hzp9Iqw51RgpjToRTLFm9w+2ug2hz+e75557XnnlleDxzy1btnTu3NnpdAohEhISOnbsqO+pkZoCAGq/GogR7tix4+WXXz5+/Pitt9566623mgucPn36hRde2LJlS3Z29p/+9KeEhAQhxDfffPPhhx8eOHAgNTV12LBh119//bnOX92h0X/+85+dO3fu2bNn8M6ff/65Xr16gc2UlJQjR45U+Xa/3+/1eqs8BACICrb8GT9y5EjPnj1TUlL69et33333zZkzx1ymoKDg66+/HjJkyNKlS0eNGqXv3LRpU1pa2i233NKkSZOBAwd+8skn57pEtXqEP/300wsvvLBo0aLdu3cfPny4vLx89+7dl156aYMGDU6dOhUoduLEidTU1CrPoGmay8UTqgAQxarzZ9zyDuGMGTOuvvrqRx99VAhRWVn53HPPFRQUBBfYtm3bkiVLDh06VK9evRtvvLFx48b79u1r3rz52LFjA2X279//ySefDBw4sMpLVKtxOnjwYHx8/G233SaEOHXq1PHjx3Nzczdv3tyyZctt27b5/X5N0yorK3fu3NmqVavqnBD2UAYJ5XVQpECOaewgRPqEMW7kjjOcPD7BsFlpTFGQQoDCPLFZiKCg8fF95fpQ5WcN9awok2OEHmPAstKrXpXJOMWafA9Vd0wI4XarVlaSc06SVAkSUlCwbrIUMpTTJ5LqGgpIc7BJqzK51cswhU6fUC7DJJUmSBhtNNPPNPRblMXXrFmTk5Ojv87JyRk5cmRZWZk++BkocOWVV+rDk6mpqZ06dSosLGzevHmgwP79+1etWvWXv/zlXJeo1tBojx49dv2f559/vmPHjrt27UpMTOzbt6/P55s9e7YQ4o033khNTb322murc0IAAKqjpKQkLS1Nf52enu73+0tKSs5VQC9TXFysv54zZ46mac2aNbv88stHjBhxrkuEPcVa3bp1s7Ky9Ncul+u9996bNGlSSkrKCy+88O6774b/sBAA4CKiRfK1fPnyrkZTpkzRz5eQkFBRUaG/Li8vF0IkJSUFXzAxMdHj8QQ2y8rKAgUKCgr8fv/OnTt//PHHRx555FxVDjtud8stt9xyyy2BzV69eu3evVvqqAIAYlNkeYGdO3d+7rnngvcEelxNmjT56aef9Nc//fRTQkJCenq6VDI4nWH//v1NmjQJLtC6devx48e/+OKLeqDRzJoHWGgFLwbhpBUKc4zQGIrzu4z5dsbMMyn9zlfpUhwVptibXDdlfp46Oa8swZhHaIoRepUxQtOqTKoYoSnz0pxHqIoRJkjrLinzCKVMQSkoKGUZClMeoXQtqSbyFGsu1SehihhhOImDiDoRrUeopaSkdOnSpcqjQ4YMmThx4sMPP5yUlDRr1qz8/HyHwyGE+OKLL5o3b962bdu8vLwxY8Z8++23Xbp0WbVqVWlpqR5TPHz4cGZmphDC5/MtXbq0bdu256oAT3ICAGqvgQMHzpgx46qrrmrSpMn27duXLl2q73/ooYcKCgratm1bv379Z555pl+/ft27dy8sLPznP/+ZmJgohOjfv7/P58vIyNixY0f9+vUXLlx4rkvQEAIArGN1/oTT6VywYMGGDRtKS0u7deumN3JCiAULFtSpU0d/PW7cuIEDB27btq1Tp06XXHKJvnPNmjWbN2/++eefGzdu3L59e0U/lYYQAGAhLdxllapT3jxnWaDB0zVp0kQKDbpcrmrOdEZDGLukf4/84aQVCvMcnsbi8tSjxqN+n9O4aYy0VV3f4LqoYm+mxYyM2XjGoGDC2RAxQnUeobwqk/KeqWdnFaHzCI1hvCTV1KOJ8qpMccaj8i++tAxTqMlFpcRBVVBQqyJGGEbiIM+h4wKgIQQAWEaLwn9faAgBANaJwoV5aQjxK9P8a/IOubw0k5lU2jjQ6vI7pMPGrRCfQ6kqpswN5exu8YZLS4N+5WWGoVGPaRkmj3GJKGk2ODnTQ7lEvWlSOtPQqHGpdylLQZrnTD3jmrSOkpRrIWVHCPNYqDynmnGgW7nQkrRpHho15UuoR0oRZSJLn6ihylRT2DPLAABwMaFHCACwTAQzy9g+DEBDCACwDjFCRDFlTDD0v3jKVAETp+KYpskfS/XKSnKqRpxqTrWKRGO+RLnhWlJEUAjh9Vy4GKGUluBWfiNxypChKaComjLNfHI5KGgMu8oLLUlBQWndJHOIUP3RsPtvImIQDSEAwDIRPCxj+38/NIQAAMtE4cgoT40CAGIbPUJUTT0Bmwg1B5sjvH+xjCGrKqZzk06uWs/ItO6SITVQCgp6KwxHPR45RigFBaUp1kyTwxmXYZJmC1NWW5iihm63FPs0xgjjVUmHppifah0lYbppUqagOigoz6kWOi+QSdQudtH2M6UhBABYJwoT6mkIAQCWiSSP0O4eJDFCAEBMo0eIajH/x+aXQ0GqHLqwQoaaJpeWpqsMK0bolUJrCYaYnzpN0LxHShz0G4uHihEaCjtMk3A65ciclMynjiCqCkuJgK5Q4Un19KEhgoLKqURF9MWPEKYIHhu1Gw0hAMAymhZ+zM/usVGGRgEAMY0eIarHPFGWNLnYeYyUmsbWQiz5pDkNI5Kmld8NR11uw9ncUjqEcYl5n1e+dKU8FqpOnzBWW9oMtXS7U8oSMQ5gOpUTy0njwyGOOk2Dz8oEFS3Ez0uIc29X8a9+tI2bISwRLMxrd4eQhhAAYB1NaFqY/+zY3Q7SEAIALBSFc6wRIwQAxDR6hIiUtNpOWCFDeQI3w5b5vzPTyj7GZ/0dhvc7jSFDnzFG6Kt0GDeNR31yjNCUL2HclIuryN+FKUaoDtSFldIQ4qjp0ur4pVxz+c3KoKDd/+zjQgs/od52NIQAAMtEsAyTZn5A7sJiaBQAENPoEQIArBOFD8vQEMIi5xMylM5kimD5ldE1hzFu53OownghYn7mGKG0w29ZjNAcSHEov031WlTSyUyJgCFWPpLHsuSgrFz83GXt/6MGe0Uy6bYpB/cCoyEEAFgmgjxC2/97IkYIAIhp9AgBAJaKsg4hDSFqiDLzzO9Xh6xM8QIpxCgHII1ljfEwvzGO5zCWDjcvUDpbeJENOfAW4rdfnXdoyq0M42gVM8eqrx3GQcQ6FuYFACDK0CMEAFgmkoT6GqpKtdEQwgbyr4ky10KESrcIMXDqV13ML6UsmK7jN1VOURH1MkyyKsYnzyOHQc7NUF4q1NpIZEQgcuQRAgBiWgRzjdrdEBIjBADENHqEAADLaNV4Itr0Fpu7hDSEqAVCPNtvivPJW+oIomo75KT3fun9YU2qpqYO+lVRPpyD4Z7b7rEpQGHz5s1Tp04tLS3Nz88vKCgwFzhx4sQzzzyzffv2Tp06TZw4sU6dOkKIwsLC+fPn79mzJzMz8+67787Ozj7X+RkaBQDUXocOHbr++utbtGgxfPjw+++//+233zaXuf3227du3Tp69Oh169aNHDlS3/nYY48lJCQMHz48NTX1mmuu2bhx47kuQY8QAGCZyNYjVBydPn16r169HnzwQSFERUXFlClT7rrrruACW7ZsWbly5eHDh+vUqXPttdc2atRo9+7dLVu2XLhwoX7mIUOGbNiwYf78+efqFNIjBABYRp9ZJqwv9cj9N998c9111+mvr7vuuu+//76srEwqcOWVV+rDoSkpKZdddtm6detEUPvq9/sPHDhwySWXnOsS9AgRBUL8fynNmhYi4uhXHTZP7iZfyrJgWugTnUfMkERA2CaCPEKlkpKStLQ0/XV6errf7y8pKWnRokWgwKFDhwIFhBBpaWnFxcXBZ3jppZdOnjwp9SOD0RACAGy2fPnyrl27Bu+57bbb9OHQpKSkQBfw7NmzQgi98xeQmJhYXl4e2CwrKwsuMHfu3KeffnrZsmWJiYnnujoNIQDAMhHFCEXnzp2fe+654J3NmjXTXzRt2nTfvn3663379iUmJgb3/6QCepkmTZror+fPnz9hwoTPP/+8Xbt2igoQIwQAWCf8GKGmiZSUlC5GGRkZ+vmGDh36wQcfnDp1Sggxc+bMIUOGOBwOIcQnn3yydetWIUReXl5xcfGaNWuEEMuWLTt9+vSNN94ohPjss8/GjRu3cOHCyy+/XF1leoSIfqHSEA2kgKLyVKEudUER9kNs6t+//6xZsy6//PKsrKyioqKlS5fq+x977LGCgoKOHTvWrVv3+eefHzhwYOfOnTds2PDSSy/Fx8cLIX73u9+Vl5fffvvtevnhw4c/9dRTVV6ChhAAUHs5HI65c+du37792LFjXbp0iYuL0/d/+umnCQkJ+uvRo0cPGDBg586d7dq1S09P13d+9913Pp8vcB69dawSDSEAwDKaiGAZptDlzUG+1NTU4M3MzMzMzMzgPcnJydWsAA0hYkxY46gAwhTBCvW2/x7ysAwAIKbRIwQAWIeFeQEAsczyuUYvAIZGAQAxjR4hAMAyETwsY/fIKA0hAMBa0RYjZGgUABDT6BECAKxUnQT5WoWGEABgmUhihHa3mzSEAADLRJI+YXcPkhghACCm0SMEAFgnCmeWqW6P8KmnnuratWtmZmaXLl3mzp0b2L958+YbbrghLS2tb9++O3furJlKAgCiQwSr8kZNQ/jzzz+/+OKLW7duffjhh0ePHv3NN98IIXw+35AhQ/r3779v376ePXsG1j8EACBaVLchfO6553r27Jmenp6fn3/11VcXFhYKIVauXHn8+PGJEyfWrVv373//+65du9atW1eTtQUA1Grar4/LhMPuLmHYD8scO3Zs06ZNV111lRDihx9+uOKKKxwOhxAiLi6uQ4cO27Zts76OAADUmPAelvF6vSNHjhw0aFDPnj2FEMeOHatXr17gaHJy8tGjR6t8o9/v93q951NRAIC9PB6P2+1Wl4nGPMIweoQ+n2/UqFEej+ff//63vic1NfXUqVOBAqWlpWlpaVW+V9M0l4snVAEgioVsBaNUdRsnv98/fvz4n3766dNPP42Pj9d3tm7desuWLX6/X9M0r9e7Y8eONm3a1FhVAQC1XvgJ9XaHCKvdIxw/fvyXX375yiuvHDp0aPfu3ceOHRNC9O7d2+Vy/fvf//b5fFOnTm3UqFH37t1rsrYAgNpNC//LbtVtCFetWlVeXp6fn5+bm5ubmzt9+nQhhNPpnDdv3quvvlq3bt133333vffes32hYQAAwlLdodHNmzdXub9bt24bN260rj4AgCgWjQ/L8AALAMAykUy6bXdLSEMIALCU3T28cLH6BAAgptEjBABYpnY8BxoeGkIAgGWiMUbI0CgAIKbRIwQAWCcKF+alIQQAWEcLe6jT9qFRGkIAgGUiSaivmZpUHzFCAEBMo0cIALAUMUIAQMzShKbZ3rKFiaFRAEBMo0cIALAMq08AAGJbDcyxtn79+ueff/748eO33nrrmDFjzAV++eWXJ554YsuWLdnZ2Q8++GD9+vWFECUlJf/5z3/Wr1+flpb2wAMPKM7P0CgAwDJaRBQnLC4u7tOnz1VXXfW73/3u8ccfnzZtmrnMsGHDioqKJk6cuHPnzhEjRug7P/vss9mzZ2/dunX+/PnqOtMjBADUXtOmTcvJyfnv//5vIcTTTz/96KOPjh07NrjA999/v2bNmiNHjiQmJnbt2rVhw4Y7d+5s06bNqFGjRo0aNWvWrJdeekl9CXqEAADL6DHCsL7UQ6nr1q3r2bOn/rpnz55btmw5e/asVODKK69MTEwUQtSvX/+yyy5bv359WHWmIQQAWEcL/0vp0KFDqamp+uu0tDQhRElJybkK6GWkAiExNAoAsNkXX3zRqlWr4D0jRoz4xz/+IYSoU6dOoAuov6hXr15wyTp16pSVlQU2z5w5U7du3bCuTkMIALBQJOsRXnvtta+99lrwzvT0dP1Fs2bN9u7dq7/es2dPUlKS3i8MaNq0aaCAEGLv3r1NmzYNqwIMjQIALBPByKgmRJ06dVoa6SkQQojhw4fPmzevtLRUCDFt2rRhw4bpDe1HH320adMmIcRNN910+PDhL7/8Ugjx+eefl5eX5+TkhFVneoQAAOtYvR5hXl5eTk5Ox44dGzVqdPLkyS+++ELfP2XKlIKCgssvvzwpKenll1/Oz8/v0KHDtm3bXn/99bi4OCHEsmXLhgwZUlFRUV5enpqaetNNN82ZM6fKS9AQAgBqL4fDMXPmzH379pWWlnbq1MnpdOr7//Of/7jdbv31nXfeOXDgwB9//LFt27aBCGKvXr127doVOE+gsBkNIQDAMiET5CPTvHlzaU+dOnWCN5OTk7t06RK8x+12N2jQoDonpyEEAFgmGuca5WEZAEBMo0cIALCMpokI0idqqDLVREMIALCU3UOd4WJoFAAQ0+gRAgCsFG0dQhpCAIB1IkifIEYIALiIWD2zzAVAjBAAENPoEQIALBONCfU0hAAAy9TQFGs1iqFRAEBMo0cIALBOFD4sQ0MIALCMFn7Mz/aRVBpCAIBlNKFpYXbxwi1vOWKEAICYRo8QAGCdCGKEdqMhBABYJ/w8QtsbToZGAQAxjR4hAMAyTLoNAIhtxAgBALFMExH0CGuoLtVFjBAAENPoEQIALBOFM6zREAIALBSFLSFDowCAmEaPEABgGU2LIB2C9AkAwEWDFeoBADGNGCEAANGFHiEAwDKRJNTb3SWkIQQAWEZj9QkAAKILDSEAIKYxNAoAsExEyzDVUF2qi4YQAGCZCGKEdreDDI0CAGIbPUIAgHWiMKGehhAAYJ3wY4S2BwkZGgUAxDR6hAAAy2jhd/Ds7hDSEAIArBNJ+oTdQcIL1BA6HI5Zzz2+4v23L8zlot2+ffuaNWsW/ppesevs2bMnTpxo2LCh3RWJJkeOHElKSqpTp47dFYkme/fubdGihd21sM13v/nNY489pi7Tp0V6uKdN3N1ob3x8pJWygOb3+y/Mlfbt21dZWXlhrhXtysvL4239WEQdv9/v8Xji4uLsrkg0qaiocLvd/L8Vlhj/3bzkkksSExMtP63f76+oqLDxxl64hhAAgFqIp0YBADGNhhAAENNoCAEAMY2GEAAQ05yTJ0+2uw4xze/3r1u37u233/7444+Lioo6duzocv2a0/LLL7+89NJLH374od/vb926tb31rJ0OHz48e/bspKSkjIwMfc+uXbteeOGFzz77LCMjo1GjRvZWr7bxer2zZ8+eOXPm2rVrMzIy9Jvm9/vnzJkzY8aMnTt3XnHFFYGPH4QQHo/n3XffnTVr1tq1a5s2bdqgQQN9/7Fjx/TfTSEEv5sXAXqENjt48ODQoUOPHTvWtGnTN954o0+fPl6vVwjh8Xiuu+667777rmXLlvfcc8+bb75pd01ro/Hjx0+cOHHFihX65k8//dStW7eKiorU1NQbbrjhm2++sbd6tYrX6+3Xr99rr73WpEkTj8ezZs0aff+DDz741FNPtWnTZuHChbfddpu9laxtfvvb37766quXXXbZyZMns7Ozd+7cKYSoqKjo1avXxo0bL7300rFjx86YMcPuauK8+WErj8fj8Xj01ydPnkxISPj222/9fv+8efM6dOhQWVnp9/sXLlzYsmVL/TUC5s6de8cdd1x//fX/+te/9D3333//nXfeqb9+5JFHhg4dal/tap1XX301Ozs78GHTlZaW1q1bd/PmzX6//8yZM8nJyevXr7epgrVOZWWl2+3euHGjvtm7d++pU6f6/f45c+ZcdtllPp/P7/cvWLCgTZs2+mtEL3qENnO5XIHBKK/X6/V669atK4RYuXJlnz59HA6HECI3N3fPnj379++3s6K1zM8//zx58uQXXngheOfKlSvz8vL017m5uYGeIoQQixcvHjFixOLFi59//vnCwkJ95/r16+vXr9+pUychRGJiYq9evbhpAQ6Ho127dt99950QorS0dNeuXfqNWrlyZd++ffWJCHJzc3fu3Hnw4EGb64rzQ0NYi0yYMGHQoEFt27YVQhQXFwfiXvHx8cnJycXFxbbWrnb54x//OHHiRGlOteCblpmZefTo0YqKCjtqVxvt2bPn3//+9/z580+ePJmfn//qq68KIUpKSgJ3TAjRsGFD/qYH++ijjx555JF27do1b958/PjxvXv3FsaPWVJSUt26dfndjHYExmuLyZMnr1+/fvny5fqmy+UKnpGO+cOCLVq0qLi4eNSoUdJ+l8ulR1iFEF6v1+l0Op3OC1252srhcHTo0EEPaHXp0mXUqFG/+93vzB+zWJ4/TFJZWTly5Mj8/Pxx48bt2bPnv/7rv7p3756TkxP8MRNCeL1efjejHQ1hrTBlypT3339/2bJlqamp+p6srKyioiL99fHjx0+fPt24cWP7Kli7zJo1a+/evVdffbUQYvv27Xv37j148OCjjz6alZUV6NAUFRU1atSIhjCgSZMmHTt21F936tTp6NGj+oequLjY5/Ppg/BFRUU9evSwtZq1yPfff//tt9+uXLnS5XK1b9/+jjvumDFjRk5OTvDH7OjRo2VlZfxuRjuGRu33wgsvTJs2bcmSJcEDfYMGDfr0009PnTolhPjggw+uvvpqkgECnnzyyXnz5r3++uuvv/56mzZt7rzzzjFjxgghBg0aNG/ePL/fL4SYN2/eoEGD7K5pLZKfnx8IDX799dctWrSoU6fO1Vdf7Xa7//Of/wghDh48uGbNmgEDBthazVokLS3N4/EcOHBA39y1a1d6eroQYtCgQYsWLTp9+rQQ4oMPPrjmmmv0/YheTLptsx9//LFNmzYtWrRIS0vT9zz77LM33nijEOK2227btWtX586dFy1aNHfu3D59+tha01rqhhtuKCgoGDdunBCitLS0V69eDRs2TElJWbNmzapVq2J5xRzJmTNnevfunZiY2LJly4ULF06bNk3/R2HmzJn333//gAEDVq5cmZ+f/9xzz9ld01rkD3/4w4IFCwYOHLh79+6tW7euXr26WbNmfr9/8ODBP/30U3Z29qJFi+bNm6f/wiJ60RDarKysbMuWLcF7WrVqlZKSIoTw+XwrVqw4dOhQr169mjRpYlMFa7vt27enpqYGHl4oKytbsmRJeXl53759k5OT7a1bbePxeJYtW3bq1KkePXoEj+Zt3759/fr1LVu27N69u43Vq502b968devWBg0a9OrVK7ACUeB387rrrsvKyrK3hjh/NIQAgJhGjBAAENNoCAEAMY2GEAAQ02gIAQAxjYYQABDTaAgBADGNhhAAENNoCAEAMY2GEAAQ02gIAQAxjYYQABDT/h8niMAfq68+CgAAAABJRU5ErkJggg==",
"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
}