From e99749183992e7c92f23a03f8a480a26b5d85f24 Mon Sep 17 00:00:00 2001 From: bartandrews Date: Thu, 7 Nov 2024 14:43:15 +0100 Subject: [PATCH] apply UCJ operator directly --- docs/how-to-guides/lucj_mps.ipynb | 29 ++--- python/ffsim/tenpy/__init__.py | 12 +- python/ffsim/tenpy/circuits/gates.py | 117 +++++++++++++++++++- python/ffsim/tenpy/circuits/lucj_circuit.py | 99 +++-------------- tests/python/tenpy/lucj_circuit_test.py | 4 +- 5 files changed, 147 insertions(+), 114 deletions(-) diff --git a/docs/how-to-guides/lucj_mps.ipynb b/docs/how-to-guides/lucj_mps.ipynb index 3b75d26b9..593d9c57b 100644 --- a/docs/how-to-guides/lucj_mps.ipynb +++ b/docs/how-to-guides/lucj_mps.ipynb @@ -34,9 +34,9 @@ "output_type": "stream", "text": [ "converged SCF energy = -77.8266321248745\n", - "Parsing /tmp/tmpibrj_a18\n", + "Parsing /tmp/tmp1e4pf6nc\n", "converged SCF energy = -77.8266321248744\n", - "CASCI E = -77.8742165643862 E(CI) = -4.02122442107773 S^2 = 0.0000000\n", + "CASCI E = -77.8742165643863 E(CI) = -4.02122442107772 S^2 = 0.0000000\n", "norb = 4\n", "nelec = (2, 2)\n" ] @@ -45,7 +45,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Overwritten attributes get_ovlp get_hcore of \n", + "Overwritten attributes get_hcore get_ovlp of \n", "/home/bart/PycharmProjects/ffsim/.ffsim_dev/lib/python3.12/site-packages/pyscf/gto/mole.py:1294: UserWarning: Function mol.dumps drops attribute energy_nuc because it is not JSON-serializable\n", " warnings.warn(msg)\n", "/home/bart/PycharmProjects/ffsim/.ffsim_dev/lib/python3.12/site-packages/pyscf/gto/mole.py:1294: UserWarning: Function mol.dumps drops attribute intor_symmetric because it is not JSON-serializable\n", @@ -120,7 +120,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "E(CCSD) = -77.87421536374032 E_corr = -0.04758323886585264\n" + "E(CCSD) = -77.87421536374035 E_corr = -0.04758323886585367\n" ] }, { @@ -184,7 +184,7 @@ "text": [ "original Hamiltonian type = \n", "converted Hamiltonian type = \n", - "maximum MPO bond dimension = 54\n" + "maximum MPO bond dimension = 58\n" ] } ], @@ -336,9 +336,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "LUCJ (MPS) energy = -77.77107802937196\n", + "LUCJ (MPS) energy = -77.77102667787551\n", "LUCJ energy = -77.84651018653345\n", - "FCI energy = -77.87421656438624\n" + "FCI energy = -77.87421656438627\n" ] } ], @@ -380,21 +380,10 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "bf98d538-c182-4ede-917f-1eed31969c9a", "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2IAAAFzCAYAAABcurqFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC2xElEQVR4nOzdeVxVdf748ddd4LIJyg6CgooIgqAiaKVpOTmkZrvtambLZMvPakZnSh2/LtO0jGXO2JRmi5XVjNa02GImVgoJQiJoLoAigqCyr/fe8/vjyk1iEfDCYXk/H4/zgHPu53zO+4Jy7vucz3l/NIqiKAghhBBCCCGE6DRatQMQQgghhBBCiN5GEjEhhBBCCCGE6GSSiAkhhBBCCCFEJ5NETAghhBBCCCE6mSRiQgghhBBCCNHJJBETQgghhBBCiE4miZgQQgghhBBCdDJJxIQQQgghhBCik+nVDqAnMJvN5OXl0adPHzQajdrhCCFEt6IoCmVlZfj7+6PVyvXBrkDOa0II0X6tPa9JImYDeXl5BAYGqh2GEEJ0aydOnCAgIEDtMARyXhNCCFu42HlNEjEb6NOnD2D5Ybu6uqocjRBCdC+lpaUEBgZa/5YK9cl5TQgh2q+15zVJxGygftiGq6urnLCE6OGqKsooeGEcAD5P7MbRWZIHW5EhcF2HnNeE6D3kvNZxLnZek0RMCCHaQFHMBJlPAFCpmFWORgghhLg0cl5TjzwVLYQQQgghhBCdTBIxIYQQQgghhOhkMjRRCNGpTCYTdXV1aofRbjU1tWhdAn/9Xl+tckRdn06nQ6/XyzNgnezTTz/liSeewGw286c//Yn77rtP7ZCEEEJcQBIxIUSnKS8vJzc3F0VR1A6l3cxmM9rLX7B8n38arbZI5Yi6BycnJ/z8/LC3t1c7lF7BaDSyYMECduzYgZubG6NHj+aGG27Aw8ND7dCEEEKcJ4mYEKJTmEwmcnNzcXJywsvLq9veHTGZjOjO1Fi+9xiITid/RluiKAq1tbUUFhaSlZVFSEiITNrcCZKSkhg+fDj9+/cHID4+nq+++orbb79d5ciEEELUk08QQohOUVdXh6IoeHl54ejoqHY47WY2mTDq7QCwd3BEq9OpHFHX5+joiJ2dHTk5OdTW1uLg4KB2SF1eQkICzz33HMnJyZw6dYotW7Zw/fXXN2izdu1annvuOfLz84mKimLNmjXExsYClgmZ65MwgP79+3Py5MnOfAtCiG5Co9FyCi8A+mrkQllnkp+2EKJTddc7YfW0Oh32/pHY+0dKEtYGchesbSoqKoiKimLt2rVNvr5582YWLFjAkiVLSElJISoqiilTpnD69OlOjlQI0d05OvfBb+kR/JYekTnEOlm3ODN+9913aDSaJpeffvoJgKVLlzb5urOzc4t9Hz9+nKlTp+Lk5IS3tzdPPfUURqOxM94WAAW5R0n/4X8U5B7ttGMKIYTo2uLj41m+fDk33HBDk6+/+OKLzJs3jzlz5hAeHs66detwcnJiw4YNAPj7+ze4A3by5En8/f2bPV5NTQ2lpaUNlvaS85oQQrROt0jELrvsMk6dOtVgue+++wgODiYmJgaAJ598slGb8PBwbrnllmb7NZlMTJ06ldraWn788UfefPNNNm7cyOLFizvlfSV99A+8XhtNxNd34fnaaJL+s7pTjiuEEKL7qq2tJTk5mcmTJ1u3abVaJk+ezO7duwGIjY0lPT2dkydPUl5ezhdffMGUKVOa7XPVqlW4ublZl8DAwHbFlvSf1Q3Oa3s2/71d/QghRG/QLZ4Rs7e3x9fX17peV1fHxx9/zCOPPGId5uTi4oKLi4u1TVpaGhkZGaxbt67Zfr/66isyMjL45ptv8PHxITo6mv/7v//jT3/6E0uXLu3Q6l4FuUeJ2f9XtBpL9TidRmHUz3+lIG46PgGDO+y4QohLYzaZqC04BIC9T6gMTxSdrqioCJPJhI+PT4PtPj4+HDx4EAC9Xs8LL7zApEmTMJvN/PGPf2yxYuKiRYtYsGCBdb20tLTNyVhB7lFG/7y0wXltbOYKzi19iSKdL+UOvtQ4+4NbIPYeA3DxDsbdfxAe3v3RyNBVIVRTXVnOiRcnAhC44DscnFxa3kHYTLdIxH7rk08+4cyZM8yZM6fZNq+//jpDhw5l/PjxzbbZvXs3kZGRDU5mU6ZM4aGHHuLAgQOMHDnSpnFfqDAnAx9NwxLeeo2ZopyDkogJ0YUpKDhwvmoi3bcMv+j5rrvuOq677rpWtTUYDBgMhks6XlPnNYB+lNPPdAQqjkAFcBo4/OvrNYodp7VelNh7U+noj8k1AF2/QJw8B9LXbzCe/YNxcGz5MYOurCD3KIU5GXgNDJfzu+iSzGYTIUbLf8pKs0nlaHqXbpmIrV+/nilTphAQENDk69XV1WzatImFCxe22E9+fn6TVxTrX2tOTU0NNTU11vX2jKX3GhiOSdGgu+CkZVS0eA4c1ua+hBCisxUXFzN58mSMRiNGo5HHHnuMefPmqR1Wr+Dp6YlOp6OgoKDB9oKCggajRzpbU+c1k6LlwKQNmI3V1JzJwVx8AvvykzhX5eNuLMBTOYdBU0egkkdgTR7UpEIxcLxh30X05azem3IHX2qd+4NbAPYeA3H1DcbDfzB9PXwa3VXriARIMZupqa6kpqqC2poqaqsrqK2upK66EmNtFaaaSow1lZjqqjDXVqHJSmB0ydf4aMCsaNgz7CnG3v4Xm8QihOj+VE3EFi5cyLPPPttim8zMTIYN+zU5yc3N5csvv+SDDz5odp8tW7ZQVlbGrFmzbBbrhVatWsVf//rXS+rDJ2AwSSOWMvrnpeg0CooCe4c9yVi5WiaE6Ab69OlDQkICTk5OVFRUEBERwY033igTBncCe3t7Ro8ezfbt260l7c1mM9u3b2f+/PmqxVV/Xhv181/Ra8wYFS0pI5YQO7HpgiMAtTXVFOVlc+7UMSoLszGeO46uNBeHylO41RbgZTqNk6YGT4rxNBZD+S9QDjTMQalS7CnUeVFs70u1ox92NWeJqtiNj0bBrGhI6jsFxX8kSl0VSl01GKvRnF+0phq0JstXvakavbkGvVKLnbkGO6UGe6UWe2oxKLU4aOpwANo0AcP5QrFajcLYQ3/nzNJ15BkGU+EWgtY7DNcBkfgPHYVrX/m/I0Rvo2oi9sQTTzB79uwW2wwaNKjB+htvvIGHh0eLwy1ef/11pk2b1uhu12/5+vqSlJTUYFv9FcaWriraYiw9QOxNj1MQNx3T61Pw1xSidZCSoUJ0VQsXLuQf//gHN954I++98Ee1w1GdTqfDyckJsIwSUBQFRZGhmrZSXl7OkSNHrOtZWVmkpqbi7u7OgAEDWLBgAbNmzSImJobY2FhWr15NRUVFi0P2O0P9ea0o5yCeA4cRe5GLi/YGB/yDh+Ef3PRoEMVspvjsaYpOHqWsIIuaMzlQfAL7ijxcqk/hbjyNJ8U4amoZYD7JgOqTUH1+5wsSoNiSbVCy7dLe3G9m3jAqWmqwp1ZjTw321GnsqdUaMGrsMWoN6M3VhBp/adSNB6V41OyD0/sswzTTgc+hAA8KHIKpdAtB5xOGW1AU/UOice7T99LiFkJ0WaomYl5eXnh5ebW6vaIovPHGG9xzzz3Y2dk12SYrK4sdO3bwySefXLS/cePGsWLFCk6fPo23tzcAX3/9Na6uroSHhze7ny3G0tfzCRjM7oE34Z+zDodDW4FHbdKvEMK2Fi1aREBAAI888gj/N/82hgQPUDukZrVmMmBoeULg1iguLubKK6/k8OHDPPfcc3h6etrwXfRue/fuZdKkSdb1+ot/s2bNYuPGjcycOZPCwkIWL15Mfn4+0dHRbNu27aIXIDuDT8Bgmw0F1Gi19PX0pa+nL3B5k21qqispysviXN5RKgtz0GT/wJiSLxq1y7CLoNLRD7PeAUVnwKx3AL0j2DmisXNAo3dAa++Izt4RncERnb0TentH9AYn7BycsXdwxM7ghMHRGQdHZ+zsDeiB5p5eK8g9ium10Y2Gav58xVqMZYWYCjJxKvkF3+osvDmLD2fwqT4D1Xstd/x+tuxzCi9OOwZT2Xcoep9w+gaNICAkWuZ7EqIH6FbPiH377bdkZWVx3333Ndtmw4YN+Pn5ER8f3+i1LVu2sGjRImtVqWuuuYbw8HDuvvtu/v73v5Ofn8/TTz/Nww8/bLNEqzUCx98NOesIr9rHmYJcPHyafvZNCKEeNzc35s6dy2OPPcb+g0e6dCJWPxnwvffey4033thkm/oJgdetW0dcXByrV69mypQpHDp0yHphKjo6usl5Fb/66iv8/f3p27cvaWlpFBQUcOONN3LzzTd3iUSgJ5g4ceJF7zDOnz9f1aGIXYXBwYn+g4bTf9BwAApyr8X02rZGz2B7zHqL8E4c/t/sUM3f3dGobcm5Ik4d3kdJzs8opzNxKT2Mb002nhTjRyF+VYVQlQSngFTL82a5Wh8KHQdR3TcEO99w+gVH039IZLOFTaRoiBBdT7dKxNavX89ll13W4JmxC5nNZjZu3Mjs2bPRNVFSuqSkhEOHDlnXdTodn376KQ899BDjxo3D2dmZWbNmsWzZsg57D00JGBLBL/qhDDX+wpHvNuEx80+denwhuptTJVVkFVUQ7OmMn5tjpx3XaDTi5OTEzwePMT3+d78dqdRlxMfHN3kx6kIXTggMsG7dOj777DM2bNhgLXSUmpraquP5+PgQFRXFrl27uPnmmy8pdiEuVbMJkArJR2uHarr188Qt9ncQ+7sG24uL8sk7vI+yE/vhfILmX5tDP00pAUo+AZX5UPkj5AEpYFI0nND6UeQ0iOp+Q7H3G457cBSn03cSc2A5PhoFk6IhacRSYm96vON/AE2QhLBrOocrAJ13G0IAaBQZ1H/JSktLcXNzo6SkBFdX13b1sWfTXxl7+EUy7YYT9pcfbRyhEOqrrq4mKyuL4OBgHBwcUBSFqrq2l8n9T3IuSz45gFkBrQb+et1wbhrdtrvIjnY66xyEbfHYY4/x8ssvc+utt7J58+Y2799WK1euZOXKlS22ycjIYMCA5u/OaTSaRkMTa2trcXJy4qOPPmqwfdasWRQXF/Pxxx9fNLaCggKcnJzo06cPJSUlXH755bz33ntERkY22f63v/8L2eJvqLCtnvA7Kcg9ak2AetoH/jMFuZw6so/y4/vRFB3EtfQI/euycaWiyfaKAhf+yVMUSHOMw2TnjKLRoWj1KBotilYP59ctX3Wg1aM5v02j1aNodWjqt+ss7TQ6vWVda/leo7VDo9OjPb+u1VleKzvwJXGn3kV7PiFMVjEhFKIjtfZvaLe6I9aTDZp4D+Zf/kFY3QHyTxzBN3CI2iEJ0aGq6kyEL/7ykvowK/DMxwd45uMDbdovY9kUnOzb9ucvOTmZdevWMXXqVNLT09u0b3s9+OCD3HrrrS228ff3b3O/rZkQ+GJycnK4//77rUU6HnnkkWaTMCHUYMtn1boaD58Ay2MMl0+3blPMZgrzj5N/ZB8VueloCzNxKzvKgLpjOGpqG+yv0UB0deKvhU060/mEUKdRGPPzEpKPbafObzTOA6LwGxqDp2/bi58J0V1JItZFePcPJsMQQXjtfrJ3vo3vXZdWHl8IYTtms5kHHniA+fPnExcXx1133UVdXV2zRYOakpeXx1NPPcWmTZtavY+7uzvu7u7tCbnDxcbGtnroohCi42m0Wrz8g/DyDwJ+nTIg//hh7NePafDMnFnRkBj8IBp7FzAbwWwCswnFbERjNqIoJjAb0ZhNoJjQmI3nv5rQnF/XKKYmF+355cLvtZhxMJXTn9MNY9bA6IoEOJIAR4BvLXPGnTIMoqJvKDq/SNwHjyYgZAQGB6fO+UEK0YkkEetCyobMgIz9eGb9D5BETPRsjnY6MpZNadM++SXVTH5xJ+YLBlRrNfDNgivxdWv9zD6Odo2fIW3JmjVrKCoqYtmyZWRnZVFXV0fqri8YfeVUtE08j9oUf3//NiVhYJuhiU3pqhMCCyFsz3dASJPPzI3r5CGBTVeR1JDkdzuG8pN4VRymv/kUnppiPGtSoCAFCt6DVKhTdGTpAjjjEoLRMxynwCj8QmPw9B3QaCJv0XbVleUcXf17AAY/vg0HJxeVI+o9JBHrQkIm3oHxwEqGmI5y4nAagSFRaockRIfRaDRtHh44yMuFVTdG8uf/pmNSFHQaDStvjGCQV8edNE6ePMkzzzzDe++9h7OzM0NChmAw2HPk4AFGXXkt2dnZzJgxg4iICJKSkpg8eTJTpkxh1apVVFRUsGXLFkJCQsjOzubmm2/mo48+YsaMGURHR5OUlMSIESN4//33m3xmraOGJnbVCYGFEB2jrfO7dYTmiqhcmBBWVZRx4lAKxVkpKPnpuJQcIrD2GK6aCoLNOQSX5kDpN3AM2Ann6EOe/SDK+g5D6xdBv+CRBA4dKYlEG5nNJobX7geg0tz2Z7dF+0ki1oW4e/fnZ8fRjKj+idxdmyQRE6IJM8cMYMJQL7KLKgnydOrwqomPPvoo8fHxTJ06FQC9Xk/YkGDSDx6hPkXKzMzkgw8+YMiQIURERODi4kJiYiKvvvoqr7zyCi+99FKDPjMzM3nvvfcICwtj0qRJfP/994wfP77Rsds7NPFikwEDXXZCYCFEx+gKz8xdLCF0dO7D0FFXwqgrrdsUs5n8k8coOLyXyhM/Y1+UiWfFYQJMufTTlNGvNg1Op8HpzZBmucuWowugyHkIdR7hOASOwHdoDD79BzW4eybVG0VXIIlYF1Mz7HpI/Qn/3M9QzH+TW+5CNMHPzbFTytZ/+umnfPvtt2RmZjbYHjlsCOmHjlrXQ0NDCQ0NBSAsLIzJkydb2kVG8vnnnzfqNzQ01Dpp/MiRI8nOzm4yEWuvi00GDHTpCYGFED1XWxNCjVaLb+CQ80XMbrNur64sJ+twKmePpWA+Zbl71r/mKP00ZQw0n2Bg2Qko2wHZwC4oxZlc+0GUuQ5FW1fO6OKvukQ5f9G7SSLWxQybdAc1+5Yy0JzL0QNJDI4cq3ZIQvRa06ZN49y5c422v/Xy/wFQP4DjwgngtVqtdV2r1WIyNR7mcWF7nU7XZJtL0ZrJgEEmBBZCdF8OTi4MiboCoq6wbquvHHnq0F4qT6ShL8rAo9xy98xVU0F47X4osgzB+231xszM96hwCsDo5AN9fNH39cPJPQBXr0DcfQNxcnFT4V2Knk4SsS6mj5s7KS5jGVWxi9O7N0kiJoQQQgjRCg0rR/46uXxNdSVHD//M2aPJaA9/yeiKnQ3300CY8SCUHoRSIL9x32WKI2d1HpTrPahy8LImbHZ9/XH0CMDNOxB3nwE4OvdpU8wyRLJ3k0SsC1KG3wRJuxh4ahuK2SzDE4UQQggh2sng4MTgyLEMjhxLQe7vm6ze+FPYH8FYA2X52FWexrGmENe6ItzNZ3HS1NBHU0Ufcy7U5kItzSZspThxTutBmZ0H1QYv6px90JxP2Jw8AnD1CrAmbEn/Wc3on5fKEMleTKO0ZvyKaFFrZ89uraqKMsx/H4KzppqDU//DsDGTbRClEOqqrq4mKyuL4OBgHBxaX2q+qzGZjJB/fkJn3wh0Orme1Rot/f5t/TdUXDr5nYieLOk/qxtVb2wuAVLMZsrLijmbf5yywuNUn83DWJJnSdiqTuNYXYib0ZKw/Xbi7JaUKo70oYoLC+YaFS1n5u3t9DtjleUl8FyIZeWpwzIM0wZa+zdUPkF0QY7OfdjrdgUxpd9QnPQ+SCImRJeh0+mhf7TaYQghhGintpTz12i19HFzp4+bO4RGN9tOMZspLT3HufzjlBWdoPrsSYzFeVBegF1lAU41hfQxnsHjfMLmqqlq1IdeY6Yo52CnJ2JOLm7w19MXbyhsThKxLsou+hZI+IYhhV9jMhrR6eVXJYQQQghhC7Yu56/RanHt64FrXw9gZLPtFLOZkpKz5GbsJuyru9FeMETSqGjxHDjMZjGJrk8ePuqiwi6/nhKc8aSYzD2Ny18LIYQQQojuRaPV4tbPk+GXT2fviKWYz+dhZgW+HrxQCnb0MpKIdVH2BgcOuV8FQGXyZpWjEULUM5tNVJ76hcpTv2A227bsvBBCiN4j9qbHKb/sjwAkmEfwpf01qsRRXVVB2rO/I+3Z31FdVaFKDL2VJGJdmNPomQCEnttBbU21ytEIIQAURcFJqcBJqWjVXF1CCCFEc1yDYwDw1pTw+f58zpTXdHoMZpORqKokoqqSMJuMnX783kwSsS4sLC6eQvrhRgUZ329VOxwhhBBCCGFL/YIACNKdptZk4oO9uerGIzqVJGJdmE6v56j37wAwpn2ocjRCCCGEEMKm+g4AwEmpoi/lbErMwWSW0Ra9hSRiXVzf2NsACC/ZRVVFmcrRCCGEEEIIm7FzBBcfAMIczpF7roqdv0gp+d5CErEuLnTUJPI03jhpasjY+YHa4QghhBBCCFs6Pzzx5kF1ALyz57iKwYjOJIlYF6fRasnx+z0A2gP/VTkaIYQQQghhU30HAjDJxzLJ845DpzlxtlLNiEQnkUSsG/C57C4AwssTKS0+o3I0Qoierri4mJiYGKKjo4mIiOC1115TOyQhhOi5+lkSMffaU4wP8URRYFOi3BXrDSQR6waCw8eQow3EoKnj0HfvqR2OEL2aTqcH/5HgP9LyfQ/Up08fEhISSE1NJTExkZUrV3LmjFwEEkKIDnF+aCLnsrlrrCUp+2DvCWqMnTNXpZOLGywtgaUllu9Fp5FErBvQaLXkBU4FwHBwi8rRCCE605kzZ/D29iY7O7vTjqnT6XBycgKgpqYGRVEazJl222238cILL3RaPEII0aOdH5pIcQ5XD/PGz82BsxW1fLE/X924RIeTRKybCBh/fnhiVQpnT59UORohRGdZsWIFM2bMICgoCICEhASmT5+Ov78/Go2GrVu3Nrnf2rVrCQoKwsHBgbi4OJKSktp03OLiYqKioggICOCpp57C09PT+trTTz/NihUrKCkpae/bEkIIUa9ffSJ2Ar1G4Y5YS0n7t/fkqBiU6AySiHUTgUMiOawbgl5j5vB3m9QOR4hey2w2UZF/mIr8w5jNHTtspLKykvXr1zN37lzrtoqKCqKioli7dm2z+23evJkFCxawZMkSUlJSiIqKYsqUKZw+/WtJ5Prnv3675OXlAdC3b1/S0tLIysri3XffpaCgwLpvREQEgwcP5p133umAdy2EEL2Ma3/Q6sFcB6V5zIwNRK/VkJxzjgN5HX/Bq7qqgpTnp5Py/HSqqyo6/HjiV5KIdSNngqcD0OfwxypHIkTvcsstt+Dl5cW///1vFEXB2VxO+t7dODg48tVXX3XYcT///HMMBgNjx461bouPj2f58uXccMMNze734osvMm/ePObMmUN4eDjr1q3DycmJDRs2WNukpqaSnp7eaPH392/Ql4+PD1FRUezatavB9unTp/P+++/b6J0KIUQvptWBW6Dl++IcvPs4MCXCF+icUvZmk5FR5QmMKk/AbDJ2+PHEryQR60aCrjw/PLEunYLcoypHI0Tv8fLLL3PTTTexbNkyAMorKrnrkad58MEHuOaaazrsuLt27WL06NFt2qe2tpbk5GQmT55s3abVapk8eTK7d+9uVR8FBQWUlVkmkC8pKSEhIYHQ0NAGbWJjY0lKSqKmpqZN8QkhhGiCtWCHZTji3eeLdnycepLS6jqVghIdTRKxbsQ3cAgZdhEAZH33tsrRCGEblbXGZpfqOpPN27aHn58fjz/+OCdPnuTMmTM8+szfMRjsWbVqVbvfd2vk5OQ0ukN1MUVFRZhMJnx8fBps9/HxIT+/dQ9+5+TkMH78eKKiohg/fjyPPPIIkZGRDdr4+/tTW1vb6j5F5zpx4gQTJ04kPDycESNG8OGHH6odkhCiJfXPiZ3LBiAu2J0Qbxcqa01sSZHaAD1Vz6y93IOVDbkOMtPxyP4UWKp2OEJcsvDFXzb72qRQL96YE2tdH/1/31BV1/RzWXHB7mx+YJx1/Ypnd3C2orZRu+y/TW1XnEOHDsXJyYmlS5eyacsXJH36Ng4ODu3qq7Wqqqo6/BhNiY2NJTU1tcU2jo6OgOU5NtH16PV6Vq9eTXR0NPn5+YwePZprr70WZ2dntUMTQjTlgsqJABqNhrvGDmTJJwd4e08O94wbiEajUTFA0RHkjlg3EzLpLoyKlhDjYU4c2a92OEL0GlqtlsjISP71r3Us/+MfiBo+tMOP6enpyblz59q8j06na1BcAyzDDX19fW0W29mzZwHw8vKyWZ/Cdvz8/IiOjgbA19cXT09P6+9MCNEFXTCXWL0bRvXHyV7HkdPlJGbJ/9+eSO6IdTPu3v352XEUI6r3krvrHQKHPKt2SEJckoxlU5p9Tfubq3/Jz0xupmXjtt//adKlBfYb9fNojRo1kiceuNu6PTs7mxkzZhAREUFSUhKTJ09mypQprFq1ioqKCrZs2UJISAgA06ZN49SpU9TU1LBo0SLuvPNOdu/ezWOPPcaPP/7ImTNnuOKKK9i1axe+vr6MHDmyzZUJ7e3tGT16NNu3b+f6668HwGw2s337dubPn2+bHwaQnp5OQEBAg7L2ovUSEhJ47rnnSE5O5tSpU2zZssX6+6q3du1annvuOfLz84mKimLNmjXExsY23WELkpOTMZlMBAYG2ih6IYTNWYcm/lqy3tXBjutH9ufdxOO8vSeHsYM8VApOdBRJxLqh6tDrIW0v/ic+QzGvQqOVG5ui+3Kyb/2foY5q2xqrV68mMTGR6OgotL/5P5eZmckHH3zAkCFDiIiIwMXFhcTERF599VVeeeUVXnrpJQDeeust3N3dqaioYMyYMdx8882MGzeOCRMm8Oyzz7Jv3z4WL15svXM1ZcoUFi1axLlz5+jXrx8A5eXlHDlyxHrsrKwsUlNTcXd3Z8AAy9wzCxYsYNasWcTExBAbG8vq1aupqKhgzpw5Nvt57Nq1q0MLlfR09dMQ3Hvvvdx4442NXq+fgmDdunXExcWxevVqpkyZwqFDh/D29gYsUxAYjY2fe/zqq6+szxaePXuWe+65h9dee61j35AQ4tL0DbJ8Lc+Huiqwswz/vituIO8mHufL9HxOl1bj7dr5w9VFB1LEJSspKVEApaSkpHOOd65IqVnsrihLXJWj+/d0yjGFuFRVVVVKRkaGUlVVpXYobfbzzz8rBoNB+cMf/qDY29sr1dVVitFYp5jNZiUrK0uJiIiwtr3hhhuUbdu2KYqiKD/88INy3XXXWV97+umnlREjRigjRoxQnJ2dlV9++UVRFMvPJjQ0VJk6dWqjY8fGxirr1q2zru/YsUMBGi2zZs1qsN+aNWuUAQMGKPb29kpsbKyyZ4/t/lZUVVUpbm5uyu7du9u0T3O//87+G9rVAMqWLVsabIuNjVUefvhh67rJZFL8/f2VVatWtbrf6upqZfz48cpbb73VqrYlJSXW5cSJE736dyJEpzObFWVFf0VZ4qoopw82eOnGf/6gDPzTp8pL3/zSMYc2mZSKsmKloqxYMZtMHXKM3qa15zW5ldINufb14IBzHAAFP8rkzkJ0pOrqau644w5mzpzJ8uXLqa2t5fDhI+h0euuD0waDwdpeq9Va17VaLSaTpbjIjh07+OGHH0hMTCQtLY1hw4ZZS7+fPn2a2tpaa8XDCy1evJiXXnoJs9kMwMSJE1EUpdGycePGBvvNnz+fnJwcampqSExMJC4uzmY/kzfeeIPY2NgG85sJ27HFFASKojB79myuuuoq7r777ou2X7VqFW5ubtZFhjEK0ck0miaHJ8KvpezfSzqO0WS2/aG1Wpxc3HBycZNRVp1MftrdlDniJgAGntqGYrb9f0ohhMXChQupqKjglVdeoV+/fgwcOJDVq1eTl5fXpn5KS0vx8PDAwcGB1NRU0tLSrK/NmzePNWvWMGbMGF544YUG+02dOpX777+fkye7TvliOzs71qxZo3YYPZYtpiD44Ycf2Lx5M1u3biU6Opro6Gj272++wNOiRYsoKSmxLidOnLik9yCEaIffVE6sFx/pi7uzPadKqtl+8LQKgYmOIs+IdVPhV95KZeIi/CngUMp3hMZcpXZIQvQ4X331FWvXrmXnzp306dMHgD//+c8sWrSQgrzjfPzpF63u6/e//z3/+te/CA8PZ/jw4daJmtevX4+3tzdTp05l4sSJxMbGMmPGjAYTKD/++OM2fV+X6r777lM7BHERV1xxhfUuamsYDIYGd3aFECpoonIigEGv49aYQNbtPMo7e3KYMtx2FXABaqor+flflmeIRzz0BgYHJ5v2L5oniVg35ejch71u44kp/YZzSe+BJGJC2Nw111xDXV1dg21z597L/dPGAGBSFIKCgti7d6/19Y8++sj6/dixY/n0008Bywfdbdu2NTpGREQEc+fOBcDZ2ZkDBw7Y/H2I7qWzpiAQQnQxv5nU+UJ3xg3g1YSj7DpcRFZRBcGetpsT0GSsY0yJ5fxUaay7SGthSzI0sRvTR90CwJDTX2FqonKWEEKI7ufCKQjq1U9BMG7cuBb2FEJ0a80MTQQIdHdiUqilYuqmPY1fF92TJGLdWPgV11OCM54Uk5nY+iFSQggh1FVeXk5qaiqpqanAr9MQHD9+HLBMQfDaa6/x5ptvkpmZyUMPPWTzKQiEEF2MdWhiDpyfu/JCd421TFHyYXIuVbWmRq+L7kcSsW7M3uDAoX6WSWsrkzerHI0QQojW2rt3LyNHjmTkyJGAJfEaOXIkixcvBmDmzJk8//zzLF68mOjoaFJTU9m2bVujAh5CiB6kryXRoqYUqs41evnKod4E9HOkpKqO//3ctoJRomuSRKybcxp1KwBDz+6gtqZa5WiEEEK0RmumIejIKQiEEF2QvRM4W4YfNjU8UafVcGecZfiiDE/sGSQR6+bCxk2liL70pZzMH7aqHY4QQgghhGivZion1rs1JgB7nZa03BLSThR3VlSig0gi1s3p9HqOelkm/axL/egirYUQQgghRJfVzKTO9TxcDFwbaame+o7cFev2JBHrAdxibwcgrGQXVRVlKkcjRM+m1eoweg3H6DUcrVandjhCCCF6khYqJ9a7e5ylzSdpeZRUXnq5eUenPpz9QwZn/5CBo1OfS+5PtJ4kYj1A6OirOIUXzppqMnbKXTEhOpJGo0FvZ4/ezh6NRqN2OEIIIXqSiwxNBBg1oB9hfq7UGM18mHzikg+p0Wpx9+6Pu3d/NFpJDTqT/LR7AI1WS7ZfvOX7A/9RORohhBBCCNEuFxmaCJYLgvWl7DclHsdsblzqXnQPkoj1EN6X3QnA8PI9lBafUTkaIXous9lE+eksyk9nYTbLPC5CCCFsqH5oYskJaOEcc310f1wMerKKKvjhaNElHbKmupLEV+aQ+MocaqorL6kv0TaSiPUQg4bHkqMNwKCp49B376sdjhA9lqIouBiLcTEWozQx4aYQQgjRbq79QasHUy2UnWq2mbNBz02j+gOXXrTDZKwjrui/xBX9F5Px0p85E60niVgPodFqyQuYCoDh4BaVoxFCCCGEEG2m04NbgOX7FoYnAtw11nL37OuMAk6VVHV0ZKIDSCLWgwSMtwxPDK9K5lxh81dRhBDCFrKyspg0aRLh4eFERkZSUVGhdkhCCNH9taJyIkCITx/igt0xK/Be4vFOCEzYmiRiPUhgSBRHdIPRa8z8suMdtcMRQvRws2fPZtmyZWRkZLBz504MBoPaIQkhRPfXisqJ9epL2b/30wnqTOaOi0l0CEnEepiioGkA9DnyscqRCNGzLFy4EIPBwF133a12KF3CgQMHsLOzY/z48QC4u7uj1+tVjkoIIXqAVlROrHdNuC9efQwUltXw1YGCDg5M2JokYj1M0JWWD4nDatIpyD2qcjRC9ByLFi3ihRde4P333+dIVtceApKQkMD06dPx9/dHo9GwdevWJtutXbuWoKAgHBwciIuLIykpqdXHOHz4MC4uLkyfPp1Ro0axcuVKG0UvhBC9XCuHJgLY67XcNiYQgLf3ZHdgUKIjSCLWw/gOCCHTbjhajULWThmeKIStuLm5MXfuXLRaLfsPHlE7nBZVVFQQFRXF2rVrm22zefNmFixYwJIlS0hJSSEqKoopU6Zw+vRpa5vo6GgiIiIaLXl5eRiNRnbt2sU///lPdu/ezddff83XX3/dGW9PCCF6tn7Blq+tGJoIcHvsALQa2HPsLIcLyjouLmFzkoj1QKVDrgPAPetTlSMRooOUnISsBMvXTmQ0GnFyciLteCm1HmFotbpOPX5rxcfHs3z5cm644YZm27z44ovMmzePOXPmEB4ezrp163BycmLDhg3WNqmpqaSnpzda/P396d+/PzExMQQGBmIwGLj22mtJTU3thHcnhBA9XP3QxLJTUFd90eb+fR2ZHOYDWCZ4bisHRxfyZieRNzsJB0eXNu8v2q9bJGLfffcdGo2myeWnn34CYOnSpU2+7uzs3Gy/aWlp3H777QQGBuLo6EhYWBgvvfRSZ72tDjNk4l2YFA1Djb+QeyRd7XCEaJqiQG1F25ek12B1BLw53fI16bW299HO+b+efvppysvLyTx4EHuDAxqNxsY/lIZWrlyJi4tLi8vx420/6dbW1pKcnMzkyZOt27RaLZMnT2b37t2t6mPMmDGcPn2ac+fOYTabSUhIICwsrM2xCCGE+A0nD7A7//m15ESrdqkv2vGf5FwqaoxtOpxWp8M/KBT/oFC0uq55gbGn6hZPVl922WWcOtWwHPszzzzD9u3biYmJAeDJJ5/kwQcfbNDm6quvZsyYMc32m5ycjLe3N++88w6BgYH8+OOP3H///eh0OubPn2/7N9JJPHwC2O8wksiaFE58/w4BQ/6mdkhCNFZXCSv9L60PxQyfP2lZ2uLPeWDf/EWapiQnJ7Nu3TqmTp1KenrnXOB48MEHufXWW1ts4+/f9p9hUVERJpMJHx+fBtt9fHw4ePBgq/rQ6/WsXLmSCRMmoCgK11xzDdOmTWtzLEIIIX5Do7FUTjx9wDI80TPkortcPtiTIA8nss9U8nFqHnfEDejwMMWl6xaJmL29Pb6+vtb1uro6Pv74Yx555BHrFen6q8P10tLSyMjIYN26dc32e++99zZYHzRoELt37+a///1vt07EAKqG3QBpKfgd/wyQREyIS2E2m3nggQeYP38+Y8aM4Z577uHcySO4+Q1Cq23dwIK8vDyeeuopNm3a1Orjuru74+7u3t6wO1x8fDzx8fFqhyGEED1Pv4G/JmKtoNVquGvsQJZ/lsnbe3K4PTaw1aM2amuqSXljAQCj5ryIvcGhvVGLNuoWidhvffLJJ5w5c4Y5c+Y02+b1119n6NCh1tLKrVVSUnLRDz41NTXU1NRY10tLS9t0jM4QOvEOalP/SpD5OFkHEgkeHqd2SEI0ZOdkuTPVFqV5sDbWciesnkYHDyeCaxvuDNk5temwa9asoaioiGXLlpGVdYy6ujpyD6bg6htEa0d4+/v7tykJA8vQxItVI8zIyGDAgLZd+fT09ESn01FQ0LDUcUFBQYOLXkIIIVTShsqJ9W4eHcBzXx4i81QpKceLGT2wX6v2M9bVMDbfcn6qrFsliVgn6hbPiP3W+vXrmTJlCgEBAU2+Xl1dzaZNm5g7d26b+v3xxx/ZvHkz999/f4vtVq1ahZubm3UJDAxs03E6g1s/Tw44xwKQ/+O7KkcjRBM0GsvwwLYsniEw/SVL8gWWr9NXW7a3pZ82PNt18uRJnnnmGdauXYuzszMhISEYDPakH7JMD5GdnU1UVBR33nknISEhPPTQQ2zdupW4uDgiIiI4fPiwtV1MTIy1/axZswgLC2PmzJkozTyz9uCDD5Kamtri0p6hifb29owePZrt27dbt5nNZrZv3864cePa3J8QQggba8OkzvX6OtlzXZTlnPDOntYncEI9qiZiCxcubLYIR/3y2+cVcnNz+fLLL1tMsrZs2UJZWRmzZs1qdSzp6enMmDGDJUuWcM0117TYdtGiRZSUlFiXEyda9yBlZzMNvxGAAXlfoJhltnXRQ4y6Bx7fD7M+tXwddU+HHu7RRx8lPj6eqVOnApZno8KGBJN+QQn7zMxMFi9ezMGDB/nuu+/44YcfSExM5JFHHuGVV15p1GdmZiZ/+tOfyMjIoKCggO+//77JY7u7uzNkyJAWl6YmUS4vL7cmagBZWVmkpqY2KOyxYMECXnvtNd58800yMzN56KGHqKioaHGkgRBCiE7ShkmdL3TXWMt+n/18ijPlNRdpLdSm6tDEJ554gtmzZ7fYZtCgQQ3W33jjDTw8PLjuuuua3ef1119n2rRpjR5Eb05GRgZXX301999/P08//fRF2xsMBgwGQ6v6VlP4lbdSmfQX+lPAL6kJDB01Ue2QhLANt/6WpYN9+umnfPvtt2RmZjbYHjlsiPWOGEBoaCihoaEAhIWFWasRRkZG8vnnnzfqNzQ0lPDwcABGjhxJdnZ2m4dRt2Tv3r1MmjTJur5ggWXs/6xZs9i4cSMAM2fOpLCwkMWLF5Ofn090dDTbtm1r9d9NIYQQHagdQxMBogL7MiLAjZ9zS/gwOZcHrxzcAcEJW1E1EfPy8sLLy6vV7RVF4Y033uCee+7Bzs6uyTZZWVns2LGDTz75pFV9HjhwgKuuuopZs2axYsWKVsfSHTi5uLHX9QpiyrZzds+7IImYEG0ybdo0zp0712j7Wy//HwCm8+sXXpjRarXWda1Wi8lk+u3uDdrrdLom21yKiRMnNjvc8ULz58/v9oWJhBCiR6q/I1ZdAlXnwLF1z3uB5a7YHz/6mU2JOcwbPwidtmOnWhHt162eEfv222/Jysrivvvua7bNhg0b8PPza7KS15YtWxg2bJh1PT09nUmTJnHNNdewYMEC8vPzyc/Pp7CwsEPiV4N+xM0ADD79FSZj2+aVUENB7lHSf/gfBblHL95YCCGEEKInsncG5/M3K9o4PHH6CH9cHfScOFtFwi895zNtT9StErH169dz2WWXNUimLmQ2m9m4cSOzZ89G18SEdCUlJRw6dMi6/tFHH1FYWMg777yDn5+fdWlp7rHuJmz89ZTijBfnOJj4pdrhtCjpP6vxfG00EV/fhedro0n6z2q1QxJCCCGEUEd9wY42Dk90tNdxS4ylkNzbUrSjS9MorRm/IlpUWlqKm5sbJSUluLq6qh1OI0mrbye2+HMS3a8j7tG31Q6nSfnHD+O9fgxaza//HI2KljPz9uITIOObe4Lq6mqysrIIDg7GwaH7lsZVFIXa6koA7B2cWj1PS2/X0u+/q/8N7Y3kdyJEF/DRXEj/CH63DC5/rE27Hiss56oXdqLRQMJTkwh0b37aFrPJxPFf9gEwYOhItE3czBBt09q/od3qjphoH8dRtwEQevZb6mq7XgWdknNFFL0zt0ESBqDXmCnKOdjMXkKoQ6PRYHB0xuDoLEmYEEKIjtPOyokAg7xcGB/iiaLAu0nHW2yr1ekICoshKCxGkrBOJolYLxA2Lp4zuNGXcjK+/1jtcBrI2LONqpfGElGbxm/vzRoVLZ4Dmx6GKoQQQgjRo7VjLrEL3RlnSeQ2/3SCGqNti0IJ25BErBfQ29lzxNNSTrs27UOVo7Goq61h9+v/j9AvbsOXQnI1vuzpPweT8usdhp8Gz5dhiaLLMZvNlBeeoLzwBGaZn08IIURHaWcJ+3qTw7zxdXXgbEUt29Lzm21XW1PN7vVPsnv9k9TWVLfrWKJ9JBHrJdxibwcgvDiB6spyVWM5eewAx/4+nnG5G9BpFH7qG0/f/7eHcfevpmheMke1QZaGinzIFV2PophxqSvCpa4IRf6NCiGE6Cj1QxOLj0M7LvzpdVruiBsAwNu7m0/mjHU1jDvxGuNOvIaxrus9wtKTSSLWSwwdfRX5eOGsqSYj4SNVYlDMZn7a+gp937yKUOMhSnEmOfZFxjz+Pi6ulvkxfAIGc2b4bAD8cj5BkTsOQgjRLpWVlQwcOJAnn3xS7VCEEO3hGgAaHZhqoexUu7q4bUwgeq2GvTnnyMgrtXGA4lJJItZLaHU6svymWFb2/6fTj19ytpCUf9zImNS/WJJB+0gq793J6GvnNmobetXd1Cp6gszHOXYgqdNjFUKInmDFihWMHTtW7TCEEO2l04NbgOX7dg5P9HZ1YMpwXwDeSZRS9l2NJGK9iNfYOwAYXr6bspKznXbcjN1fUPXyWEaX7aBO0bEn6GFC//gdvgNCmmzv1s+TdJdxABT+8FanxSmEED3F4cOHOXjwIPHx8WqHIoS4FJdQObHeXWMtfWzdd5Ky6jpbRCVsRBKxXmRw5DiOa/tj0NRx8Lv3O/x4dbU17H7tMYZtux1fisjV+JE147+Mnb0SnV7f8s4jZgIwKP8LTEZjh8cqhBCdJSEhgenTp+Pv749Go2Hr1q2N2qxdu5agoCAcHByIi4sjKaltowOefPJJVq1aZaOIhRCqucTKiQBjB7kzxNuFyloTW/adtElYwjYkEetFNFotJ/tfC4B95n879Fi5R9LJ+vsVjDu5Ea1GIanvtfRbsIehoya2av/hV95EKc54c5bMPZ93aKxCCNGZKioqiIqKYu3atU2+vnnzZhYsWMCSJUtISUkhKiqKKVOmcPr0aWub6OhoIiIiGi15eXl8/PHHDB06lKFDh3bWWxJCdJRLrJwIlvkv77qgaIfy2/mChGoucltC9DT9x98F775GeFUK5wpP0c/Lz6b9K2YzP338ChGpywnQ1FCKM4djVxB77Zw29WNwcCLV/Wrizn5C5d734IrrbBqnEEKoJT4+vsUhgy+++CLz5s1jzhzL381169bx2WefsWHDBhYuXAhAampqs/vv2bOH999/nw8//JDy8nLq6upwdXVl8eLFze5TU1NDTc2v1dJKS+WhfiG6BOsdsUt7vuvG0QE8u+0Qh0+Xk5R1lrhBHpcem7hkckeslxkwNJojusHYaUz88t0mm/ZdcraQfS9eT2zaMzhpajhgP4LKuQmMbmMSVq/PGMszbeHndqhecl+Ielqtjmq3IVS7DUGr1akdjuhhamtrSU5OZvLkydZtWq2WyZMns3v37lb1sWrVKk6cOEF2djbPP/888+bNazEJq9/Hzc3NugQGBl7S+xBC2IgNhiYCuDrYcf3I/gC8vadhUmdwcOaX6z7hl+s+weDgfEnHEW0jiVgvVDRwKgAuhz+2WZ8Hfvyc6pfjGFW+kzpFx+7g+Qz74w58A4e0u89hsdeQjxcumioOfLfZZrEKcSk0Gg0Ozn1wcO6DRqO5+A7dVFZWFpMmTSI8PJzIyEgqKirUDqlXKCoqwmQy4ePj02C7j48P+fnNT8h6qRYtWkRJSYl1OXHiRIcdSwjRBvVDE8tOgfHS5vi6a6xleOK29HxOl/06cbNOr2foqCsZOurKiz/DL2xKErFeaOCVdwEQVrOf0yezLqmvutoadv/7UcK+vAMfznBC40/WjC2Mm7Xikv8za3U6svwtz7Rp09WZ+0yI3mr27NksW7aMjIwMdu7cicFgUDsk0Q6zZ8/m+eefv2g7g8GAq6trg0UI0QU4e4KdE6BA8aVdIBnu78aoAX0xmhU2J8nFlq5AErFeyG9gKAftwtFqFI7tfKfd/Zw4sp/sv1/OuLw3LQU5+k3FfcFuho660max+o+fBUBERSLnCts3maEQtmQ2myk/c5LyMycxd8KE42fOnMHb25vs7OwOP1a9AwcOYGdnx/jx4wFwd3dHf8GFldtuu40XXnih0+LpTTw9PdHpdBQUFDTYXlBQgK+vr0pRCSFUo9HYbHgiwN3jLHfY3ks6jtFkOYfV1lSz5+3F7Hl7MbU11S3tLmxMErFeqmSwpfiF+7FP2ryvYjbz039fwuPtqwkxHqYEZ1LGvkTsY+/i3KevTeMcGDaao7pBlmfadrxt076FaA9FMeNScxqXmtMoSscnYitWrGDGjBkEBQUBrSt9DpdW/vzw4cO4uLgwffp0Ro0axcqVKxu8/vTTT7NixQpKSkra+7ZEM+zt7Rk9ejTbt2+3bjObzWzfvp1x48apGJkQQjXWyonZl9xVfIQf7s725JVU8+1BSyVWY10NY4++xNijL2Gsu7Thj6JtJBHrpQZPvAuTomGo8RdOHjvQ6v1KzhSw78UZjPl58fmCHFFU37eLUb+f3WGxFgbPAMDt8JYOO4YQXVFlZSXr169n7ty51m0XK30Ol17+3Gg0smvXLv75z3+ye/duvv76a77++mvrvhEREQwePJh33mn/HfXerLy8nNTUVGvlw6ysLFJTUzl+/DgACxYs4LXXXuPNN98kMzOThx56iIqKCmsVRSFEL2ODSZ3rOdjpuCUmAGhctEN0PknEeilP30AyHKIBOJHQuuqJB374jOo14xhVnmApyDHoUYb98Vt8AgZ3YKQw+KrZmBUNw+oyOHkss0OPJURTbrnlFry8vPj3v/9t3ZaYsh9HRye++uqrDjvu559/jsFgYOzYsdZt8fHxLF++nBtuuKHZ/S4sfx4eHs66detwcnJiw4YN1japqamkp6c3Wvz9/enfvz8xMTEEBgZiMBi49tprG5VLnz59Ou+/3/ETw/dEe/fuZeTIkYwcORKwJF4jR460VjacOXMmzz//PIsXLyY6OprU1FS2bdvWqICHEKKXsOHQRIA7Ywei0cCuw0VkFUkhJjVJItaLVQ29HgCf45+22K62pprdrz5C2Fd3WgtyZN/wMePu+b9Oqa7j5R/EgfNJ4/GdGzv8eKKT1VY0v9RVt6FtVevatsPLL7/MTTfdxLJlywAor6jkrkee5sEHH+Caa65pV5+tsWvXLkaPHt2mfWxR/nzMmDGcPn2ac+fOYTabSUhIICwsrEGb2NhYkpKSGsw9JVpn4sSJKIrSaNm4caO1zfz588nJyaGmpobExETi4uLUC1gIoS4bTOp8oQEeTkwc6gXAu4lyV0xNUqOyFwuddCe1Py8j2JxDVsZPBIePadTmxOE0qt+/l3GmI6CBpH7TiJj7T5xc3Do11qphN0HaPvqf+B+KeRUarVxD6DFW+jf/Wsg1cOeHv64/NwTqKptuO/AKmPPZr+urI6HyTON2S9v+XJOfnx+PP/44r776KmfOnGHhM3/HYLBn1apVbe6rLXJycvD3b+Hn04SWyp8fPHiwVX3o9XpWrlzJhAkTUBSFa665hmnTpjVo4+/vT21tLfn5+QwcOLBNMQohhGgDGw5NrHf3uIHsOFTIB3tzeegyX5xs1rNoC/k024u5uXtxwNlylTX/h4bDExWzmaT/rMbjnd8RYjpCMS7sG/cysY9t6vQkDCDsqjupVuwYYD7JkbTvO/34QgwdOhQnJyeWLl3Kpi1fsGnNChwcHDr0mFVVVR1+jObEx8ezf/9+0tPTefHFFxu97ujoCFieYxNCCNGB6u+IVRdDVbFNurxyqDf9+zpSUlXHVwcKLr6D6BByR6yXM4XfAHt/JDDvCxTzi2i0WkrOFHB0w33EViSABtIN0Xjfs5GR/YNVi7OPmzvJrlcwumwHZ3a/Q8jICarFImzsz3nNv6bRNVx/6kgLbX9zXenx/e2PqQlarZbIyEj+9a91/P3px4gaPhSTTY/QmKenJ+fOnWvzPp1R/vzs2bMAeHl52axPIYQQTTC4gJMnVBZZhic69r3kLnVaDXeOHcDftx3iw5ST3HbpUYp2kDtivVz4xJlUKgYClHwSP3yOvZ+9Rs2asYyqSKBW0bFn8GOE//FbvFVMwurpoi1/Joac/hJjXa3K0QibsXdufrFzaENbx9a1bSdFUQAYNWoUf3hqMVWuwWi1OrKzs4mKiuLOO+8kJCSEhx56iK1btxIXF0dERASHDx+29jFt2jRGjx5NREQEmzZZ7kLv3r2b2NhYjEYjBQUFhISEkJ+fD8DIkSPJyMhoU5ydVf48PT2dgIAAPD09bdanEEKIZnTA8MRbYwKx12nZl1fF12PWc+B372JwaP95UrSdJGK9nJOLG7n6AQCMzVzJ6KQn8eYsJzT+5NzwMWPvXoZWp7tIL51j+PgbOEcfPCkm44f/qR2O6GVWr15NYmIiZrMZZ1d3HF36otFoAMjMzGTx4sUcPHiQ7777jh9++IHExEQeeeQRXnnlFWsfb731FsnJySQmJrJixQpqamoYN24cEyZM4Nlnn+Xhhx9m8eLF1jtXU6ZM4cCBAw3uil2s9Dl0TvnzXbt2dWihEiGEEBewceVEAE8XA9dG+mJGy9dVIQy/fGqnFGETv5JErJcryD3KEOOvw700GjAroNz+PiHR41WMrDE7ewO/eP4OgJp9UjZbdJ79+/ezaNEi/vCHP5CRkYHRaGzwemhoKKGhoeh0OsLCwqwVCyMjI8nOzra2+8c//kFUVBSXXXYZx48ftyZPy5cv5+2336a6upq7777b2j4yMpJRo0bxwQcfWLddrPQ5dHz58+rqarZu3cq8efNs0p8QQoiLsHHlxHp3jbX0+0laHiWVdTbtW1ycJGK9XGFOBlqN0mCbVgNlhbkqRdQyt7g7ARhevJPK8rZXvxOiraqrq7njjjuYOXMmy5cvp7a2lpTEXZSfPYXZbAbAYDBY22u1Wuu6VqvFZLI8SbZjxw7rnbK0tDSGDRtmLf1++vRpamtrrRUPL7R48WJeeukl67FaU/ocOrb8+RtvvEFsbGyD+c2EEEJ0oA4YmggwemA/wr0duMW8jY/+tZgThcU27b+tTpVU8ePRIk6VVF28cQ+IQ+4/9nJeA8MxKRp0FyRjRkWL58BhKkbVvNDRV3Hycx/6U8DeHe8TM/0BtUMSPdzChQupqKjglVdeoU+fPgwcOJB/r32ZZU8+hGNU6wtVlJaW4uHhgYODA6mpqaSlpVlfmzdvHmvWrGHbtm288MIL/PGPf7S+NnXqVA4fPszJkycJDAy06XtrLzs7O9asWaN2GEII0Xt0wNBEAI1GQ5i3gf8r3QhlEPZCDLOuDOeqYZ0/gfy3Bwt4NeEYimIZofXAhEGqx6HVwKobI5k5ZkCHHEsSsV7OJ2AwSSOWMurnv6LXmDEqWlJGLCE2YLDaoTVJo9Vyov80+ueux+7AhyCJmOhAX331FWvXrmXnzp306dMHgD//eRF/XriQM+dK+OjTq1rd1+9//3v+9a9/ER4ezvDhw60TNa9fvx5vb2+mTp3KxIkTiY2NZcaMGYSGhlr3ffzxx236vi7Vfffdp3YIQgjRu1iHJh4HsxlsNJ/qqZIqvkgv4IULamOt23mMdTuP2aT/9lKUrhGHWYE//zedCUO98HNzvPgObaRR6kuBiXYrLS3Fzc2NkpISXF1d1Q6nXQpyj1KUcxDPgcPw6aJJWL3jv6Qy4N0rMSpaSv6wHw+fALVDEq1QXV1NVlYWwcHBqs2NZQsmkxFdgaU0vsknEp1Orme1Rku//57wN7Snkd+JEF2MqQ6We4NihgUHwdXPJt3+eLSIua/tJNPhXgDCqjdQhQN+bg442ndesbaqWhOnSqobbe8qcbw3byzjBnu0up/W/g2VTxACsNwZ6+oJWL0BQ6P5RT+UocZfOPztW3jc/me1QxJCCCGE6Dg6O3ALsNwRO5dts0Qs2NMZzW8PpdHw3z9c1iF3gJpzqqSKy//2LeYLbg91pTiCPJ065HhSrEN0S2cHXw+A+9Et6gYihBBCCNEZOqByop+bI89MDbOua4GVN0Z0avJTH8eqGyPRnZ8WRqfR9Io45I6Y6JaGTLoH48HnGWr8hRNH9hM4JFLtkIQQQgghOk6/gZC9y+aVE2eM7A/bLd9/+uh4gv29bdp/a80cM4AJQ73ILqokyNOp05MwNeKQO2KiW/L0DSTD0VLsIHfnmypHI4QQQgjRwTqocuKFfFwNF2/UgfzcHBk32EO1JKyz45BETHRbteE3AxCY+z+U83MsCdHRtBotlS4DqXQZiFYjf0KFEEJ0kr5Blq82ntTZ3uBI2oRXSZvwKvYGdROg3kY+RYhuK/yq26lUDAQo+RxK2aF2OKKX0Gi1OLm64+TqjsZG5YOFEEKIi+qgSZ31dvZEXXUbUVfdht7O3qZ9i5bJpwjRbTm5uJHhNh6Akj3vqByNEEIIIUQHqh+aWHoSjDWqhiJsQ4p1iG7NfuTtsPMbhhZ9Q11tDXb26o5tFj2f2WymqqQQAEc3L7RyV6zbCQ4ORqP5bcHmi3v88cd59NFHOyAiIYRoBWcvsHOCukooyQUP20w7VFdbw77P/g3AyKn3y2epTiSJmOjWwq+4jjM73fCghLTvtxB11W1qhyR6OEUx41yVB4DJ1QMZWND9bNy4sV37BQUF2TQOIYRoE43GUsK+MNNSsMNmiVg1sWlPA1D5u3skEetEkoiJbk1vZ89h7yl4nP4A477NIImYEOIirrzySrVDEEKI9ul3QSImuj25lCu6PfexdwEQXvo95aXnVI5GCCGEEKKDdMCkzkI9koiJbi8kejwnNP44amrJ/PZdtcMRPdSVV16JRqNBr7dD038Umv6j0OvtuOeee9QOTdiAl5cX3t7eTS6BgYFMmDCBHTukOqsQQmXWucQkEesJZGii6PY0Wi25A64jMGcdDpkfwfUPqx2S6GEURWHfvn08//zz3HbbTHSFmQCYvMJwc+urbnDCJgoLC5t9zWQykZ6ezl133cX+/fs7MSohhPgNawn7bFXDELYhd8REjzDgytkAhFfvoyhPrhIJ2zp8+DBlZWVMmDABX19ffL09LYuvLy4uLmqHJ2ygsLCQjIyMRtszMjI4e/YsUVFRPPHEEypEJoQQF5ChiT2KJGKiR+g/KIyD+jB0GoUjO95UOxzRBhUVFc0u1dXVrW5bVVXVqrbtkZycjF6vZ8SIEe1+n6Jrmz9/PufONX7G9Ny5c9aS9bNnz+7kqIQQ4jfq74hVnYPqEnVjEZdMEjHRY5SE3ACA57Gt6gYi2sTFxaXZ5aabbmrQ1tvbu9m28fHxDdoGBQU12a49UlJSMJlMeHh44ObWF5ehV+Ay9AoeeugPAHh6ejbaJzs7m5iYmAbbZs+ezaeffgpAbm4uN954I4MHDyYmJoZbbrmFgoKCZvsTHSsrK4vLL7+80fbLL7+c9PR0FSISQogmGPqAk4flexs9J2ZvcCQ5djXJsauxNzjapE/ROvKMmOgxhl51D3UZzzLEdJSczGQGho1WOyTRQ6SkpHD77bfz17/+tcF2d3f3dvWnKAozZszgD3/4A//9738B2LVrF4WFhfj4+FxyvKLtmrobVu+3d1uFEEJVfQdC5RnL8ES/Sx+pobezZ/S1c2wQmGgrScREj9HPy49U51iiK3eT9/1bkoh1E+Xl5c2+ptPpGqyfPn262bZabcMb/NnZ2ZcU14VSUlJYsWIFQ4YMsUl/27dvx8XFhblz51q3jR8/3iZ9i/YZMWIEGzdubDT88K233iIyMlKdoIQQoin9giAvRQp29ACSiIkexTT8ZvhpNwPzPsdsehHtbz7Ii67H2dlZ9bYtOXbsGMXFxURFRQGgmM1Ulp0FwKmPOxpt20d4Z2RkMGrUKJvEJ2zj5ZdfZsaMGbz55pvW301KSgplZWVs3bpV3eCEEOJC1sqJthmaaKyrJe3rTQBE/e5O9Hb2NulXXJwkYqJHCZ84k/Kkp/HnNJl7vyEsboraIYluLjk5GQAfHx/y8/MxmYzW8vX2w6/ETtv0CUuj0bRpu1BX//792bt3L9u3b7dWT4yPj2fy5MkqR9Z+WVlZ3HvvvRQUFKDT6dizZ4/NLlAIIVRk48qJtTVVjE56HIDKCTdKItaJJBETPYqjcx9+6jeRMcVfUJq0CSQRE5coJSUFgJCQkAbbDQZ7zp07h10zJywPD49Gzx2dPXsWT09P7O3trc+Gia7h888/t34/ePBgNBoNbm5uVFZW4uTkpGJk7Td79myWL1/O+PHjOXv2LAaDQe2QhBC2YJ3UOVvNKIQNSNVE0eM4jr4dgGFnvqG2pvoirYVo2apVq1AUxboYjXUoJ1OoPrYHe/vmrxq6uLjQt29ffvzxR8BSJXH//v0MHz6cyZMnU1paysaNG63tv//+e6nOp6IPP/ywwfLBBx+wfPlyIiMj2b59u9rhtdmBAwews7OzPnvo7u6OXi/XXoXoEeqHJhYfB0VRNxZxSSQREz1O2LipFNIPNyo4sPMjtcMRvcC5c+cICAiwLu+99x4Ab775JgsXLiQ6Oprrr7+eV199FRcXFzQaDVu3bmXr1q0MHjyY4cOHs2bNGry8vDAajXLnQgVvvPFGo+Xjjz9m586dLFq0yObHS0hIYPr06fj7+1v/PfzW2rVrCQoKwsHBgbi4OJKSklrd/+HDh3FxcWH69OmMGjWKlStX2jB6IYSq3AJBowVjNZQXqB2NuARyeUz0ODq9nqM+v8er4D2UnzfDNXepHZLo4UwmU5PbIyIiSEhIaPK1AQMGNPnhOy0tjeDgYFuGJy5BQEAAdXV1Nu+3oqKCqKgo7r33Xm688cZGr2/evJkFCxawbt064uLiWL16NVOmTOHQoUN4e3sDEB0djdFobLTvV199hdFoZNeuXaSmpuLt7c3vf/97xowZw+9+9zubvxchRCfT2YFrAJQctwxP7OOrdkSinSQREz2S1+X3wH/fY3jZbkrOFeHWTybIFV3fG2+8wfPPP89LL72kdijivN27d+Pq6mrzfuPj4xtNQn6hF198kXnz5jFnjmVun3Xr1vHZZ5+xYcMGFi5cCEBqamqz+/fv35+YmBgCAwMBuPbaa0lNTW02EaupqaGmpsa6Xlpa2ta3JIToTP0Gnk/EcmDAWLWjEe0kiZjokQZFjCV76wCCzMdJ+/YdYm96XO2QhLioOXPmWD94i841ZsyYRhUtz549S79+/XjzzTc7NZba2lqSk5MbDInUarVMnjyZ3bt3t6qPMWPGcPr0ac6dO4ebmxsJCQk88MADzbZftWpVownLhRBdWN+BwC6bVU4U6pBETPRIGq2WUwOnE5S1FudD/wUeVzsk0UNoNFoqHP0BcNTIY7Y9xUcfNXyeVKPR4OHhgbOzM5s3byY8PLzTYikqKsJkMuHj49Ngu4+PDwcPHmxVH3q9npUrVzJhwgQUReGaa65h2rRpzbZftGgRCxYssK6XlpZa76YJIbogG1ZOtLN3IClqOQAj7R0uuT/Ret3iU8R3332HRqNpcvnpp58AWLp0aZOvt3bOlDNnzhAQEIBGo6G4uLgD343oLEETZwEQVvMz+SeOqByN6Cm0Wi3O/Xxw7ueDth2TOYuuaeDAgQ2WAQMGWM8fTz31lMrRtU98fDz79+8nPT2dF198scW2BoMBV1fXBosQoguz4aTOdvYGYm94hNgbHsHOXopFdaZu8Snisssu49SpUw2W++67j+DgYGJiYgB48sknG7UJDw/nlltuadUx5s6dy4gRIzrybYhO5jcwlAz7SLQahewdnTu0SAjRcyidXB7a09MTnU5HQUHDamgFBQX4+spD+UIIbD6ps1BHt0jE7O3t8fX1tS4eHh58/PHHzJkzxzqm38XFpUGbgoICMjIymDt37kX7/9e//kVxcTFPPvlkR78V0cnKh1qqkXlnf6JyJKKnUMxmKkvPUll6FsVsVjsc0Ql+++xYR7O3t2f06NEN5i8zm81s376dcePGdWosQoguqn5oYkkuGGsvqStjXS1p375P2rfvY6y7tL5E23TLZ8Q++eQTzpw50+JD7a+//jpDhw61TmbZnIyMDJYtW0ZiYiLHjh1r1fGlulT3EXrV3dTuX8EgczbH0hMZFBGndkiimzMrZpzKLVcgTc6u6LrH9SxxEV5eXk0mXIqidMhw9fLyco4c+XXIdFZWFqmpqbi7uzNgwAAWLFjArFmziImJITY2ltWrV1NRUSHFXIQQFi7eoHcEYxWUnACPwe3uqramiqgESzGfyth49Hb2topSXES3TMTWr1/PlClTCAgIaPL16upqNm3aZC3x25yamhpuv/12nnvuOQYMGNDqREyqS3Ufbu5e7HMZy8iK7yn44S1JxIQQTSosLOzU4+3du5dJkyZZ1+sLZcyaNYuNGzcyc+ZMCgsLWbx4Mfn5+URHR7Nt27ZGBTyEEL2URgN9B0DRIcvwxEtIxIR6VL2Uu3DhwmaLcNQvv60QlZuby5dfftnikMMtW7ZQVlbGrFmzWjz+okWLCAsL46672jbh76JFiygpKbEuJ06caNP+opONuBWAQae+wNzMxLtCiN5l9uzZVFZWqnb8iRMnoihKo2Xjxo3WNvPnzycnJ4eamhoSExOJi5MLSUKIC9iwcqJQh6qJ2BNPPEFmZmaLy6BBgxrs88Ybb+Dh4cF1113XbL+vv/4606ZNu+iVw2+//ZYPP/wQvV6PXq/n6quvBiwPSi9ZsqTZ/aS6VPcSNuFmSnHGhzNk7tmmdjhCiC7g7bffpry83Lr+0EMPNRqCaDQaOzkqIYRoAxtWThTqaPXQxNmzZ/PPf/4TJycnmx3cy8sLLy+vVrdXFIU33niDe+65Bzs7uybbZGVlsWPHDj755OLFGf7zn/9QVVVlXf/pp5+499572bVrF4MHyy3ensLB0Zmf+00i9tynVOzdBJdPVTskIYTKflsJcdOmTTz11FP07dsXsFQoDA4OVvWumRBCtEgqJ3Z7rb4j1hWuHn777bdkZWVx3333Ndtmw4YN+Pn5ER8f3+i1LVu2MGzYMOv64MGDiYiIsC7BwcEAhIWF4e3tbfs3IFTjPOZOAMLOfkt1VYXK0QghupqmStRXV1erEIkQQrSSDE3s9lqdiDV19fDs2bPW9YKCgg4ford+/Xouu+yyBsnUhcxmMxs3bmT27NnodLpGr5eUlHDo0KEOjVF0TWFxU8jHkz6aKjK++0DtcEQ3dOWVV6LRaNDr7dD0H4Wm/yj0ejvuuecetUMTHaSzy9YLIUSbyNDEbq/dVRPVuHr47rvvtvi6VqttsXDG7NmzmT17drOv1z88LXoerU5Hlv+1+Oa9hWb/hxAvJaBF6ymKwr59+3j++ee5/fbbqSqxVNhzdPOSZ0S7sXfffZcJEyYQGRmpdihCCNF29UMTq85CdSk4tO98ZGfvQGLYIgBG2TvYKjrRCjYtXy9XD0VX5nvFPfDBWwyv2EPJmQLcPKQMtGidw4cPU1ZWxoQJE/D39wd/f7VDEpdo/PjxLFmyhLKyMuzs7DAajSxZsoTLL7+c6OjoNj2/LIQQqnBwBUd3SyJWnAO+7buoZGdvIG5my1M+iY7RpqqJ7777LikpKdTV1XVUPEJ0mODwMRzVBWOvMXHw27fVDkd0I8nJyej1ekaMGKF2KMJGdu7caR2u/uabb/LEE09w6tQp/vznP3PZZZcxdOhQtUMUQoiLk+GJ3Vqr74jJ1UPRExQGz2DwkdW4/vJf4Em1wxFARYWleIqTk5P1rnptbS11dXXo9XoMBkOjto6Ojmi1lutIdXV11NbWotPpcHBwuGjb5iqutiQlJQWTyYSHh0eD7XfccQf//ve/0ev1REREWLfv3r0bR0dHcnNzefTRR0lLS6Nfv34EBwfzyiuv4OPjg6enJ0VFRW2ORdhWSEgIISEh3HbbbdZtWVlZ7N27l3379qkYmRBCtELfgZC375IqJ5qMRg4mfgnAsLgp6PQ2HTAnWtDqO2Jy9VD0BIMmzcasaAirO0BethRu6QpcXFxwcXFpkJQ899xzuLi4MH/+/AZtvb29cXFx4fjx49Zta9euxcXFpdEk70FBQbi4uJCZmWndduFkuW2RkpLC7bffTmpqKsnJe0n9chOpX25ixYrlAPTt25fU1FTr4ujoiKIozJgxg6lTp3L06FH27t3Lo48+SmFhYbtiEJ0nODiYW265hZUrV6odihBCtMwGlRNrqisY/vUdDP/6DmqqpbJ0Z2pzyitXD0V35t0/mHSHKCJqUjn+3Zv4z5YPWuLiUlJSWLFiBUOGDMFkMqLrY5l/0OTu3uw+27dvb5Qgjh8/vsNjFUII0YvI0MRuzSb3HoODg61XEIXo6iqH3QRpqfgd/xjFvByNtk2PSgobq5+f8MLJ4p966ikef/xx9L8ZHnH69GnAMtyw3sMPP8y8efMaTVmRnZ3dqG1LVVObc+zYMYqLi4mKimq2TXFxMdHR0QDExMTw+uuvk5GRwahRo9p8PCGEEKLVZFLnbk0GgYpeZ9ikO6lOXcZAcy5H9v/IkKgr1A6pV3N2dm60zd7eHnt7+1a1tbOza/K5r+batlVycjIAPj4+5OfnW+6IFVqGUXp4mdHpfh2aKLq3X375hUGDBjW6ACCEEF2WdWhiDigKSAXzbkVuBYhex7WvBwf6XA5A0Y/vqByN6OpSUlIAy7BsPz8/AgIC8Rt5DUFjp2E0GpvdLywsTIZrdzNhYWEcO3ZM7TCEEKL13AIBDRiroPy02tGINpJETPRKuqiZAAwp+AJTCx+mhVi1ahWKolgXo7EO5WQK1cf2NHnXrt7kyZMpLS1tUCDk+++/Jz09vROiFu2hKIraIQghRNvo7cEtwPK9DE/sdiQRE71S+IQbKcYFT4rJ+PF/aocjeiCNRsPWrVvZunUrgwcPZvjw4axZs0am+hBCCGFb9c+JXULlRKEOGQgveiV7gwP7PCYTd2Yr1cnvw4Qb1A5JdBMajZZygzcAThrLtazm5gMbMGAAW7dubfI1mUNMCCGETfQbCDnft7tyot7OwJ7BjwEwys5wkdbCliQRE72WW9xd8PlWwou/o6qiDEfnPmqHJLoBrVaLi0d/tcMQQgghLOoLdhRnt2t3e4MDY+9eZrNwROvJ0ETRa4XGXE2exgdnTTUHvntf7XCEEEIIIdqur8wl1l1JIiZ6LY1WS07/qQDYpX+ocjSiu1AUheqKMqoryqS4gxBCCPVd4qTOJqORX1J28kvKTilg1skkERO9mv/4WQCEV+7l7OmTKkcjugOz2YRDyREcSo5gNpvUDkcIIURvVz80sTQXTHVt3r2muoKhn1zH0E+uo6a6wraxiRZJIiZ6tYGh0RzWh2CnMXF4x9tqhyOEUNGf/vQnPDw81A5DCCHaxsUH9A6gmKHkhNrRiDaQREz0emcGzQCg7+EtKkfSO8hwvt7JbDarHcJFrVq1ShIxIUT3o9FA3wGW7+U5sW5FqiaKXm/IVbMwHXqBUONBco+kEzAkQu2QeiQ7Ozs0Gg2FhYV4eXmh0WjUDqldTCYjOqMlmTRVV6PTyZ/RliiKQm1tLYWFhWi12hYnwRZCCNFO/YKg6BeZS6ybkU8Qotfz9B3Az46jGVG9lxMJbxEw5O9qh9Qj6XQ6AgICyM3NJTs7W+1w2s1sNqMtLbR8X2ZAq5WBBa3h5OTEgAED5OclhBAdob5yYrHcEetOJBETAqgNvxlS9hJw4n8o5r+hkQ+LHcLFxYWQkBDq6tr+MHFXUVVRhuMXMy3fz9kh88+1gk6nQ6/Xd9u7oEII0eVdYuVEoQ5JxIQAwibdTmXyEgLJ45fUBIaOmqh2SD2WTqdDp9OpHUa7mY01OJRbHoY2G+xxcHBQOSLRHrNnz+af//wnTk5OaocihBCXrr5yogxN7Fbksr8QgHOfvmS4jQfg3G6pniiap7czsDtwHrsD56G3M6gdjmint99+m/Lycuv6Qw89RHFxcYM2RplPRwjRXVzC0EQ5r6lHEjEhzrMbeRsAIYVfUVdbo3I0oquyNzgwbu7zjJv7PPYGuRvWXf22euemTZs4e/asdb2goABXV9fODksIIdqnfmhi5RmoKWvTrnJeU48kYkKcN/yKGZzFFXdKyfjhE7XDEUJ0oqamVaiurlYhEiGEaAcHN3DsZ/lenhPrNiQRE+I8vZ09h72uAaBu3/sqRyO6KrPJRHbmXrIz92I2mdQOR3Sg7l5c5B//+AfDhw8nPDycRx99VObwE6Kna+fwRDmvqUcSMSEu0G/c3QCEl+yioqxY3WBEl1RdVU7Q5qsJ2nw11VXlF99BdFnvvvsuKSkp3bqKZ3MKCwt55ZVXSE5OZv/+/SQnJ7Nnzx61wxJCdKR2Vk6U85p6JBET4gIh0RM4ofHHSVPDz+8tpiD3qNohCSE6wPjx41myZAkxMTG4uLhQWVnJkiVLWLduHXv27GlQyKO7MhqNVFdXU1dXR11dHd7e3mqHJIToSFI5sduRREyIC2i0WoocLFeUxuW9iedro0n6z2p1gxJC2NzOnTspKSnh0KFDvPnmmzzxxBOcOnWKP//5z1x22WUMHTq0Q4+fkJDA9OnT8ff3R6PRsHXr1kZt1q5dS1BQEA4ODsTFxZGUlNTq/r28vHjyyScZMGAA/v7+TJ48mcGDB9vwHQghuhyZ1LnbkXnEhLhAQe5RRlTugfOPhug0CmN+XsLPv3xMhecI7HzD6RcURf8hkTg4OqsbrBDikoWEhBASEsJtt91m3Xbs2DGSk5PZt29fhx23oqKCqKgo7r33Xm688cZGr2/evJkFCxawbt064uLiWL16NVOmTOHQoUPWO1vR0dFNltj/6quvcHR05NNPPyU7OxtHR0fi4+NJSEhgwoQJHfaehBAqk0mdux1JxIS4QGFOBj6ahg+0azQwoiYFTqbASSAZTIqGE1o/ipwGUdMvFDu/cNyDo+g/OFJKvwrRzQ0aNIhBgwZxyy23dNgx4uPjiY+Pb/b1F198kXnz5jFnzhwA1q1bx2effcaGDRtYuHAhAKmpqc3u/+GHHzJkyBDc3d0BmDp1Knv27Gk2EaupqaGm5tdpO0pLSy/6HkwmU498vk40T6fTodfru30hmx6rX7Dla3EOKIrlA4zo0iQRE+ICXgPDMSkadBckYyZFQ9LA+9FW5ONaepj+dTm4aioIVPIIrMiDiu8hF/gJ6hQdOTp/zjgFU+M+DINfOB7BUfgPGo6dvUySKERXEBwc3K4Pko8//jiPPvpoB0TUUG1tLcnJySxatMi6TavVMnnyZHbv3t2qPgIDA/nxxx+prq7Gzs6O7777jvvvv7/Z9qtWreKvf/1rq2MsLy8nNzdXKjH2Qk5OTvj5+WFvb692KOK33AIADdRVQkUhuMhzoV2dJGJCXMAnYDBJI5Yy6ue/oteYMSpaUkYsYdxNj1vbKGYzhfnHyT+yj4rcdLSFB3ErO0L/uhxcNFUMNJ9gYPkJKE+A40Ai1Co6snQBnHUeRK17KAb/CLwGjcA/eDg6vfw3FKIzbdy4sV37BQUF2TSO5hQVFWEymfDx8Wmw3cfHh4MHD7aqj7Fjx3LttdcycuRItFotV199Ndddd12z7RctWsSCBQus66WlpQQGBjbZ1mQykZubi5OTE15eXnJ3pJdQFIXa2loKCwvJysoiJCQErVZKDXQpegO4+kPpScvwREnEujz5BCjEb8Te9DgFcdMpyjmI58BhxAY0fMBdo9Xi5R+El38QcIN1u2I2k3/yGKePpFJ5Mh1t0SH6lh8hoC4HJ00NweYcgstyoGwH5AC7oUaxI1sfwDnnwdR5hOLoH4HX4Cj8Bg5Dq9NRkHuUwpwMvAaG4xMgD9p3BXo7A3t87wRglJ3c5eyOrrzySrVD6BQrVqxgxYoVrWprMBgwGFr377murg5FUfDy8sLR0fFSQhTdjKOjI3Z2duTk5FBbW4uDgwzF73L6BZ1PxLIhcEyrdpHzmnokEROiCT4Bg9uc+Gi0WnwDh+AbOAS42brdbDKRd+IIhcf2UZV7AN2ZQ/QrP0J/4wkcNbUMNmVBaRaUfgNZwA9QpdhTrHHFVynCR3N+eOSIpcRecGdOqMPe4MDYB/+pdhiiB/P09ESn01FQUNBge0FBAb6+vipF1ZjcCeud5C5YF9d3IOT8AMXZrd5FzmvqkURMiA6m1enwDwrFPyi0wXaT0cjJ44coPJpGVV46dmcO0a/iKAHGXBw1tThS1KB646if/0pB3HS5MyZED2dvb8/o0aPZvn07119/PQBms5nt27czf/58dYMTQnRtUjmxW5FETAiV6PR6+g8aTv9BwxtsN9bV8tNnrzEm9c8Ntus1ZopyDkoipjKzyUT+iSMA+AYOQavTqRyR6I7Ky8s5cuSIdT0rK4vU1FTc3d0ZMGAACxYsYNasWcTExBAbG8vq1aupqKiwVlEUQogmtWNSZzmvqUcSMSG6GL2dPQNifo9p318aVG9UFPAI7NhJZsXFVVeV478xFoDKJ4/j5OKmckSiO9q7dy+TJk2yrtcXypg1axYbN25k5syZFBYWsnjxYvLz84mOjmbbtm2NCngIIUQD7ZjUWc5r6pGBvkJ0QT4Bg0kesRSjYvkvWj8dSG7atypHJoSwhYkTJ6IoSqPlwoqO8+fPJycnh5qaGhITE4mLi1MvYNFjFRcXExMTQ3R0NBEREbz22mtqhyQuRf3QxJKTYJJ5/ro6uSMmRBd1YfXG0p//x7iC9xiSvIyiMfF4+g5QOzwhhBA9QJ8+fUhISMDJyYmKigoiIiK48cYb8fDwUDs00R4uvqAzgKkGSnLBPVjtiEQL5I6YEF2YT8Bghl8+lZi5L3FEN5i+lHPirQdQzGa1QxNCiG5p4cKFGAwG7rjjDrVD6RJ0Oh1OTk4A1NTUWO/Oim5Kq4W+5y/WtmF4olCHJGJCdAN29ga0N/yLWkXHyMofSf7032qHJIQQ3dKiRYt44YUXeO+99xoUTOmKEhISmD59Ov7+/mg0GrZu3dpku7Vr1xIUFISDgwNxcXEkJSW16TjFxcVERUUREBDAU089haenpw2iF6qRyondhiRiQnQTgyLiSA6aB0BIyjKK8uQPrBBCtJWbmxtz585Fq9Wyf/9+tcNpUUVFBVFRUaxdu7bZNps3b2bBggUsWbKElJQUoqKimDJlCqdPn7a2qX/+67dLXl4eAH379iUtLY2srCzefffdRnPYiW6mHZUThTokEROiG4m5cxlHdINxo4Lcd2SIohBCtIfRaMTJyYn09HS1Q2lRfHw8y5cv54Ybbmi2zYsvvsi8efOYM2cO4eHhrFu3DicnJzZs2GBtk5qaSnp6eqPF39+/QV8+Pj5ERUWxa9euDntPohO0o3KiUIckYkJ0I3b2BnQ3WoYoRlfuZu//1qkdUq+j09uR6HkjiZ43otPbqR2OEN3aqZIqfjxaxKmSqk497tNPP015eXmnJWIrV67ExcWlxeX48eNt7re2tpbk5GQmT55s3abVapk8eTK7d+9uVR8FBQWUlZUBUFJSQkJCAqGhoW2ORXQhbRyaKOc19UjVRCG6meDhcewOvp9x2f8idN9yCsdci5d/kNph9RoGByfi5r+hdhhCdBmKolBVZ2rzfv9JzmXJJwcwK6DVwF+vG85NowPa1IejnQ6NRtOmfZKTk1m3bh1Tp07ttETswQcf5NZbb22xzW/vTrVGUVERJpOp0fxyPj4+HDx4sFV95OTkcP/991uLdDzyyCNERka2ORbRhbRxaKKc19QjiZgQ3dCYO5dx+G9fE2I6QurbD+D51BdotHKDWwjR+arqTIQv/vKS+jAr8MzHB3jm4wNt2i9j2RSc7Fv/UcZsNvPAAw8wf/584uLiuOuuu6irq8POrvV3AfLy8njqqafYtGlTq/dxd3fH3d291e07U2xsLKmpqWqHIWypfmhiZRHUlIPBRd14RLPkk5sQ3ZDezh79TeuoVfREV+1h7yf/UjukXkMxmzl7+iRnT5+UZ/SE6GbWrFlDUVERy5YtIzIykrq6ulbfOarn7+/fpiQMOm5ooqenJzqdrlFxjYKCAnx9fdvcn+ghHPuCQ1/L98UX/3cl5zX1yB0xIbqp4PAx7B70AOOy1hKaupzTY67Fu79M3NjRqirLcP9nOACVTx7HycVN5YiEUJejnY6MZVPatE9+STWTX9yJ+YLpqrQa+GbBlfi6ObTp2K118uRJnnnmGd577z2cnZ0JCQnBYDCQnp5OZGQk2dnZzJgxg4iICJKSkpg8eTJTpkxh1apVVFRUsGXLFkJCQsjOzubmm2/mo48+YsaMGURHR5OUlMSIESN4//33mxwq2VFDE+3t7Rk9ejTbt2/n+uuvByx3/bZv3878+fPb3J/oQfoNhFPFluGJPuEtNpXzmnokEROiGxtzx1J+efZrhhp/Ieud+/F66ksZoiiE6FQajaZNwwMBBnm5sOrGSP7833RMioJOo2HljREM8uq4IVSPPvoo8fHxTJ06FQC9Xk9YWFiD58QyMzP54IMPGDJkCBEREbi4uJCYmMirr77KK6+8wksvvdSgz8zMTN577z3CwsKYNGkS33//PePHj2907PYOTSwvL28w11lWVhapqam4u7szYIBl0t4FCxYwa9YsYmJiiI2NZfXq1VRUVDBnzpw2H0/0IH0Hwqk0qZzYxUkiJkQ3prezx3DTOmrfv4aoqiSSPl5L7A2PqB2WEEJc1MwxA5gw1IvsokqCPJ3wc3PssGN9+umnfPvtt2RmZjbYHhkZ2SARCw0NtVYMDAsLs1YjjIyM5PPPP2/Ub2hoKOHhljsJI0eOJDs7u8lErL327t3LpEmTrOsLFiwAYNasWWzcuBGAmTNnUlhYyOLFi8nPzyc6Oppt27Y1KuAhehlrwQ5JxLoyScSE6OYGho1m96AHGZf1CsPSVnI6dpoMURRCdAt+bo4dmoDVmzZtGufOnWu0/a233mqwbjAYrN9rtVrrularxWRqXBnywvY6na7JNpdi4sSJKIpy0Xbz58+XoYiiIWsJ+2xVwxAtkzFMQvQAY+5Ywi/6obhSyam358nDtkIIIURv1jfI8lWGJnZpkogJ0QNYhygqeqKqf+Knj19ROyQhhBBCqOXCoYmtuKsq1CFDE4XoIQaGjWbP4D8w9tjLhKWuJD/mWnwDh6gdlhBCdAtBQUHs3bvXuv7RRx9Zvx87diyffvppo3YXtn/++ec7KVIhWqFvIKCBugqoKAIXL7UjEk2QO2JC9CBj7ljCIX0ofTRVnN70gAxR7AA6vR0/uf2en9x+j07f+klghRBCiE6jN0AfP8v3FxmeKOc19XSLROy7775Do9E0ufz0008ALF26tMnXnZ2dL9r/xo0bGTFiBA4ODnh7e/Pwww939FsSokPo9Hocbvk3NYodI6r38tOWl9UOqccxODgx5v9tZsz/24zBwUntcIQQQoimWYcnZrfYTM5r6ukWidhll13GqVOnGiz33XcfwcHBxMTEAPDkk082ahMeHs4tt9zSYt8vvvgif/nLX1i4cCEHDhzgm2++YcqUtk1MKURXMjA0mn1D/gBA2M9/I//EkYvsIYQQQogeRyondnnd4hkxe3t7fH19ret1dXV8/PHHPPLII9YZ7F1cXHBx+XUiyLS0NDIyMli3bl2z/Z47d46nn36a//3vf1x99dXW7SNGjOiAdyFE5xlz+2IO/e1LQo0Hydr0AD5//FomerYRxWymqrIMAEenPvJzFUII0TX1PZ+IXWRoopzX1NMtf9KffPIJZ86caXHW+Ndff52hQ4e2OLHi119/jdls5uTJk4SFhREQEMCtt97KiRMnWjx+TU0NpaWlDRYhuhLLEMVXLxii+JLaIfUYVZVlOD0/AKfnB1hPXEIIIUSX08qhiXJeU0+3TMTWr1/PlClTCAgIaPL16upqNm3axNy5c1vs59ixY5jNZlauXMnq1av56KOPOHv2LL/73e+ora1tdr9Vq1bh5uZmXQIDAy/p/QjREQaGRrMvxDLBZ/jPz5J//LDKEQkhhBCi01iHJspcYl2VqonYwoULmy3CUb8cPHiwwT65ubl8+eWXLSZZW7ZsoaysjFmzZrV4fLPZTF1dHS+//DJTpkxh7NixvPfeexw+fJgdO3Y0u9+iRYsoKSmxLhe7gyaEWsbc9jQH7cJx0VRRuOl+qaIohBBC9Bb1QxNLcsFkVDcW0SRVnxF74oknmD17dottBg0a1GD9jTfewMPDg+uuu67ZfV5//XWmTZuGj49Pi337+VnKeoaHh1u3eXl54enpyfHjx5vdz2AwYDAYWuxbiK5Ap9fjfOs6qt+5msiaFBL/8w/ibnlC7bCEEEII0dH6+IHOHky1UJr761BF0WWomoh5eXnh5dX6CeYUReGNN97gnnvuwc6u6XkOsrKy2LFjB5988slF+7v88ssBOHTokHWY49mzZykqKmLgwIGtjkuIriwwJIo9Qx9h7OEXiUj/O6dip+E3MFTtsIQQQgjRkbRa6DsAzhyxDE+URKzL6VbPiH377bdkZWVx3333Ndtmw4YN+Pn5ER8f3+i1LVu2MGzYMOv60KFDmTFjBo899hg//vgj6enpzJo1i2HDhjFp0qQOeQ9CqGHMzL+QaReOs6aaondlomchhBCiV2hl5UShjm6ViK1fv57LLrusQTJ1IbPZzMaNG5k9ezY6na7R6yUlJRw6dKjBtrfeeou4uDimTp3KlVdeiZ2dHdu2bWv2jpsQ3ZFOr8fl1nVUK3ZE1uwj6T8vqh2SEEIIITpaKysnCnV0i3nE6r377rstvq7ValssnDF79uxGz6S5urqyfv161q9fb4sQheiyAkOi2BP6GGN/eZ6I9OfIGzMd/yAZothWWp2eFJcJAITrutWfUCFEL1NcXMzkyZMxGo0YjUYee+wx5s2bp3ZYojO1onKinNfU063uiAkhLk3szD+TaTccZ001Z9+7H7PJpHZI3Y6DozOjnvwfo578Hw6OzmqHI3q5G264gX79+nHzzTc3eu3TTz8lNDSUkJAQXn/9dRWiE2rr06cPCQkJpKamkpiYyMqVKzlz5ozaYYnO1IqhiXJeU48kYkL0Ilqdjj4z/02VYk9ETSo//ecFtUMSQlyCxx57jLfeeqvRdqPRyIIFC/j222/Zt28fzz33nHwAV9mZM2fw9vYmOzu7046p0+lwcnICoKamBkVRUBTF+vptt93GCy/IeaBHk6GJXZokYkL0MgFDIkgLfQyAyAPPk5d18CJ7CCG6qokTJ9KnT59G25OSkhg+fDj9+/fHxcWF+Ph4vvrqKxUiFPVWrFjBjBkzCAoKAiAhIYHp06fj7++PRqNh69atTe63du1agoKCcHBwIC4ujqSkpDYdt7i4mKioKAICAnjqqafw9PS0vvb000+zYsUKSkpK2vu2RFdXPzSxohBqK9SNRTQiiZgQvVDszEVk2EXgpKnhnAxRbJPK8hJY6gZL3SzfC9GM1nzQvtQP2c3Jy8ujf//+1vX+/ftz8uRJm/Qt2q6yspL169czd+5c67aKigqioqJYu3Zts/tt3ryZBQsWsGTJElJSUoiKimLKlCmcPn3a2iY6OpqIiIhGS15eHgB9+/YlLS2NrKws3n33XQoKCqz7RkREMHjwYN55550OeNeiS3DsBwY3y/fFTc+RK+c19UgiJkQvpNXpcJ35KlWKPcNr0/jpo+fVDkmIHudiH7Rt8SFbtM0tt9yCl5cX//73v63bEhMTsbe379A7hp9//jkGg4GxY8dat8XHx7N8+XJuuOGGZvd78cUXmTdvHnPmzCE8PJx169bh5OTEhg0brG1SU1NJT09vtPj7+zfoy8fHh6ioKHbt2tVg+/Tp03n//fdt9E5Fl2Qt2JGtahiiMUnEhOilAoZEkDbs/wEQmfECJ49lqhyRED3LxT5o2/JD9m/5+/s3uAN28uTJFvepqamhtLS0wdITvfzyy9x0000sW7YMgPLycu666y4eeughrrnmmg477q5duxg9enSb9qmtrSU5OZnJkydbt2m1WiZPnszu3btb1UdBQQFlZWWAZQqfhIQEQkMbVsuNjY0lKSmJmpqaNsUnupFWVE4U6pBETIheLPbWP3HAPhInTQ3F78+TIYpCdBJbfMhuSWxsLOnp6Zw8eZLy8nK++OILpkyZ0mz7VatW4ebmZl0CAwPbfMzKWmOzS3WdyeZt28PPz4/HH3+ckydPcubMGR599FEMBgPPPvtsu/prrZycnIsmz79VVFSEyWTCx8enwXYfHx/y8/Nbfdzx48cTFRXF+PHjeeSRR4iMjGzQxt/fn9ra2lb3KbohmdS5y5LJAoToxbQ6HX1ve43KN69keO1+Ej/8O3G3LVI7LCF6vJY+ZB882PoCOpMnTyYtLY2KigoCAgL48MMPGTduHHq9nhdeeIFJkyZhNpv54x//iIeHR7P9LFq0iAULFljXS0tL25yMhS/+stnXJoV68cacWOv66P/7hqq6pi/8xAW7s/mBcdb1K57dwdmK2kbtsv82tU3x1Rs6dChOTk4sXryYTZs2kZSUhIODQ7v6aq2qqqoOP0ZTYmNjSU1NbbGNo6MjYHmOTfRQUjmxy5JETIherv+gMBLD/h9xB/9GZOY/OHnsOvoPGq52WEKIVvjmm2+afe26667juuuua1U/BoMBg8Fgq7C6NK1WS2RkJP/85z/5+9//TlRUVIcf09PTk3PnzrV5H51O16C4BliGG/r6+tostrNnzwLg5eVlsz5FF2NNxOSOWFcjiZgQgjG3/JEDz37O8NqfyX7/Afz+tBOtTqd2WEL0WJ31IbszZSxrfuijVqNpsJ78zORmWjZu+/2fJl1aYL9RP4/WqFGjeOKJJ6zbs7OzmTFjBhERESQlJTF58mSmTJnCqlWrqKioYMuWLYSEhAAwbdo0Tp06RU1NDYsWLeLOO+9k9+7dPPbYY/z444+cOXOGK664gl27duHr68vIkSPbXJnQ3t6e0aNHs337dq6//noAzGYz27dvZ/78+bb5YQDp6ekEBAQ0KGsvepgLhyYqCvzm/5hQjzwjJoQ4P0Tx31QqBsJr95P0Qcc+L9GdaXV60hxjSXOMRauTa1mifS78kF2v/kP2uHHjWtiz63Ky1ze7ONjpbN62vVavXk1iYiJmsxmttuHHoMzMTBYvXszBgwf57rvv+OGHH0hMTOSRRx7hlVdesbZ76623SE5OJjExkRUrVlBTU8O4ceOYMGECzz77LA8//DCLFy+2JtVTpkzhwIEDDe6KlZeXk5qaah06mJWVRWpqKseP/1pifMGCBbz22mu8+eabZGZm8tBDD1FRUcGcOXPa/f5/a9euXR1aqER0AX0HWL7WlkNl44nd5bymHvlpCyGA80MUw58gLnMlIw6uJvfIdQQMiVA7rC7HwdGZqD99rXYYohsoLy/nyJEj1vX6D9ru7u4MGDCABQsWMGvWLGJiYoiNjWX16tU2/5AtGtq/fz+LFi3iD3/4A6+//jpGoxG9/tePQqGhodaqgmFhYdZiKpGRkXz++efWdv/4xz/45JNPADh+/DjHjx8nJCSE5cuXEx0dzZAhQ7j77rut7SMjIxk1ahQffPABDzzwAAB79+5l0qRf7/bVP6M3a9YsNm7cCMDMmTMpLCxk8eLF5OfnEx0dzbZt2xo9W9he1dXVbN26lW3bttmkP9FF2TlAHz8oO2UZnujc8O6nnNfUI3fEhBBWY25+kgP2UThpaijd/IBUURTiEuzdu5eRI0cycuRIwPJBe+TIkSxevBiwfMh+/vnnWbx4MdHR0aSmptr0Q7ZoqLq6mjvuuIOZM2eyfPlyamtrGxVGufA5Oa1Wa13XarWYzv893LFjh/VOWVpaGsOGDbOWfj99+jS1tbXWYiwXWrx4MS+99BJmsxmAiRMnoihKo6U+Cas3f/58cnJyqKmpITExkbi4OJv9TN544w1iY2MbzG8meijr8MRsVcMQDUkiJoSw0up09Lv9/BDFunSSNq9SOyQhuq3WfNDuyA/ZoqGFCxdSUVHBK6+8Qr9+/Rg4cCCrV69u8+TYpaWleHh44ODgQGpqKmlpadbX5s2bx5o1axgzZgwvvPBCg/2mTp3K/fff32B+N7XZ2dmxZs0atcMQnUEqJ3ZJkogJIRrwDx7G/uFPAhB16CVOHNmvckRdS2V5CZVLvC1LeYna4QghWuGrr75i7dq1vPPOO/Tp0weAp59+mq1bt/Lwww+3qa/f//73lJWVER4ezooVK6wTNa9fvx5vb2+mTp3K3/72N958800OHTrUYN/HH3+8XXO0dZT77ruv0QTPoodqYVJnOa+pR6PUlw8S7VZaWoqbmxslJSW4urqqHY4Ql8xsMpHx96uIqEkl0244oQt3SRXF8yrLS3B63vLgc+WTx3FycVM5ou5P/oZ2PS39Tqqrq8nKyiI4OFiVubGEuuT3303t2wQf/wEGTYR7Pm7wkpzXbK+15zW5IyaEaESr0+F++7+pUBwIqztA4sY/kv7D/yjIPap2aEIIIYRoKxma2CVJIiaEaJJ/UCjpEU8BMPb460R8fReer40m6T+r1Q1MCCGEEG1TPzSxJBfMUoirq5BETAjRrIFxMxrM/ajTKIz6+a9yZ0wIIYToTvr4gc4ezEYo7ToFY3o7ScSEEM0qOnHQmoTV02vM5KbtVCcgIYQQQrSdVgdu5wvFyPDELkMSMSFEs7wGhmNSNI22RyY9SeKae8jLPtTEXkIIIYToclqonCjUIYmYEKJZPgGDSR6xFKNi+VNhUjQc1/hjrzERd+ZjvN4YR9Lq28k9ekDlSDuPVqvjgH0kB+wj0WqlkqQQQohuor5gR3HDREzOa+rRqx2AEKJri73pcQriplOUcxDPgcMYEDCYjN1fYPruWSJr9hFb/DnGt7bxU9/JeF/7FwaGRqsdcodycHJh+J+/VzsMIYQQom361t8Ry26wWc5r6pFETAhxUT4Bg/EJGGxdDx8XD+PiOfjTN9Rs/xtR1T8xpuQrzO9+zV7Xq/CI/zPB4TEqRiyEEEKIBmRoYpcjQxOFEO02bMxkohZ+w+EZn7DP6TK0GoWYsu0Ef3A1Kc9N5+jPP6odohBCCCGg2aGJQj2SiAkhLlnIyCsZ+ccvOHrjF6Q4TwBgVEUCg/8bz76//55fUhJUjtB2KstLOLc0kHNLA6ksL1E7HCGEEKJ16ocmlhdAbaV1s5zX1COJmBDCZgaPuIxRT/2P7Fu/IbnPVZgVDSMrdzP0k+mk/W0yB3/6Ru0QbaIfpfSjVO0whBBCiNZz7AcGV8v3xccbvCTnNXVIIiaEsLmg8DGMfmILuXd+x09u12BSNERV/8Swz25i/6qJHNj9hdohCiGEEL2LRvPrc2IyPLFLkERMCNFhBgyNZsz/+5D8e77np75TqVN0RNbsY/iXt5G+cjz7d32CYjarHaYQQvRqWVlZTJo0ifDwcCIjI6moqFA7JNFRmqmcKNQhiZgQosP1HxzBmMffpeje3SR5zKBW0RFR+zOR2+/m4KrLSPvuP5KQCSGESmbPns2yZcvIyMhg586dGAwGtUMSHaW+YIdUTuwSJBETQnQav4GhxD7yFufmJZHodTM1ih1hdZlEfXcvh1fGse+b9yUhE0J0qIULF2IwGLjjjjvUDqVLOHDgAHZ2dowfPx4Ad3d39HqZ3ajHksqJXYokYkKITucTMIS4h9dT9mAyiT63UaXYM9T4CyO/f4CjK2JI3vYWZpNJ7TCFED3QokWLeOGFF3jvvfc4cuSI2uG0KCEhgenTp+Pv749Go2Hr1q1Ntlu7di1BQUE4ODgQFxdHUlJSq49x+PBhXFxcmD59OqNGjWLlypU2il50STI0sUuRREwIoRpPv4HEPfQqVX/Yxx6/u6lUDAwxHWX0nkfIWTGSnz57HZPRqHaYDWi1Og7rQzisD0Gr1akdjhCijdzc3Jg7dy5arZb9+/erHU6LKioqiIqKYu3atc222bx5MwsWLGDJkiWkpKQQFRXFlClTOH36tLVNdHQ0ERERjZa8vDyMRiO7du3in//8J7t37+brr7/m66+/7oy3J9Rw4aTOigLIeU1Ncu9ZCKE6d58Axj7wCsWFf+Hnj59l+In3CTbnEPzTE+Qk/4OCqEcYde29nCk4QWFOBl4Dw/EJGKxKrA5OLoQ8vVeVYwvR45SchLNHwX0wuPXvtMMajUacnJxIT0/nhhtu6LTjtlV8fDzx8fEttnnxxReZN28ec+bMAWDdunV89tlnbNiwgYULFwKQmpra7P79+/cnJiaGwMBAAK699lpSU1P53e9+Z5s3IbqWvgMsX2vLoPIsOHvIeU1FckdMCNFl9PXyY+x9qzE/vp89A+6nFGcGmnOJ3fcnSlYMxuu10UR8fReer40m6T+r1Q5XCAGWq+q1FW1fkl6D1RHw5nTL16TX2t7H+Sv6bfX0009TXl5Oenq6jX8YTVu5ciUuLi4tLsePH794R79RW1tLcnIykydPtm7TarVMnjyZ3bt3t6qPMWPGcPr0ac6dO4fZbCYhIYGwsLA2xyK6CTtHcPG1fF+crWooQu6ICSG6ILd+Xoy99znKShaxZ+vzDMvaiAeloLG8rtMojPr5rxTETVftzpgQ4ry6Sljpf2l9KGb4/EnL0hZ/zgN75zbtkpyczLp165g6dWqnJWIPPvggt956a4tt/P3b/jMsKirCZDLh4+PTYLuPjw8HDx5sVR96vZ6VK1cyYcIEFEXhmmuuYdq0aW2ORXQj/QZCeb5leGL/0WpH06tJIiaE6LL6uLkzdtZK0rZH0XfX/Q1e02vMFOUc7PRErKqijOLnRgLQ96l9ODr36dTjCyHaz2w288ADDzB//nzi4uK46667qKurw87OrtV95OXl8dRTT7Fp06ZW7+Pu7o67u3t7Qu4UrRkCKXqQfkFwItFasEPOa+qRREwI0eX5hsZgStCg0/w6DElRwNkzoNNjURQzfhQCUKlIqX0hsHOy3Jlqi9I8WBtruRNWT6ODhxPBtQ13huyc2nTYNWvWUFRUxLJlyzh+/Dh1dXUcPHiQyMjIVvfh7+/fpiQMLEMTL1aNMCMjgwEDBrSpX09PT3Q6HQUFBQ22FxQU4Ovr26a+RC9SXznxfAl7Oa+pR54RE0J0eT4Bg0kesRSjYvmTpSig0UDpf/8ftTXVKkcnRC+n0ViGB7Zl8QyB6S9Zki+wfJ2+2rK9Lf1oNK0O8+TJkzzzzDOsXbsWZ2dnQkJCMBgM1uGJ2dnZREVFceeddxISEsJDDz3E1q1biYuLIyIigsOHD1vbxcTEWNvPmjWLsLAwZs6cidLMM2sPPvggqampLS7tGZpob2/P6NGj2b59u3Wb2Wxm+/btjBs3rs39iV7iwsqJQlVyR0wI0S3E3vQ4BXHTKco5SF1VKUN3PcaImmT2vnInox7/AK1OSu4K0a2MugcGXw1nj4H7oA6vmvjoo48SHx/P1KlTAcuzUWFhYQ2eE8vMzOSDDz5gyJAhRERE4OLiQmJiIq+++iqvvPIKL730UoM+MzMzee+99wgLC2PSpEl8//331omRL9TeoYnl5eUN5jrLysoiNTUVd3d3692zBQsWMGvWLGJiYoiNjWX16tVUVFRYqygK0Uj9pM4yl5jqJBETQnQbPgGDrc+E7dfbMWzH/cSUfcPuVx9i7IPr0GjlJr8Q3Ypb/04pW//pp5/y7bffkpmZ2WB7ZGRkg0QsNDSU0NBQAMLCwqzVCCMjI/n8888b9RsaGkp4eDgAI0eOJDs7u8lErL327t3LpEmTrOsLFiwAYNasWWzcuBGAmTNnUlhYyOLFi8nPzyc6Oppt27Y1KuAhhFX90MSSXDCb1I2ll5NETAjRLUVOvJm9pYXEpCxk3OnN7H7Hl3H3LFM7LCFEFzRt2jTOnTvXaPtbb73VYN1gMFi/12q11nWtVovJ1PgD64XtdTpdk20uxcSJE5sd7nih+fPnM3/+fJseW/Rgrv6gtQNzneV5Tb2r2hH1WnL5WAjRbcVc9xCJQ/4fAOOOvUTSlldUjkgIIYTo4rQ66GuZwFuGJ6pLEjEhRLcWd9dS9vjdCcCo1GdI+3Zzhx5Po9GSrQ0kWxuIRiN/QoX6brjhBvr168fNN9/cYPuJEyeYOHEi4eHhjBgxgg8//FClCIUQXc4FlRPlvKYejdKae96iRaWlpbi5uVFSUoKrq9zeFaKzKWYTe1+6jTElX1GpGDg+7T2Gjbla7bBEK8nf0Evz3XffUVb2/9u787AmrvUP4N+ENSyCgCwRWVwRRUBAKnpdKhWtRa2tYrUoar2tgJWqVGsrUuuK+4KlbmirXqvtD9yqYr2A1ovghgviUouIVcAF2SFAzu8PykgkYdEwYXk/z5OnzZk3M2+ScV5O5syZfOzevRu//PIL1/748WNkZWXByckJmZmZcHFxwZ07d6CrW/cNkGv7TkpKSpCWlgZbW1toa2sr/f2Qpo2+/xbiyCzg0i5gwJfA21+rOpsWp751jbq9hJBmTyBUg1PAHlzTdoOOoBTmxyYh/dZlVadFCC8GDRoEff2aN2C1sLCAk5MTAMDc3BwmJiZ4/vw5z9kRQpokmjmxSaCOGCGkRdDQ1ELnwF9wR70rDFEArf1jkfXwnqrTIq3cmTNn4O3tDbFYDIFAgOjo6Box4eHhsLGxgba2Ntzd3ZGUlKT0PC5duoSKigp06NBB6esmhDRDr9zUmagGdcQIIS2Gjp4hTD89jAfC9jDHUxTtHI3c59lK3UZxYT7uL+6J+4t7orgwX6nrJi1PYWEhHB0dER4eLnf5zz//jNmzZ2PRokW4fPkyHB0d4eXlhezsl/utk5MTevbsWePx6NGjeuXw/PlzTJo0CVu3blXKeyKEtADVbupMdU11aPp6QkiLYtjOAiV+0XiycyhspQ+Q+v1oaAbFQKSrp5T1MyaFjTQDAFDEpEpZJ2m5hg8fjuHDhytcvnbtWkyfPp27+W5ERASOHTuGnTt3Yv78+QCA5OTk195+aWkpRo8ejfnz58PDw6PWuNLSUu55Xl7ea2+TENIMtLWt/G9BJpikiOqaitAZMUJIi2Nu1RWFY39GHnTQvSwFt8I/RHmZRNVpESJDIpHg0qVL3E2Dgcr7VXl6eiIhIeGN188Yg5+fH95++234+vrWGrt8+XIYGBhwDxrCSEgLJ2oLaFZeWyrIy1BxMq0XdcQIIS2STQ93/O21E6VMA85FCbi8ZQqYlH7pI03H06dPUVFRATMzM5l2MzMzZGZm1ns9np6eGDt2LH777TdYWlpynbhz587h559/RnR0NJycnODk5ITr16/LXcdXX32F3Nxc7pGRQX+YEdKiCQTc8ERBLv17VxUamkgIabG69x2OK/nr0etcIPrkHEXCjjnoO32dqtMiRKl+//13ue39+/eHtJ4/PmhpaUFLS0uZaRFCmrq2NkDWDQhfPFB1Jq0WnREjhLRozkM/xiWHhQCAvn/vxPn9y1WcESGVTExMoKamhqysLJn2rKwsmJubqygrQkirYUhnxFSNOmKEkBavz4dzcN7q08r/T12JS7/tVHFGhACamppwcXHB6dOnuTapVIrTp0+jb9++KsyMENIqcEMT6YyYqjSLjlhcXBwEAoHcx4ULFwAAoaGhcpfr6urWuu4LFy5gyJAhMDQ0RNu2beHl5YWrV6/y8bYIITxy91uBJJP3IRQwOCQG48YfR15rPQKBEI/RDo/RDgJBsziEEhUqKChAcnIyN/NhWloakpOT8eBB5R8+s2fPxrZt27B7926kpqZixowZKCws5GZRJISQRvPPTZ2FuRlU11SkWVwj5uHhgcePH8u0LVy4EKdPn4arqysAYO7cufjss89kYoYMGQI3NzeF6y0oKMCwYcMwcuRIbNmyBeXl5Vi0aBG8vLyQkZEBDQ0N5b8ZQohKCIRCuHy2HZfXP0PvgjOwOTUdf+r/is6O/Rq0HpGuPkShfzZSlqSluXjxIgYPHsw9nz17NgBg8uTJ2LVrF3x8fPDkyROEhIQgMzMTTk5OOHHiRI0JPAghROn+GZoozHsIi0UPKifwILxqFh0xTU1NmfHyZWVlOHToEGbOnAnBPzuNnp4e9PRe3ifo6tWruHnzJiIiIhSu99atW3j+/DkWL17MTdW7aNEi9OrVC+np6ejcuXMjvSNCiCqoqaujR+B+3Fw7DPaSazCMmoC/9Y+jfUd7VadGWqhBgwaBMVZrTGBgIAIDA3nKiLQ2aWlpmDp1KrKysqCmpobz58/XOVqItBKGVpX/Lc0DinMAHSPV5tMKNcvzj4cPH8azZ89qHbqxfft2dO3aFf/6178UxnTr1g3GxsbYsWMHJBIJiouLsWPHDnTv3h02NjYKX1daWoq8vDyZByGkedDS1oWlfzTuqdnCBC/AfnofTzPpQmVCSMvk5+eHxYsX4+bNm4iPj6fZMclLmjqA3j9n33PuqzSV1qpZdsR27NgBLy8vWFpayl1eUlKCvXv3Ytq0abWuR19fH3FxcdizZw9EIhH09PRw4sQJHD9+HOrqik8W0o0vCWne2hgaw+CTQ3gkMIUly8SLbaOQn/u8Xq8tKSrA3SWuuLvEFSVFBY2cKSGkpXj27BlMTU1x//593raZkpICDQ0N7kdpIyMjmb9vxo8fjzVr1vCWD2mC/hme+HjHR1TXVEClHbH58+crnISj6nHr1i2Z1zx8+BAnT56stZMVFRWF/Px8TJ48udbtFxcXY9q0aejXrx/Onz+Pc+fOoWfPnhgxYgSKi4sVvo5ufElI82diYQ3pxP/Dc7RB54p7uL9lDEpLFf+7ryKVVqBL+V10Kb8LqbSCh0wJIS3B0qVLMWrUKG7EzZkzZ+Dt7Q2xWAyBQIDo6Gi5rwsPD4eNjQ20tbXh7u6OpKSkem/z7t270NPTg7e3N3r37o1ly5bJLP/mm2+wdOlS5Obmvu7bIs3dPxN2WEgzqa6pgEqvEZszZw78/PxqjenYsaPM88jISBgbG2PkyJEKX7N9+3a89957dV7svG/fPty/fx8JCQkQCoVcW9u2bXHo0CGMHz9e7uvoxpeEtAyWnR3w56i90I7+EA6lV3Bx8wT0DvoFQjU1VadGCGlBioqKsGPHDpw8eZJrKywshKOjI6ZOnYoxY8bIfd3PP/+M2bNnIyIiAu7u7li/fj28vLxw+/ZtmJqaAgCcnJxQXl5e47UxMTEoLy/H2bNnkZycDFNTUwwbNgxubm545513AAA9e/ZEp06dsGfPHgQEBDTCOydN3j9T2BPVUOkZsXbt2sHOzq7Wh6amJhfPGENkZCQmTZqkcEbDtLQ0xMbG1jksEag8MAqFQm7CDwDcc6lU+uZvkBDS5HV2HoB7QyJQxtTgmv9fJEV8Ckb//glpkcaOHYt27dph69atXFtiYiI0NTURExPTaNv97bffoKWlhbfeeotrGz58OJYsWYL3339f4evWrl2L6dOnY8qUKbC3t0dERAR0dHSwc+fLeyEmJyfjxo0bNR5isRjt27eHq6srOnToAC0tLbz77rvcrRSqeHt7Y//+/Up/z6SZMKSOmCo1q2vE/vvf/yItLQ2ffPKJwpidO3fCwsICw4cPr7EsKioKdnZ23PN33nkHOTk5CAgIQGpqKlJSUjBlyhSoq6vLTDdMCGnZHAaMwTW3FQCAt54cxPmfFqo4I0KaIUmh4kdZSQNii+sX+xo2btyIDz74AIsXLwZQeRubjz/+GDNmzMDQoUNfa531cfbsWbi4uDToNRKJBJcuXYKnpyfXJhQK4enpiYSEhHqtw83NDdnZ2cjJyYFUKsWZM2fQvXt3mZg+ffogKSkJpaWlDcqPtBD/DE0kqtEspq+vsmPHDnh4eMh0pqqTSqXYtWsX/Pz8oCZnaFFubi5u377NPbezs8ORI0fw7bffom/fvhAKhXB2dsaJEydgYWHRaO+DENL0uLz3byTmZsL97hr0TduMpP8zQ58xn6s6LUKaj2Vixcu6DAUmHnz5fFVnoKxIfqx1f2DKsZfP1zsARc9qxoU2/LomCwsLBAUF4YcffsCzZ88QHBwMLS0trFy5ssHraoj09HSIxbV8PnI8ffoUFRUVNS6zMDMzq3H9vCLq6upYtmwZBgwYAMYYhg4divfee08mRiwWQyKRIDMzE9bWdHak1aGhiSrVrDpi+/btq3W5UCisdeIMPz+/GtekvfPOO9xYaUJI6+Y+MQTnf8jGW49/Qu+ri5DcxhROnvKvFSWENE9du3aFjo4OQkJCsHfvXiQlJUFbW7tRt1lcXNzo21Bk+PDhckcJVRGJRAAqL9cgrVCb9mACNQgYTdKhCs2qI0YIIY3NffpGXNz4BK4vTqDb2Zm41cYYdn1kf6zJQRsAAE3ZQ0g1Cx4pXiZ4ZZRK8J+1xL5y1UTQ9dfPSQ6hUAgHBwds2bIFYWFhcHR0VOr65TExMUFOTk6DX6OmpoasrCyZ9qysLJibmystt+fPK2/d0a5dO6WtkzQjQjUwA0sIXqQjHzqgqar41ayuESOEkMYmEArh5P8jror6QCSQwOK3ybifeolbrqNngLahGWgbmgEdPQMVZkpIE6Opq/ihod2AWFH9Yl8TYwwA0Lt3b8yZM4drv3//PhwdHTFx4kR06dIFM2bMQHR0NNzd3dGzZ0/cvXuXi33vvffg4uKCnj17Yu/evQCAhIQE9OnTB+Xl5cjKykKXLl2QmZkJAHB2dsbNmzcblKempiZcXFxw+vRprk0qleL06dPo27fva7//V924cQOWlpYwMTFR2jpJ8yI0sgUA6I9aTXWNZ9QRI4SQV6hraqFrwK+4rd4NBiiE6OexyMyo5Rd8QkizsX79eiQmJkIqlXK3rqmSmpqKkJAQ3Lp1C3FxcTh37hwSExMxc+ZMbN68mYv78ccfcenSJSQmJmLp0qUoLS1F3759MWDAAKxcuRIBAQEICQnhzlx5eXkhJSVF5qxYQUEBkpOTuVkM09LSkJycjAcPHnAxs2fPxrZt27B7926kpqZixowZKCwsxJQpU5T2eZw9e7ZRJyohzUDVzIkv0lWbRytEHTFCCJFDpNcGZp8dRrrQEmZ4hpLI0ch9llX3CwkhTdb169fx1Vdfwd/fHzdv3qxx/61u3bqhW7duUFNTQ/fu3bkZCx0cHHD//n0ubt26dXB0dISHhwcePHjAdZ6WLFmCn376CSUlJfD19eXiHRwc0Lt3bxw4cIBru3jxIpydneHs7AygstPl7OyMkJAQLsbHxwerV69GSEgInJyckJycjBMnTtR5n9T6KikpQXR0NKZPn66U9ZFmqmrmxBzqiPGNOmKEEKKAoYk5tP2ikQ0j2Egz8DhiFHKeZiJlWX+kLOuPkqICVadICKmnkpISTJgwAT4+PliyZAkkEkmN2Qe1tF5e+SkUCrnnQqEQFRWVkxnExsZyZ8quXr0KOzs7bur37OxsSCQSbsbD6kJCQrBhwwbuPqWDBg0CY6zGY9euXTKvCwwMRHp6OkpLS5GYmAh3d3elfSaRkZHo06ePzP3NSOsj0a3s2BddP0x1jWfUESOEkFqYWXVB0bgDyIMu7MpSkbHtY/SQXEcPyXVIpTTLFCHNxfz581FYWIjNmzejbdu2sLa2xvr16/HoUS2TjMiRl5cHY2NjaGtrIzk5GVevXuWWTZ8+HZs2bYKbmxvWrFkj87oRI0bg3//+N/7++2+lvB9l0NDQwKZNm1SdBlExqX57AIAOK6K6xjPqiBFCSB1s7N3waHgkSpgGepVe4NqfPrqvuqQAZD28hxvnjiDr4T3Kg5BaxMTEIDw8HHv27IG+vj4A4JtvvkF0dDQCAgIatK5hw4YhPz8f9vb2WLp0KXej5h07dsDU1BQjRozAihUrsHv3bpl7lwJAUFAQOnTooJw3pQSffPIJunXrpuo0iIpJDV/eS0zw4r7qEgGA3L+BtDOV/20FeQhY1fRB5LXl5eXBwMAAjx49Qps2bWosV1NTk7l/SGFhocJ1CYVC7p4eDY0tKiqCoq9TIBBAR0fntWKLi4u5oRTy6OrqvlZsSUlJjaEbrxuro6MDgUAAACgtLa0x7v91Y0UiEXcxt0QiQVlZmVJitbW1uZuONyS2rKwMEolEYayWlhbU1dUbHFteXs4NrZFHU1MTGhoaDY6tqKhASUmJwlgNDQ1oamo2OFYqlaK4uFgpserq6tzwI8ZYrffSuRl3EC5JsyAUVMYWSICLRt7Q7DKwRqxQKIDWPzkAQHGJ4s+sIbECgQDaWpqQ3IlFn2fRKC6TQsoEuGj0Xo08qmLrs14AEGm/HJZV31jJnVg4ZEYBYKhgAlzpuQCuowNlYus6RuTl5UEsFiM3N1fuMZTwr6quyftOSkpKkJaWBltbW5XdG4uoDn3/LU/p/yKgFTMPAMAACFynAR0H8Z/IX3HAxZ1VWQCuU1Wfh0AIeG8Aek9q0CpqO4ZWRx0xJaj6sBV59913cezYMe65rq6uwj/2Bg4ciLi4OO55u3bt8PTpU7mxrq6uuHDh5a/zNjY2SE+Xf6Glvb09UlJSuOc9evRQOJWutbW1zEXJbm5uuHjxotxYExMTPHnyhHs+aNAgxMfHy43V0dGR6ViOGDECv/32m9xYADIdxbFjx+KXX35RGFtQUMD9sefn54fdu3crjM3OzubulxIQEIAtW7YojE1LS4ONjQ0AIDg4GKtXr1YYe+PGDfTo0QMAEBoaim+//VZhbFJSEtzc3AAAq1atwpdffqkwNjY2FoMGDQIAhIeHIzAwUGHs0aNHMWLECADArl27ap1Z68CBAxg7diwA4ODBgxg3bpzC2MjISO5m6MeOHcN7772nMHbz5s3cL8xxcXEYPHiwwtiwsDAEBwcDAC5cuIA+ffoojF20aBFCQ0MBACkpKejZs6fC2Llz52LVqlUAKqektrW1VRjr7++P8PBwAMCTJ09gamqqMHbch2PwH/vfIRQAhRIGveX5CmM/tFfHwbEvf9AQfJunMPbdLuo4NuFlrO6yPBQp6JsPtFZDnN/Ljk27Vfl4WiT/MO4qFuLCdD3uuc36fKTnyo+1bydEiv/L2B5bCnDzifwfVawNBLgfpM89d9tWgIuP5Mc25BhBHbGmgzpiRBH6/luY3L/B1vWAANQdUEigVnk/Q4P29X5JfTtidENnQgipp5LCFxAK6hebL2iDWxpdq7WcVxhbINDDLQ077jlDEgD5HZsigQ7ShFawlT6Qu1wmX4EItzTsuedlgssA5J8lLRVoycSWCq4CkH8msUyggVsa9tCqKKhXHoQQQpqo5/fkd8LadQe0ebynWEku8CS1aebBKoDnfzWoI1ZfdEZMCWhoIg1NbGgsDU2s1NyGJj7LTEf7H/tC7Z+hiUVlQDkT4tmUszBr31EmtiH/7ht6jMh79ggm21ygJmAolFT+O5aXR2MfI7Ie3oPJNhdIyqWQ/hMqLw8amtj80Bkxogh9/y2MvDNir3EGSBl5YH1PgFWrD804DzojpgK6uroyf3DUFteQddZX9T+ilBlb/Q85ZcY25ADekFgtLS2ZKYiVFaupqcn9ca+qWA0NDa6To8xYdXV1rlOmzFg1NbV678MNiRUKhY0SKxAIao3V7WSPRPuv0OfmcggEAmhpCJDSKxR9ujrUuW5l/rsXWXZCUq9Q9L72LXQ1pShnQqT0WlRnHso+RphVy0NdUL885B0javuRhRBCSCMyaA/J0JXQPPklBAKAQQiB93p+Oz//5AHvDcCRoMozUAI1oBXkQWfElKC+vV5CSMuQ9fAenqbfgom1HcwsO1Eeb5gHHUObHjojRhSh77+Fyv27cvidUUf+Oz8tMA86I0YIIY3EzLKTSjs+lAdpCuh33NaJvvcWyqC9ajs+rTQPuo8YIYQQQuqt+vWtpPWpuo62vkPfCSGK0RkxQghpgJLiQtzeOBoA0O3zaGiL6n/tFyEtgbq6OnR0dPDkyRNoaGhwkxSRlq1qMqPs7GwYGhpyHXLS/FFdUx3qiClRYWEh9PX1uRn5qmbDqz4rW1UcIDvLXtUMd4pmT6tPbNUsZ9Vn2aua4U7R7Gn1ia2a5az6LHtVM9w1JPbVWdmqZkKUN8tefWKrz4ZXfXKDqpkQFc2cV1ds9Znzqs+wWPV9NiS2Pt+9MvYTed+nMvaTqu/zTfcTRTNxvu5+ouj7fNP9pPr3qShWWlEOx+IkMMbwNC8XFdKGffevu5+09GMEaT4EAgEsLCyQlpam8N6VpOUyNDSEubm5qtMgSlRV1wCgqELxTNKkETDyxnJzcxkqbwPOsrOzufYlS5YwAOyTTz6RidfR0WEAWFpaGte2bt06BoBNmDBBJtbExIQBYDdu3ODatm7dygCwUaNGycRaW1szACwpKYlr27NnDwPAPD09ZWLt7e0ZABYbG8u1RUVFMQDMw8NDJtbV1ZUBYEePHuXaYmJiGADm6OgoEztw4EAGgB04cIBr++OPPxgA1rlzZ5nYd999lwFgkZGRXNuVK1cYACYWi2ViP/zwQwaAbd68mWu7c+cOA8AMDAxkYidPnswAsLCwMK7t4cOHDABTV1eXifX392cA2KJFi7i2nJwc7vuUSCRc+9y5cxkANnfuXK5NIpFwsTk5OVz7okWLGADm7+8vsz11dXUGgD18+JBrCwsLYwDY5MmTZWINDAwYAHbnzh2ubfPmzQwA+/DDD2VixWIxA8CuXLnCtUVGRjIA7N1335WJ7dy5MwPA/vjjD67twIEDDAAbOHCgTKyjoyMDwGJiYri2o0ePMgDM1dVVJtbDw4MBYFFRUVxbbGwsA8Ds7e1lYj09PRkAtmfPHq4tKSmJAWDW1tYysaNGjWIA2NatW7m2GzduMADMxMREJnbChAkMAFu3bh3XlpaWxgAwHR0dmdhPPvmEAWBLlizh2rKzs7nvs7pZs2YxAGzBggWsMP8FY4vasIKv9LnYgoICLnbBggUMAJs1a5bMOugYUUneMeLkyZMMAMvNzWWkaaiqa7V9JxUVFay4uJgerehRXl7O415I+FJV19iiNpX/T95YfY6hjDFGZ8QIIYQQ0mBCoZBmzSOEkDdA09crQfUbOpubm9OwIxqaSEMTW/DQxHJJMXRWW1UOTfS/BR09Axqa+IbHiJycHBgZGdH09U0I3VKAkNajqCAXOqutKv9/7gPo6BmoOKPmr77HUOqIKQEVLEJaDypYykfH0KaHvhNCWg+qa8pX32MoTXVECCGEEEIIITyja8SUoOqkYl5enoozIYQ0tqKCPJSXVv6bL8rLQ7lUoOKMmr+qYycN0Gg6qK4R0npQXVO++tY1GpqoBA8fPkSHDh1UnQYhhDRrGRkZsLS0VHUaBFTXCCFEGeqqa9QRUwKpVIpHjx7J3EOsvvLy8tChQwdkZGSodBw+5UF5UB6Uh6ryYIwhPz8fYrGYbg7cRFBdozwoD8qD8mj8ukZDE5VAKBS+8a+4bdq0aRIXRFMelAflQXmoIg8DA7o4vCmhukZ5UB6UB+XR+HWNfnokhBBCCCGEEJ5RR4wQQgghhBBCeEYdMRXT0tLCokWLZG7mSnlQHpQH5UF5kOaqqewPlAflQXlQHk09D5qsgxBCCCGEEEJ4RmfECCGEEEIIIYRn1BEjhBBCCCGEEJ5RR4wQQgghhBBCeEYdMUIIIYQQQgjhGXXEVGT58uVwc3ODvr4+TE1NMXr0aNy+fVulOa1YsQICgQBBQUEq2f7ff/+Njz/+GMbGxhCJRHBwcMDFixd5zaGiogILFy6Era0tRCIROnXqhO+++w6NPafNmTNn4O3tDbFYDIFAgOjoaJnljDGEhITAwsICIpEInp6euHv3Lq95lJWVYd68eXBwcICuri7EYjEmTZqER48e8ZrHqz777DMIBAKsX79eJXmkpqZi5MiRMDAwgK6uLtzc3PDgwQNe8ygoKEBgYCAsLS0hEolgb2+PiIgIpeYA1O+4VVJSgoCAABgbG0NPTw8ffPABsrKylJ4LaXqortVEdY3qWn3yeBXVtdZT16gjpiLx8fEICAjA+fPncerUKZSVlWHo0KEoLCxUST4XLlzADz/8gF69eqlk+zk5OejXrx80NDRw/Phx3Lx5E2vWrEHbtm15zWPlypX4/vvvsXnzZqSmpmLlypUICwvDpk2bGnW7hYWFcHR0RHh4uNzlYWFh2LhxIyIiIpCYmAhdXV14eXmhpKSEtzyKiopw+fJlLFy4EJcvX8b//d//4fbt2xg5cqRSc6grj+qioqJw/vx5iMVipedQnzzu3buH/v37w87ODnFxcbh27RoWLlwIbW1tXvOYPXs2Tpw4gT179iA1NRVBQUEIDAzE4cOHlZpHfY5bX3zxBY4cOYKDBw8iPj4ejx49wpgxY5SaB2maqK7JorpGda2+eVRHda1Sq6lrjDQJ2dnZDACLj4/nfdv5+fmsS5cu7NSpU2zgwIFs1qxZvOcwb9481r9/f963+6oRI0awqVOnyrSNGTOGTZw4kbccALCoqCjuuVQqZebm5mzVqlVc24sXL5iWlhb7z3/+w1se8iQlJTEALD09nfc8Hj58yNq3b89u3LjBrK2t2bp16xotB0V5+Pj4sI8//rhRt1ufPHr06MEWL14s09a7d2/29ddfN2ourx63Xrx4wTQ0NNjBgwe5mNTUVAaAJSQkNGoupOmhukZ1rQrVtfrlQXXtpdZS1+iMWBORm5sLADAyMuJ92wEBARgxYgQ8PT1533aVw4cPw9XVFWPHjoWpqSmcnZ2xbds23vPw8PDA6dOncefOHQDA1atX8ccff2D48OG851IlLS0NmZmZMt+PgYEB3N3dkZCQoLK8gMr9ViAQwNDQkNftSqVS+Pr6Ijg4GD169OB129VzOHbsGLp27QovLy+YmprC3d291uEmjcXDwwOHDx/G33//DcYYYmNjcefOHQwdOrRRt/vqcevSpUsoKyuT2Vft7OxgZWWl8n2V8I/qGtU1Raiu1UR1TVZrqWvUEWsCpFIpgoKC0K9fP/Ts2ZPXbe/fvx+XL1/G8uXLed3uq/766y98//336NKlC06ePIkZM2bg888/x+7du3nNY/78+Rg/fjzs7OygoaEBZ2dnBAUFYeLEibzmUV1mZiYAwMzMTKbdzMyMW6YKJSUlmDdvHj766CO0adOG122vXLkS6urq+Pzzz3ndbnXZ2dkoKCjAihUrMGzYMMTExOD999/HmDFjEB8fz2sumzZtgr29PSwtLaGpqYlhw4YhPDwcAwYMaLRtyjtuZWZmQlNTs8YfMKreVwn/qK5RXasN1bWaqK7Jai11Tf2N10DeWEBAAG7cuIE//viD1+1mZGRg1qxZOHXqlNLH/jaUVCqFq6srli1bBgBwdnbGjRs3EBERgcmTJ/OWx4EDB7B3717s27cPPXr0QHJyMoKCgiAWi3nNo6krKyvDuHHjwBjD999/z+u2L126hA0bNuDy5csQCAS8brs6qVQKABg1ahS++OILAICTkxP+97//ISIiAgMHDuQtl02bNuH8+fM4fPgwrK2tcebMGQQEBEAsFjfaGQFVHbdI80B1jepac0N1jeqaKo5bdEZMxQIDA3H06FHExsbC0tKS121funQJ2dnZ6N27N9TV1aGuro74+Hhs3LgR6urqqKio4C0XCwsL2Nvby7R1795d6bP01CU4OJj79dDBwQG+vr744osvVPrLqrm5OQDUmKEnKyuLW8anqmKVnp6OU6dO8f6r4dmzZ5GdnQ0rKytuv01PT8ecOXNgY2PDWx4mJiZQV1dX+X5bXFyMBQsWYO3atfD29kavXr0QGBgIHx8frF69ulG2qei4ZW5uDolEghcvXsjEq2pfJapBda0S1TXFqK7JoromqzXVNeqIqQhjDIGBgYiKisJ///tf2Nra8p7DkCFDcP36dSQnJ3MPV1dXTJw4EcnJyVBTU+Mtl379+tWYLvTOnTuwtrbmLQegcgYloVD2n4Wamhr3K5Eq2NrawtzcHKdPn+ba8vLykJiYiL59+/KaS1Wxunv3Ln7//XcYGxvzun0A8PX1xbVr12T2W7FYjODgYJw8eZK3PDQ1NeHm5qby/basrAxlZWW87Ld1HbdcXFygoaEhs6/evn0bDx484H1fJfyjuiaL6ppiVNdkUV2T1ZrqGg1NVJGAgADs27cPhw4dgr6+PjfO1MDAACKRiJcc9PX1a4zd19XVhbGxMe9j+r/44gt4eHhg2bJlGDduHJKSkrB161Zs3bqV1zy8vb2xdOlSWFlZoUePHrhy5QrWrl2LqVOnNup2CwoK8Oeff3LP09LSkJycDCMjI1hZWSEoKAhLlixBly5dYGtri4ULF0IsFmP06NG85WFhYYEPP/wQly9fxtGjR1FRUcHtt0ZGRtDU1OQlDysrqxqFUkNDA+bm5ujWrZvScqhPHsHBwfDx8cGAAQMwePBgnDhxAkeOHEFcXByveQwcOBDBwcEQiUSwtrZGfHw8fvzxR6xdu1apedR13DIwMMC0adMwe/ZsGBkZoU2bNpg5cyb69u2Lt956S6m5kKaH6posqmtU1+qbB9W1VlzX3njeRfJaAMh9REZGqjQvVU3zyxhjR44cYT179mRaWlrMzs6Obd26lfcc8vLy2KxZs5iVlRXT1tZmHTt2ZF9//TUrLS1t1O3GxsbK3R8mT57MGKuc6nfhwoXMzMyMaWlpsSFDhrDbt2/zmkdaWprC/TY2Npa3PORprGl+65PHjh07WOfOnZm2tjZzdHRk0dHRvOfx+PFj5ufnx8RiMdPW1mbdunVja9asYVKpVKl51Oe4VVxczPz9/Vnbtm2Zjo4Oe//999njx4+Vmgdpmqiu1UR1jepaffKQh+pa66hrgn+SIIQQQgghhBDCE7pGjBBCCCGEEEJ4Rh0xQgghhBBCCOEZdcQIIYQQQgghhGfUESOEEEIIIYQQnlFHjBBCCCGEEEJ4Rh0xQgghhBBCCOEZdcQIIYQQQgghhGfUESPNyq5du2BoaKjqNJodGxsbrF+/XiXbHjRoEIKCghr0mtDQUDg5OXHP/fz8MHr0aKXm1RhU+TkTQponqmuvh+oaP6iuNS7qiJFmxcfHB3fu3FF1GvUWFxcHgUCAtm3boqSkRGbZhQsXIBAIIBAIasRXPczMzPDBBx/gr7/+4mKuXr2KkSNHwtTUFNra2rCxsYGPjw+ys7N5e19827BhA3bt2qXqNOp04cIF/Pvf/1Z1GoSQZoTqGtW1pozqWuOijhhpVkQiEUxNTVWdRoPp6+sjKipKpm3Hjh2wsrKSG3/79m08evQIBw8eREpKCry9vVFRUYEnT55gyJAhMDIywsmTJ5GamorIyEiIxWIUFhby8VZUwsDAoFn8YtyuXTvo6OioOg1CSDNCdY3qWlNGda1xUUeMvJZBgwZh5syZCAoKQtu2bWFmZoZt27ahsLAQU6ZMgb6+Pjp37ozjx49zr6moqMC0adNga2sLkUiEbt26YcOGDdzykpIS9OjRQ+aXl3v37kFfXx87d+4EUHMIR9Wp/p07d8LKygp6enrw9/dHRUUFwsLCYG5uDlNTUyxdupR7zf379yEQCJCcnMy1vXjxAgKBAHFxcQBe/oJ38uRJODs7QyQS4e2330Z2djaOHz+O7t27o02bNpgwYQKKiorq/LwmT57MvQcAKC4uxv79+zF58mS58aamprCwsMCAAQMQEhKCmzdv4s8//8S5c+eQm5uL7du3w9nZGba2thg8eDDWrVsHW1vbWnPIz8/HRx99BF1dXbRv3x7h4eEyyx88eIBRo0ZBT08Pbdq0wbhx45CVlVXjs/7pp59gY2MDAwMDjB8/Hvn5+VxMYWEhJk2aBD09PVhYWGDNmjV1fjYAsGLFCpiZmUFfXx/Tpk2r8Svrq0M4Xmf/A4AbN25g+PDh0NPTg5mZGXx9ffH06VOZ9X7++ef48ssvYWRkBHNzc4SGhnLLGWMIDQ2FlZUVtLS0IBaL8fnnn3PLXx3CoYzPlBDCD6prVNeorlFd4xt1xMhr2717N0xMTJCUlISZM2dixowZGDt2LDw8PHD58mUMHToUvr6+3AFdKpXC0tISBw8exM2bNxESEoIFCxbgwIEDAABtbW3s3bsXu3fvxqFDh1BRUYGPP/4Y77zzDqZOnaowj3v37uH48eM4ceIE/vOf/2DHjh0YMWIEHj58iPj4eKxcuRLffPMNEhMTG/weQ0NDsXnzZvzvf/9DRkYGxo0bh/Xr12Pfvn04duwYYmJisGnTpjrX4+vri7Nnz+LBgwcAgF9//RU2Njbo3bt3na8ViUQAAIlEAnNzc5SXlyMqKgqMsQa9l1WrVsHR0RFXrlzB/PnzMWvWLJw6dQpA5XczatQoPH/+HPHx8Th16hT++usv+Pj4yKzj3r17iI6OxtGjR3H06FHEx8djxYoV3PLg4GDEx8fj0KFDiImJQVxcHC5fvlxrXgcOHEBoaCiWLVuGixcvwsLCAlu2bKnz/TR0/3vx4gXefvttODs74+LFizhx4gSysrIwbty4GuvV1dVFYmIiwsLCsHjxYu5z+vXXX7Fu3Tr88MMPuHv3LqKjo+Hg4CA3P2V9poQQ/lBdo7pGdY3qGq8YIa9h4MCBrH///tzz8vJypqury3x9fbm2x48fMwAsISFB4XoCAgLYBx98INMWFhbGTExMWGBgILOwsGBPnz7llkVGRjIDAwPu+aJFi5iOjg7Ly8vj2ry8vJiNjQ2rqKjg2rp168aWL1/OGGMsLS2NAWBXrlzhlufk5DAALDY2ljHGWGxsLAPAfv/9dy5m+fLlDAC7d+8e1/bpp58yLy8vhe+vaj05OTls9OjR7Ntvv2WMMTZ48GC2YcMGFhUVxar/M6wezxhjjx49Yh4eHqx9+/astLSUMcbYggULmLq6OjMyMmLDhg1jYWFhLDMzU2EOjDFmbW3Nhg0bJtPm4+PDhg8fzhhjLCYmhqmpqbEHDx5wy1NSUhgAlpSUxBiT/1kHBwczd3d3xhhj+fn5TFNTkx04cIBb/uzZMyYSidisWbMU5ta3b1/m7+8v0+bu7s4cHR2555MnT2ajRo3inr/O/vfdd9+xoUOHymwnIyODAWC3b9+Wu17GGHNzc2Pz5s1jjDG2Zs0a1rVrVyaRSOS+F2tra7Zu3TrGmHI+U0IIf6iuVaK6RnWtOqprjYvOiJHX1qtXL+7/1dTUYGxsLPMripmZGQDIXGwbHh4OFxcXtGvXDnp6eti6dSv3a1qVOXPmoGvXrti8eTN27twJY2PjWvOwsbGBvr6+zHbt7e0hFApl2l7not/q79HMzAw6Ojro2LHja6136tSp2LVrF/766y8kJCRg4sSJCmMtLS2hq6vLjZH/9ddfoampCQBYunQpMjMzERERgR49eiAiIgJ2dna4fv16rdvv27dvjeepqakAgNTUVHTo0AEdOnTgltvb28PQ0JCLAWp+1hYWFtz7v3fvHiQSCdzd3bnlRkZG6NatW615paamyrxGXq7yNHT/u3r1KmJjY6Gnp8c97OzsuNzlrffV9zh27FgUFxejY8eOmD59OqKiolBeXq7wfb3pZ0oI4RfVNaprVNeorvGJOmLktWloaMg8FwgEMm1VsyZJpVIAwP79+zF37lxMmzYNMTExSE5OxpQpUyCRSGTWk52djTt37kBNTQ1379594zyq2qryqCpkrNoQiLKysjrXXdd66zJ8+HAUFxdj2rRp8Pb2rrUQnz17FteuXUNeXh6Sk5NrHNCNjY0xduxYrF69GqmpqRCLxVi9enW98ngTb/L++ciltv2voKAA3t7eSE5OlnncvXsXAwYMqHW9Vevo0KEDbt++jS1btkAkEsHf3x8DBgxQuP+87vtQ1WdKSGtHdY3qGtU1qmt8oo4Y4c25c+fg4eEBf39/ODs7o3PnzjK/2FSZOnUqHBwcsHv3bsybN0/mVxZlaNeuHQDg8ePHXFv1C5wbi7q6OiZNmoS4uLharw0AAFtbW3Tq1EnmFyVFNDU10alTpzpnlzp//nyN5927dwcAdO/eHRkZGcjIyOCW37x5Ey9evIC9vX2dOQBAp06doKGhIXPNQk5OTp3TMnfv3r3GdQ6v5qoMvXv3RkpKCmxsbNC5c2eZh66ubr3XIxKJ4O3tjY0bNyIuLg4JCQlyf7VVxmdKCGnaqK5RXZOH6hqpL3VVJ0Bajy5duuDHH3/EyZMnYWtri59++gkXLlyQmRUpPDwcCQkJuHbtGjp06IBjx45h4sSJOH/+PDeE4U2JRCK89dZbWLFiBWxtbZGdnY1vvvlGKeuuy3fffYfg4OA6h6UocvToUezfvx/jx49H165dwRjDkSNH8NtvvyEyMrLW1547dw5hYWEYPXo0Tp06hYMHD+LYsWMAAE9PTzg4OGDixIlYv349ysvL4e/vj4EDB8LV1bVeuenp6WHatGnc+zM1NcXXX38tM5RGnlmzZsHPzw+urq7o168f9u7di5SUFJmhMsoQEBCAbdu24aOPPuJmj/rzzz+xf/9+bN++HWpqanWuY9euXaioqIC7uzt0dHSwZ88eiEQiWFtb14hVxmdKCGnaqK5RXZOH6hqpLzojRnjz6aefYsyYMfDx8YG7uzuePXsGf39/bvmtW7cQHByMLVu2cOOPt2zZgqdPn2LhwoVKzWXnzp0oLy+Hi4sLgoKCsGTJEqWuXxFNTU2YmJjI3OyyIezt7aGjo4M5c+bAyckJb731Fg4cOIDt27fD19e31tfOmTMHFy9ehLOzM5YsWYK1a9fCy8sLQOWwgUOHDqFt27YYMGAAPD090bFjR/z8888Nym/VqlX417/+BW9vb3h6eqJ///5wcXGp9TU+Pj5YuHAhvvzyS7i4uCA9PR0zZsxo0HbrQywW49y5c6ioqMDQoUPh4OCAoKAgGBoa1llUqxgaGmLbtm3o168fevXqhd9//x1HjhyR+weIsj5TQkjTRXWN6po8VNdIfQkYa+BcoYQQQgghhBBC3gidESOEEEIIIYQQnlFHjBBCCCGEEEJ4Rh0xQgghhBBCCOEZdcQIIYQQQgghhGfUESOEEEIIIYQQnlFHjBBCCCGEEEJ4Rh0xQgghhBBCCOEZdcQIIYQQQgghhGfUESOEEEIIIYQQnlFHjBBCCCGEEEJ4Rh0xQgghhBBCCOEZdcQIIYQQQgghhGf/D2graAxzgfhfAAAAAElFTkSuQmCC", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "import matplotlib.gridspec as gridspec\n", "\n", diff --git a/python/ffsim/tenpy/__init__.py b/python/ffsim/tenpy/__init__.py index 622bfdbaf..bb0f692cf 100644 --- a/python/ffsim/tenpy/__init__.py +++ b/python/ffsim/tenpy/__init__.py @@ -11,8 +11,10 @@ """Code that uses TeNPy, e.g. for emulating quantum circuits.""" from ffsim.tenpy.circuits.gates import ( - gate1, - gate2, + apply_diag_coulomb_evolution, + apply_gate1, + apply_gate2, + apply_orbital_rotation, givens_rotation, num_interaction, num_num_interaction, @@ -25,8 +27,10 @@ from ffsim.tenpy.util import product_state_as_mps __all__ = [ - "gate1", - "gate2", + "apply_diag_coulomb_evolution", + "apply_gate1", + "apply_gate2", + "apply_orbital_rotation", "givens_rotation", "lucj_circuit_as_mps", "MolecularChain", diff --git a/python/ffsim/tenpy/circuits/gates.py b/python/ffsim/tenpy/circuits/gates.py index 87cfe97a9..9175dbd24 100644 --- a/python/ffsim/tenpy/circuits/gates.py +++ b/python/ffsim/tenpy/circuits/gates.py @@ -6,6 +6,7 @@ from tenpy.networks.mps import MPS from tenpy.networks.site import SpinHalfFermionSite +from ffsim.linalg import givens_decomposition from ffsim.spin import Spin # ignore lowercase argument and variable checks to maintain TeNPy naming conventions @@ -247,7 +248,7 @@ def num_num_interaction(theta: float, spin: Spin) -> np.ndarray: return NNgate_sym -def gate1(U1: np.ndarray, site: int, psi: MPS) -> None: +def apply_gate1(U1: np.ndarray, site: int, psi: MPS) -> None: r"""Apply a single-site gate to a `TeNPy MPS `__ wavefunction. @@ -257,7 +258,8 @@ def gate1(U1: np.ndarray, site: int, psi: MPS) -> None: site: The gate will be applied to `site` on the `TeNPy MPS `__ wavefunction. - psi: The wavefunction MPS. + psi: The `TeNPy MPS `__ + wavefunction. Returns: None @@ -268,13 +270,13 @@ def gate1(U1: np.ndarray, site: int, psi: MPS) -> None: psi.apply_local_op(site, U1_npc) -def gate2( +def apply_gate2( U2: np.ndarray, site: int, psi: MPS, eng: TEBDEngine, chi_list: list, - norm_tol: float, + norm_tol: float = 1e-5, ) -> None: r"""Apply a two-site gate to a `TeNPy MPS `__ wavefunction. @@ -308,3 +310,110 @@ def gate2( # recanonicalize psi if below error threshold if np.linalg.norm(psi.norm_test()) > norm_tol: psi.canonical_form_finite() + + +def apply_orbital_rotation( + mat: np.ndarray, + psi: MPS, + eng: TEBDEngine, + chi_list: list, + norm_tol: float = 1e-5, +) -> None: + r"""Apply an orbital rotation gate to a + `TeNPy MPS `__ + wavefunction. + + The orbital rotation gate is defined in + `apply_orbital_rotation `__. + + Args: + mat: The orbital rotation matrix of dimension `(norb, norb)`. + psi: The `TeNPy MPS `__ + wavefunction. + eng: The + `TeNPy TEBDEngine `__. + chi_list: The list to which to append the MPS bond dimensions as the circuit is + evaluated. + norm_tol: The norm error above which we recanonicalize the wavefunction, as + defined in the + `TeNPy documentation `__. + + Returns: + None + """ + + # Givens decomposition + givens_list, diag_mat = givens_decomposition(mat) + + # apply the Givens rotation gates + for gate in givens_list: + theta = np.arccos(gate.c) + s = np.conj(gate.s) + phi = np.real(1j * np.log(-s / np.sin(theta))) if theta else 0 + conj = True if gate.j < gate.i else False + apply_gate2( + givens_rotation(theta, Spin.ALPHA_AND_BETA, conj, phi=phi), + max(gate.i, gate.j), + psi, + eng, + chi_list, + norm_tol, + ) + + # apply the number interaction gates + for i, z in enumerate(diag_mat): + theta = float(np.angle(z)) + apply_gate1( + np.exp(1j * theta) * num_interaction(-theta, Spin.ALPHA_AND_BETA), i, psi + ) + + +def apply_diag_coulomb_evolution( + mat: np.ndarray, psi: MPS, eng: TEBDEngine, chi_list: list, norm_tol: float = 1e-5 +) -> None: + r"""Apply a diagonal Coulomb evolution gate to a + `TeNPy MPS `__ + wavefunction. + + The diagonal Coulomb evolution gate is defined in + `apply_diag_coulomb_evolution `__. + + Args: + mat: The diagonal Coulomb matrices of dimension `(2, norb, norb)`. + psi: The `TeNPy MPS `__ + wavefunction. + eng: The + `TeNPy TEBDEngine `__. + chi_list: The list to which to append the MPS bond dimensions as the circuit is + evaluated. + norm_tol: The norm error above which we recanonicalize the wavefunction, as + defined in the + `TeNPy documentation `__. + + Returns: + None + """ + + # extract norb + assert np.shape(mat)[1] == np.shape(mat)[2] + norb = np.shape(mat)[1] + + # unpack alpha-alpha and alpha-beta matrices + mat_aa, mat_ab = mat + + # apply alpha-alpha gates + for i in range(norb): + for j in range(norb): + if j > i and mat_aa[i, j]: + apply_gate2( + num_num_interaction(-mat_aa[i, j], Spin.ALPHA_AND_BETA), + j, + psi, + eng, + chi_list, + norm_tol, + ) + + # apply alpha-beta gates + for i in range(norb): + apply_gate1(on_site_interaction(-mat_ab[i, i]), i, psi) diff --git a/python/ffsim/tenpy/circuits/lucj_circuit.py b/python/ffsim/tenpy/circuits/lucj_circuit.py index c3b57013c..1a826a9ca 100644 --- a/python/ffsim/tenpy/circuits/lucj_circuit.py +++ b/python/ffsim/tenpy/circuits/lucj_circuit.py @@ -1,19 +1,12 @@ from __future__ import annotations import numpy as np -from qiskit.circuit import QuantumCircuit, QuantumRegister from tenpy.algorithms.tebd import TEBDEngine from tenpy.networks.mps import MPS -import ffsim -from ffsim.spin import Spin from ffsim.tenpy.circuits.gates import ( - gate1, - gate2, - givens_rotation, - num_interaction, - num_num_interaction, - on_site_interaction, + apply_diag_coulomb_evolution, + apply_orbital_rotation, ) from ffsim.tenpy.util import product_state_as_mps from ffsim.variational.ucj_spin_balanced import UCJOpSpinBalanced @@ -21,7 +14,7 @@ def lucj_circuit_as_mps( norb: int, - nelec: tuple, + nelec: int | tuple[int, int], ucj_op: UCJOpSpinBalanced, options: dict, norm_tol: float = 1e-5, @@ -52,82 +45,20 @@ def lucj_circuit_as_mps( # prepare initial Hartree-Fock state psi = product_state_as_mps(norb, nelec, 0) - # construct the qiskit circuit - qubits = QuantumRegister(2 * norb) - circuit = QuantumCircuit(qubits) - circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits) - # define the TEBD engine eng = TEBDEngine(psi, None, options) - # execute the tenpy circuit - for ins in circuit.decompose(reps=2): - if ins.operation.name == "p": - qubit = ins.qubits[0] - idx = qubit._index - spin_flag = Spin.ALPHA if idx < norb else Spin.BETA - lmbda = ins.operation.params[0] - gate1( - np.exp(1j * lmbda) * num_interaction(-lmbda, spin_flag), idx % norb, psi - ) - elif ins.operation.name == "xx_plus_yy": - qubit0 = ins.qubits[0] - qubit1 = ins.qubits[1] - idx0, idx1 = qubit0._index, qubit1._index - if idx0 < norb and idx1 < norb: - spin_flag = Spin.ALPHA - elif idx0 >= norb and idx1 >= norb: - spin_flag = Spin.BETA - else: - raise ValueError("XXPlusYY gate not allowed across spin sectors") - theta_val = ins.operation.params[0] - beta_val = ins.operation.params[1] - # directionality important when beta!=0 - conj_flag = True if idx0 > idx1 else False - gate2( - givens_rotation( - theta_val / 2, spin_flag, conj_flag, phi=beta_val - np.pi / 2 - ), - max(idx0 % norb, idx1 % norb), - psi, - eng, - chi_list, - norm_tol, - ) - elif ins.operation.name == "cp": - qubit0 = ins.qubits[0] - qubit1 = ins.qubits[1] - idx0, idx1 = qubit0._index, qubit1._index - lmbda = ins.operation.params[0] - # onsite (different spins) - if np.abs(idx0 - idx1) == norb: - gate1(on_site_interaction(-lmbda), min(idx0, idx1), psi) - # NN (up spins) - elif np.abs(idx0 - idx1) == 1 and idx0 < norb and idx1 < norb: - gate2( - num_num_interaction(-lmbda, Spin.ALPHA), - max(idx0, idx1), - psi, - eng, - chi_list, - norm_tol, - ) - # NN (down spins) - elif np.abs(idx0 - idx1) == 1 and idx0 >= norb and idx1 >= norb: - gate2( - num_num_interaction(-lmbda, Spin.BETA), - max(idx0 % norb, idx1 % norb), - psi, - eng, - chi_list, - norm_tol, - ) - else: - raise ValueError( - "CPhase only implemented onsite (different spins) " - "and NN (same spins)" - ) - else: - raise ValueError(f"gate {ins.operation.name} not implemented.") + # construct the LUCJ MPS + n_reps = np.shape(ucj_op.orbital_rotations)[0] + for i in range(n_reps): + apply_orbital_rotation( + np.conj(ucj_op.orbital_rotations[i]).T, psi, eng, chi_list, norm_tol + ) + apply_diag_coulomb_evolution( + ucj_op.diag_coulomb_mats[i], psi, eng, chi_list, norm_tol + ) + apply_orbital_rotation( + ucj_op.orbital_rotations[i], psi, eng, chi_list, norm_tol + ) return psi, chi_list diff --git a/tests/python/tenpy/lucj_circuit_test.py b/tests/python/tenpy/lucj_circuit_test.py index 6b7fc2888..c683566d5 100644 --- a/tests/python/tenpy/lucj_circuit_test.py +++ b/tests/python/tenpy/lucj_circuit_test.py @@ -75,7 +75,7 @@ def test_lucj_circuit_as_mps(norb: int, nelec: tuple[int, int], connectivity: st interaction_pairs=_interaction_pairs_spin_balanced_( connectivity=connectivity, norb=norb ), - with_final_orbital_rotation=True, + with_final_orbital_rotation=False, ) params = rng.uniform(-10, 10, size=n_params) lucj_op = ffsim.UCJOpSpinBalanced.from_parameters( @@ -85,7 +85,7 @@ def test_lucj_circuit_as_mps(norb: int, nelec: tuple[int, int], connectivity: st interaction_pairs=_interaction_pairs_spin_balanced_( connectivity=connectivity, norb=norb ), - with_final_orbital_rotation=True, + with_final_orbital_rotation=False, ) # generate the corresponding LUCJ circuit