diff --git a/docs/source/tutorials/ssh_model_myqlm.ipynb b/docs/source/tutorials/ssh_model_myqlm.ipynb new file mode 100644 index 0000000..c76c3e6 --- /dev/null +++ b/docs/source/tutorials/ssh_model_myqlm.ipynb @@ -0,0 +1,170 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Simulating the SSH model\n", + "\n", + "Here we simulate the SSH model using the MyQLM backend.\n", + "\n", + "For the background on the SSH model see the [previous notebook](https://ichec.github.io/qse/tutorials/ssh_model.html)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "import qse\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Angle is: 54.74 degrees\n" + ] + } + ], + "source": [ + "angle = np.arccos(np.sqrt(1.0 / 3.0)) * 180 / np.pi\n", + "print(f\"Angle is: {angle:.2f} degrees\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create two chains of qubits, A and B.\n", + "By rotating them by the angle above, there will be no interactions between qubits on the same chain." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAGwCAYAAAAqkitTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAARJpJREFUeJzt3QlclNX+P/DPDMsAgoADqKS4gApuWGpqm5XmnkvSZouWWZTWLa9d81aaZel1abGbem/Xsl83l6tlqX+tLFMrN1zIXMstXFhc2GQZkHn+r+9hhgYURAVm5uHzfr2e1zAzzyyPg/PhnPM95zFomqaBiIhIR4zOfgNERERVjeFGRES6w3AjIiLdYbgREZHuMNyIiEh3GG5ERKQ7DDciItIdT+ic1WrFqVOnEBAQAIPB4Oy3Q0RUa2mahuzsbISHh8NorN62le7DTYKtcePGzn4bRERkc/z4cTRq1AjVSffhJi02+z9m3bp1nf12iIhqraysLNXYsH8vVyfdh5u9K1KCjeFGROR8NTFExIISIiLSHYYbERHpDsONiIh0h+FGRES6w3AjIiLdYbgREVGFLhRZsWnPMSxZl4gFaxLQvfdAVfFo38xmM/r06YPdu3eXPGbgwIGIiIiAj48PGjZsiEceeQTJyckVvs6IESNKPa99a9OmDa6UQe9n4pZ5FYGBgcjMzORUACKiK3A2Mwert+zHso2/4NDZMyg0WCGBkbpuGax553HzkKcwoGtrtGsUiJn/eFOFW1JSknrsO++8g27duqlgO3nyJMaNG4eioiJs27at3O9juT0vL6/k+oULFxAbG4tnn30Wr7322pW8dYYbERFdbPfhUxg3byWOZp+DxVQIq58GeElqAOlfLoc1Lx+hfYfBZPFCs4B6GNapPh57YAjS0tIQGhp60fOtWLECgwcPVktwVfb7+Msvv8Q999yDo0ePokmTJrgS7JYkIqKLgm307M9xMD8VeeYCFNXXoIUBWgNAawhovgBMULfL/fuzkjBhykxENGmquijLOnfuHD777DN06dIFV2L+/Pno2bPnFQdbrVihhIiIrqwrUlpsSdZ0XAjWgCBZxxCAj22TJpEPkJ/wG5JffxPST6kVFMLDLwCdhz6L9Ow8mAPrqOcaP348/vnPfyI3Nxddu3bFokWL0KxZs0qvC7xmzRosXLgQV4MtNyIiKiFjbNIVWVTXFmzBts3ssJkA7/bNEDojHqGT4xHyl1EwRUUiYdl7+HTFDyXP9eKLL2LXrl349ttv4eHhgaeeeqrkPn9//5ItPj4eZX3yyScICgpSXZlXgy03IiIqqYqU4hEZY9Pq2FpssgUCqGNruXkA8AaMAV7wjDYDmcWB5xXSECmvT8PcefMwZlg/eHoYERISoraWLVsiJiam1BlaEhMTS34uO/4m43IfffSRqrD09vbG1WC4ERGRsm1/kqqKtEp3pK8tzPxswRZg+9kWbio95PZCAFLgKPsbgDPnM5FwIAnd2jRF2XNrOoqKikJ5NmzYgEOHDmHkyJG4Wgw3IiJSjqdlqHJ/VRVpchhns4ecP4rv8wS0oiIUnc8GcgHrmXzkrNgGraAAvk2isW7Dj9jxwyrccsstCA4OxuHDh/Hqq6+q8TapfKxMIYkUn7Rt2xZXi+FGRERKnqVQzWOTFpiqyLBvHrbNFmxym2XjIaRunKUeZ/D1hmdYCOrdex98zc1ghSe++OILTJo0CTk5OWqum0zyltCS7smKyDSBzz//HO+99x6uBcONiIgUX5OXyjWVcFaHrci2SRek1Je8OwTBU4YUj7dlSIklgNOAIQUwpBnQKiYGr65bh0vNO74cmZcs1ZXXiuFGRERK47AgeGnG4hCzAMh32KTVBlsLToIut8z9sn+B7GZUz+NsnApARETKjTERiDKHwJhrKC4SybeFWA6AbOkzdNiybbfbQy4PMOYZ0MIcis7REXA2hhsRESlSvh93W6xaUstgDzR7qGXYtnSHn+0hlw21v8niibjusep5nM2p72Du3Llo3769muMgmyyyKTPS7fLz8zF69Gi1nItM9Bs6dChSU1Od+ZaJiHStX9cYtVakR5bhzzBLt42rld3s92VA7d88wIy+XaLhCpwabo0aNcK0adOwY8cObN++HXfeeScGDRqEvXv3qvtfeOEFrFy5EkuXLlXzHmQ5FllEk4iIqoc5sA5mxt+NCGMwPNMNMNiKRXDGdnm69HW5X/aT/WfGDyxZesvZXO6sAPXq1cOMGTMQFxenVpaWdcXkZ3HgwAFVRrp582a1TtmlWCwWtTlW58iseJ4VgIjoyhZPfnHeShy5xFkBVDVlQfEYm3RFSotNgq1dZEOXOUuLy1RLynl+pIUmcyKke1Jac4WFhWpFaLvo6Gh18ruKwm3q1KmYPHlyDb5zIiL9aR8ZjoWvPIw1Ww9g6YbE4vO5wQrNoMGgGVRVpBSPxPWPVV2RrtJic5lw+/XXX1WYyfiajKstX74crVu3VuuOyZpisnCmo/r16yMlJaXc55swYQLGjh17UcuNiIiujATWw7064oEe16sltWQFk9z8Qvj5eKlyf6mKdIXiEZcMt1atWqkgk2bqsmXLMHz4cDW+drVMJpPaiIioakiAyVqR3drAbTg93KR1Zl9As2PHjkhISFDLrtx///0oKChARkZGqdabVEs2aNDAie+YiIhcncu1J2XlaCkIkaDz8vLC999/X3LfwYMHkZSUpLoxiYiIXLLlJuNjffv2VUUi2dnZqjJy/fr1+Oabb1RFjZzuQMbPpIJSKmueffZZFWzlFZMQERE5PdzS0tLw6KOPIjk5WYWZTOiWYLvrrrvU/e+88w6MRqOavC2tud69e2POnDn85IiIyL3muVW1mpxXQURErvF97HJjbkRERNeK4UZERLrDcCMiIt1huBERke4w3IiISHcYbkREpDsMNyIi0h2GGxER6Q7DjYiIdIfhRkREusNwIyIi3WG4ERGR7jDciIhIdxhuRESkOww3IiLSHYYbERHpDsONiIh0h+FGRES6w3AjIiLdYbgREZHuMNyIiEh3GG5ERKQ7DDciItIdhhsREekOw42IiHSH4UZERLrDcCMiIt1huBERke4w3IiISHcYbkREpDsMNyIi0h2GGxER6Q7DjYiIdIfhRkREusNwIyIi3WG4ERGR7jDciIhIdxhuRESkOww3IiLSHYYbERHpDsONiIh0h+FGRES6w3AjIiLdcWq4TZ06FZ07d0ZAQADCwsIwePBgHDx4sNQ+t99+OwwGQ6ktPj7eae+ZiIhcn1PDbcOGDRg9ejS2bNmCtWvXorCwEL169UJOTk6p/UaNGoXk5OSSbfr06U57z0RE5Po8nfniX3/9danrCxYsUC24HTt24Lbbbiu53c/PDw0aNHDCOyQiInfkUmNumZmZ6rJevXqlbv/ss88QEhKCtm3bYsKECcjNzS33OSwWC7KyskptRERUuzi15ebIarXi+eefx80336xCzG7YsGFo0qQJwsPDsXv3bowfP16Ny33xxRfljuNNnjy5Bt85ERG5GoOmaRpcwNNPP401a9bgp59+QqNGjcrdb926dejRowcOHTqEyMjIS7bcZLOTllvjxo1Vq7Bu3brV9v6JiKhi8n0cGBhYI9/HLtFyGzNmDFatWoWNGzdWGGyiS5cu6rK8cDOZTGojIqLay6ljbtJolGBbvny5apE1a9bsso9JTExUlw0bNqyBd0hENWXEiBGlpvyYzWb06dNHDUeIY8eOYeTIkep7wtfXV/1xO2nSJBQUFFz2uaU35+WXX1ZDHPLHb9OmTfHRRx/VwFGRszi15SbTABYuXIivvvpKzXVLSUlRt0uzVX55Dx8+rO7v16+f+kWXX/IXXnhBVVK2b9/emW+diK7RhSIrtu1PwvG0DORZCnE0+Ry63twdS5cshKeHUX0fvPLKKxgwYACSkpJw4MABNTb/r3/9C1FRUdizZ4+aJiRTh2bOnFnha913331ITU3F/Pnz1WNlSpE8F+mXU8fc5K+zS/n444/VX3HHjx/Hww8/rH6J5RdYxs6GDBmifuEr219bk328RHR5ZzNzsHrLfizb+AsOnT2DQoMV8iWUum4ZUJCPO4aNRdxtsejXNQb7f92FW2+9FWlpaQgNDb3ouWbMmIG5c+fiyJEjFU45euCBB9Q+ZSuxqWbVmjG3y+WqhJlM9CYifdh9+BTGzVuJo9nnYDEVwhqsAV7yly5Q5FMEq1aEXedPYt/yNHy8+meYT29XLS3pubkU+ZK8XGCtWLECnTp1Uos/fPrpp6hTpw4GDhyIN954Q/UQkT65REEJEdWOYBs9+3MkWdNRZNag1QEg2SL1X0ZA8wXy9/yGE/OmQJpyhwoL4V2nLv7vs//BaLy4PECKyt5///3LdklKi02qsH18fNT4/pkzZ/DMM8/g7NmzqpeI9InhRkQ10hUpLTYJtgvSWgsCEADAx7ZJdvkA3jHNEPRgfyAHsKblI+/nBDw67F602rIVHdrFlDzfyZMnVbHJvffeq8bd7Pz9/Ut+liGNefPmqbE1GQKRxSCkS0y8/fbbiIuLw5w5c9h60ymGGxFVOxljk65IabGpYAu2hZufQ7iZAGOAFzxbmYFsACGAd2hDJM+chklTZ+KrhfPVc506dQp33HEHbrrpJvz73/++ZDW1sI/pSGX1ddddVxJsIiYmRg2LnDhxAi1atKixfweqOQw3Iqr2qkgpHpExNtUVGWDbJGvq2MLNA4A3isff7Bl0AdBkpT2DAXsPn1DPk5qSrIKtY8eOqkuxbHeljM+VJaseLV26FOfPny9p2f3222/qsZebV0vuy6XWliQi/ZFyf6mKtPppxWNsPrYWmz3oAm2bN6BZi1CUn622wrOnkfn1GmiFFlyo3wSrN2xTp8CKiIhQ42ynT59W0wXsU4jKI0v4SUHKY489hn379qnFIl588UU8/vjj7JLUMbbciKhayTw2KfdXrTKTwzibPeSkMSX3eQKWnw8hdeAs9TiDrzc8w0JQL+4+eJub4etvvlVFJLKVbXFVVHktrTU5pdazzz6rqiYl6GTe25QpU6r92Ml5GG5EVK1kgraKHoOtr8i+edg2W7AFzx2C4H8MAeTkIOlShQLgNGBIBrQ0DTfe3hdz/jHxqt5DdHS0CjiqPRhuRFStfE1eKtdUwlkdtiLbVmjbsdDhNsf9NFltwgA/H0lBosphuBFRtWocFgQvzVgcXnLCjnyHzZ5XHrZQyy1zv+xfILsZ1fMQVRYLSoioWt0YE4EocwiMuQYgzxZaEmI5KC75z3TYsm2320MuDzDmGdDCHIrO0RHOPhRyIww3IqpWsgiyrBVpsnjBYA80e6hl2LZ0h5/tIZcNtb/J4om47rHqeYgqi78tRFTtZBHkZgH14JFl+DPM7EUjZTf7fRlQ+zcPMKNvl2hnHwK5GYYbEVU7c2AdzIy/GxHGYHimG2CwVULijO3ydOnrcr/sJ/vPjB+oHk90JRhuRFQj2keG44PnhiLapz58z3rDI9UAQxpgSCku91eXqVC3+571UvvNeS4O7SJ5YmJys/O51QSez43I9RZRXrP1AJZuSCw+nxus0AyaKveXqkgpHpExNumKZItNX7Jq8PuY4UZETiFrRSYcKD4Td25+oZrHJuX+UhXJ4hF9yqotJyslotpLAqxbm6bo1sbZ74T0iH8eERGR7jDciIhIdxhuRESkOww3IiLSHYYbERHpDsONiIh0h+FGRES6w3AjIiLdYbgREZHuMNyIiEh3GG5ERKQ7DDciItIdhhsREekOw42IiC55SqJNe45hybpELFiToC7l+o8//QwPDw/079+/1P5nz55Fnz59EB4eDpPJhMaNG2PMmDHqNDcV+eKLL9CpUycEBQWhTp066NChAz799FNcK57yhoiISp1MdvWW/Vi28Zfik8karJCTfhoAeGlGnN+yCj0H3ocN363CqVOnVJgJo9GIQYMGYcqUKQgNDcWhQ4cwevRonDt3DgsXLkR56tWrh5dffhnR0dHw9vbGqlWr8NhjjyEsLAy9e/fG1eLJSomISNl9+BTGzVuJo9nnYDEVwuqnAV62ZNMAa44FKe/PQuN7RsOS+CPu698Ds2dNQ3lmz56NGTNm4Pjx41f0fXzDDTeoluEbb7yBq8VuSSIiggTb6Nmf42B+KvLMBSiqr0ELA7QGgNaw+DLv5F54hoWgKDIQ1qiW+HD+fPxy6OQln09addLl2L1790q/B2lrff/99zh48CBuu+22azoehhsRUS13NjNHtdiSrOm4EKxBMwMIBRBiu7RtOTt3wfeW9up+zw6RKLDk4omXZ6nH2z344IPw8/PDddddp1pn//nPfy77+tKS8/f3V92S0mJ7//33cdddd13TMTHciIhqudVb9quuyKK6GhAEINi2mf/cLuSeQeGRk/C9q626z1DPA75t2uDgrh+wZuuBkud65513sHPnTnz11Vc4fPgwxo4dq25PSkoqGZ+Ty7feeqvkMQEBAUhMTERCQgLefPNN9Zj169df0zFxzI2IqJZXRQ6d+DF2nT+puiJVK02CLRBAHQA+ADyAzLe/Rc7/bQKMMgBnowEGD0/0fe5tfDX9GXh6lG4v/fTTT7j11ltVF6UUmfz6669qPE3Cr0mTJqqY5FKeeOIJNU73zTffXPVxsVqSiKgW27Y/SVVFWoM1wNcWZn62YAso/lnTipC3+hfUHdcLpthIQKr7MwGkA+c+XIydCd8j4cAAdGvTtNRzW61WdWmxWODp6YnIyEh1XS4ramzI4+Qx14LhRkRUix1Py1Dl/qoq0mQLNx+HkPMH8r/9DdasfPg9dAOMVh8gQwbqAJwGfGNaI/1AAlasWIkD2/zRuXNnNX62d+9evPjii7j55pvRtGnp0HM0depUNc9NAk8CbfXq1Wqe29y5c6/puBhuRES1WJ6lUM1jU+X+RofNw7Z5AbkLd8HUvTmMQT7FLTaH/Xxbx+D85p/xx7Fj2LByM1544QUVUjKJ+5577sFLL71U4evn5OTgmWeewYkTJ+Dr66vmu/33v//F/ffff03HxTE3IqJabMm6RLy6dA1yQwtUub8ac5MikiDbuJt/ccChEMB5W3ekQ8vNkAL4pXljyv19cd8dHVzm+5gtNyKiWqxxWJBaeUSFlwxz5TtsEmqwteCKAOSWuV/2L5DdjOp5XAmnAhAR1WI3xkQgyhwCY64ByLOFloSYTF3LtrXU7Fu27XZ7yOUBxjwDWphD0Tk6Aq6E4UZEVIt5ehgRd1ssTBYvGOyBlu3Q/ShbusPP9pDLhtrfZPFEXPfYi6YBOJtT341UyUhljUzgk0UyBw8erJZdcZSfn68W3zSbzaoCZ+jQoUhNTXXaeyYi0pt+XWPQLKAePLIMf4ZZum1crexmvy8Dav/mAWb07RINV+PUcNuwYYMKri1btmDt2rUoLCxEr169VPWMnVTerFy5EkuXLlX7y2RAqcAhIqKqYQ6sg5nxdyPCGAzPdAMMtmIRnLFdni59Xe6X/WT/mfED1eNdjUtVS54+fVq14CTEZNFMqaiRWe1yuoS4uDi1z4EDBxATE4PNmzeja9eul31OVksSEVV+8eQX563EkXLOCiDFIzLGJl2R0mKTYGsX2bCSz16LqyXlgIV9SZYdO3ao1lzPnj1L9pE5EBEREeWGm8yvcJzZfrkT5RERUbH2keFY+MrDaq3IpRsSi8/nBis0gwaDZlBVkVI8Etc/VnVFumKLzeXCTZZbef7559Vs9rZt26rbUlJS1CrRcoZWR/Xr11f3lTeON3ny5Bp5z0REemMOrIOHe3XEAz2uR8KBJLWCSW5+Ifx8vFS5v1RFulrxiEuHm4y97dmzRy20eS0mTJhQsgq1veUmM+WJiKjyJMBkrchubeCWXCLcxowZo04tvnHjRjRq1Kjk9gYNGqCgoAAZGRmlWm9SLSn3XYrJZFIbERHVXk5tW0otiwTb8uXLsW7dOjRr1qzU/R07doSXl5c6M6udTBWQ8wJ169bNCe+YiIjcgaezuyKlElJOaidz3ezjaFJNIwtoyuXIkSNVN6MUmUh1zbPPPquCrTKVkkREVDs5dSqAweBw0jsHH3/8MUaMGFEyifuvf/0rFi1apKoge/fujTlz5pTbLVkWpwIQEbmGmvw+dql5btWB4UZEVPu+j12/npOIiOgKMdyIiEh3GG5ERKQ7DDciItIdhhsREekOw42IiHSH4UZERLrDcCMiIt1huBERke4w3IiISHcYbkREpDsMNyIi0h2GGxER6Q7DjYiIdIfhRkREusNwIyIi3WG4ERGR7jDcSFcuFFmxac8xLFmXiAVrEtC990AYDIaSzWw2o0+fPti9e3fJY958803cdNNN8PPzQ1BQUKVe59ixY6We175t2bKlGo+OiCrLs9J7Ermws5k5WL1lP5Zt/AWHzp5BocEKDUDqsZOo07glbh7yFAZ0bY12jQIx8x9vYsCAAUhKSlKPLSgowL333otu3bph/vz5V/S63333Hdq0aVNyXcKTiJyP4UZub/fhUxg3byWOZp+DxVQIa7AGeAEwAEU+RbBqBuzHeRzZsAvNAuph2CNP4P89MASnT59GaGgoJk+erJ5nwYIFV/zaEmYNGjSohqMiomvBbkly+2AbPftzHMxPRZ65AEX1NWhhgNYA0BoCmi8AE9Ttcv/+rCRMmDITEU2aVkkra+DAgQgLC8Mtt9yCFStWVMkxEdG1Y8uN3LorUlpsSdZ0XJDWmgyXBQDwsW3yp5sPkJ/wG5JffxPST6kVFMLDLwCdhz6L9Ow8mAPrXNVr+/v7Y9asWbj55pthNBrx+eefY/Dgwfjyyy9V4BGRczHcyG3JGJt0RRaZbcEWbAs3P4dwMwHe7Zsh6PH+QA5gTclH7g8JSFj2Hj69qxuef2TAZV9HxtT++OMP9fOtt96KNWvWICQkBGPHji3Zp3Pnzjh16hRmzJjBcCNyx3A7evQofvzxR/WfPTc3V41ZXH/99Wow3sdHvlGIaqYqUopHZIxNq2MLNdkCAch1+VX0AOANGAO84BltBjJlkAzwCmmIlNenYe68eRgzrB88PSrunV+9ejUKCwvVz76+0s95aV26dMHatWur/FiJqBrD7bPPPsN7772H7du3o379+ggPD1f/0c+dO4fDhw+rYHvooYcwfvx4NGnS5CreClHlbdufpKoiVfGIry3MpMVmDzq/P8NN/ZbL7ZJPeZJQxcUmZ85nIuFAErq1aVrha1X29zkxMRENGzasoiMkomoPN2mZeXt7Y8SIEWpsoXHjxqXut1gs2Lx5MxYvXoxOnTphzpw5qrSaqLocT8tQ5f6qKtLkMM5mDzl/FN/nCWhFRSg6nw3kAtYz+chZsQ1aQQF8m0Sr57kuIEn9kSZTA4qKilRIiaioKDW2dimffPKJ+j8h/zfEF198gY8++gj/+c9/avTfgYiuIdymTZuG3r17l3u/yWTC7bffrjaZECsTXImqU56lUM1jkxaYGluzbx62zRZscptl4yGkbpylHmfw9YZnWAjq3XsffM3NkJtfiIkTJ6qwsrMH1g8//KB+p8vzxhtvqO55T09PREdHY8mSJYiLi6v+gyeiyzJomqa+I/QqKysLgYGByMzMRN26dZ39dqiKyAokry5dg9zQAlX2j9Di8TRVWBLo0HKTrsjzKB5vy5ASSwCnAUMK4JfmjSn398V9d3Rw9uEQ1QpZNfh9fNXVkmlpaWqzWq2lbm/fvn1VvC+iCjUOC4KXZiwOLwuAfIdNQg22FlwRVHdkqftl/wLZzaieh4j054rDbceOHRg+fDj2798Pe6NP1tSTn+VSxiyIqtuNMRGIModg1/mTKMrTikMr1yHYCh3CTe7LcQi5PMCYZ0ALcyg6R0c4+UiIyCXC7fHHH0fLli3VGnxSNSmBRlTTpHw/7rZY7FuehrycAmg+Dr/NhQ7z3KRjwR582cWbIQcwWTwR1z/2stMAiKiWhNuRI0dUxaRUkhE5U7+uMfj0ux04mJWKCx62oeMLtnL/suGWbwu3DMAjy4DmAWb07RLt7EMgompyxX+29ujRA7/88kv1vBuiKyBLZ82MvxsRxmB4phtgsBWL4Izt8nTp63K/7Cf7z4wfeNVLbxGRDqslz5w5o8bcbrzxRrRt2xZeXvZBjmKutvQQqyVrx+LJL85biSP2swL4/XlWADVfoKB4jE26IqXFJsHWLpKTrYn0/H18xeG2cuVKPPLII+pNXvRkLlhQwnCrPYsor9l6AEs3JBafzw1WaAYNBs2gqiKleCSue6zqimSLjcg5XDrcmjZtqk70+Oqrr6qCElfHcKtdZM1JWVJLVh6RCdp+Pl6q3F+qIlk8QuRcLh1uAQEBanmiyMhIuAOGGxFR7fs+vuI/Ze+55x61LBEREZFupgLIHLcJEybgp59+Qrt27S4qKHnuueeq8v0RERFVf7dks2bNyn8yg0HNg3Ml7JYkInINLr22pJyslIiIyJVVWflYcnIypk+fXlVPR0REVLNrS16KnNdq27Zt+Nvf/nb174aIiMgZ4Zaenl7qukzalnE2OUuAnIGbiIjI7bolly9fXmpbsWIF9uzZg9dffx1ffvnlFT3Xxo0bcffddyM8PFwVo5R9/IgRI9TtjlufPn2u9C0TEVEtU2Vjbg8++CDWr19/RY/JyclBbGwsPvjgg3L3kTCT8Tz7tmjRoip4t0REpGdXfSbusuRMAddff/0VPaZv375qq4jJZEKDBg0q/ZwWi0VtdpdaA5OIiPTtisNt7NixF92WmpqKr776Cv379y91/9tvv33Nb1Bag2FhYQgODsadd96JKVOmwGw2l7v/1KlTMXny5Gt+XSIiqkWTuO+4447KPbHBgHXr1lX+jRgMagxv8ODBJbctXrwYfn5+auL44cOH8fe//x3+/v7YvHkzPDw8Kt1ya9y4MSdxExE5mUtP4q7JdSUfeOCBkp9lqa/27durBZulNScnTS2vG1M2IiKqvdzqHCDNmzdHSEgIDh065Oy3QkRE7h5uUrG4ZcuWy+6XnZ2Nf/zjHxVWP16LEydO4OzZs2jYkGdRJnLX8+1t2nMMS9YlYsGaBHUp13/86Wc11CDj9mUL1aQSW4YWfH19ERMTg/fee++yr/Phhx/i1ltvVWP1svXs2VMtMkG1R6W6Je+9914MHTpU9ZXKvLROnTqpuWk+Pj5qUve+ffvUWQJWr16tfjlnzJhRqRc/f/58qVaYrFsp54qrV6+e2qQwRF5XqiVlzE1WP4mKikLv3r2v/oiJyClnSl+9ZT+Wbfyl+EzpBitksN8AwEsz4vyWVeg58D5s+G4VTp06pb5fxI4dO1RB2X//+18VcJs2bcKTTz6pgnDMmDHlvp4MXUgo3nTTTep7Sv7o7tWrF/bu3YvrrruuBo+cXL6gRIo0li5diiVLlqggkwFB9QQGA1q3bq0CZ+TIkeovq8qSX8BLFagMHz4cc+fOVcUlu3btQkZGhvpll1/ON95444rOAM6zAhA51+7DpzBu3koczT4Hi6kQVj8N8LIlmwZYcyxIeX8WGt8zGpbEH3Ff/x6YPWtauc83evRotSLSlRSsyUpK0oL75z//iUcffbSKjox0UVAiRRoPP/yw2oS8uby8PFWWX/acbpV1++23o6Js/eabb67qeYnIdYJt9OzPkWRNR5FZg1YHgK98odgGRaxA3sa98AwLQVFkIKznW+LD+fMx8ulnERt16RaWfPdIz86VyM3NRWFh4RU/jmphQYmkr3QXXm2wEZH+uyKlxSbBdiFYgybTU0MBhNgubVvOzl3wvaW9ut+zQyQKLLl44uVZ6vFlSbek9B5J1+SVGD9+vOr9kbE3qh3cqlqSiNyHjLFJV2RRXQ0IAhBs28x/bhdyz6DwyEn43tVW3Weo5wHfNm1wcNcPWLP1QKnnkzVsBw0ahEmTJqkhCpGUlKTmvtq3t95666L3MW3aNDVnVubRyvgb1Q5VtvwWEZFjVaQUj8gYm+qKDLBtgQDkumSMB5Dz8U6gyIrUJ2f9+WANMHh4YtE3m/FAj+vh6WFURWsyt1VabK+88krJrtIakyI0u7LdjjNnzlTh9t1336l5slR7MNyIqMpt25+kqiKtwVrxGJuEmZ8t2CTk/ABNK0Le6l9Qd1wvmGIjAVkGVurU0oFzHy7GzoTvkXBgAOoiRy29J4Vmb775ZqnX8fT0VBXUlyInT5b9ZexeKrypdmG4EVGVO56Wocr9VVWkyRZuPg4h5w/kf/sbrFn58HvoBhitPkCGDNQBOA34xrRG+oEE/Lx1O6a/9IyqxpZ1a1NSUtTzy1SA0FAZtLs0Kf2fOHEiFi5ciKZNm5Y8zt59Sfp3xWNu8teTnIeNiKg8eZZCNY9NlfsbHTYP2+YF5C7cBVP35jAG+RTf5rCfb+sYWE6fxKL5H+D06dNqnpss3mDfOnfuXOHry1SigoICxMXFlXqcdFNS7XDFLTcpw5WKoyZNmuCxxx5TYcdJkUTkyNfkpXJNJZzVYSuybYWA+bNh6hLnbbc57Ocd3ggtnnoL4+/vi/vu6HDFr3/s2LHqOCzSc8tNzpZ98uRJPP3006okV5r8ck62ZcuWqXkkRESNw4LUyiMqvOQkHfkOW64t0LJsl7ll7pf9C6RxZ1TPQ1RjUwGkr1v6v2Xdt61bt6oB3UceeURVLr3wwgv4/fffr+rNEJE+3BgTgShzCIy5BiDPIdRk6lq2rXDEvmXbbreHXB5gzDOghTkUnaMjnH0oVBvnuSUnJ2Pt2rVqkwHefv364ddff1XLcb3zzjtV9y6JyK1I+X7cbbEwWbxgsAeaPdQybFu6w8/2kMuG2t9k8URc91j1PERX44p/c6Tr8fPPP8eAAQPUuJusN/n888+rxU4/+eQTNZ/kf//7H15//fWrekNEpA/9usagWUA9eGQZ/gyzdFtFZNnNfl8G1P7NA8zo2yXa2YdAtamgRCqOrFarWnFbTiHRocPFg72yGHJQEPvKiWozc2AdzIy/u3htyfR0FBVp0PIvXltSjbHlFbfYJNgijMGYGT9QPZ6o2s8KYPfpp5+qU+C4yzI2PCsAkfMXT35x3kocKeesAFI8ImNs0hUpLTYJtnaRPGejHmXV4PfxFYebu2G4ETmfLIIsa0Uu3ZBYfD43WKEZNBg0g6qKlOIRGWOTrki22PQri+FWdRhuRK615mTCgSS1gklufiH8fLxUub9URbJ4RP+yXPF8bkRE10oCrFubpujWxtnvhPSOfyoREZHuMNyIiEh3GG5ERKQ7DDciItIdhhsREekOw42IiHSH4UZERLrDcCMiIt1huBERke4w3IiISHcYbkREpDsMNyIi0h2GGxER6Q7DjYiIdIfhRkREusNwIyIi3WG4ERGR7jDciIhIdxhuRESkOww3IiLSHYYbERHpDsONiIh0h+FWCReKrNi05xiWrEvEgjUJ6rLfoDgYDIaSzWw2o0+fPti9e3fJ486dO4eHHnoIdevWRVBQEEaOHInz589X+nV//vlneHp6okOHDtV0ZERE+uTp7Dfgys5m5mD1lv1YtvEXHDp7BoUGKzQABgBpB44irHlbvDJlBnp0bImC3Cy88sorGDBgAJKSktTjJdiSk5Oxdu1aFBYW4rHHHsOTTz6JhQsXXva1MzIy8Oijj6JHjx5ITU2tgaMlItIPg6Zp8n2tW1lZWQgMDERmZqZqQVXW7sOnMG7eShzNPgeLqRBWPw3wsiWbBqQvXw4t14Lr7noUzQLqYWb83chKPoJbb70VaWlpOHPmDFq3bo2EhAR06tRJPefXX3+Nfv364cSJEwgPD6/w9R944AG0aNECHh4e+PLLL5GYmHjN/xZERO74fXw12C1ZTrCNnv05DuanIs9cgKL6GrQwQGsAaA1tl76AZtLU/bJf/MyFeG/OvxAVFaW6KDdv3qy6Iu3BJnr27Amj0YitW7dW+Poff/wxjhw5gkmTJtXA0RIR6Q+7JS/RFSkttiRrOi4Ea0AQgAAAPrZN/hywAjAB+Tt+w6mZb6qWXFJhIUx1ArH6//0/FWApKSkICwsr9dwyflavXj11X3l+//13vPTSS/jxxx/V/kRE5GYtt40bN+Luu+9WXXRSlCHdb46kx3TixIlo2LAhfH19VctHvvyrk4yxSVdkUV1bsAXbNnOZzQR4t2mG0MnxCB0Xj9ARo+B9XSSGDB6IP/74o1Kv5e/vX7LFx8ejqKgIw4YNw+TJk9GyZctqPU4iIj1zatMgJycHsbGxePzxx3HPPfdcdP/06dMxe/ZsfPLJJ2jWrBleffVV9O7dG/v27YOPjzSjqr4qUopHZIxNq2NrsckWCECuy0t6ACgC4A0Y/b3gGWUG0iWpAJN/OJL/NQ3/+te/ERUVqcbeSj3/hQuqgrJBgwbquuM4mvQ/Z2dnY/v27di1axfGjBmjbrdarSrkpRX37bff4s4776zy4yYi0hunhlvfvn3Vdinyhf7uu++qCsRBgwap2/7v//4P9evXVy08Kbioatv2J6mqSKt0R/rawszPFmwBtp/t4eZl+9eT2/Kg9peiE6nOOXoqFY888rCqeNyxYwc6duyonn/dunUqrLp06aKuy/icI7nv119/LXXbnDlz1OOWLVumAp6IiC7PZQd1jh49qsampCvSTqpsJBikWKO8cLNYLGpzrM6prONpGarcXwWXyWGczR5y/rZQKyz+l9OKilCUmw3kA9b0fOSu2wZrYQFaxnZFTEyMmvc2atQozJs3T00FkNaYvO/yKiVlrK5t27albpNxO2mllr2diIjcMNzsRRfSUnMk1ysqyJg6daoas7oaeZZC1fJS5f5Gh83Dttlba7Z9LJsOIXXTrOKrPt7wrBeCBr0eQLPoWHXbZ599pgJN5qpJcA0dOlR1sxIRUS0Nt6s1YcIEjB07tlTLrXHjxpV6rK/JS+WaSjirw1Zk26TFJgqB4H8MQfDfhwAZUmIJ4DRgSAH80rzh5yMpCFUZWZkJ2xV57bXX1EZERDoIN3vRhazOIdWSdnK9ouWoTCaT2q5G47AgeGnG4hCTns18h604r/4cc8stc7/sXyC7GdXzEBGR87jsJG4pnpCA+/7770u1wmQCdLdu3arlNW+MiUCUOQTGXENxkUi+LcRyAGQDyHTYsm2320MuDzDmGdDCHIrO0RHV8v6IiMgNWm6yiPChQ4dKFZFIebx050VEROD555/HlClT1DJU9qkAUowxePDgank/nh5GxN0Wi33L05CXUwDNx+FfqLDMJG578EnIZQOGHMBk8URc/1j1PEREVEvDTeZ03XHHHSXX7WNlw4cPx4IFC/C3v/1NzYWTxYalrP6WW25R6zNWxxw3u35dY/DpdztwMCsVFzxsy25esJX7lw23fFu4ZQAeWQY0DzCjb5foantvRERUOVw4uYK1JWUJLlmpRE3olnlvJodwkzG2vOIWmwRbhDEYc56LQ7vIP8cHiYjoT1w42cnaR4bjg+eGItqnPnzPesMj1QBDWnE1pCHZdpkKdbvvWS+1H4ONiMh1sOV2mUWU12w9gKUbEovP5wYrNIMGg2ZQVZFSPBLXPVZ1RZoDpXlHRESu0HJjuFVyzcmEA0lqBZPc/EI1j03K/aUqksUjRESuF24uO8/NlUiAdWvTFN3aOPudEBFRZbDZQUREusNwIyIi3WG4ERGR7jDciIhIdxhuRESkOww3IiLSHYYbERHpDsONiIh0h+FGRES6w3AjIiLdYbgREZHuMNyIiEh3GG5ERKQ7DDdye3JKok17jmHJukQsWJOgLuX6jz/9DA8PD/Tv3/+ixzz33HPo2LEjTCYTOnToUKnXSU5OxrBhw9CyZUsYjUY8//zz1XA0RFQVeMobcltyMtnVW/Zj2cZfik8ma7BCTk5oAOClGXF+yyr0HHgfNny3CqdOnUJ4eHipxz/++OPYunUrdu/eXanXs1gsCA0NxSuvvIJ33nmnmo6KiKoCw43c0u7DpzBu3koczT4Hi6kQ1mAN8LIlmwZYcyxI2bcNiO4M/+uiMW3WbMyeNa3k8bNnz1aXp0+frnS4NW3aFO+99576+aOPPqqmIyOiqsBuSXLLYBs9+3MczE9FnrkARfU1aGGA1gDQGhZf5p3cC8+wEBRFBsIa1RIfzp+PXw6ddPZbJ6IawnAjt+uKlBZbkjUdF4I1aGYAoQBCbJe2LWfnLvje0l7d79khEgWWXDzx8iz1eCLSP4YbuRUZY5OuyKK6GhAEINi2mf/cLuSeQeGRk/C9q626z1DPA75t2uDgrh+wZuuBSr2Ov79/yRYfH1/tx0VEVYtjbuRWVZFSPCJjbFodAAG2LRCAXPcB4AHkfLwTKLIi9clZfz5YAwwenlj0zWY80ON6eHpU/HddYmJiyc9169atzsMiomrAcCO3sW1/kqqKVMUjvrYw87MFm4ScH6BpRchb/QvqjusFU2wkkAUgE0A6cO7DxdiZ8D0SDgxAtzZNK3ytqKioGjsuIqp6DDdyG8fTMlS5v6qKNNnCzcch5PyB/G9/gzUrH34P3QCj1QfIkIE6KYsEfGNaI/1AgnqeUNMhnD9/HikpKcjLyytpqbVu3Rre3t7lvgf7fvJYqbSU67K/PI6IXAfDjdxGnqVQzWNT5f5Gh83DtnkBuQt3wdS9OYxBPsUtNof9fFvH4Pzmn3Fw/z7MmfxvbNiwoeS5r7/+enV59OhRVfJfHvt+YseOHVi4cCGaNGmCY8eOVe/BE9EVYbiR2/A1ealcUwlnddiKbFshYP5smLrEedttDvt5hzdCi6feQquY1nh1/fqreg+apuKViFwcw43cRuOwILXyiAovC4B8h026KmFrwUmo5Za5X/YvkN2M6nmISN84FYDcxo0xEYgyh8CYawDybKElISZT17JthSP2Ldt2uz3k8gBjngEtzKHoHB3h7EMhomrGcCO3IeX7cbfFwmTxgsEeaPZQy7Bt6Q4/20MuG2p/k8UTcd1jLzsNgIjcH/+Xk1vp1zUGzQLqwSPL8GeYpdsqIstu9vsyoPZvHmBG3y7Rzj4EIqoBDDdyK+bAOpgZfzcijMHwTDfAYCvzxxnb5enS1+V+2U/2nxk/UD2eiPSP4UZup31kOD54biiiferD96w3PFINMKQBhhTAkGy7TIW63fesl9pvznNxaBfZ0NlvnYhqiEHTeW1zVlYWAgMDkZmZyWWUdEYWQZa1IpduSCw+nxus0AwaDJpBVUVK8YiMsUlXJFtsRLXr+5jhRrpYczLhQJJaeSQ3vxB+Pl6q3F+qIlk8QlQ7v485z43cngSYrBXZrY2z3wkRuQr+WUtERLrDcCMiIt1huBERke4w3IiISHcYbkREpDsMNyIi0h2GGxER6Y5Lh9trr70Gg8FQaouO5sK3RETk5pO427Rpg++++67kuqeny79lIiJyMpdPCgmzBg0aVHp/i8WiNsflXoiIqHZx6W5J8fvvvyM8PBzNmzfHQw89hKSkpAr3nzp1qlq7zL41bty4xt4rERG5BpdeOHnNmjU4f/48WrVqheTkZEyePBknT57Enj17EBAQUOmWmwQcF04mInIunhWgHBkZGWjSpAnefvttjBw5slKP4VkBiIhcQ01+H7t8t6SjoKAgtGzZEocOHXL2WyEiIhfmVuEmXZSHDx9Gw4Y8ozIREblpuI0bNw4bNmzAsWPHsGnTJgwZMgQeHh548MEHnf3WiIjIhbn0VIATJ06oIDt79ixCQ0Nxyy23YMuWLepnIiIitwy3xYsXO/stEBGRG3LpbkkiIqKrwXAjIiLdYbgREZHuMNyIiEh3GG5ERKQ7DDciItIdhhsREekOw42IiHSH4UZERLrDcCMiIt1huBERke4w3IiISHcYbkREpDsMNyIi0h2GGxER6Q7DjYiIdIfhRkREusNwIyIi3WG4ERGR7jDciIhIdxhuRESkOww3IiLSHYYbERHpDsONiIh0h+FGRES6w3AjIiLdYbgREZHuMNyIiEh3GG5ERKQ7DDciItIdhhsREekOw42IiHSH4UZERLrDcCMiIt1huBERke4w3IiISHcYbkREpDsMNyIi0h2GWzkuFFmxac8xLFmXiAVrEtSlXP/xp5/h4eGB/v37X/SYpKQkdbufnx/CwsLw4osv4sKFC5d9raVLlyI6Oho+Pj5o164dVq9eXU1HRURUO3g6+w24mrOZOVi9ZT+WbfwFh86eQaHBCg2AAYCXZsT5LavQc+B92PDdKpw6dQrh4eHqcUVFRSrYGjRogE2bNiE5ORmPPvoovLy88NZbb5X7erLvgw8+iKlTp2LAgAFYuHAhBg8ejJ07d6Jt27Y1eORERPph0DRNvrt1KysrC4GBgcjMzETdunUr3Hf34VMYN28ljmafg8VUCKufBnjZkk0DrDkWpLw/C43vGQ1L4o+4r38PzJ41TT12zZo1Kpwk8OrXr69umzdvHsaPH4/Tp0/D29v7kq95//33IycnB6tWrSq5rWvXrujQoYN6PBFRbfw+vlbslnQIttGzP8fB/FTkmQtQVF+DFgZoDQCtYfFl3sm98AwLQVFkIKxRLfHh/Pn45dBJ9fjNmzerLkV7sInevXurD3Pv3r3lvq48rmfPnqVuk8fJ7UREdHUYbrauSGmxJVnTcSFYg2YGEAogxHZp23J27oLvLe3V/Z4dIlFgycUTL89Sj09JSSkVbMJ+Xe4rT3mPq+gxRESkg3D74IMP0LRpU1Vw0aVLF2zbtq1Kn1/G2KQrsqiuBgQBCLZt5j+3C7lnUHjkJHzvaqvuM9TzgG+bNji46wes2Xrgsq8hxSb+/v4lW0XjcEREpPOCkiVLlmDs2LFq/EmC7d1331XddgcPHlQViVVRFSnFIzLGptUBEGDbAgHIdR8AHkDOxzuBIitSn5z154M1wODhiUXfbMb1YfUvCt3U1FR1KUUmUniSmJhYcl+9evVK7rPv5/g4uZ2IiHTacnv77bcxatQoPPbYY2jdurUKOSm1/+ijj6rk+bftT1JVkap4xNcWZn62YLOFnFanCHmrf0Hdcb0Q+kk8Qt+PR+iUeISOi4fR3x87E75HvUaR+PXXX5GWllby3GvXrlWDpvK+PT09ERUVVbLZw61bt274/vvvS70neZzcTkREOmy5FRQUYMeOHZgwYULJbUajURVglFdwYbFY1GYnBR0VOZ6Wocr9VVWkyRZuPg4h5w/kf/sbrFn58HvoBhitPkCGDNQBOA34xrRG+oEEhEe9okLskUcewfTp09WY2SuvvILRo0fDZJInvrS//OUv6N69O2bNmqWmEixevBjbt2/Hv//976v6NyMiIhdvuZ05c0bNH7uSgguZLyalpvatcePGFb5GnqVQzWNT5f5Gh83DtnkBuQt3wdS9OYxBPsW3Oezn2zoGltMn8ftvB1U5v0zwllbXww8/rOa5vf766xW+/k033aTmtkmYxcbGYtmyZfjyyy85x42ISK8tt6shrTwZo3NsuVUUcL4mL5VrKuGsDluRbSsEzJ8NU5c4b7vNYT/v8EZo8dRbaBXTGk2aNLmq1UXuvfdetRERUS0It5CQENUSupKCC+kCrKgbsKzGYUFq5REVXtKbme+wSVel8LCFWm6Z+2X/AtnNqJ6HiIhcg0t3S8qqHh07dixVcGG1WtX1qiq4uDEmAlHmEBhzDUCeLbQkxHIAZAPIdNiybbfbQy4PMOYZ0MIcis7REVXyfoiISOctNyFdjMOHD0enTp1w4403qqkAslyVVE9WBU8PI+Jui8W+5WnIyymA5uPwryKtOR/bnwBWh+CTkMsGDDmAyeKJuP6x6nmIiMg1uHy4ydqLsjbjxIkTVRGJrLn49ddfX1Rkci36dY3Bp9/twMGsVFzwsC21KYv5510i3PJt4ZYBeGQZ0DzAjL5doqvsvRAR0bXjwsll1paUJbhkpRI1oVvmvZkcwk3G2PKKW2wSbBHGYMx5Lg7tIhvW6DEREbmjLC6cXPPaR4bjg+eGItqnPnzPesMj1QBDGmBIAQzJtstUqNt9z3qp/RhsRESuiS23MmQRZFkrcumGxOLzucEKzaDBoBlUVaQUj8R1j1VdkeZAad4REZGrtdwYbhWsOZlwIEmtYJKbXwg/Hy9V7i9VkSweISJy7XBz+YISZ5EA69amKbq1cfY7ISKiK8UmCBER6Q7DjYiIdIfhRkREusNwIyIi3WG4ERGR7ui+WtI+0+FyJy0lIqLqZf8erokZaLoPt+xsWQgSlz1pKRER1dz3ssx3q066n8Qtp8g5deoUAgICYDCo05JeEfvJTo8fP17tkw5rml6PTa/Hpedj0+txCR7bnyRuJNjCw8NhNFbvqJjuW27yD9ioUaNrfh754PT2i6n3Y9Prcen52PR6XILHVqy6W2x2LCghIiLdYbgREZHuMNwuw2QyYdKkSepSb/R6bHo9Lj0fm16PS/DYnEP3BSVERFT7sOVGRES6w3AjIiLdYbgREZHuMNyIiEh3GG6X8cEHH6Bp06bw8fFBly5dsG3bNriz1157Ta3U4rhFR0fDHW3cuBF33323Wu1AjuPLL78sdb/USk2cOBENGzaEr68vevbsid9//x16OLYRI0Zc9Dn26dMHrm7q1Kno3LmzWjEoLCwMgwcPxsGDB0vtk5+fj9GjR8NsNsPf3x9Dhw5Famoq3P24br/99os+s/j4eLi6uXPnon379iUTtbt164Y1a9a4/OfFcKvAkiVLMHbsWFXqunPnTsTGxqJ3795IS0uDO2vTpg2Sk5NLtp9++gnuKCcnR30m8gfIpUyfPh2zZ8/GvHnzsHXrVtSpU0d9fvKf0dVd7tiEhJnj57ho0SK4ug0bNqgvwi1btmDt2rUoLCxEr1691PHavfDCC1i5ciWWLl2q9pfl8+655x64+3GJUaNGlfrM5HfU1TVq1AjTpk3Djh07sH37dtx5550YNGgQ9u7d69qfl0wFoEu78cYbtdGjR5dcLyoq0sLDw7WpU6dq7mrSpElabGyspjfyq7x8+fKS61arVWvQoIE2Y8aMktsyMjI0k8mkLVq0SHPnYxPDhw/XBg0apLm7tLQ0dXwbNmwo+Yy8vLy0pUuXluyzf/9+tc/mzZs1dz0u0b17d+0vf/mLpgfBwcHaf/7zH5f+vNhyK0dBQYH6S0W6shzXqZTrmzdvhjuTrjnp7mrevDkeeughJCUlQW+OHj2KlJSUUp+frGknXcvu/vnZrV+/XnWBtWrVCk8//TTOnj0Ld5OZmaku69Wrpy7l/5y0ehw/N+k2j4iIcKvPrexx2X322WcICQlB27ZtMWHCBOTm5sKdFBUVYfHixapFKt2Trvx56X7h5Kt15swZ9UHWr1+/1O1y/cCBA3BX8uW+YMEC9YUo3SKTJ0/Grbfeij179qjxAr2QYBOX+vzs97kz6ZKUrp9mzZrh8OHD+Pvf/46+ffuqLxQPDw+4yxk7nn/+edx8883qy17IZ+Pt7Y2goCC3/dwudVxi2LBhaNKkifrDcvfu3Rg/frwal/viiy/g6n799VcVZtKlL+Nqy5cvR+vWrZGYmOiynxfDrZaRL0A7GSSWsJP/cP/73/8wcuRIp743qrwHHnig5Od27dqpzzIyMlK15nr06AF3IGNU8keVu475XulxPfnkk6U+Myl0ks9K/jiRz86VtWrVSgWZtEiXLVuG4cOHq/E1V8ZuyXJI14H8BVy26keuN2jQAHohf3G1bNkShw4dgp7YPyO9f3520sUsv7Pu8jmOGTMGq1atwg8//FDqlFTy2ciQQEZGhlt+buUd16XIH5bCHT4zb29vREVFoWPHjqoyVIqd3nvvPZf+vBhuFXyY8kF+//33pbob5Lo0z/Xi/Pnz6i9H+StST6S7Tv5zOX5+cmJFqZrU0+dnd+LECTXm5uqfo9THSABIt9a6devU5+RI/s95eXmV+tyk607GhV35c7vccV2KtISEq39mlyLfhRaLxbU/L6eWs7i4xYsXq+q6BQsWaPv27dOefPJJLSgoSEtJSdHc1V//+ldt/fr12tGjR7Wff/5Z69mzpxYSEqKqu9xNdna2tmvXLrXJr/Lbb7+tfv7jjz/U/dOmTVOf11dffaXt3r1bVRc2a9ZMy8vL09z52OS+cePGqWo0+Ry/++477YYbbtBatGih5efna67s6aef1gIDA9XvYHJycsmWm5tbsk98fLwWERGhrVu3Ttu+fbvWrVs3tbnzcR06dEh7/fXX1fHIZya/k82bN9duu+02zdW99NJLqupT3rf8P5LrBoNB+/bbb13682K4Xcb777+vPjhvb281NWDLli2aO7v//vu1hg0bquO57rrr1HX5j+eOfvjhB/XFX3aTMnn7dIBXX31Vq1+/vvojpUePHtrBgwc1dz82+cLs1auXFhoaqsqwmzRpoo0aNcot/ui61DHJ9vHHH5fsI398PPPMM6rc3M/PTxsyZIgKCnc+rqSkJBVk9erVU7+LUVFR2osvvqhlZmZqru7xxx9Xv2PynSG/c/L/yB5srvx58ZQ3RESkOxxzIyIi3WG4ERGR7jDciIhIdxhuRESkOww3IiLSHYYbERHpDsONiIh0h+FGRES6w3AjchHz589XZ2+ubl9//TU6dOig1gck0iuGG5ELkPNkvfrqq5g0aVKNnAtOFruVE2cS6RXDjcgFyDmy6tatq05wWRNGjBiB2bNn18hrETkDw42oCp0+fVqdauett94quW3Tpk3qFEqOpwUpa/Hixbj77rtL3Xb77berMzo7Gjx4sAomu6ZNm2LKlCl49NFH1RmS5cSzK1asUO9j0KBB6jY5ken27dtLPY+8ltwmpzsi0iOGG1EVCg0NxUcffYTXXntNhUd2djYeeeQRda6vis6QLWdt7tSp01W95jvvvKNafLt27UL//v3V60nYPfzww9i5c6c6y7Ncd1wjPSIiAvXr18ePP/54Va9J5OoYbkRVrF+/fhg1ahQeeughxMfHo06dOursxeWRsxhnZmYiPDz8ql/vqaeeQosWLTBx4kR1UtbOnTvj3nvvVWdZHz9+PPbv33/RWcnl9f7444+rek0iV8dwI6oGM2fOxIULF7B06VJVuGEymcrdNy8vT136+Phc1WtJt6OdtMZEu3btLrotLS2t1ON8fX2Rm5t7Va9J5OoYbkTVQMayTp06pcrtjx07VuG+ZrMZBoMB6enpl33eoqKii26Tykc7eZ7ybitb+n/u3DnVjUqkRww3oipWUFCgxrvuv/9+vPHGG3jiiScuajU5kmKT1q1bY9++fRfdV7Yr8ciRI1U29UAC+Prrr6+S5yNyNQw3oir28ssvqzE0KbWX8S4Z93r88ccrfEzv3r1VUUlZX331Fb744gsVRG+++aYKQBknO3ny5DW9xy1btqiu0m7dul3T8xC5KoYbURVav3493n33XXz66adq3prRaFQ/S1Xi3Llzy33cyJEjsXr1ahWKjqT6cfr06aplt3HjRsyZMwfbtm1Tz3ktFi1apApe/Pz8rul5iFyVQXOsDyYip5HqxhtuuAETJkwomecmy2RJWFalM2fOoFWrVmqqQrNmzar0uYlcBVtuRC5ixowZatJ1dZMCF2kBMthIz9hyI3JR1dVyI6oNGG5ERKQ77JYkIiLdYbgREZHuMNyIiEh3GG5ERKQ7DDciItIdhhsREekOw42IiHSH4UZERNCb/w/M9bXPmhOOYQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "lattice_spacing = 12\n", + "repeats = 4\n", + "\n", + "qbits1 = qse.lattices.chain(lattice_spacing=lattice_spacing, repeats=repeats)\n", + "qbits1.labels = [f\"A{i}-{i}\" for i in range(repeats)]\n", + "\n", + "qbits2 = qse.lattices.chain(lattice_spacing=lattice_spacing, repeats=repeats)\n", + "qbits2.labels = [f\"B{i}-{i+repeats}\" for i in range(repeats)]\n", + "qbits2.translate((lattice_spacing * 0.5, 8, 0))\n", + "\n", + "qbits_ssh = qbits1 + qbits2\n", + "qbits_ssh.rotate(\n", + " 90 - angle\n", + ") # by rotating to this angle interactions in the A & B chains will cancel.\n", + "\n", + "qbits_ssh.draw(show_labels=True, units=\"µm\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can verify that there are no interactions in the same chain by computing the couplings between qubits." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "ename": "TypeError", + "evalue": "'NoneType' object is not callable", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mTypeError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[4]\u001b[39m\u001b[32m, line 3\u001b[39m\n\u001b[32m 1\u001b[39m magnetic_field = np.array([\u001b[32m0.0\u001b[39m, \u001b[32m1.0\u001b[39m, \u001b[32m0.0\u001b[39m])\n\u001b[32m----> \u001b[39m\u001b[32m3\u001b[39m pcalc_ssh = \u001b[43mqse\u001b[49m\u001b[43m.\u001b[49m\u001b[43mcalc\u001b[49m\u001b[43m.\u001b[49m\u001b[43mMyqlm\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 4\u001b[39m \u001b[43m \u001b[49m\u001b[43mqbits\u001b[49m\u001b[43m=\u001b[49m\u001b[43mqbits_ssh\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 5\u001b[39m \u001b[43m \u001b[49m\u001b[38;5;66;43;03m# amplitude=simple_pulse.amplitude,\u001b[39;49;00m\n\u001b[32m 6\u001b[39m \u001b[43m \u001b[49m\u001b[38;5;66;43;03m# detuning=simple_pulse.detuning,\u001b[39;49;00m\n\u001b[32m 7\u001b[39m \u001b[43m \u001b[49m\u001b[43mchannel\u001b[49m\u001b[43m=\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mmw_global\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[32m 8\u001b[39m \u001b[43m \u001b[49m\u001b[43mmagnetic_field\u001b[49m\u001b[43m=\u001b[49m\u001b[43mmagnetic_field\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 9\u001b[39m \u001b[43m)\u001b[49m\n\u001b[32m 10\u001b[39m pcalc_ssh.build_sequence()\n\u001b[32m 11\u001b[39m pcalc_ssh.calculate()\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/quantum/qse/qse/calc/myqlm.py:123\u001b[39m, in \u001b[36mMyqlm.__init__\u001b[39m\u001b[34m(self, qbits, amplitude, detuning, channel, magnetic_field, qpu, label, wtimes)\u001b[39m\n\u001b[32m 111\u001b[39m installation_message = (\n\u001b[32m 112\u001b[39m \u001b[33m\"\u001b[39m\u001b[33mmyQLM is not installed. To install, \u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 113\u001b[39m \u001b[33m\"\u001b[39m\u001b[33msee https://myqlm.github.io/01_getting_started/:myqlm:01_install.html.\u001b[39m\u001b[33m\"\u001b[39m\n\u001b[32m 114\u001b[39m )\n\u001b[32m 116\u001b[39m \u001b[38;5;28msuper\u001b[39m().\u001b[34m__init__\u001b[39m(\n\u001b[32m 117\u001b[39m qbits=qbits,\n\u001b[32m 118\u001b[39m label=label,\n\u001b[32m 119\u001b[39m is_calculator_available=CALCULATOR_AVAILABLE,\n\u001b[32m 120\u001b[39m installation_message=installation_message,\n\u001b[32m 121\u001b[39m )\n\u001b[32m--> \u001b[39m\u001b[32m123\u001b[39m \u001b[38;5;28mself\u001b[39m.qpu = \u001b[43mAQPU\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mif\u001b[39;00m qpu \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;28;01melse\u001b[39;00m qpu\n\u001b[32m 124\u001b[39m \u001b[38;5;28mself\u001b[39m.wtimes = wtimes\n\u001b[32m 125\u001b[39m \u001b[38;5;28mself\u001b[39m.results = \u001b[38;5;28;01mNone\u001b[39;00m\n", + "\u001b[31mTypeError\u001b[39m: 'NoneType' object is not callable" + ] + } + ], + "source": [ + "magnetic_field = np.array([0.0, 1.0, 0.0])\n", + "\n", + "pcalc_ssh = qse.calc.Myqlm(\n", + " qbits=qbits_ssh,\n", + " # amplitude=simple_pulse.amplitude,\n", + " # detuning=simple_pulse.detuning,\n", + " channel=\"mw_global\",\n", + " magnetic_field=magnetic_field,\n", + ")\n", + "pcalc_ssh.build_sequence()\n", + "pcalc_ssh.calculate()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Version" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "qse.utils.print_environment()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "qse-myqlm", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/qse/calc/myqlm.py b/qse/calc/myqlm.py index 3629543..2600e28 100644 --- a/qse/calc/myqlm.py +++ b/qse/calc/myqlm.py @@ -53,7 +53,7 @@ default_params = { - "rydberg": { + "rydberg_global": { "C6": 5420158.53, # unit (rad/µs)(µm)**6 "min_omega": None, "max_omega": 12.57, # rad/µs from pulser @@ -62,8 +62,9 @@ "min_atom_distance": 4, # µm from pulser "max_duration": 4000, # ns from pulser (max sequence duration), }, - "ssh": {}, - "fermion": {}, + "mw_global": { + "C3": 3700, # unit (rad/µs)(µm)**3 + }, } @@ -79,6 +80,13 @@ class Myqlm(Calculator): The amplitude pulse. detuning : qse.Signal The detuning pulse. + channel : str + Which channel to use. For example "rydberg_global" for Rydberg or + "mw_global" for microwave. + Defaults to "rydberg_global". + magnetic_field : np.ndarray | list + A magnetic field. Must be a 3-component array or list. + Can only be passed when using the Microwave channel ("mw_global"). qpu : qat.qpus The Quantum Processing Unit for executing the job. label: str @@ -94,6 +102,8 @@ def __init__( qbits=None, amplitude=None, detuning=None, + channel="rydberg_global", + magnetic_field=None, qpu=None, label="myqlm-run", wtimes=True, @@ -111,56 +121,30 @@ def __init__( ) self.qpu = AQPU() if qpu is None else qpu - self.label = label self.wtimes = wtimes self.results = None - self.system = "rydberg" - self.params = dict(default_params[self.system]) + self.channel = channel + self.magnetic_field = magnetic_field + self.params = dict(default_params[self.channel]) self.amplitude = amplitude self.detuning = detuning - self.C6 = self.params["C6"] - - @property - def Hamiltonian(self): - return self._get_hamiltonian() - - def _get_hamiltonian(self): - if self.system == "rydberg": - ham = self._generate_rydberg_hamiltonian() - else: - ham = None - return ham - - def _generate_rydberg_hamiltonian(self): - nqbits = self.qbits.nqbits - # > add checks for ensuring whether amplitudes/detunings etc - # > are within allowed limits compatible with pulser virtual device - H_amplitude = qat.core.Observable( - nqbits, pauli_terms=[qat.core.Term(0.5, "X", [i]) for i in range(nqbits)] - ) - amplitude = _waveform(self.amplitude) - - H_detuning = qat.core.Observable( - nqbits, pauli_terms=[qat.core.Term(0.5, "Z", [i]) for i in range(nqbits)] - ) - detuning = _waveform(self.detuning) - - rij = self.qbits.get_all_distances() - H_interact = 0 - for i in range(nqbits): - for j in range(i + 1, nqbits): - H_interact += ( - (self.C6 / rij[i, j] ** 6) * _occ_op(nqbits, i) * _occ_op(nqbits, j) - ) - - return [ - (amplitude, H_amplitude), - (detuning, H_detuning), - (1, H_interact), - ] + def get_hamiltonian(self): + if self.system == "rydberg_global": + return self._generate_rydberg_hamiltonian( + self.qbits, self.amplitude, self.detuning, self.params["C6"] + ) + elif self.system == "mw_global": + return self._generate_microwave_hamiltonian( + self.qbits, + self.amplitude, + self.detuning, + self.params["C3"], + self.magnetic_field, + ) + return None def calculate(self): """ @@ -172,7 +156,7 @@ def calculate(self): tmax = ( max(self.amplitude.duration, self.detuning.duration) / 1000 ) # convert to microseconds. - self.schedule = qat.core.Schedule(drive=self.Hamiltonian, tmax=tmax) + self.schedule = qat.core.Schedule(drive=self.get_hamiltonian(), tmax=tmax) self.job = self.schedule.to_job() self.async_result = self.qpu.submit(self.job) self.results = self.async_result.join() @@ -225,3 +209,87 @@ def _waveform(pulse): def _occ_op(nqbits, qi): ti = qat.core.Term(1.0, "Z", [qi]) return (1 - qat.core.Observable(nqbits, pauli_terms=[ti])) / 2 + + +def _plus_op(nqbits, qi): + return ( + qat.core.Observable(nqbits, pauli_terms=[qat.core.Term(1.0, "X", [qi])]) + - 1j * qat.core.Observable(nqbits, pauli_terms=[qat.core.Term(1.0, "Y", [qi])]) + ) / 2 + + +def _minus_op(nqbits, qi): + return ( + qat.core.Observable(nqbits, pauli_terms=[qat.core.Term(1.0, "X", [qi])]) + + 1j * qat.core.Observable(nqbits, pauli_terms=[qat.core.Term(1.0, "Y", [qi])]) + ) / 2 + + +def _hopping_op(nqbits, qi, qj): + return _plus_op(nqbits, qi) * _minus_op(nqbits, qj) + _plus_op( + nqbits, qj + ) * _minus_op(nqbits, qi) + + +def _generate_rydberg_hamiltonian(qbits, amplitude, detuning, c6): + nqbits = qbits.nqbits + H_amplitude = qat.core.Observable( + nqbits, pauli_terms=[qat.core.Term(0.5, "X", [i]) for i in range(nqbits)] + ) + amplitude = _waveform(amplitude) + + H_detuning = qat.core.Observable( + nqbits, pauli_terms=[qat.core.Term(0.5, "Z", [i]) for i in range(nqbits)] + ) + detuning = _waveform(detuning) + + rij = qbits.get_all_distances() + H_interact = 0 + for i in range(nqbits): + for j in range(i + 1, nqbits): + H_interact += ( + (c6 / rij[i, j] ** 6) * _occ_op(nqbits, i) * _occ_op(nqbits, j) + ) + + return [ + (amplitude, H_amplitude), + (detuning, H_detuning), + (1, H_interact), + ] + + +def _generate_microwave_hamiltonian(qbits, amplitude, detuning, c3, magnetic_field): + nqbits = qbits.nqbits + H_amplitude = qat.core.Observable( + nqbits, pauli_terms=[qat.core.Term(0.5, "X", [i]) for i in range(nqbits)] + ) + amplitude = _waveform(amplitude) + + H_detuning = qat.core.Observable( + nqbits, pauli_terms=[qat.core.Term(0.5, "Z", [i]) for i in range(nqbits)] + ) + detuning = _waveform(detuning) + + rij = qbits.get_all_distances() + + if magnetic_field is None: + H_interact = 0 + for i in range(nqbits): + for j in range(i + 1, nqbits): + H_interact += (c3 / rij[i, j] ** 3) * _hopping_op(nqbits, i, j) + else: + H_interact = 0 + magnetic_field = np.array(magnetic_field) + magnetic_field /= np.linalg.norm(magnetic_field) # normalize + for i in range(nqbits): + for j in range(i + 1, nqbits): + qubit_vec = qbits.positions[i] - qbits.positions[j] + qubit_vec /= np.linalg.norm(qubit_vec) + m_term = 1 - 3 * np.dot(magnetic_field, qubit_vec) ** 2 + H_interact += (m_term * c3 / rij[i, j] ** 3) * _hopping_op(nqbits, i, j) + + return [ + (amplitude, H_amplitude), + (detuning, H_detuning), + (1, H_interact), + ]