{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Chains and PT-TEBD\n", "\n", "An introduction on how to use the OQuPy package to compute the dynamics of a chain of open quantum systems using the process tensor approach to time evolving block decimation (PT-TEBD). We illustrate this by applying PT-TEBD to a 5-site XYZ Heisenberg spin chain. This method and example is explained in detail in [Fux2022] ([arXiv:2201.05529](https://arxiv.org/abs/2201.05529)).\n", "\n", "- [launch binder](https://mybinder.org/v2/gh/tempoCollaboration/OQuPy/HEAD?labpath=tutorials%2Fpt_tebd.ipynb) (runs in browser),\n", "- [download the jupyter file](https://raw.githubusercontent.com/tempoCollaboration/OQuPy/main/tutorials/pt_tebd.ipynb), or\n", "- read through the text below and code along.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Contents:**\n", "\n", "* Example - Heisenberg spin chain\n", " * 1. Closed Heisenberg spin chain\n", " * 2. Open Heisenberg spin chain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's import OQuPy and some other packages we are going to use" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import sys\n", "sys.path.insert(0,'..')\n", "\n", "import oqupy\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and check what version of tempo we are using." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'0.4.0'" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "oqupy.__version__" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also import some shorthands for the spin-1/2 operators and and density matrices." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "sx = 0.5 * oqupy.operators.sigma(\"x\")\n", "sy = 0.5 * oqupy.operators.sigma(\"y\")\n", "sz = 0.5 * oqupy.operators.sigma(\"z\")\n", "up_dm = oqupy.operators.spin_dm(\"z+\")\n", "down_dm = oqupy.operators.spin_dm(\"z-\")\n", "mixed_dm = oqupy.operators.spin_dm(\"mixed\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "------------------------------------\n", "## Example - Heisenberg spin chain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Closed Heisenberg spin chain\n", "\n", "Let's calculate the dynamics of a short XYZ Heisenberg spin chain with the same parameters as in [Fux2022]. Before we include any environment coupling we first consider the closed chain with the Hamiltonian\n", "\n", "$$ H_\\mathrm{chain} = \\sum_{n=1}^N \\epsilon s_n^z\n", " + \\sum_{n=1}^{N-1} \\sum_{\\gamma \\in \\{x,y,z\\}}\n", " J^\\gamma s_{n}^\\gamma s_{n+1}^\\gamma \\,\\mathrm{,}$$\n", "\n", "with $N=5$, $\\epsilon=1.0$, $J^x = 1.3$, $J^y = 0.7$, and $J^z = 1.2$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "N = 5\n", "epsilon = 1.0\n", "J = [1.3, 0.7, 1.2]\n", "\n", "system_chain = oqupy.SystemChain(hilbert_space_dimensions=[2]*N)\n", "\n", "for n in range(N):\n", " system_chain.add_site_hamiltonian(\n", " site=n,\n", " hamiltonian=epsilon*sz)\n", "\n", "for n in range(N-1):\n", " for J_gamma, s_gamma in zip(J, [sx, sy, sz]):\n", " system_chain.add_nn_hamiltonian(\n", " site=n,\n", " hamiltonian_l=J_gamma*s_gamma,\n", " hamiltonian_r=s_gamma)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We begin with an initial chain state where the first spin is 'up' (in the $z$-basis) and all the others are 'down'.\n", "$$\\tilde{\\rho}^A_\\mathrm{chain} = | \\uparrow \\rangle \\langle \\uparrow |\n", " \\otimes | \\downarrow \\rangle \\langle \\downarrow |\n", " \\otimes | \\downarrow \\rangle \\langle \\downarrow |\n", " \\otimes | \\downarrow \\rangle \\langle \\downarrow |\n", " \\otimes | \\downarrow \\rangle \\langle \\downarrow | $$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "initial_augmented_mps_A = oqupy.AugmentedMPS([up_dm, down_dm, down_dm, down_dm, down_dm])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To run a PT-TEBD simulation we need to choose three computation parameters:\n", "\n", "1. The time step length `dt`,\n", "2. The Trotterization order `order` (currently only `1` and `2` are implemented), and\n", "3. The relative singular value truncation tolerance `epsrel`.\n", "\n", "We describe details of the computation parameters in the supplemental material of [Fux2022]." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "pt_tebd_params = oqupy.PtTebdParameters(\n", " dt=0.2,\n", " order=2,\n", " epsrel=1.0e-6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we insert all this information into a new `PtTebd` object to prepare the simulation. We also specify a list of process tensors that represent the environment of each site of the chain. Because in this first example we don't have any environment, we specify all to be `None`. The `dynamics_sites` parameter allows us to specify a list of sites whose reduced density matrix shall be recorded during the propagation." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "pt_tebd_closed = oqupy.PtTebd(\n", " initial_augmented_mps=initial_augmented_mps_A,\n", " system_chain=system_chain,\n", " process_tensors=[None, None, None, None, None],\n", " parameters=pt_tebd_params,\n", " dynamics_sites=[0, 1, 2, 3, 4])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then go ahead and set off the actual computation. Let's (for the sake of computation time) be modest and say we are only interested in 20 time steps." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PT-TEBD computation (closed spin chain):\n", "--> PT-TEBD computation:\n", "100.0% 20 of 20 [########################################] 00:00:00\n", "Elapsed time: 0.6s\n" ] } ], "source": [ "num_steps = 20\n", "print(\"PT-TEBD computation (closed spin chain):\")\n", "results_closed = pt_tebd_closed.compute(num_steps, progress_type=\"bar\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The computation returns a results dictionary which in addition to the the dynamics of the before specified sites carries all sorts of other information about the propagation (such as the total spin chain norm, the bond dimensions of the augmented MPS, etc.)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Keys of the results dictionary: dict_keys(['time', 'norm', 'bond_dimensions', 'dynamics', 'pt_bond_dimensions'])\n", "Keys of the dynamics results: dict_keys([0, 1, 2, 3, 4])\n" ] } ], "source": [ "print(f\"Keys of the results dictionary: {results_closed.keys()}\")\n", "print(f\"Keys of the dynamics results: {results_closed['dynamics'].keys()}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use the dynamics results to then compute the evolution of the $s^z_n$ observables." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABJCElEQVR4nO3deXhU1f3H8feZLZOZ7AvZQ9j3RQgEBBVwA1FQCgrigqK4tbXVttraX6ttrbWtdlUR3CiioKgIiIpsrux7wr4nIfsymWyznt8fN4QAYU8yM8l5Pc88M5O5mfvNwHzmzLnnniOklCiKoiitn87XBSiKoigtQwW+oihKG6ECX1EUpY1Qga8oitJGqMBXFEVpIwy+LuBsYmJiZFpamq/LUBRFCSibN28ullLGNvaY3wZ+WloamzZt8nUZiqIoAUUIcfRsj6kuHUVRlDZCBb6iKEoboQJfURSljfDbPnxFUZSm4HK5yMnJoba21telNCmz2UxycjJGo/GCf0cFvqIorVpOTg6hoaGkpaUhhPB1OU1CSklJSQk5OTl06NDhgn9PdekoitKq1dbWEh0d3WrCHkAIQXR09EV/a1GBryhKq9eawv6ES/mbWl2Xjtvj5W/L95IcaSElMpiUKAtJEcGYjXpfl6YoiuJTrS7wiyudvP3dEZwe7yk/jwsLIiXSQmqUheSokx8GKVEW4sPM6HWtrwWgKIrSUKsL/PhwM3v+OJpCu4PssmqyS6s5VlpNdmkN2WXVrDtUQt62XBqu+2LUCxIjgkmJtJASFcygtCjG9E4g2KS+FSiK0nq0usAH0OkE8eFm4sPNDEqLOuNxp9vL8fKaug+EmvoPhuyyGr7IzOf9Ddn87tMsbu6bwKT0ZAakRrbKPkBFUVrWqFGjWL58OQbDmdF77NgxZs2axeHDhwH473//y/33388nn3zSZPtvlYF/PiaDjrQYK2kx1jMe83olG4+U8uHmHBZvP878jdl0jLUycWAyPxqQTFyY2QcVK4oS6LKysoiOjm407AFSU1OZPn06v//973nttdewWq2UlpZSUlJCdHR0k9SgRumcRqcTZHSM5u+T+rHhmev468S+RFtN/PWLvQx9YSXT3t7AZzvycLg9vi5VUZQA8umnn3LrrbcCMGfOHAYOHEjfvn0ZPnw4AEeOHOHZZ5+tD3uAsWPHsmTJkiarQfjrIubp6enSn2bLPFxcxcLN2Xy8JZc8Wy0RFiPj+yUyKT2FXolhqstHUfzU7t276dGjBwDPLcli1/GKJn3+nolh/P6WXufdbvjw4SxduhS9Xk9GRgbbtm3DZDJRXl5OREQEPXv25MYbb8RqtfLEE08QFRXF3r17eeqpp1i0aFGjz9nwbztBCLFZSpne2PZtskvnUnSIsfLLG7vzxPXd+P5AMR9syub9jdnMWXuU7vGhTEpP4db+iUSHBPm6VEVR/Ex1dTVOp5OIiAiqq6upqanhySef5N577yU9XcvmXbt2nfF73bp1Y+/evU1Whwr8i6TXCa7uGsvVXWOxVbtYvD2XDzfn8Melu3hh2W7G9k3g6THdSQgP9nWpiqKc5kJa4s3BYrEghKCyspKQkBAyMzNZsmQJM2bM4IEHHuDRRx9t9PeOHj16UVMnnI8K/MsQbjFy99A07h6axt58Ox9syubddUf5alcBj1/bhfuGdcBkUIdJFEWBG2+8kS+++IJ+/frRpUsXJk+ezK5du+qnR6iqquLRRx/FZDIxYsQIpk6dyqeffsr48eObrAaVRk2kW3wo/3dzT1Y8cQ1Xdormhc/3cNO/v2XtwRJfl6Yoih8YP348ixYt4vnnn6dbt24MGDCAw4cP17fuP/74YyZOnMjs2bNZvHgxAEuWLGnSwFct/CaWEmXhjXsHsWJXAc8uyWLK7HWM75/IMzf1oJ0a0qkobdbAgQPZsWMHW7ZsaXRoZk5ODn369AFAr9dTVlaGw+EgPj6+yWpQLfxmcl3POFY8cQ0/HdWZz3fmM+qlr3nzu8O4T5vyQVGUtmPHjh1nHYefnJxMTk4OAF6vl8jISL755psm3b8K/GZkNup54oZufPnzqxnQPpI/Lt3Fzf/5jo1HSn1dmqIofmbChAl89NFHPPLII9xyyy3Nsg/VpdMCOsRYmXPfIL7MyucPS3YxaeZafjQgmafHdCc2VA3jVBQFrFYrb7/9drPuQ7XwW4gQgtG9E1jx5DU8MqITi7fnMuqlNfxv7RE8Xv88+U1RlNZFBX4Ls5gMPDW6O58/fjV9ksL53adZjPvvd2w5Vubr0hRFaeVU4PtI53YhzHsgg/9MuYLiSgcTXv2BX3+8g2qn29elKYrSSjVJ4AshRgsh9gohDgghnj7Hdj8SQkghRKPzPLQ1Qghu6ZfIyidH8OBVHZi/MZsJr/7A0ZIqX5emKEordNmBL4TQA68AY4CewBQhRM9GtgsFHgfWX+4+W5uQIAPPjO3JO/cNJs9Wyy3/+Y7Vewp9XZaiKK1MU7TwBwMHpJSHpJROYD7Q2KlhfwReBC5umfU25JqusSz9yXCSIy3cP2cj/165H686oKsoShNpisBPArIb3M+p+1k9IcQAIEVK+dm5nkgIMUMIsUkIsamoqKgJSgs8KVEWPnrkSm7rn8TLX+1jxtxNVNS6fF2WoiitQLMftBVC6ICXgSfPt62UcpaUMl1KmR4bG9vcpfmtYJOel27vx3PjerFmbxHj//s9+wrsvi5LUZQA1xSBnwukNLifXPezE0KB3sAaIcQRYAiwWB24PTchBPdemcb7M4ZQ6XBz6yvf89mOPF+XpSjKZRg1ahRud+Mj8Y4dO8Zvf/tbpk6dytSpUykrK+O2225r0v03ReBvBLoIIToIIUzAZGDxiQellDYpZYyUMk1KmQasA8ZJKf1nOSs/NigtiqU/GU6PhDAee28LLyzbrebjUZQAdKFr2ur1embNmkVkZGT9mrZN5bIDX0rpBn4MfAnsBj6QUmYJIf4ghBh3uc+vQFyYmfcfHMLdQ9rz+jeHuOetDZRUOnxdlqIoF0GtaXsO/ramrb/4cFM2zyzKJDYkiNfuGkDf5Ahfl6Qofq2xdV99Qa1pq1y0SekpdI8P4+F3NzNx5lr+dGtvbk9POf8vKooCnz8N+Tub9jnj+8CYv5xzE39Z01ZNrRCA+iSHs+QnwxmcFsWvFu7gmU924nB7fF2Woihn0XBNW4vFQmZmJsOGDWPGjBm8+uqrZ/09taatAkCU1cSc+wfzty/3MvPrg+zKq2DW3elqumVFOZfztMSb0/nWtD106BDPP/88NpuNhQsXAqg1bZWT9DrB02O689rUAezJszNp5g9kl1b7uixFURpxvjVtO3bsyJtvvnnK76g1bZUzjOmTQFy4mfve3sjEmT8wd3oGXeNCfV2WoigNnG9N29OpNW2VsxqQGskHDw1FSrj99bVsVfPrK4rfOdeatqdTa9oq59QtPpSPHrmS8GAjU99Yz7f72+Z8RIoSiEpKSnj44YfZunUrL7zwQrPsQ3XptDIpURY+fHgo97y5gfvf2ci/Jl/BTX0SfF2WoijnER0dzcyZM5t1H6qF3wq1CzWz4KGh9EuO4LH3tvD+hmO+LklRFD+gAr+VCg82Mnd6Btd0jeXXH+/ktTUHfV2Soig+pgK/FQs26Zl9Tzrj+yfy4hd7eGHZbvx1Kg1FUZqf6sNv5Yx6Hf+4vT/hwUZe/+YQZdVO/nxbHwx69VmvKG2NCvw2QKcTPDeuFxEWE/9euZ+KGjf/nNwfs1Hv69IURWlBqpnXRggheOL6rvzu5p58kZXP/e9spNLR+EIMiqK0Tirw25j7h3fg5dv7sf5wKVNnr6O0yunrkhRFaSEq8NugCQOSef2ugezOt3P762vJs9X4uiRFUVqACvw26rqecfzv/sEU2GqZ+NpaDhdX+bokRWn1WsOatkqAGtIxmvdnDKHW5WHyLBX6itKcWsWatkpg650UznsPDsHtkSr0FaUZ+cOatirwFbrFh/Leg0NweSRTZq3jiAp9RWlyy5YtY+zYsdjtdl588UXWrl3Ljh07WLp0KQA33XQTUVFRvPDCC5SWlgIn59BvKmocvgKcCP0M7py9nimz1zF/xhDaR1t9XZaiNKkXN7zIntI9Tfqc3aO689Tgp865jVrTVvE73ePDmPdARl2f/jqOlqiWvqI0BbWmreKXeiSEMe+BIUx9Yx1TZq1j/oyhpEZbfF2WojSJ87XEm9P51rRdtGgRn332GRUVFUyfPp0bbrhBrWmrNL+eiVroV9eN3jlWotbJVZTLdb41bW+99VZmz57NzJkzWbBgAaDWtFVaiBb6GUx942SffkqUaukryqW60DVt//SnP/HYY4+pNW2VltUrMZx3p2dQ6XAzedY6sktVS19RLse51rSVUvLUU08xZswYBgwYoNa0VVpe76Rw5j2gQl9Rmtt//vMfVqxYwcKFC5ttqUPhrwtipKeny02bNvm6DKVOZq6NO2evIyzYyPwZQ0iOVN07SmDYvXs3PXr08HUZzaKxv00IsVlKmd7Y9qqFr1wQraU/hIoaF5NnrSO3XE24piiBRgW+csH6JIfz7gMZ2GpcTJ61VoW+ogSYJgl8IcRoIcReIcQBIcTTjTz+hBBilxBihxBipRCifVPsV2l5fZMjeHd6BuXVLqbMWsdxFfqKEjAuO/CFEHrgFWAM0BOYIoToedpmW4F0KWVfYCHw18vdr+I7/VIimDs9g7IqJ5NV6CtKwGiKFv5g4ICU8pCU0gnMB045U0BKuVpKeWJ4xzoguQn2q/hQ/5QI5j6ghf6U2esoqKj1dUmKopxHUwR+EpDd4H5O3c/OZjrweRPsV/Gx/ikRzJk+mGK7gymz11Fkd/i6JEVRzqFFD9oKIe4C0oG/neXxGUKITUKITUVFRS1ZmnKJBqRG8s79g8krr2XqG+soqVShryj+qikCPxdIaXA/ue5npxBCXAc8A4yTUjaaClLKWVLKdCllemxsbBOUprSEQWlRvDktnaMl1dz15gbKq9XC6Irij5oi8DcCXYQQHYQQJmAysLjhBkKIK4DX0cK+sAn2qfiZKzvFMPuedA4WVXL3mxuw1bh8XZKiKKe57MCXUrqBHwNfAruBD6SUWUKIPwghxtVt9jcgBPhQCLFNCLH4LE+nBLCru8by+l0D2ZNfwb1vbcBeq0JfURpqFYuYSymXSSm7Sik7SSmfr/vZ76SUi+tuXyeljJNS9q+7jDv3MyqBamT3drxy5wAyc23c9/ZGqhyN/+dWlLZGLWKutEo39Irn31OuYGt2OdPnbKTG6fF1SYric/6wiLmaPE1pNp9uy+XnC7ZxZacY3rg3HbNR7+uSlDao4QRj+X/+M47dTbumbVCP7sT/5jfn3W748OEsXboUvV5PRkYG27Ztw2QyUV5eTkREBD179uTGG2/EarXyxBNPEBUVxd69e3nqqafOupD5xU6ephZAUZrN+P5JuDySXy7czsPvbub1uwcSZFChr7Q9/rKIuQp8pVlNHJiM2+Pl6Y938ti8rbw6dQAmg+pJVHzjQlrizaHhIuYhISFkZmayZMkSZsyYwQMPPFC/zOHpmnoRc/XOU5rd5MGp/HF8L1bsLuDx+Vtxe7y+LklRWtyJRcz379+P1Wpl8uTJ3HzzzfWLmO/evZuHH36YiRMn8tprrwGoRcyVwHT30DT+7+aefJ6Zz88/2I7H65/HjhSluZxvEfMePXowc+ZMPvjgA77//ntALWKuBLDpwzvg8nj5y+d7MOoFf5/YD51O+LosRWkRF7KI+eLFi3nttde4++671SLmSuB7+JpOPHl9Vz7ekstvPtmJV7X0lTbkXIuYA4wbN47PP/+cefPmNcsi5qqFr7S4n1zbBZfHy79XHcCgF/xxfG+EUC19pW1bs2YNH3/8MQ6Hg5tuuqlZ9qECX/GJn1/fFadHMvPrgwD8YVxv1b2jtGkjRoxgxIgRzboPFfiKTwgheGp0NwBmfn0Qr4Q/jVehryjNSQW+4jMnQl+vg1dWH0RKyfO39lGhryjNRAW+4lNCCH5xQzd0QvCfVQfweuGFCSr0FaU5qMBXfE4IwRPXd0UIwb9X7scjJS/+qC96FfqK0qRU4Ct+4UTo6wT8c8V+pIS/TlShryhNSQW+4ld+dl1XdELw8lf7kFLyt0n9VOgrShNRga/4nZ9e2wWdgL8v34dHSl6a1A+DXp0jqCiXSwW+4pd+PKoLQgj+9uVepISXb1ehryiXS72DFL/12MjOPDW6O4u3H+dnC7apWTaVgNcq1rRVlObyyIhO/Oam7izdkcfj87fhUqGvBCi1pq2iXIAZV3fit2N78NnOPH76/lYV+kpA8oc1bVXgKwHhgas61s+n/+P3tuB0q9BXAsuyZcsYO3YsdrudF198kbVr17Jjxw6WLl0KwE033URUVBQvvPACpaWlwMk59JuKOmirBIzpwzugE/Dckl089t4WXrlTLZeoXJxvP9hHcXZlkz5nTEoIV93e9Zzb+MuaturdogSU+4Z14LlxvfhqVwGPzttMrcvj65IU5bwarmlrsVjIzMxk2LBhzJgxg1dfffWsv9fUa9qqFr4ScO69Mg2dTvB/izKZ9vYGZt+TTqjZ6OuylABwvpZ4czqxpm2/fv3o0qULkydPZteuXfVr2gJUVVVxzTXX8Oyzz3LzzTerNW0VBeDuIe35xx392HSkjMmz1lFkd/i6JEU5p/OtaQvw4osvcvvtt9ffV2vaKkqd265IJtJi4pF3tzBx5g/MvT+D1GiLr8tSlEadb03br776ip49e9a3+NWatopymhHd2vHegxnYalxMeO0Hso7bfF2SopzVuda0XbNmDevWreO9995j9uzZhIeHqzVtFeV0V6RGsvDhodzz5gYmv76OWfekM7RTtK/LUpSL8vzzzwPwzjvvEBMTg07X9O1x1cJXWoXO7UJZ+MiVxIWbufftDXyRmefrkhTlkkybNo2bb765WZ5bBb7SaiRGBPPhQ0PplRjGo/O28P6GY74uSVH8SpMEvhBitBBirxDigBDi6UYeDxJCLKh7fL0QIq0p9qsop4u0mpj3QAZXd43l1x/v5D8r9yOl9HVZiuIXLjvwhRB64BVgDNATmCKE6HnaZtOBMillZ+AfwIuXu19FORuLycDse9KZcEUSL321j2cXZ+H1qtBXlKZo4Q8GDkgpD0kpncB84PSBo+OBOXW3FwLXCiHUMkZKszHqdfx9Uj8evKoDc9Ye5fEF29T8O21Ya/yWdyl/U1MEfhKQ3eB+Tt3PGt1GSukGbMAZwyiEEDOEEJuEEJuKioqaoDSlLdPpBM+M7cmvx3RnyfbjTJ+zkUpH43ORK62X2WympKSkVYW+lJKSkhLMZvNF/Z5fDcuUUs4CZgGkp6e3nn8dxaceuqYTUVYTT3+8kztnr+PtaYOIDgnydVlKC0lOTiYnJ4fW1og0m80kJydf1O80ReDnAikN7ifX/ayxbXKEEAYgHGi6Wf0V5TwmpacQaTHx2HtbmDRzLXPuH0xKlDorty0wGo1NOgFZIGuKLp2NQBchRAchhAmYDCw+bZvFwL11tycCq2Rr+n6lBITresYx74EMiisdTJz5A5m56qxcpW257MCv65P/MfAlsBv4QEqZJYT4gxBiXN1mbwLRQogDwBPAGUM3FaUlpKdF8eHDV6ITgokzf2DR1tO/jDYjrxcqi8CeD141rbPS8oS/NrTT09Plpk2bfF2G0koVVzp4bN4W1h8u5f5hHfj1Td0x6i+x/SMl1JSBPa/ukl93XdDgfj5U5oO37qCxzgChCdolLBHCkuquE07eDokHg6np/milTRBCbJZSpjf2mF8dtFWUlhITEsS7D2Tw52W7eev7w2Qdt/HK1AHEnOtgrtcLx7fA3mVQcqBBsOeDx3nm9uaIulCPh5iu2nVoAgih/V7FcajIhYJM2L8cXNVnPoe1XYMPhASI7ACdr4XY7trzKMpFUC18pc37ZGsOT3+0kyiridfvHkjf5IiTD7qdcPQ72L1UC3p7Hgg9RHXUAvhEoJ9+HRIHxuALL0JKqLXVfQgcB/vxkx8IFQ0+HGrLte3DU6HrDdB1NKQNv7h9Ka3auVr4KvAVBcjMtfHQ3M0UVTr4y80dmRC2Rwv5fV+CwwaGYK1l3eMW6HIDWKJ8U6gtV/s2sO9LOLQG3DVabR1HaB8AXW6E8NNPg1HaEhX4inI+VSVU7lzC3jXv06tmM2bhQgZHIrrdBN3HQseRYPKzYZyuGjjynRb++7+E8rrJ4uL6nGz9Jw0End63dSotSgW+ojSm/Bjs+UxryR/7AaQXGZ7MJvOVvJTdFW/yEP579yDahV7c2Yw+ISUU7akL/+VwbB1ID1iiofN10PVG6HQtBEf4ulKlmanAV5QTpIQ9S+Hbl7UDsADtemqt+O43Q0I/EILF24/z1MIdhJoNvHbXQAa2j/Rt3RerpgwOrNTCf/9XUFOqjQzqeSsMfVRr+Sutkgp8RZFSC77Vf4K87RDdGQbco4V8dKdGf2V3XgUPzd1Mnq2G58b15s6M1BYuuol4PZCzCbI+ga3vgtMOKRkw5FHt79erwXqtiQp8pW079DWs+hPkbICI9jDi19Bn0gUFXXm1k8fnb+PrfUVMGZzCs+N6EWQI4D7x2grYNg/Wz4SyI9pon4wZcMXdqrunlVCBr7RNx9ZpQX/kW20c+9W/hCvuAr3xop7G45W8/NVeXll9kP4pEcy8ayDx4QHQr38uXg/s/RzWvaYNOzVatdcm46GzfuNRAoMKfKVtOb4VVj0PB77STly66kkYOA2MlxfSX2Tm8eQH2wk2Gfj35P5c2Tmmaer1tbztWvDvXKidCdx1tNbPn3aVOrkrAKnAV9qGgixY/WftoGxwJAz7GQx+EEzWJtvF/gI7D83dzKHiKu4aksrTY3oQEtRK+sDtBbDpTdj4JlQXQ1xvGPII9J542R+WSstRga+0bsUHYM0LkPkRBIXC0B9rQWUOa5bd1Tg9vLR8L29+f5jE8GBemNCHq7vGNsu+fMJVCzs/1Fr9hVlgjYX06TB4BljPWLdI8TMq8JXWqewofP1X2P4eGMyQ8TBc+ZMWOwt289EyfrVwOweLqrg9PZlnxvYkPPjijg/4NSnh8Nda8O/7AkyhcOWPYehj2ger4pdU4Cuti8cF374E3/wdhA4GPQDDfw4hLd/KrnV5+NfK/cz65hAxISb+fFsfru0R1+J1NLvCPdqQ1t1LtJO5rvoFpN+vunr8kAp8pfUo2gsfz4C8bdrQyuv/oM0m6WM7csr51cId7Mm3c9sVSfzu5p5EWlvh1Ma5m2HlH7R5fMKSYcTT0G+KGsvvR1TgK4HP64X1r8GK57SDsLf8E3qO93VVp3C6vfx39QFeXX2ACIuJP93ai9G9E3xdVvM49DWsfE77AIjuAqN+q/17qFE9PqcCXwlsZUdh0aPaePGuY+CWf0Go/3ab7DpewS8XbifreAVj+yTw3Phe555nP1BJqc1FtOqP2jw+Cf3h2t9Bp1Eq+H1IBb4SmKTUpgL44tfa/TF/gf5TAyJMXB4vs745xL9W7McapOfZcb0Y1y8REQC1XzSvB3Z8oA2JtR3Txu9f+3tIGeTrytokFfhK4LEXwJKfaqND0q6C8a9AZHtfV3XR9hXY+eXCHWzPLuf6nnE8f2tv2oW10gOdbgdsngPf/A2qCqHbTTDq/yCup68ra1NU4CuBJWsRLP25tuTftb/XhlvqLnG9WT/g8Ure/O4QLy3fR5BBx69v6sGkgckYLnUNXX/nrNKGcn7/b3BUQN/bYeRvIDKt2XYppaTGXUOVq4pKV2X9tdPjJCIogihzFFHmKIINwa3zW1YDKvCVwFBTBst+BTs/gMQr4LbXIbabr6tqMoeKKnn6o51sOFJKhxgrj1/bhVv6JaLXtdIAqi6F7/8J61/Xun0Gz4Crf3FB50m4vC4O2w6zp3QPOfackyHurDwl1Bve9krveZ/XrDcTaY4kyhxVfx1tjj7lZw3vmw2B921MBb7i/w6shE9/rHUFXP1Lbf6bi5zkLBBIKflqVwEvf7WPPfl2urQL4WfXdWVM73h0rTX4K/Jg9fPaLJ1BodoY/sEz6sfwV7uq2Ve2jz2le9hTuofdpbs5UHYAp/fkwvAWg4UQYwhWk1W7Nja4Np16v+Fto96IzWGjtLaUstoySmtLz7zUlJ6yr4bah7Wnd0xv+sT0oXdMb7pHdSdI798H4FXgK/7LWQXL/0+bwyWmG0x4XWvdt3Jer+TzzHz+sWIfBwor6ZEQxs+v68L1PeNab5dDQRaly3/DnuPr2B3Wjr2JvdjtreJoxVEkWg5FBEXQPap7/aVHVA9Sw1Ix6JpvnL+Ukmp39SkfAGWOMgqqC9hTsoedxTspqikCwCAMdInsUv8B0CemDx3CO6D3o2UkVeAr/ilnM3w0XZuXfehj2lhuY7Cvq2pRHq9kyfbj/HPFPo6UVNM3OZyfX9+VEV1jAz74pZTsL9/Pmuw1bC/azp6SPRTWFNY/nuhy010XTPeON9K94w30iO5BnMU/P/AKqgrILM5kZ/FOMksyySrOotJVCWjfPnrF9Dr5TSC6N/HWeJ/9HSrwFf+zdR4s/RmExMNtr0HacF9X5FNuj5ePt+Tyr5X7yS2vYWD7SJ68vmvATcHs9rrZWriV1dmrWX1sNTmVOQB0juh8Squ9W0RXwvd+ro3hr8jVzq+4/rmAOWbjlV6OVBzRPgSKdpJZnMnesr24vC4Aos3RZCRkMCp1FMOThmM1Nt2MreejAl/xHx43LP+tdtZsh6th0pwWm+wsEDjdXj7YlM1/Vx0gv6KWIR2jePKGbgxK89/XqNpVzdrja1mVvYpvcr6h3FGOUWdkSMIQRqaOZETyCGItZ5nnyFWjjej57h9a996Ae7QVyfz4xLqzcXqc7Cvbx87inewo2sH3ud9T5ijDqDMyOGEwo1JGMTJl5NlfizpSSrw2G/qIiEuqQwW+4h+qS+HDe+HwN9p6qtf/Uc3Bcha1Lg/vbzjGK6sPUlzp4KouMTx5Qzf6p0T4ujQASmpK+Drna1YdW8W6vHU4PA7CTGFcnXw1I1NGMixp2MW1aquKtZlPN70J+iAY9rg2M2cTrmXQ0jxeD9uKtrH62GpWZa8i254NQJ+YPoxK1cK/Y3jH+q4f6fViX7GC4pkz0YeH0/7tty9pvyrwFd/Lz4T5d4I9X5sHp/+dvq4oINQ4Pcxdd4TX1hykrNrFiG6xTB6Uyqju7TAZWnYc/xHbEVZlr2L1sdVsL9qORJJoTWRk6khGpoxkQNwAjLrLHFlVchBWPAu7F2vdfSN/oy296EcHRS+FlJKD5QdZnb2aVcdWkVmSCUBqaCojk67h+kNWQt77Euf+/RjbpxLz0MOE33brJR0HUIGv+FbWIlj0CJjD4Y55kDzQ1xUFnEqHm3e+P8w7PxyluNJBhMXIuH6JTBiQTL/k8GY7QFheW87ig4v55MAnHCg/AECPqB6MTB3JqJRRdI3s2jz7PrZe6/rL2QCx3bVunh7jAvoEvIYKqgr4+ugqchd9QJ/P9pJUIsmL1XP41oF0nHA3Q1KGEWy4tAEMKvAV3/B6tfHX3/4dkgfBHe9CaLyvqwpobo+Xbw8U8/GWXJZn5eNwe+kYa+VHA5K59YokkiIuf5STlJIthVv4cN+HfHXkK5xeJ31j+zK2w1hGpowkIaSFZgCVEnZ9qv0fKt4HcX20Fn+3MQExn9LZSJcL25KlFL8+E9fRYxi7dOb4pCv5LLWUb49/j91lp3NEZz4Z/8klPX+zBb4QIgpYAKQBR4DbpZRlp23TH3gNCAM8wPNSygXne24V+AGu1qbNW7/vC7jibhj7Ehj8+4SVQFNR6+LznXl8tCWXDYdLEQKGdoxmwoBkRveOv+i1dm0OG4sPLmbhvoUcsh0ixBjCzR1vZlK3SXSN7NpMf8UF8Hq0Bda//guUHoLEATDyGeh8bUAFv3Q6Kf9kESWzZuHKzSWoZw9iHnmE0GuvRdR9c3F5XGwq2ESVq4rr2l93SftpzsD/K1AqpfyLEOJpIFJK+dRp23QFpJRyvxAiEdgM9JBSlp/ruVXgB7DiAzB/itYfO/ov2kLiAfTGDETHSqr5ZGsuH2/N4WhJNcFGPaN7xzNhQBJXdoo56/QNUkq2Fm7lw30fsvzIcq01H9OXiV0ncmPajViMlhb+S87B44bt72sHd23HICVDa/F3uMav/395HQ7KFy6k5I03ceflYe7bl5hHHiZkxIhm6Q5rzsDfC4yQUuYJIRKANVLKcw6kFUJsByZKKfefazsV+AFq/1ewcLo2+mbSHOhwla8ralOklGw5VsbCzbks3XEce62b+DAzt16RxG1XJNE1LgQhBDaHjSUHl7Bw30IO2g4SYgxhbMexTOo6iW5Rfj4W3u2ErXO1JS7tx6H9cBj1DLS/0teVncJbU0P5Bx9oQV9URPAVVxDz6KNYhw9r1pOymjPwy6WUEXW3BVB24v5Zth8MzAF6SXnumY5U4AcYKbWJslY8B3G9YfK8gJzOuDWpdXlYubuQj7fksGZfER6vl+jo40TEbaGEjXikk97Rvbm92+3+15q/EK5a2DJHW9+4sgA6jtS6enw8D7+3tpay996n5M038ZSUYBk8mJhHH8GSkXHWoPd6JWXVTgrtDgrtDgw6wbBLPOnusgJfCLECaOxI2zPAnIYBL4Qok1JGnuV5EoA1wL1SynVn2WYGMAMgNTV14NGjR89Zm+InnNWw+MeQ+RH0uk2buz6Ax0+3Nk6Pk/d3fcI7mXMpdh4FbxBO2xW4ygYTLFMY0D6S9PZRDEqLpH9qBBZTgJ0b4azWxu9/9w+oLoEuN2hdPS08J5N0uyn/5BOK//sK7oICrFcOJfKhh6nu0ZfCCgeF9lot0E+5rV0XVzpweU5mcd/kcBb/+NLOPvd5l44QIgwt7P8spVx4Ic+tWvgBovwYzJ8K+Tvh2v+D4U/4dX9qW1LhrODDvR8yb/c8imqK6BbZjSndpzCmwxjKqgSbjpSy6UgZG4+UsrfAjpRg0Al6JYUzqH0k6WlRpKdFBs7yjI5K2DALfvi3NtV2t7Ew8tcQ36fJd+X1SkqrneTbaimsqKF21UpiF7xNSGEuxxM78emgW/khNI2SKgeNRWyU1US70CBiQ4NoF2qmXVgQ7RrcTgg3kxx5ad+4mjPw/waUNDhoGyWl/NVp25iAz4ElUsp/Xuhzq8APALlb4L3btZWOfvQGdL3R1xUpQH5VPu/uepeF+xdS5apiaMJQpvWextCEoWftUrDVuNhyrIxNR0rZeKSMbdnlON1ar2vHGCvpaZF0bhdCUoSFxAgzSZHBxFiD/HNK59oKbbqGta+Aw6ZN4THoAW0FrvNMue3yeCmrclJS5aS4UmuN51fUUlhRS0GD24V2B26vpH/RfqZlLaNbeTZHQ+P4eMA4cnum0y7MTFyYWQvx065jQoKa9aS55gz8aOADIBU4ijYss1QIkQ48LKV8QAhxF/A2kNXgV6dJKbed67lV4Pu5vV/AwvvAEgN3LQyYSa9as/1l+3kn6x2WHVqGRHJD2g3c1+s+ekT3uOjncrg9ZOZW1H8AbD5aSlm165RtTHodiRFmEiOCSYoI1q4jtdtJEcHEh5sxG314hmxNGd6NbyE3vYW+IgdncByH025nZ9x4sl3hlFQ5KKnUwr2k0kFJlZPy0/7GE8LMBuLqQjwuzExXWw79vphHeNZWZGwcloceIfH2CZhMvl/DQZ14pTStjW/Asl9CfF+484OAnOiqtZBSsqlgE29nvs23ud8SbAhmQpcJ3N3zbpJCkpp0XxW1LnLLajheXkNueQ25ZXXX5drPCu1ndl/EhgaRGBFMpMWIUa/DZNARpNfV3zYZTt4OMugw6gUmvQ6jQYep7uegTTFR4/JQ7fRQ46y7drmpdp782cnH3fU/q3K6QXoZqdvKPfqvuEa/A5fU86V3EIuMYzhi7U90iNbqjg4xEWU1afet2u0TAR9s0j64HIcPU/Svf2P/4gv0ERFEP/wQkVOmoAvyn24vFfhK0/B6YeWz8P2/oMuNMPEtCArxdVVtksfrYeWxlbyd+TaZJZlEmaO4s/ud3NHtDiLMET6pyen2km+rJae8muPltad8OFTUunC6vTg9XpxuL676a1n/8wulE2AxGQg26bGY9AQb665NeoKNBiymk/etJgPRIVqIR1tNxHtySdj3HsFZ8xG15RDbAwZNh753gDnsrPt0FRRQ/MqrlH/0ESIoiOhp04i6/z70If73/18FvnL53A5tPpzMjyD9fhjzNzXTpQ/Uumv59MCnzNk1h2x7Nqmhqdzb617GdRoXkOuvniCl1ML/tA8Ep8eLlNQHe7BJT5BBd/nj2J3VkPUxbJgNedvAFAL9JkP6dIjrWb+Zx2ajZPZsSue+i/R6ibzjDmIefghDjP+uU6ACX7k81aWw4C44+j1c95w2da0aidOi7E478/fM593d71JaW0qfmD7c1/s+RqWM8qvl9QKOlNrgg42zIfNj8Dig/TC8fe+hdGMZJW+9g9duJ3zcLcT85CeYkpN9XfF5qcBXLl3ZEZg3Sbu+9TXoM9HXFbUp5bXlvLv7Xd7b/R52l53hScOZ3ns6A+MG+uVSgAGtqgS5aQ7l896keIMLd62ekB6xxM64C/OoKdoC7AHgXIGvvpMrZ5e7Bd67Q2v13L0I0ob5uqI2o7immP/t+h8L9iyg2l3NdanX8WDfB+kZ3fP8v6xcNCkl9m83UvSP5TiPegnu3pmkK41Y3Jtg/S9h428gdQh0vg66XA/tegbkt1wV+ErjGg67nLZUDbtsIflV+byT9Q4L9y3E5XUxOm00D/Z5kM6RnX1dWqtV9cMPFL70MrVZWQR16Uzyq68QMnKk9g3K7YTs9XBghXZZ8XvtEpqozdbZ+TroOAKCI3z9Z1wQ1aWjnGnjm7DsF2rYZQvKtmfzVuZbLDqwCCTc0ukWpveZTvswNR9Rc6nJzKLo5Zeo+mEthsQEYn/yU8LH3YLQn+OYSMVxOLASDnwFB9doJ3YJPaQM1sK/83Xa+8aHC7WoPnzlwni9sPI5bRI0NeyyRRy2HeaNnW/w2aHP0AkdE7pM4P7e95MYkujr0lot55EjFP7rX9g/v8yx9B435Gysa/1/BXnbtZ9b22mt/8QBENUBIjtARCoYTE3/xzRCBb5yfmrYZYvaW7qXN3a+wZdHviRIH8SkbpOY1msa7SztfF1aq+UqKKT41VcpX7iwbiz9vUTdf3/TjaWvLKxr/a+Ag6ugpvTkY0IHYckQlaZ9AJz4IDhxfY5zAC6WOmirnFtNmTYB2tHv4bpnYdjPAvKAVCDIKs7i9R2vszp7NVajlel9pnN3z7uJMkf5urRWy1NRQcnsNyidOxfpdmtj6R99pOnH0oe0g/5TtIuU2gdA2WEoPXzq9Z7PoLr41N+1RJ/6ARDfG3qOb9r6UIGvlB2tG3Z5GH70php22Uw25W9i9s7Z/HD8B8JMYTza/1Hu7H4n4UHhvi6t1fLW1lI2bx7Fs2bjtdkIGzuW2Md/iik1tfl3LoR27Cs0Thvdc7raCm2o8+kfCNnrtW/ZKRkq8JUmlru5btilE+7+BNIubf5tpXFSSr4//j2zd8xmS+EWosxR/Hzgz7mj2x1YjWq9gObiramhbMECSt54E09xMdarrqLdEz/H3OPiJ5FrNuYwSOirXU7ndmprQjcDFfht1Z7PtKUIQ2Jh2mdq2GUT8kovq46tYvbO2ewq2UW8NZ5fD/41E7pMCOjpD/ydt6aGsvkLtJWmiouxZGQQ+4+XsQzy7QpYF81g0t6XzfHUzfKsin9bNxO+eBqSBsCU+Vrfo3LZ3F43nx/+nDd3vslB20FSQ1N57srnuKXjLRjPMw+7cum81dUng76kBMvQIcT+8x9Y0hs9btmmqcBvS7weWP5bWPeqthrQj94AU4CtY+qHnB4nnx78lLd2vkVOZQ6dIzrz4lUvckPaDRh06i3WXLzV1ZS9P18L+tJSrFcOJeaxx7AMHOjr0vyW+t/YVjir4eMHYc9SyHgEbnwe1KRbl6XGXcPCfQt5J+sdCqsL6R3dm18O+iUjUkagE7478aa104L+fUrefKsu6K8k5sePYRkwwNel+T0V+G1BZRG8f4c2N87ov8CQR3xdUUCzO+0s2LuAubvmUlpbSnpcOn8c9sdzLiGoXD5vVdXJoC8rwzpsmNaiH9Cyi5UHMhX4rV3RPpg3URsTfMe70ONmX1cUsPKr8nlv93t8uO9DKl2VDEsaxow+MxgQp1qWzclbVUXpe+9R+tbbWtAPH07MY49iuUIF/cVSgd+aHfke5t+pLdw87TNIVn2bl2Jv6V7mZM3h88Of48XLDe1vYFrvafSK7uXr0lo1d3ExZQsWUDb3XTzl5VivuorYxx4luH9/X5fWrKRX4qh2Yw5p+gP9KvBbqx0fwqePQkR7bZHxyDRfVxRQpJSszVvLnKw5/HD8B4INwUzuPpmpPaaSHOr/i2AEsprMLMrmzqVi2TKky0XINdcQ8+gjBPfr5+vSmoWz1k3hUTv5B23kH9Iu0Ukh3PZk039zVIHf2kgJ374Eq/4I7YfDHXPBok7bv1Auj4svjnzBO1nvsK9sHzHBMTw+4HEmdZ2kzoptRtLlwv7VV5TOfZearVvRWSxE3HEHkVPvJKhDB1+X12SklNhLa7VgP1hB/iEbxTmVSK82p1lkgpVOV8SS2DWyWfavAr818bjgsydgy/+gzyQY/woYLnIGwDbK7rTz0b6PmLt7LoXVhXQK78QfrvwDYzuOxaRvmVkO2yJ3aSnlH3xI2fvv4y4owJiaStxvfk34bbehDw2MFabOxeP2UpTdoPV+0EaVzQmAIUhPXFoYA0e3J75jOHEdwjBbm/d8DRX4rUVtBXx4rzZL39W/hJHPqAnQLkB+VT7v7nqXhfsXUuWqIiM+g2eHPsvwpOFqxE0zqt29m9K571KxdCnS6cQ6bBjxzz1LyNVXI3w4l/zlkF6JraiGomN2Co/ZKThso/CIHY/bC0BotJnErpEkdAonvmM40UlWdPqW/VtV4LcGZUe1g7OFu2Hcf2DAPb6uyO/tKd3DO1nv8OXhL5FIbki7gXt73asOxDYj6XZjX7mKsrlzqd60CREcTPiPJhA1dSpBnQNrRS+vV1KeX01Rtp2iY3WXbDuuWg8AOoMgNiWU3iOSSOioBbw1wvfftlXgB7oDK+Gj6driJVM/1BZeUBrl9DhZfnQ5H+z9gK2FW7EYLEzpMYW7etylFhxpRu6yMmwffUTpe+/hPp6HMSmJdk89RcSE29CH+/9xEa/HS1l+NYVHtVAvOmqnOMeO26m13PVGHTHJIXTLiCc2NZTY1FCiEq3oW7j1fiFU4Acqrxe+exlW/Qna9dDG2Ed38nVVfinHnsOH+z5k0YFFlNaWkhqayi/Sf8FtXW4jzNR0C08oJ3kqq6hcvYqKz5ZR+f334HJhGTKE+GeeIWTEiHMvI+hDtVUuSo9XUpJbRenxKoqy7RTnVOJxaeFuCNITmxJCz2GJ9eEeGW9p8a6ZS6UCPxDVVmirU+1ZCr0nwrh/g0lNt9uQx+vh++PfM3/PfL7L/Q4hBCNTRnJ7t9sZkjBETX3QDLy1tVR+/Q0Vy5ZRuWYN0uHAkJBA1D13Ez5uPOZuXX1dYj1nrZvSPC3US3OrKM2rpOR4FdV1B1QBTGY9MSmh9L4midiUUNq1DyW8nQWdLnCP7ajADzSFe2DBVG2xhNF/gYyH1cHZBkpqSvjkwCd8uPdDjlcdJyY4hof6PcSPuvyIeGu8r8trdaTTSdXatdg++4zKFSvxVlejj44mYuJEwsbeRHD//j49COt2eSjLr9aC/bgW6qXHq7CX1NZvYzDqiEq0ktojiqjEEKKSrEQnWrFGBLW6A/cq8ANJ1iew6DGtNX/vEkgb5uuK/IKUkq2FW1mwdwHLjy7H7XUzOH4wT6Y/ycjUkRh1amripiQ9Hqo3bqTis2XYly/HY7OhCw8nbOxNhN10E5ZBgxCGlosW6ZVUljsoL6g+eSnUru0ltZxYtlunF0TGW4jvGE7PYYlEJVqJTrISGh0c0K32i6ECPxB43LDyWfjhP5A8GG6fA2HqIGOVq4qlB5cyf+98DpQfINQYyuRuk5nUbRIdwzv6urxWRUpJzbZtVCz7nIovPsdTVIywWAi99lrCbhpDyLBhCFPzna8gpaS2ykV5QU19oNtOBHthTX0fO2j97BHtgmmXFkbXjHiiEqxEJ4YQHhfslwdSW5IKfH9XWQQL74Mj38KgB+HGP2sr4rRRXullY/5Glh5ayvIjy6l2V9MjqgfPXfkco9NGYzGq+f2biis3l6r1G6hev56q9etx5+cjTCZCrrmGsLFjCbnmanTBwU22P69XUlXuoKK4pu5SS0VxDbYiLeQd1e76bXU6QVhsMBFxFlJ6RBERZyGinYWIOAuWcFOr64ppKpcV+EKIKGABkAYcAW6XUpadZdswYBewSEr548vZb5uRswk+uAeqS+DWmdB/iq8r8pn9ZftZemgpnx36jILqAqxGKzek3cCkrpPoE9NHvcGbgKugoD7cq9dvwJWTA4A+MhLL4MGE/vxnhFx7LfqQkEveR22VC3tJLbaiulAvqa0PeHtJLV6PrN9WCAiJNBMWa6ZzehwR7YLrgz00xtzmW+uX4nJb+E8DK6WUfxFCPF13/6mzbPtH4JvL3F/bICVsfgc+/xWExsP05ZDQOieOOpei6iKWHV7G0kNL2VO6B73QMyxpGL9I/wUjUkao9WEvk7u4uD7cq9evx3n0KAC68HAsg9KJuuceLBkZBHXpfEEHXqWU1NhdVJbVYi+txV5SS2WpA3tZbX2gN2ylAwRZDYTHBBOTHEqnK2IJiwkmLDqYsFgzIZFm9AYV6k3pcgN/PDCi7vYcYA2NBL4QYiAQB3wBqIUmz8VVC8uehK3vQufrYMLsNjX5WbWrmpXHVrL00FLW5a3DK730ju7N04OfZnTaaKKDo31dYsByl5Vp4b5hPVXrN+A8eBAAXUgIlvR0IiZPxpoxmKBu3RodJ+9xeaks14LcXuo4NdjLHNhLa0/pSwcwmHSERpkJjTYT3zGcsJhgwmOCCY0xExYTTFCw6lVuSZf7asdJKfPqbuejhfophBA64CXgLuC6cz2ZEGIGMAMgNTX1MksLQOXHYMHdkLcNrv4VjHi6TSxD6PF6WJ+3nqWHlrLi2Apq3DUkWhOZ3ns6N3e6WR2AvUjemhocBw/h2L//lIs7Px8AYbFgGTiQiNtuxZKRgblHD7zoqCp3UFrmoHJLEZVlDirLHFSVacFeWeagusJ5xr4sYSZCosxEJ4WQ1ieakCizFvB1lyCrQXW3+ZHzBr4QYgXQ2ADmZxrekVJKIYRsZLtHgWVSypzz/cNLKWcBswDS09Mbe67Wa/cSWPxTbaHxKfOh2xhfV9SsvNJLVnEWXx75kmWHl1FUU0SoMZSbOtzELZ1u4Yp2V6iTo85Dulw4jx6tD/Tafftw7N+P61g2J8YiCpMJfaeuyPSrIakL7qQuOENiKaxwa0G+pJbKueuoaSTMTWY91kgzIZFBRCeFNAjzIO12pBm9Uf0bBZLzBr6U8qytciFEgRAiQUqZJ4RIAAob2WwocJUQ4lEgBDAJISqllE9fctWtiS0Hlv0K9n4G8X1g0pxWO0WCy+NiQ/4GVh1bxZrsNRTWFGLQGbgq6Spu6XQLVydfTZDe9xNM+RPpdOIqLMJdkI8rPx9Xds7JgD98GLc04DBF4AiOxJPYGXfn8TgHxOEwhVMjzVRXeamtckMlsBfYWwVUYQo2EBIZREhkELEpIfXBHhIZREiEdtukultancv9F10M3Av8pe7609M3kFJOPXFbCDENSFdhj9aS3zBLmwvH64Hr/wBDHtWWI2xF7E473+Z8y+rs1Xyb+y1VriqCDcEMTxrOyJSRXJ18dZtdWMTrdOIuKMCVl6dd5+fjzi/AVZCPKy+fqpIqqis9OIIicQRF1F+cYT1wJF5FbZIFjzyty88DwW4j1pAgwiKCSIjQQtwaoV1C6oLdZFZh3hZd7r/6X4APhBDTgaPA7QBCiHTgYSnlA5f5/K1T3nat+yZvm3ZgduxLrWoJwvyqfNZkr2HVsVVsLNiI2+smyhzF6LTRjEodRUZCRqtryUsp8VZV4Skvx1NWjsdm026Xl+OxleMpt9XfdhaXUlVSTXWNOCXIHUEROCyJOMy9ccSEIGNP7S4ROrDW9ZlHRpgJiQjCGhl06nV4kOpmUc5KSOmfXeXp6ely06ZNvi6jaTkqYc0LsO5VsMTAmL9ArwkBPxeOlJID5QdYdWwVq7NXk1WSBUBaWBojU0cyKmUUfWL6oPfzA9DS5cJTWYnXbsdTYcdbacdTUYHXXonHftr1iUBvcI3bjUdX18USFHFqy9waiyM4ilpTOE7dmSeHGYxCa4FHmQmJMNe1xoNOubaEmhBtZAoA5dIJITZLKRsdDam+17WUfV/CZ0+CLRsGToPrnoXg5lm3siVUOCvYnL+Z9fnr+Tr7a3IqtZN0+sb25fEBjzMqdZRPRtd4nU68NpsWxDYbHlsFHpsNb8Wp9z0VNrwVdjyVdu3abkfW1JzzuT06I67IRJzhCTjD4nBaO1AbEUltpzBqRTA1niAc7jM/1IKCDVgjg4g40QqvC3ZrxMn7QRY1mkVpfirwm5s9XzuBatenENsd7vsC2g/1dVUXrcZdw9bCrWzI28CG/A1klWThlV7MejPp8enc3+d+RiSPINYS22T7lFLitdlwl5TgLinBU1KCu6QUd0kxnpJS3KUleMtteCpOhHjFuUNbCHRhYejDw9GHhaEPC8XQrh26sFBESDi15ggchnBqdBZqpYVqt5Fqh56qakmV3aMd/DyNOdiINTKI8Mggkupb4ycPgFojVH+54j/U/8Tm4vXC5rdgxXPgdsCo38KVjwfMPDgur4vM4kzW5a1jQ94Gthdtx+V1YRAG+sb2ZUbfGWTEZ9A3tu9FL/LtrarCVVCAOz8fV2FhfZB7SorrAr0u3EtLwX1myKLToY+KwhAVhT48HFP71Logj9CCPEILdF14OPqwcERYGA4RTLXLQGWZE3uDseWVpXVjzI864bTezSCLJCRST0iMmfgudQc8o7RW+YmDnwaTf3dTKUpDKvCbQ0EWLHkccjZCh2vg5n/4/VBLj9fD3rK9bMjbwLr8dWwp2EKNuwaBoHtUd+7qcReDEwYzoN2Ac05Q5qms0oYQ5uXXDyV0558YgZKPq6AAb0XFGb8ngoIwREejj4nBGBeHuVdPDFHRGGKi0dddG6Kj0UdHo4+IqD/VX0qJo8qN/ZQAr8We56ByVy2VpRVUlRfh9Z6a5oYgPaGRWtdKdHJIfYBrFzWSRWmd1P/opuSshm/+qk1jbA6H216Hvnf45UFZu9POrpJd7Czeyc6inWwq2ESFUwvijuEdGd9pPBkJGQyKH1Q/bFJKiaekhJqcfThzcnHl5uLKyakbVqiFvLey8ox96aOjMcbHY0xNxTJoEIaEeO1+fDyGdu3QR8egs1oa7cN21rqpLHNgK6ulssCBfU85lWUF9S3zyrLa+rVFT9DpRX1wJ3QJJyRSO2HoxM9Co7Qx5qrPXGlrVOA3heIDsPlt2DYPasqg/1S4/o9g9Y95X5weJ/vK9rGzeCeZxZnsLN7JEdsRZF0fRkpoSv1wyXRrDyJKHDhzcnDtyqYmZy0Vubk4c3Nw5R4/o49cHxWFMSEBY/v2WAZnYEyIxxAXjzE+DkNCgtZHfpZ50t0uD1XlDipzHVSWVdTPz1LZoMvl9Mm2ECeHJkYnWWnfO/qUMA+JUqNZFOVsVOBfKrdTOzt201tw+BvQGaD7zTDkEUgd4rOyvNLL0Yqj9cGeWZzJntI9uLwuAGINkQzVdeEu0ZuOVRbalYPueCHOnCxcucspq6ig4fzWupAQjMnJmNLSCBk2HGNyMsakJIzJSZiSktBZG19L1+X0aPOwHKqksrzBvCzlWphXlTuosbvO+D1ziJGQyCBCo4NJ7ByhDVOs63o5cRBUTYurKJdGBf7FKjsCm+fA1rlQVQQRqXDt76D/XRB6xtxxzarGXcMR2xEO2g5yoOwAmSWZZBVnISvsxJVDSoWRq52xTKtMJbbci6XQjje/ELxF9c9RHRSkBXhKMpYrrjg10JOT0YefehbsiZWHKsudVB2pparcprXSyx1UljqoKj9Ly5yTYR4SEURch/D6IYmhUeogqKK0BBX4F8Ljhv3Ltdb8gRVan3zXMZB+H3Qa1ewzWlY6KzlkO8TB8oMcth3maOE+yo/uQx4vJNYmibVJ4m2CKXYTMWUeTNWeE4UD2ehjYjAlJ2Mc2B1TSgrGlBRMKckYU1IxxMYgdDqklDhr3FSVO7HZHFTZHFStLaPKVkB1ed39cidVFQ687jNP1gsOMxESEURYjNYyt57oYmlwFqgKc0XxLRX451JxHLb8T2vR249DaAJc8xQMuBvCk5t8d+W15RyyHeJQ4W4KDmRScfQAjtxszIUVtLNBrE1ypQ3GVJ/2iyYjxsRETJ1SMSWnYExNwZSSgiEpGREbT63bSLXdSU2Fk9IKp3Z7r5PqjfnU2I9RXeGk2ubEfdpc5gCmYIM2D0u4icSuEVjDg7BGmOqug7CEm7CGqdP5FSUQqMA/ndcLB1dpB2H3fg7SA52uhZv+Bl1Hg/7SXjKX10VBVQH5+Qcpyd5HxfGjVOfl4C4sRJSUYSytItzmItYGvaqhV8OSjHq8cdEYU5MJSukG7VLwRsXhCY3BbYnAKYKoqPZowV3hpDrfSc0+JzX2PNyu3DOLEWC2GrGEmQgONRHXIRxruKku2LVAt4Rrt41BqlWuKK1F2wx8KbV1YksPnXkpOQi15dpcN8N+CgPuhagOZ30qt9dNhb0YW1Eu9tI8KgpyqDh+lNr843gKixAl5QSVVRFicxJZCRa3wGQw005vxm0w49GbqbYm4wiLxJ0cSXHveMrC4hHBUbj0FpxeIw6npLbKjaPKhSwGik/svbLuAkInCA4xEhxmwhJmIjLOot0ONWEJO/nz4FATwSFGdOrAp6K0Oa038KXUpjU4LdC9pQdxlxzFWWvHJQUurw631OOytMNpicMdeiXV8e2pMiRSs6cC5/rXcdmq8VY5kNVORI0X4ZTonTr0Lh1GjwGBCY/eiFdnxKM34dGH4Db0xWUw4wyyUJZiwZsWjFcE4eU80x+7wGDXYfYaMYcYMVuNhLbTrk/cN1sNBJ1y30hQsEENRVQU5ZxaXeDnZW7gs7/vR6IDdHUnPQUj6QOiHyCQQocUOhA65In7pdr9cwqqu5yVRG+QmMx6zNYggiwmwoINGM0GTMF6TGYDpmADJrO+7vq028F6zFajOripKEqzaHWBb4qIQ893IACd0Fq9eh1Cr6u/L3QCoRcInQ6hF+j0OoRBr13rdRjNQZjDQrCEhWGNiCAkIhKT2YTepMNg1GMw6jCYdOgNOgwm7b7epFPjwxVF8WutLvCjk9tz3ztP+LoMRVEUv6OapIqiKG2ECnxFUZQ2QgW+oihKG6ECX1EUpY1Qga8oitJGqMBXFEVpI1TgK4qitBEq8BVFUdoIIeWZc5v7AyFEEXD0Mp4ihgbTjPkRVdfFUXVdHFXXxWmNdbWXUsY29oDfBv7lEkJsklKm+7qO06m6Lo6q6+Koui5OW6tLdekoiqK0ESrwFUVR2ojWHPizfF3AWai6Lo6q6+Koui5Om6qr1fbhK4qiKKdqzS18RVEUpQEV+IqiKG1EQAe+EGK0EGKvEOKAEOLpRh4PEkIsqHt8vRAizU/qmiaEKBJCbKu7PNBCdb0lhCgUQmSe5XEhhPh3Xd07hBAD/KSuEUIIW4PX63ctVFeKEGK1EGKXECJLCPF4I9u0+Gt2gXW1+GsmhDALITYIIbbX1fVcI9u0+HvyAuvyyXuybt96IcRWIcTSRh5r2tdLShmQF0APHAQ6AiZgO9DztG0eBWbW3Z4MLPCTuqYB//XBa3Y1MADIPMvjNwGfoy0QOQRY7yd1jQCW+uD1SgAG1N0OBfY18m/Z4q/ZBdbV4q9Z3WsQUnfbCKwHhpy2jS/ekxdSl0/ek3X7fgJ4r7F/r6Z+vQK5hT8YOCClPCSldALzgfGnbTMemFN3eyFwrRBC+EFdPiGl/AYoPccm44H/Sc06IEIIkeAHdfmElDJPSrml7rYd2A0knbZZi79mF1hXi6t7DSrr7hrrLqePCmnx9+QF1uUTQohkYCzwxlk2adLXK5ADPwnIbnA/hzP/09dvI6V0AzYg2g/qAvhRXRfAQiFESjPXdKEutHZfGFr3lfxzIUSvlt553VfpK9Bahw359DU7R13gg9esrntiG1AIfCWlPOvr1YLvyQupC3zznvwn8CvAe5bHm/T1CuTAD2RLgDQpZV/gK05+giuN24I2P0g/4D/AopbcuRAiBPgI+JmUsqIl930u56nLJ6+ZlNIjpewPJAODhRC9W2K/53MBdbX4e1IIcTNQKKXc3Nz7OiGQAz8XaPgpnFz3s0a3EUIYgHCgxNd1SSlLpJSOurtvAAObuaYLdSGvaYuTUlac+EoupVwGGIUQMS2xbyGEES1U50kpP25kE5+8Zuery5evWd0+y4HVwOjTHvLFe/K8dfnoPTkMGCeEOILW9TtKCPHuads06esVyIG/EegihOgghDChHdBYfNo2i4F7625PBFbJuqMfvqzrtD7ecWh9sP5gMXBP3ciTIYBNSpnn66KEEPEn+i2FEIPR/t82e0jU7fNNYLeU8uWzbNbir9mF1OWL10wIESuEiKi7HQxcD+w5bbMWf09eSF2+eE9KKX8tpUyWUqah5cQqKeVdp23WpK+X4VJ/0deklG4hxI+BL9FGxrwlpcwSQvwB2CSlXIz2ppgrhDiAdlBwsp/U9VMhxDjAXVfXtOauC0AI8T7a6I0YIUQO8Hu0A1hIKWcCy9BGnRwAqoH7/KSuicAjQgg3UANMboEPbtBaYHcDO+v6fwF+A6Q2qM0Xr9mF1OWL1ywBmCOE0KN9wHwgpVzq6/fkBdblk/dkY5rz9VJTKyiKorQRgdyloyiKolwEFfiKoihthAp8RVGUNkIFvqIoShuhAl9RFKWNUIGvKIrSRqjAVxRFaSP+H0WXG2H4f4H9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "for site, dynamics in results_closed['dynamics'].items():\n", " plt.plot(\n", " *dynamics.expectations(sz, real=True),\n", " color=f\"C{site}\", linestyle=\"solid\",\n", " label=f\"$\\\\langle s_{site}^z \\\\rangle$\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yay! We can observe how the excitation travels along the chain." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Open Heisenberg spin chain\n", "\n", "Now, let's add a strongly coupled and structured environment to the chain. We will couple a ohmic bosonic environment to the first spin through the environment Hamiltonian\n", "$$ H^E = \\sum_{k}^{\\infty} s_1^y \\left( g_k b_k^{\\dagger} + h.c \\right) + \\omega_k b_k^{\\dagger} b_k \\,\\mathrm{,}$$\n", "where $b_k^{(\\dagger)}$ are the bosonic lowering (raising) operators, and $s_1^y$ is the $y$ spin operator of the first site.\n", "The coupling constants $g_k$ and frequencies $\\omega_k$ are determined by the spectral density\n", "$$ J(\\omega) = \\sum_{k}^{\\infty} |g_k|^2 \\delta(\\omega - \\omega_k) = 2 \\alpha \\, \\omega \\, e^{-\\frac{\\omega}{\\omega_c}} \\,\\mathrm{.}$$\n", "\n", "We choose the values $\\alpha=0.32$ and $\\omega_c=4.0$, and specify that the bosonic environment is initially at temperature $T=1.6$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "alpha = 0.32\n", "omega_cutoff = 4.0\n", "temperature = 1.6\n", "\n", "correlations = oqupy.PowerLawSD(\n", " alpha=alpha,\n", " zeta=1,\n", " cutoff=omega_cutoff,\n", " cutoff_type='exponential',\n", " temperature=temperature)\n", "bath = oqupy.Bath(sy, correlations)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the process tensor approach to TEBD, we first need to compute the process tensors of the environments we wish to add to the TEBD evolution. For this we choose suitable parameters and carry out the computation (see Tutorial 02 - Time dependence and PT-TEMPO)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "tempo_parameters = oqupy.TempoParameters(\n", " dt=pt_tebd_params.dt,\n", " dkmax=40,\n", " epsrel=1.0e-5)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Process tensor (PT) computation:\n", "--> PT-TEMPO computation:\n", "100.0% 20 of 20 [########################################] 00:00:00\n", "Elapsed time: 0.2s\n" ] } ], "source": [ "print(\"Process tensor (PT) computation:\")\n", "pt = oqupy.pt_tempo_compute(\n", " bath=bath,\n", " start_time=0.0,\n", " end_time=num_steps * pt_tebd_params.dt,\n", " parameters=tempo_parameters,\n", " progress_type='bar')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To see the effect of the environment clearly we start in a fully mixed chain state." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "initial_augmented_mps_B = oqupy.AugmentedMPS([mixed_dm, mixed_dm, mixed_dm, mixed_dm, mixed_dm])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we prepare a `PtTebd` object to collect all necessary information. We can reuse the before created `SystemChain` and `PtTebdParameters` objects. This time we set the first item in the list of process tensors to the process tensor we just computed. This attaches the above specified environment to the first chain site. Because we don't want to couple the other sites to any environment we keep them free by setting them to `None`." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "pt_tebd_open = oqupy.PtTebd(\n", " initial_augmented_mps=initial_augmented_mps_B,\n", " system_chain=system_chain,\n", " process_tensors=[pt, None, None, None, None],\n", " parameters=pt_tebd_params,\n", " dynamics_sites=[0, 1, 2, 3, 4],\n", " chain_control=None)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PT-TEBD computation (open spin chain):\n", "--> PT-TEBD computation:\n", "100.0% 20 of 20 [########################################] 00:00:01\n", "Elapsed time: 1.5s\n" ] } ], "source": [ "print(\"PT-TEBD computation (open spin chain):\")\n", "results_open = pt_tebd_open.compute(num_steps, progress_type=\"bar\")" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABAxElEQVR4nO3deXxU1fn48c+ZLRtJyE7ICsquLCEKCMgOblWL1tpaS1utrdra1vZXbWurrVq1dbet1Lp+1VatUhU3BAQEBZQdZN8CCUtCAglJyKzn98e9WQkkZJk7mTzv1+u+7nZm7jMXJs+ce+49R2mtEUII0X3ZrA5ACCGEtSQRCCFENyeJQAghujlJBEII0c1JIhBCiG7OYXUAbZGcnKxzc3OtDkMIIbqU1atXH9FapzTd3iUTQW5uLqtWrbI6DCGE6FKUUgXNbZdLQ0II0c1JIhBCiG5OEoEQQnRzkgiEEKKbk0QghBDdXIckAqXURUqpbUqpnUqpO5vZH6GUet3cv1Ipldtg32/M7duUUjM6Ih4hhBCt1+5EoJSyA38HLgYGA99SSg1uUuwG4KjW+mzgMeAh87WDgWuBIcBFwD/M9xNCCBEkHfEcwfnATq31bgCl1GvAFcDmBmWuAO4xl98E/qaUUub217TWbmCPUmqn+X7LOyCuk7zwk9+ia9qeZ2w2hdNuA5QxKcxlQKkmyzQuoxQoR4N9p9NMGVW/VZ20v8EW1XS/Mg6NQill7DPnytxH3XLjufFRbNiVDaUUNmXDpuyNj35SqI03KNVkW+3nb3TuGi6rBkVVg/JNz6W5n+bLKJsCZTO31a/XbbedYrtSKLsNbDaUzYay28Fuzm3m9kYxNv381O9v+G/RKMyG5+Pk7cpWe+5rz0HT5dpy5vub24zwjc+kFI3mtpP2N1m3gd1unAebXWGzqUbLNrtC1W5v1f9h0ZV0RCLIAPY3WC8ERp2qjNbap5QqB5LM7SuavDajuYMopW4CbgLIzs5uU6D+ysG4I3q16bV1Au17ebeiQqkJSgN+q4MIC0pRlxRsNoXNYcNuN+cOG3aHwmY35naHrW5/w2Wb04bdLONw2bE7bTjNucNlw+G0m3MbDpe9yTY7dpcNh8NWlzRF+3SZJ4u11s8AzwDk5+e3aTSdG1/8TpuPv/dIFZMeWcRtE/rwi8m54PeA32fMA17we81t5nKgybqnCiqKjKm8CCoKjXn1kcafEyA6GeJ7Q2wmxGeg47PQmSPRvfPQdicaTUAH0Gi01nXzAAFqBxqq3R/QAQI6gD/gx6d9+AN+/NqPL+Brdt5wvy/gwxPw4Pa5qfHX4PbV1C/73dR4a8ztxjaP3123XOOr4YS3mipPJTX+GgCUNn+w187NbZH2SOKcscS5YunhiiXOGUusswdJUUmkRCSRHJFEclQSSa5EkqIScSknaA2BADqgjTcMBEBrtDnH7zeW/X60P4D2+9D+ANTOA/4m8wDa50cH/Gi/H+3xor0+tM+L9nrBZywHvH7wetE+HwGv19zvM/Z7vQQ8HgIeL7rGTcDjRrs9xrLXQ6DGjfZ4CLg9EPBTV7MEdF2tSKFVw+3mssuFLSoaFR2NijImW1SUsR5pzG0xMagesajoHthieqB69EBFx5hTNNgdEICA1uiARmuMeUATaDAP+M1lvzkFGqw3XPYH6rYFfAH8Po3fHyDg0/h9Afw+Y7/nhK9u2e8N1JfxB/B7A/i8AfM//plzRthxRtpxRTpwRthxRdpxNll21e6PtJvbjeWIaAeRMU4ioh04nN37inRHJIIiIKvBeqa5rbkyhUopBxAPlLbytSEhNzmGqYN68fKXB7h5ykAio2I65o29J6DiAJQXQkURqmGSOLoHCpah3BVGWUcUZI+C3PHQ50LoPQIczo6JoxO5/W6Oe45T4akw5u6K+uWmc3cFhZ4yyo/voaS4BG/Ae9L7xUfEkxKVQmp0av08OoXUqFRSetRvt9tC98utvV4Cbg/aXYN2u40kUXOCQHU1gaoqYzKX/VVV6Opq/E22B44fJXDIXK+sxFdZaSTBU1BRUdjj4rDHxxvznvHY4+KN9cQEHEnJOJISsScl40hOwpGYiHK5Ov9caCNx+DxGYvB6/EaC8ATwefz4vAF8Xn/jdY+x7q3x43H7jHmNH6/bx/GyGrw1PrxuY5vf23I13uG0ERHjJDLGQUS0sz5JxDgbJYzIaCcRMQ6iYl1Ex7qwO0Op1tt2HZEIvgT6KaX6YPwRvxb4dpMy7wKzMK79Xw18orXWSql3gX8rpR4FegP9gC86IKZOceO4PszffJj/rS3iW+e37fLUSZxRkHSWMZ1KVSns+xz2LIW9S+GTe43trh6QPdpMDOMhfTiE4B+/CHsEEVERJEcln9HrtNaUu8spPlFMSXUJxdXFlJww59UllJwoYdexXRw5cQS/bnzZx2lzkhWbRXZcNtmx2eTE5ZAdl01ObA5pMWnYLL5spZxO7E4n9OigHxSA9vsJHD+Ov6ICf3k5/vIK/OXHCDRaL8dfUU7gWDneffupKd+Ev7wcXVPT7Hva4uNxJCXhSErCbs4dyeZycrKxnpKCIyUF5WzbjxKllHHZp5N+lQf8ATNJ+PHU1CYNH+7q2slLTZUPd5WXmiov7mofFUdOUFxwHHeV16ixnEJEtJkU4owpKs5IEI3W40I/aaiOGLNYKXUJ8DhgB57XWt+vlPoTsEpr/a5SKhJ4GRgBlAHXNmhc/h3wA8AH/Fxr/WFLx8vPz9dWdDqnteZrf1tGjTfA/F9caF2jWWUJFCyrTwxHthvbI+Ig54L6xJB2rtHAGeYCOkBZTVldcjhcfZj9x/ezr2IfBRUF7D++H7ffXVc+wh5hJImGCSIuh+zYbFKjU7tlY2iguhpfaSm+I0fwl5biO1KKr6wU/5FSY3vpkbrlwPHjJ7+BzWYkhl69cKalGfNeaTjSzHmvdJypKUGpYXQ0n9ePu9pnJIkqY37iuIcTxz1UV3iprvBQXeHmxHFj2XPC1+z71CaNmJ4ueiREEpsYSY+ECHokRhKbEEmPxAhckZ17tV4ptVprnX/S9q44eL1ViQDgf2sL+cXr63nx++cxcUCqJTGc5Pgh2LsM9nxqJIay3cb2yJ6QOw4GfQ0GXQ6uaEvDtEpAByiuLqagooCCigIjQRw35vuP7290+SnSHknfnn0ZlDiIAYkDGJQ4iP4J/Yl2ds9z15yAx1OfLEqP4CsuxnfoMN7Dh+rnBw8RqKo66bX25OT6RJGWhrN3Os7MTJwZmTgzM7D37NnlE7HP66e6wsOJCi/Vxz1Ul7sbJA03Vcc8VB6toeqY+6QreRHRjrrkYCSLiEbzmIQI7Pa2/7iTRNBBPL4A4//yCf3TYnn5hqY3R4WI8kIzMSyF3YuNNoeIeDj3asj7LvQebnWEIcMf8HOo+lB9gqgoYMexHWwt20q5uxwwbtHMicthQOIABiYOrJvO9FJXd+OvrMR36BDeQ4fxHTpozA83Xm9au7DFxODMysKZmYErI9NIEpkZuLKycGZkYIuKsujTdLyAP0BVuYfKshqOH62hssxtLrupPFrD8bIa3FVNahcKbnz0QiKi2lZzkETQgf6+aCd/nbeNeT+/kAG9Yi2Lo1UCASj4DNb8H2x5F3w10OtcyJtlJIaoBKsjDElaaw5XH2Zr2Va2lG1hW9k2tpZtpaiy/l6G5KjkulpD7TwrNsvy9oeuxF9ZibewEG9hIZ79hfXLRYV4C4tOaruwJyfjysgwkkVWJq7sHFw5Obhyc7AnJHT52kRTXrefSjNJHDdrEfmX5Lb5c0oi6EDHqj2MfmAhVwzL4KGrh1oWxxk7cRQ2vglrXoJDG8ERaVwyyrsecsZ1i/aE9qrwVNQlhdpp97Hd+LTxyy3WGcuw1GGMSB3BiNQRnJN8DlGO8PkVG0xaa/ylpfVJoqgQT6GRILyFhXgPHgR//U0Ctrg4IynUTrm5uHKNZXtcnIWfJHRIIuhgd729kTdWFfL5nZNJ7hFhaSxtcmAdrH0ZNvwX3OWQ0MdICMO+DXHpVkfXpXj8HnYe28nWsq1sPLKRdcXr2HlsJwAOm4PBiYPrEsPw1OEkRSVZHHF40F4v3qIi3Hv34i0owFNQgGdvAZ69e40k0eBvmz0hwUgMZu3BlZODq29fXDk52CIjLfwUwSWJoIPtKqlkyiNL+PnUfvx8an9LY2kXTzVsmWtcOipYBsoO/aYbSaHfdLCH/nMKoajcXc664nWsLV7L2uK1bDqyCU/AA0BuXC7DU4eTl5rHiNQR5MTlhN0lDasF3G68+/fj2bu3PkEUGEnCV1xcX1ApnBkZuPr2IaJPX1x9+xLRtw+uvn2xJyaG3b+LJIJOcMOLX7Ju/zE+u3MykeHwZGLpLqOWsO7fUHkYeqQZbQmjb4boRKuj69I8fg+bSzeztngta4rXsK54HcfcxwBIjExkeMpw8tLyGJ0+mv4J/cPuD1AoCVRXG0lhzx7cu/fg2b0b9549ePbsadQmYYuPJ6JPn0bJwdWnD66sLJSjy3TK0Igkgk7w+c4jfPvZlfzlqqFcc15Wyy/oKvw+2PGx0Zaw/SPjwbXzboQxP4EeKVZHFxa01uyp2MPaw2vrag37ju8DICkyiTG9xxhT+hhSouWcB4MOBPAdPGgkhz27ce/ejWf3Htx7duMvadAVjNOJKzubiL59cZ19FhF9zyLiLCNJhPpdTZIIOoHWmoufWIrW8NHPx4fnr7jDm2Hpw7BpjtG4nP8DGHsbxLaz8z5xksNVh1l+cDnLDyxnxcEVlNWUAXB2z7O5oPcFjOk9hpFpI6Xx2QL+iooGNYhduHftxrNrF579+40788C4zNS7d+PkYM7t8fHWfgCTJIJO8ubqQn713/W8fMP5jO8Xxr/cjuyApY/AhjfA5jCeRxj3c4jPtDqysBTQAbYf3c7nBz5n+YHlrDm8Bk/Ag9PmJC81r67GMDBxoNyuaqGAx2O0Q+zejXvXLjy7zCSxZw/a46krZ09JPik5uPr2xZEa3CfZJRF0ErfPz9gHF3FORhwvfv98q8PpfGW7YemjsP4/gIIR18G4X0BCrtWRhbUTvhOsPbzWSAwHl7P9qNGtSEJEAqPTRzOm9xjGZ46Xh9xChPb7jTuadu0yk8Ru3Lt24tm1m0BlZV05W0yM2QZhNlSbicKVldnmvptORxJBJ3pq4Q4emb+dBbdfyNmpIf6AWUc5tg+WPW40Lgf8MOxaGP/L03eeJzrMkRNHWH7AuIy0/OByjpwwrmGfm3wuE7MmMiFzgjQ6hyCtNb6SkvoaxO49uHcbc9/hw/UF69oh+hg1CHPu6tMHezs6KpRE0InKqjyMeWAhM/MyeWDmuVaHE1wVB+CzJ2H1C8b4C+dcBeN/BakDrY6s29Bas/3odpYULmHx/sVsPLIRgPSYdCZmTWRi5kTye+Xjsne9Dt+6E39lpZEgdu/Gs2s37j3G3LNvX6MH5/qvXNHmNgdJBJ3sN3M2MGdNEct/M4XEmG74hTt+GJY/BV8+D95qGHwFXPj/oNc5VkfW7ZRUl/Bp4acsLlzMigMrqPHXEO2IZmzGWCZmTWR8xngSIqVrka5Cezx49u/HvWsX3sIikn7w/Ta/lySCTrbj8HGmPfYpv5zWn59O6Wd1ONapKoUVf4eVz4Cn0rhkNPkuaVS2SI2vhpUHV7K4cDFL9i+h5EQJNmVjWMqwutpCn/g+cgmpm5BEEASznv+CzQcrWHbHJCIcYfCAWXucOArLHoMVs41BbkffbDQqR4bGbXTdUUAH2FK6hcWFi1m8fzFby7YCkB2bzaSsSUzNmcrQlKFyF1IYk0QQBEt3lHD9c1/w8DeGcfVI+QUMGI3Kn9wHG16HqESYcIfxLIKjG14+CzGHqg6xZP8SFhUuYuXBlfgCPlKjU5mSPYVpOdPIS80L6eE+xZmTRBAEWmsuenwpNpvig9vGSXW7oQPrYP7vjcFzEvvClLuNdgQ5RyGhwlPBkv1LWFCwgM8OfIbb7yYxMrGupjCq1yic0u9UlyeJIEhe/3Ifd7y1kX//cBQXnCX3dDeiNexcAPP/AMWbIfM8mH6fMe6yCBnV3mqWFS1jQcEClhQuodpXTawrlomZE5maM5ULel9ApKP79NgZTiQRBEmN18/YBz9heFZPnvveeVaHE5oCflj3KnxyP1QegoGXwdQ/QvLZVkcmmnD73aw4sIL5BfNZtH8RFZ4KohxRjM8Yz7ScaYzPHE+Ms+33tYvgkkQQRI/N384TC3fwyS8n0Delh9XhhC5PFSz/O3z2hDFy2sjvG20I0rFdSPIGvKw6tIoFBQtYuG8hpTWluGwuLsi4gBm5M5iYOZEeLvn/HsokEQRRyXE3Yx/8hGvOy+S+K7vZA2ZtUVkMix+E1S+CMxrG/QxG3wouGTA+VPkDftaXrGd+wXzmF8zncPVhXDYXYzPGMj13uiSFECWJIMh+/eZ65q4/yPLfTKZntNwh0yol22HBPbDtfYjtDVPvgXO/IUNohriADrChZAPz9s6TpBDiJBEE2dZDFVz0+FJ+fdEAbpko177PyN7P4OPfwYG1RoPyRQ9B5kiroxKt0DApfFzwMcXVxZIUQogkAgtc/9xKth8+ztJfT8blkF+1ZyQQMHo4XXAPVBUbYylPvVvGQehCTpcUZuTOYGLWRGloDjJJBBZYtK2Y77/wJY9/czhXjsiwOpyuqabCGAdhxT/A7oLxtxvtB065fbErCegA60vW8/HejyUpWEgSgQUCAc20x5YQ5bIz9yfygFm7lO6Cj39vtB/0zIEZ9xu3nco57XJqk8K8vfOYv3c+xSeMpDAuY5xx+UiSQqeRRGCRf6/cx2//t5HXbxrNqL5JVofT9e1aBB/9Bkq2QJ8L4aIHIW2I1VGJNpKkEFySCCxS4/Uz5oGFnJebyDPfPen8i7bw+4zxDz65D9wVRt9FE38LMZJou7LTJYUZuTOYkDVBkkI7SSKw0MPztvH3xTv5/M7JpMfLwOMdproMFj8AXz4HET2MZHDeDSB94nR5AR1gXfE6Pi74uC4pRNgjjJpCznRJCm0kicBCu0oqmfLIEv54+RBmXZBrdTjh5/BmmPcb2L0YkgfARQ/A2VOsjkp0kNMlhak5U5mQOYFYVzcZIradJBFYbMoji+kVH8mrN0oHa51Ca9j2Icz7LRzdA/0vNhqUZQzlsNJcUnDanIzpPYap2VOZlDWJnpE9rQ4zZEkisNhfPtrKPz/dzZq7phEfLZcuOo3Pbdxq+unDxvLom40hMyPjrI5MdLDa5xQWFCxgwb4FFFUWYVd2zut1HtNypjE5ezLJUdIDcEOdkgiUUonA60AusBe4Rmt9tJlys4C7zNX7tNYvKaWigf8CZwF+YK7W+s7WHLcrJoJ1+49x5d8/49FrhjEzTwat6XTHD8PCP8G6VyAm1XgYbdi3pbuKMKW1ZnPZZiMpFCxgb8VeFIoRqSOYljONqTlT6RUjDyN2ViL4C1CmtX5QKXUnkKC1vqNJmURgFZAPaGA1MBJwA6O01ouUUi5gIfBnrfWHLR23KyaCQEAz5sGFjMhKYPb10l1C0BSthg/vgMIvIX04XPwXyB5ldVSiE2mt2XlsJwsKFjB/33x2HN0BwLnJ59YlhazYLIujtEZnJYJtwESt9UGlVDqwWGs9oEmZb5llfmSu/9Ms958m5Z4ANmmt/9XScbtiIgD4/dubeHN1IWv/MI1IpwwBGDRaw8b/GgPiHD9odGQ39Y8QL097dwd7y/eyYN8C5hfMZ3PpZgAGJg5kctZkJmVPYkDCgG7zsGdnJYJjWuue5rICjtauNyjzKyBSa32fuf574ITW+uEGZXoCa4CpWuvdpzjWTcBNANnZ2SMLCgraHLdVasc0fva7+UwdnGZ1ON2PuxI+exw+exJsdhj3C7jgp+CUW3q7i6LKorrLR+tL1qPR9I7pzcSsiUzKnsTItJE4beHbhtfmRKCUWgA0d3Htd8BLDf/wK6WOaq0Tmrz+tIlAKeUA5gLztNaPt+bDdNUagccXYOR987loSC/++o1hVofTfR3da9QONr8D8dkw/V4ZP7kbOnLiCJ8WfsqifYtYfnA5br+bWFcs4zPGMyl7EuN6jwu7nlJPlQgcLb1Qaz31NG96WCmV3uDSUHEzxYqAiQ3WM4HFDdafAXa0Ngl0ZS6HjSkDU1mw5TA+fwCHXRouLZGQC9f8H+xZCh/dCf+dBTnj4OIHoZcMJNRdJEclM7PfTGb2m0m1t5oVB1ewaP8iluxfwgd7PsBhczCq1ygmZk1kYtbEsG5sbu+lob8CpQ0aixO11r9uUiYRo4E4z9y0BhiptS5TSt0HDAK+obUOtPa4XbVGAPDBxoPc8uoa6XsoVPh9sOYlo7uKmmOQNwsm3wUxctthd1U7+tqi/YtYtH8RBRXGZejBSYOZlDWJSVmT6J/Qv0u2K3RWG0ES8AaQDRRg3D5appTKB36stb7RLPcD4Lfmy+7XWr+glMoE9gNbMe4gAvib1vrZlo7blRNBldvHiHvn851ROfzha4OtDkfUOnEUFj8EXzwDrhgY/0sY9WPp7rqb01qzp3wPn+z/hEX7F7GxZCMaTXpMOuMzxjMuYxyj0kcR7ewaw6rKA2Uh5IYXv2Tb4eMs/fWkLvmrIqyVbDO6u94xz2g/mHo3nHOVtB8IwGhXWLx/MUsLl7Li4AqqfdU4bU7y0/IZlzGO8ZnjyY3LDdnvtSSCEPL6l/u4462NfHDbeAb3lideQ9LuxTDvLji8ETLyYcaf5fkD0YjX72VN8RqWFi5lWdEydpXvAiCzR2ZdUjiv13lEOULnrjRJBCHkSKWb8+5fwG2T+/GLaf2tDkecSsBvDJe58F6oPGTcWTT1Hkjsa3VkIgQVVRaxrHAZy4qWsfLQSk74ThBhjyC/Vz7jM8ZzYcaFZMVZ+yCbJIIQ843Zn1Pp9vPhz8ZbHYpoiacKPn8KPnsC/F4Y9SO48FcQldDya0W35Pa7WX1oNUuLjNrC3oq9AOTE5TC291jG9B5Dflp+0G9PlUQQYp5dupv73t/C0l9PIiuxazQ0dXsVB427i9a9ClE9YcKdMv6BaJX9FftZWrSUpUVLWXVoFTX+GuzKzjnJ5zA6fTSj00czLGUYzk7+vySJIMQUlFYx4a+LuevSQdw4Xi41dCkHN8DHd8GeJZB4Fkz7Ewy8VBqURat4/B7Wl6xn+YHlrDy4kk2lmwjoAFGOKEamjaxLDP0S+mFTHfuskSSCEHTR458SF+XkjR+NsToUcaa0hh0fG3cYHdlmPJA24z7oPcLqyEQXU+Gp4MtDX7LiwApWHlrJnvI9ACRGJjIqfVRdYujdo3e7jyWJIAQ9On87f/tkB1/+bipJPSKsDke0hd8Ha16ERX+G6lI49xqY9FtI7GN1ZKKLOlR1iJUHV7Li4ApWHFzBkRNHAMiOzWZ0+mhuy7uN+Ij4Nr23JIIQtKmonMueWsZfrh7KNfnds1vcsFFTDksfhZWzIeCDvO8aA+LEtf9XnOi+tNbsOrarLilsKd3CvKvn4bC12DtQsyQRhCCtNeMeWsSg9FienXWe1eGIjlBxEJY+DKtfBJsDzv8hjP0FxEh3IqL9AjrQrnaDUyUC6fXMQkoppg9J49MdR6hy+6wOR3SEuHS49BH4ySoY8nX4/G/wxDBY/CDUVFgdnejiOrrxuO59O+VdRatNH9wLjy/A0h0lVociOlJiH/j6bLhlOZw1ERY/YCSEz58C7wmroxOiEUkEFjsvN4GEaCfzvjpsdSiiM6QOgm++Aj/8BHoPN247fXIErHreeDhNiBAgicBiDruNKYPSWLjlMF5/q3viFl1Nxki4/n8w6z3omQ3v/QL+dh5seMPoykIIC0kiCAHTB6dRUePjiz1lVociOluf8fCDefDtN8DVA+b8EGaPg63vG88mCGEBSQQhYHy/FCKdNuZ9dcjqUEQwKAX9Z8CPPoWrnwefG177Njw7BbZ+AAGpGYrgkkQQAqJcdi7sl8LHXx2mK97OK9rIZjPGOrj1C/jak1BVAq99C2aPhQ3/NR5WEyIIJBGEiBlDenGoooaNReVWhyKCze6AkbPgp2vg68+ADsCcG+FvI41GZW+N1RGKMCeJIERMHpiK3ab4WO4e6r7sThj2Tbh5OVz7b4hOMhqVnxgKnz0J7uNWRyjClCSCEJEQ4+L83ERpJxDGJaOBl8KNC+G770DKQJj/e3jsHLNPI7mpQHQsSQQhZMaQNHYUV7K7pNLqUEQoUAr6ToRZ78KNn0DuOFjykJEQPvotVBywOkIRJiQRhJBpQ3oBMH+zXB4STWSOhGtfhVtWwKDLjM7tHh8K7/4USndZHZ3o4iQRhJCMnlGckxEnl4fEqaUOgpnPwG1rjB5O178Of8uH/34fClfJswiiTSQRhJgZg3uxdv8xiivkThFxGgm5cNmj8PONcMFPYcd84zmEZybC2lekPyNxRiQRhJjpQ3qhNSzYUmx1KKIriE0zhsr85Ra45GEjAbxzKzw6yBg97eheqyMUXYAkghDTP60HOUnRcnlInJmIWGPsg1tXwqy5kDselv8dnhgO//4m7FggTyyLU2rbMDei0yilmDGkFy98tofjNV5iI51WhyS6EqWgz4XGVF4Eq1+A1S/B9qsgsS/k3wAjroOoBKsjFSFEagQhaPrgNLx+zeJtMkaBaIf4DJh8F/ziK7jqOYhJhY9/B48MMu42OrjB6ghFiJBEEIJGZCeQ3MPFx3IbqegIDhecezXcMA9+tBSGfsPoy+if4+G56bDxTaPjO9FtSSIIQXabYuqgNBZtLcbtk77qRQdKHwqXP2U0Ls/4M1QWw1s3wMP9YO7PoOBzaUvohiQRhKgZQ3pR6faxfFep1aGIcBSVAGNuNTq6+84c6H+RMUjOCxfDk8Ng4b1Qss3qKEWQSGNxiBpzVhIxLjsfbz7MxAGpVocjwpXNBmdPMSZ3pTFAzobXYdmjsPRhSB8OQ79pdJcdm2Z1tKKTSI0gREU67UwckMr8zYcJBORpUREEET2M3k+vnwO3b4UZDxjb5/0GHh0IL880nmR2S19Y4aZdiUAplaiUmq+U2mHOm70nTSk1yyyzQyk1q5n97yqlNrUnlnA0fUgaJcfdrN1/zOpQRHcTmwZjboEfLTEGzhl3OxzZAf+7yWhPeOuHxrMJMnhOWGhvjeBOYKHWuh+w0FxvRCmVCNwNjALOB+5umDCUUjMB+YnRjEkDU3HaFR9vlofLhIVSBsCU38PP1sP3PzIuFe34GF69yniC+f1fwc4FcudRF9beRHAF8JK5/BJwZTNlZgDztdZlWuujwHzgIgClVA/gduC+dsYRluIinYzumyRDWIrQYLNBzhj42uPwq+3wzVcge5TRt9ErV8Ff+sLr18O6f0PVEaujFWegvY3FaVrrg+byIaC51qQMYH+D9UJzG8C9wCNAdTvjCFszhvTirrc3saO4kv5psVaHI4TBEQGDvmZM3hOw51PY9iFs/wi2vAsoyDrfuBtpwMXG4DpKWR21OIUWawRKqQVKqU3NTFc0LKeNn6yt/tmqlBoOnKW1/l8ry9+klFqllFpVUtJ9nridPjgNm4K562UQEhGinFHQf4ZRU7h9C9y0GCbcAb4aWPhH+MdoeHI4fHgn7F4Mfq+18YqTqPZcclBKbQMmaq0PKqXSgcVa6wFNynzLLPMjc/2fwGKgJ/B7wINRM0kFPtdaT2zpuPn5+XrVqlVtjrur+e7zX7CruJKlv56EzSa/qkQXUl5k1BK2fwS7l4DfDRFxxu2q/S+GftMgOtHqKLsNpdRqrXX+SdvbmQj+CpRqrR9USt0JJGqtf92kTCKwGsgzN60BRmqtyxqUyQXe01qf05rjdrdE8M66In722jr+88PRjDkryepwhGgbT5VRI9j2AWyfB1UlgIK0c4xhOHPHQc4Fkhg60akSQXvbCB4E3lBK3QAUANeYB8sHfqy1vlFrXaaUuhf40nzNnxomAdGy6YN70SPCwVtrCiURiK7LFQMDLzWmQAAOrIFdn8DepUYvqSufNso1SgxjJTEEQbtqBFbpbjUCgDve3MB7Gw7w5V1TiXbJA+EizPjcULQG9i6DgmWwbyX4zFHWUoc0Tgwx8mOorTqrRiCCZGZeBq+v2s+8rw7x9RGZVocjRMdyRBi3puaMAf4f+DxGjWHvUtj7Gax9Gb74p1E2dXB9UsgYCfGZckdSO0mNoIsIBDQTHl5EblIML98wyupwhAgunwcOrDVqC3vNGoO3ytgXnQy9RzSe4tKtjTdESY2gi7PZFDNHZPLkJzs4WH6C9Pgoq0MSIngcLuPhtexRMP6Xxi2ohzYYyeHAWihaC7sWgja70I5NPzk5xCRb+xlCmCSCLmRmXgZPLNzB22sPcPPEs6wORwjr2J3GZaGMkfXbPNVwaGN9cjiw1njIrfbxpvgs6D28PjGknQMxKXJZCUkEXUpOUgzn5Sbw1ppCfjyhL0r+AwtRzxVdX2uo5T5uDMl5YE19ctgyt35/VCKkDjKefG4472a1B0kEXczMvEx+M2cjG4vKGZrZ0+pwhAhtEbGQO9aYap04aiSH4i1QsgWKt8KmN6GmvL5MdLKZFAY2SBKDwvaOJUkEXcylQ9O5+92veGt1oSQCIdoiKgH6TjCmWlrD8UP1iaF2vuENcFfUl4tJMRJDcn9I7Fs/JeSCMzLoH6WjSCLoYuIinUwfnMa76w/wu0sH43LI2EJCtJtSxp1Gcelw1uT67VpDxYGTE8RXc4yaRf0bQFwGJPapTw5JZ9UnCVdMsD/RGZFE0AVdNTKT9zYcZNG2YmYM6WV1OEKEL6UgPsOYzp7aeF91GZTtgbLdjaet70N1k264Y9PNBNHHSAzx2dAzG3pmGfts9qB9pOZIIuiCxp+dTEpsBG+tLpREIIRVohONKXPkyftqyptPEts/hqrixmVtDqM20dNMDvFZ9UkiPst4YM7u7NSPEjaJwOv1UlhYSE1NjdWhdKjIyEgyMzNxOuv/IzjsNq4c3psXP99LWZWHxBiXhREKIU4SGW/eqjr85H2eaigvhPJ9cGwfHNsP5fuN5V2L4PhBGvXor2xGraE2QVz2qNEI3oHCJhEUFhYSGxtLbm5u2NxWqbWmtLSUwsJC+vTp02jfzLxM/rV0D3PXH2DWBbnWBCiEOHOuaEjpb0zN8XmgotBIEMf2mUnCXC5aBc7oDg8pbBJBTU1NWCUBAKUUSUlJNDcQz6D0OAanxzFnTaEkAiHCicNV3+AcJGF1y0k4JYFap/tMM/MyWF9Yzs7i40GMSAgRbsIqEXQ3VwzPwG5TvLWmyOpQhBBdmCSCLiwlNoIJ/VP435oi/IGu14usECI0SCLoBJMnT8bn8zW7b9++fdx1111cd911XHfddRw9epSvf/3rbT7WVXmZHKqoYfmu0ja/hxCie5NE0MG++uorkpKScDiab4fPzs7mhhtuwG6388wzz5CQkEBZWRmlpW37Qz5lUCpxkcYwlkII0RaSCDrYO++8w5VXXgnASy+9xMiRIxk6dCjjxo0DYO/evdxzzz08/fTTxMQYj51feumlzJ0791RveVqRTjuXDevNR5sOUeluvhYihBCnI4mgg33wwQdceumlHD9+nIceeojly5ezYcMG3nvvPQAuueQSEhMTeeCBBygrKwPgiiuu4O23327zMa/Ky+CE18+HGw92xEcQQnQzYfMcQUN/nPsVmw9UtFzwDAzuHcfdXxty2jLV1dV4PB569uxJdXU1J06c4Je//CWzZs0iP98YHW7z5s0nvW7AgAFs27atzbHlZSeQmxTNnDVFfCM/q83vI4TonqRG0IGio6NRSlFZWUl0dDSbNm1i7Nix3HTTTfzjH/845esKCgpOenL4TCilmJmXyfLdpRQerW7z+wghuqewrBG09Mu9M82YMYOPPvqIYcOG0a9fP6699lo2b95c1wdSVVUVt9xyCy6Xi4kTJ3LdddfxzjvvcMUVV7TruF8fkcGj87fz9toifjK5X0d8FCFENyE1gg5We73//vvvZ8CAAeTl5bFnzx5uueUWAObMmcPVV1/Nv/71L959910A5s6d2+5EkJUYzag+iby1pgit5ZkCIUTrhWWNwEojR45kw4YNrFmzptlbSAsLCzn33HMBsNvtHD16FLfbTa9e7e9O+qqRmfz6zQ2s3X+MvOyEdr+fEKJ7kBpBJ9iwYcMpnyPIzMyksNC45z8QCJCQkMCnn37aIce9+JxeRDptvLVanikQQrSeJIIgmzlzJm+99RY333wzX/va1zr0vWMjnVw0pBdz1x+gxuvv0PcWQoQvuTQUZDExMbzwwgud9v4z8zJ5e90BPtlazCXnpnfacYQQ4UNqBGFm7NnJpMVFyOUhIUSrSSIIM3ab4soRGSzeXsKRSrfV4QghugBJBGHo6rxM/AHNO+sOWB2KEKILkEQQhvqlxTI0M5450iOpEKIVJBGEqZkjMvjqQAVbD3Vsn0tCiPDTrkSglEpUSs1XSu0w580+xaSUmmWW2aGUmtVgu0sp9YxSartSaqtS6qr2xCPqXT48A4dNMUeGsRRCtKC9NYI7gYVa637AQnO9EaVUInA3MAo4H7i7QcL4HVCste4PDAaWtDMeYUqMcTFpYCr/W1uEzx+wOhwhRAhrbyK4AnjJXH4JuLKZMjOA+VrrMq31UWA+cJG57wfAAwBa64DW+kg74xENXJWXSclxN8t2ymkVQpxaexNBmta6djSUQ0BaM2UygP0N1guBDKVUT3P9XqXUGqXUf5VSzb0eAKXUTUqpVUqpVSUlJe0Mu3MFc8zi05k0MIWkGBfPLdvTKe8vhAgPLSYCpdQCpdSmZqZG3WVqo8vLM+n20gFkAp9rrfOA5cDDpyqstX5Ga52vtc5PSUk5g8MEV7DHLD6dCIedH084i6U7jrBytwxuL4RoXouJQGs9VWt9TjPTO8BhpVQ6gDkvbuYtioCGw2ZlmttKgWpgjrn9v0BeOz5LSAj2mMUtuX5MDmlxETz88TbpnloI0az2Xhp6F6i9C2gW8E4zZeYB05VSCWYj8XRgnlmDmAtMNMtNAU4ex7GLsWLM4tOJdNr5yeR+fLn3KEu2h/YlNSGENdrb6dyDwBtKqRuAAuAaAKVUPvBjrfWNWusypdS9wJfma/6ktS4zl+8AXlZKPQ6UAN9vZzyGD++EQxs75K3q9DoXLn7wtEWsGrO4Jd/Mz+KfS3bxyMfbmdA/BaVUpx1LCNH1tKtGoLUu1VpP0Vr3My8hlZnbV2mtb2xQ7nmt9dnm9EKD7QVa6wu11kPN99nXnnisZtWYxS1xOWz8bEo/NhaVM++rw512HCFE1xSe3VC38Mu9M7U0ZvHu3bu5//77KS8v58033wTokDGLW/L1ERk8vWQXj87fxrTBadhtUisQQhiki4kO1tKYxX379uW5555r9JqOGLO4JQ67jdun9Wf74UrmrpfO6IQQ9cKzRmChlsYsbqojxyxuySXnpDMofRePLdjOpUPTcdrld4AQQmoEneJ0YxY31ZFjFrfEZlP8anp/CkqrZeAaIUQdSQRBVlpayo9//GPWrl3LAw88EPTjTx6YyvCsnjy5cAdun4xrLISQS0NBl5SUxOzZsy07vlKK/zdjANc9u5J/r9zH98d23t1KQoiuQWoE3dDYs5MZ0zeJvy/aSbWn+T6RhBDdhySCbupXMwZwpNLDi5/vtToUIYTFJBF0UyNzEpg8MJV/LtlNRY3X6nCEEBaSRNCN3T6tP+UnvDy7VLqpFqI7k0TQjZ2TEc8l5/biuaW7KavyWB2OEMIikgi6udun9eeE18/sJbusDkUIYRFJBN3c2amxXDkig5c+38vhihqrwxFCWEASgeDnU/rjD2j+9slOq0MRQlhAEoEgOymaa87L4rUv97G/rNrqcIQQQSaJoBOEyuD1Z+Knk89GKcWTC3dYHYoQIsgkEXSwUBq8/kykx0dx/egc3lpTyK6SSktjEUIElySCDhZqg9efiZsnnkWk085j87dbHYoQIogkEXSwUBu8/kwk94jgB2P78N6Gg2w+UGF1OEKIIAnL3kcf+uIhtpZt7dD3HJg4kDvOv+O0ZUJ18Poz8cPxfXlp+V4enb+dZ2flWx2OECIIpEbQgUJ18PozER/t5EcX9mXBlsOs3XfU6nCEEEEQljWCln65d6aWBq9/++23ef/996moqOCGG25g+vTpQRm8/kx8f2wfXvhsL498vJ1XbhxldThCiE4mNYIO1tLg9VdeeSX/+te/mD17Nq+//joQnMHrz0RMhIObJ57Fsp1HWL7L2ruZhBCdLyxrBFZq7eD19913H7feemtQB68/E98ZncOzS/fwp/c2879bLiDSabc6JCFEJ5EaQSc43eD1WmvuuOMOLr74YvLy8oI6eP2ZiHTa+fPMc9hysILfztmI1trqkIQQnUQSQZA99dRTLFiwgDfffNPSsYtbY/LANH4xtT9z1hbxkoxkJkTYkktDQXbbbbdx2223WR1Gq/108tlsLCrn3ve3MDA9jtF9k6wOSQjRwaRGIE7LZlM8+s1h5CRFc+urazhw7ITVIQkhOpgkAtGiuEgnz1w/khqvn5tfWU2N1291SEKIDiSJQLTK2amxPHLNcNYXlvOHdzZJ47EQYUQSgWi1i87pxU8nn80bqwp5deU+q8MRQnQQSQTijPx8an8mDUjhj3O/YtXeMqvDEUJ0gHYlAqVUolJqvlJqhzlPOEW5WWaZHUqpWQ22f0sptVEptUEp9ZFSKrk98YjOZ7cpHr92BBk9o7j51TUyzrEQYaC9NYI7gYVa637AQnO9EaVUInA3MAo4H7hbKZWglHIATwCTtNZDgQ3AT9oZjwiC+Cgn/7w+nyq3j5tfWY3HF7A6JCFEO7Q3EVwBvGQuvwRc2UyZGcB8rXWZ1vooMB+4CFDmFKOUUkAccKCd8YggGdArlr9ePYw1+47xx7lfWR2OEKId2psI0rTWB83lQ0BaM2UygP0N1guBDK21F7gZ2IiRAAYDz53qQEqpm5RSq5RSq0pKStoZdufqimMWt8WlQ9P58YSzeHXlPl77QhqPheiqWkwESqkFSqlNzUyNusvUxv2Erb6nUCnlxEgEI4DeGJeGfnOq8lrrZ7TW+Vrr/JSUlNYeJui66pjFbfX/ZgxgfL9k/vDOVzJ+gRBdVIuJQGs9VWt9TjPTO8BhpVQ6gDkvbuYtioCsBuuZ5rbh5vvvMpPIG8AF7fs41uvKYxa3hd2meOpbI0iLj+DmV9ZQctxtdUhCiDPU3ktD7wK1dwHNAt5ppsw8YLrZQJwATDe3FQGDlVK1P++nAVvaGY/luvKYxW3VM9rFP7+Tz7ETHm59dQ1evzQeC9GVtLfTuQeBN5RSNwAFwDUASql84Mda6xu11mVKqXuBL83X/ElrXWaW+yPwqVLKa77+e+2MB4BDf/4z7i0dO2ZxxKCB9Prtb09bJhzGLG6rwb3jeOiqofzstXXc//4W7rl8iNUhCSFaqV01Aq11qdZ6ita6n3kJqczcvkprfWODcs9rrc82pxcabJ+ttR6ktR6qtf6a1rprXig3hcOYxe1xxfAMbhzXhxc/38ubqwutDkcI0Uph2Q11S7/cO1NLYxZv2bKFJ554giNHjjBlyhRuvvnmkBuzuD3uvHggmw9W8Nv/bWRAWiznZsZbHZIQogXSxUQHa2nM4kGDBjF79mzeeOMNPvvsMyD0xixuD4fdxlPfGkFKjwhu/L8vWbf/mNUhCSFaEJY1Aiu1Zszid999l6effprrr78+ZMcsbo+kHhE89718bnhxFd+Y/Tm/v2ww14/OwXhuUAgRaqRG0AlON2YxwOWXX86HH37Iq6++GrJjFrfXwF5xvH/bOMb3S+EP73zFz15bR5W7+YfshBDWkhpBkC1evJg5c+bgdru55JJLrA6nU/WMdvHsd/N5eskuHvl4G5sPVjD7O3mcnRprdWhCiAYkEQTZxIkTmThxotVhBI3Nprh10tkMz+rJbf9Zy+V/+4wHrxrK5cN6Wx2aEMIkl4ZEUIw9O5n3bxvP4PQ4bvvPWu5+Z5P0WipEiJBEIIKmV3wk/7lpNDeO68NLywu45p/LKTp2wuqwhOj2JBGIoHLabdx12WCevi6PncWVXPbkUpZsD+3eZIUId5IIhCUuPjedd38ylrS4SL73whc8Nn87/kCrO68VQnQgSQTCMn1TevC/W8by9REZPLFwB9974QvKqjxWhyVEtyOJQFgqymXnkW8M44GZ57JyTxmXPblUxjUQIsgkEQjLKaX41vnZzLn5Aux2xTX/XM4Ln+2RS0VCBIkkAhEyzsmI572fjGdC/xT+OHczkx5ezPPL9nC8xmt1aEKENUkEIqTERzt55vp8nr4uj9TYCP703mYueOAT7n1vM/vLqq0OT4iwJImgE3SXwes7i82muPjcdN68+QLeuXUskwam8tLne5nw10Xc/MpqvtxbhjG6qRCiI0gi6GDdbfD6zjYsqydPfmsES++YxE0XnsXnu0r5xuzlXPH3z3hnXZEMiylEB5BE0MG62+D1wZIeH8WdFw9k+W8mc++V51BZ4+Nnr61j/EOL+MfinRyrlttOhWgr1RWr2Pn5+XrVqlWNtm3ZsoVBgwYBsPSN7RzZX9mhx0zO6sH4a/q3WG7cuHG899572O12Ro0axbp163C5XBw7doyePXsyePBgZsyYQUxMDLfffjuJiYls27aNO+6445QD2Df8bMIQCGgWby/muWV7+GxnKZFOG1flZfKDcX04K6WH1eEJEZKUUqu11vlNt0vvox2oOw9eH2w2m2LywDQmD0xjy8EKnl+2h/+uKuTVlfuYNCCF6UN6kZ+TwFkpPbDZZEAcIU4nLBNBa365d4aGg9f36NGDTZs2MXfuXG666SZuvPHGuuEqmwqXweutMig9jr9+Yxi/vmggr64s4N8r97Fom9F/UXyUk5E5CYzMSSA/J4FhWT2JdNotjliI0BKWicBKLQ1eD1BVVcWECRO45557uOyyy8Jq8HorpcRG8POp/fnZlH7sOVLFqoKjrN57lFUFZXyytRgAh00xJCOefDMxjMxNIDU20uLIhbCWJIIOdsUVV/DYY4/x3nvvsXz5cmJiYhgyZAj/+te/6so89NBDXHPNNXXrc+fO5eWXX7Yi3LCklKJvSg/6pvTgmvwsAI5WeVhdcJTV+4zk8MqKAp5btgeA7MTouqQwMieBs1N64LDLfRSi+5BE0MFaGrx+/vz5DB48uK6GEI6D14eihBgXUwenMXVwGgAeX4BNB8rragyf7ihhztoiAJSC5B4RpMVFkBYbSWpcpLFszlNjI0mLiyQpxiXtDyIsSCLoBBs2bDjlvsWLF1NVVcXmzZuJiorikksuCcvB60Ody2EjLzuBvOwEfkhftNbsK6tm1d6jFJRVU1xRw+GKGg6W17C+8BhHKk++PdVhU6TERhiJIjaCXvGR9Ix2EeGwGZPTXr/ssBPhbLDssBHprF+OcNhxOhRagwa01uYc0GCs0ex+bRTArzU+vyagNf6AOTW3reG+gCZgrhtljPcJNCgTaDgPmK8xywa0NmIxbz6svQexfl03WW9SAKMGZ7cpbMpYtilj2W5T5jrGNluDZXPusCscNhsOm8Jht5nr5jZz2Wm3YbcpnHaF3SzrtNtw2hVOhw2X3Zi6c1KXRBBk999/PwAvvvgiycnJ2GxyCSIUKKXISYohJymm2f0eX4AjlW4OVdSYScLNYXNefLyGvaVVrNxTRvkJ6ReptZRqlA8sV5ssXHYbLocNZ4O5027DZVeNtteWi3DYzbkxNd5nw9Vgf+P5yT8UXHab+YPBjj2IiUkSgUW+973vWR2COAMuh43ePaPo3TPqtOUCAY3HH8DtC+D2+XF7Gyz7Ari9AWO/11xvUM7jD6Aw/kAqFMr8O6CUqtsOmMvKLGfsUBh/yOw2hd38pWxTxi9im7nNbjf31W5rUN74RV67jbrl+m0Nlhu8l80sWxtXbbyN1xtvb0prTUAbtYuAWbuoraEEtLHfWK8vW1tD8QU0/kAAr9+o+fgCAXyBBst+ba6b25uU9fiM8+71aTx+P16/xuML4PUH6uZev8Zdt2xsr3T7jNfW/RsG8Jj/xh5/oEMSnMOmGiSY+mQy96fjOvzON0kEQnQgm00RabObX1Sn1eF0CUop7ArshMelGW1ecvPUJYjahOFvkDSM9caJpOkPhvofEg33OTqhpiCJQAghOpBSxiUmp91GTITV0bSOXKAWQohuLqwSQVfsN6kl4fiZhBChJWwSQWRkJKWlpWH1h1NrTWlpKZGR8uSrEKLztKuNQCmVCLwO5AJ7gWu01ieNPK6U+ggYDSzTWl/WYHsf4DUgCVgNXK+1blN/wpmZmRQWFlJSUtKWl4esyMhIMjMzrQ5DCBHG2ttYfCewUGv9oFLqTnP9jmbK/RWIBn7UZPtDwGNa69eUUrOBG4Cn2xKI0+mUjtuEEKIN2ntp6ArgJXP5JeDK5gpprRcCxxtuU8ZNxZOBN1t6vRBCiM7T3kSQprU+aC4fAtLO4LVJwDGtde3gvoVAxqkKK6VuUkqtUkqtCrfLP0IIYaUWLw0ppRYAzfWI9ruGK1prrZTqtJZarfUzwDNgjFDWWccRQojupsVEoLWeeqp9SqnDSql0rfVBpVQ6UHwGxy4FeiqlHGatIBMoas0LV69efUQpVXAGx2ooGTjSxtd2JonrzEhcZ0biOjPhGldOcxvb21j8LjALeNCcv9PaF5o1iEXA1Rh3DrX69VrrlDMP1aCUWtXcmJ1Wk7jOjMR1ZiSuM9Pd4mpvG8GDwDSl1A5gqrmOUipfKfVsbSGl1FLgv8AUpVShUmqGuesO4Hal1E6MNoPn2hmPEEKIM9SuGoHWuhSY0sz2VcCNDdbHn+L1u4Hz2xODEEKI9gmbJ4vPwDNWB3AKEteZkbjOjMR1ZrpVXCqcumQQQghx5rpjjUAIIUQDkgiEEKKbC9tEoJS6SCm1TSm10+wHqen+CKXU6+b+lUqp3BCJ63tKqRKl1DpzurG59+ngmJ5XShUrpTadYr9SSj1pxrxBKZXX2TG1Mq6JSqnyBufqD0GKK0sptUgptVkp9ZVS6mfNlAn6OWtlXEE/Z0qpSKXUF0qp9WZcf2ymTNC/j62MK+jfxwbHtiul1iql3mtmX8eeL6112E2AHdgF9AVcwHpgcJMytwCzzeVrgddDJK7vAX8L8vm6EMgDNp1i/yXAhxjD0I4GVoZIXBOB9yz4/5UO5JnLscD2Zv4dg37OWhlX0M+ZeQ56mMtOYCUwukkZK76PrYkr6N/HBse+Hfh3c/9eHX2+wrVGcD6wU2u9WxvdWr+G0UFeQw07zHsT4xmHzh40tTVxBZ3W+lOg7DRFrgD+TxtWYDwRnh4CcVlCa31Qa73GXD4ObOHkfrKCfs5aGVfQmeeg0lx1mlPTu1SC/n1sZVyWUEplApcCz56iSIeer3BNBBnA/gbrzXVoV1dGG11clGM81GZ1XABXmZcT3lRKZXVyTK3R2ritMMas2n+olBoS7IObVfIRGL8mG7L0nJ0mLrDgnJmXOdZhdEMzX2t9yvMVxO9ja+ICa76PjwO/BgKn2N+h5ytcE0FXNhfI1VoPBeZTn/XFydYAOVrrYcBTwNvBPLhSqgfwFvBzrXVFMI99Oi3EZck501r7tdbDMfoUO18pdU4wjtuSVsQV9O+jUuoyoFhrvbqzj1UrXBNBEdAwczfXoV1dGaWUA4jH6AjP0ri01qVaa7e5+iwwspNjao3WnM+g01pX1FbttdYfAE6lVHIwjq2UcmL8sX1Vaz2nmSKWnLOW4rLynJnHPAYsAi5qssuK72OLcVn0fRwLXK6U2otx+XiyUuqVJmU69HyFayL4EuinlOqjlHJhNKa826RMbYd5YHR894k2W16sjKvJdeTLMa7zWu1d4LvmnTCjgXJdPw6FZZRSvWqviyqlzsf4/9zpfzzMYz4HbNFaP3qKYkE/Z62Jy4pzppRKUUr1NJejgGnA1ibFgv59bE1cVnwftda/0Vpnaq1zMf5GfKK1/k6TYh16vtrb+2hI0lr7lFI/AeZh3KnzvNb6K6XUn4BVWut3Mb4wLyujw7syjBMeCnHdppS6HPCZcX2vs+NSSv0H426SZKVUIXA3RsMZWuvZwAcYd8HsBKqB73d2TK2M62rgZqWUDzgBXBuEZA7GL7brgY3m9WWA3wLZDWKz4py1Ji4rzlk68JJSyo6ReN7QWr9n9fexlXEF/ft4Kp15vqSLCSGE6ObC9dKQEEKIVpJEIIQQ3ZwkAiGE6OYkEQghRDcniUAIIbo5SQRCCNHNSSIQQohu7v8D4WzKF6X5KlcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "for site, dyn in results_open['dynamics'].items():\n", " plt.plot(*dyn.expectations(sz, real=True),\n", " color=f\"C{site}\", linestyle=\"solid\",\n", " label=f\"$\\\\langle s_{site}^z \\\\rangle$\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the environment starts to hybridize with the first spin and then the other spins start being affected too." ] } ], "metadata": { "interpreter": { "hash": "3306e98808c0871e8a1685f50cc307ae5b4a4a013844b10634a4efe89132c3fe" }, "kernelspec": { "display_name": "oqupyPR", "language": "python", "name": "oqupypr" }, "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.6.15" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 1, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }