{ "cells": [ { "cell_type": "markdown", "id": "915e9fb5-3434-4b36-8b7b-03f464d97893", "metadata": {}, "source": [ "2D Antiferromagnetic Model Example\n", "========" ] }, { "cell_type": "markdown", "id": "40b92c6e-f77e-4d70-9226-971fc56171aa", "metadata": {}, "source": [ "This is an example of how ``ikron`` can be used to deal with a large, ``ndim > 1``, Hilbert space.\n", "We'll define a simpler version of ``ham_heis_2D`` first:" ] }, { "cell_type": "code", "execution_count": 1, "id": "6e739766-50d6-49ea-aed3-49ae827799ef", "metadata": {}, "outputs": [], "source": [ "import itertools\n", "from operator import add\n", "import numpy as np\n", "from quimb import *\n", "\n", "def ham_heis_2D(n, m, j=1.0, bz=0.0, cyclic=False,\n", " sparse=True):\n", "\n", " dims = [[2] * m] * n # shape (n, m)\n", "\n", " # generate tuple of all site coordinates\n", " sites = tuple(itertools.product(range(n), range(m)))\n", "\n", " # generate neighbouring pairs of coordinates\n", " def gen_pairs():\n", " for i, j, in sites:\n", " above, right = (i + 1) % n, (j + 1) % m\n", " # ignore wraparound coordinates if not cyclic\n", " if cyclic or above != 0:\n", " yield ((i, j), (above, j))\n", " if cyclic or right != 0:\n", " yield ((i, j), (i, right))\n", "\n", " # generate all pairs of coordinates and directions\n", " pairs_ss = tuple(itertools.product(gen_pairs(), 'xyz'))\n", "\n", " # make XX, YY and ZZ interaction from pair_s\n", " # e.g. arg ([(3, 4), (3, 5)], 'z')\n", " def interactions(pair_s):\n", " pair, s = pair_s\n", " Sxyz = spin_operator(s, sparse=True)\n", " return ikron([j * Sxyz, Sxyz], dims, inds=pair)\n", "\n", " # function to make Z field at ``site``\n", " def fields(site):\n", " Sz = spin_operator('z', sparse=True)\n", " return ikron(bz * Sz, dims, inds=[site])\n", "\n", " # combine all terms\n", " all_terms = itertools.chain(map(interactions, pairs_ss),\n", " map(fields, sites) if bz != 0.0 else ())\n", " H = sum(all_terms)\n", "\n", " # can improve speed of e.g. eigensolving if known to be real\n", " if isreal(H):\n", " H = H.real\n", "\n", " if not sparse:\n", " H = qarray(H.toarray())\n", "\n", " return H" ] }, { "cell_type": "markdown", "id": "2087e16d-b0bc-4aed-b603-01f72b93e614", "metadata": {}, "source": [ "Note that in general, for 1D or flat problems, ``dims`` is just the list of subsystem hilbert space sizes and indices are just specified as single integers. For 2D+ problems ``dims`` should be nested lists, or equivalently, an array of shape ``(n, m, ...)``, with coordinates specifying subsystems then given as ``inds=[(i1, j1, ...), (i2, j2, ...), ...]``.\n", "\n", "We can now set up our parameters and generate the hamiltonian:" ] }, { "cell_type": "code", "execution_count": 2, "id": "d8c73b41-df30-48c5-9072-7cfccee76271", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2, 2, 2, 2, 2]\n", "[2, 2, 2, 2, 2]\n", "[2, 2, 2, 2, 2]\n", "[2, 2, 2, 2, 2]\n" ] } ], "source": [ "n, m = 4, 5\n", "dims = [[2] * m] * n\n", "\n", "for row in dims:\n", " print(row)" ] }, { "cell_type": "code", "execution_count": 3, "id": "10868158-99b4-4d21-a90a-b68b79d09fa1", "metadata": {}, "outputs": [], "source": [ "H = ham_heis_2D(n, m, cyclic=False)" ] }, { "cell_type": "markdown", "id": "0c9dcddd-3b83-4b19-bb07-bc6b6a151d56", "metadata": {}, "source": [ "Let's also break the symmetry, so that we get a nice alternating pattern, by adding a small field on the site ``(1, 2)``:" ] }, { "cell_type": "code", "execution_count": 4, "id": "b189979d-04c1-4eed-9a04-dc9ae824d961", "metadata": {}, "outputs": [], "source": [ "H = H + 0.2 * ikron(spin_operator('Z', sparse=True), dims, [(1, 2)])" ] }, { "cell_type": "markdown", "id": "78b7aa5f-f7ab-4aef-ab0a-33ef91125831", "metadata": {}, "source": [ "Note that ``ikron`` automatically performs a sparse kronecker product when given a sparse operator(s).\n", "\n", "Next we find its ground state and energy, which should be reasonably quick for 20 qubits\n", "(we could also make use of symmetries here to project the hamiltonian first):" ] }, { "cell_type": "code", "execution_count": 5, "id": "edf1a2d0-0c81-42ee-a3da-3dcc1aeca23e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 52.3 s, sys: 35.5 s, total: 1min 27s\n", "Wall time: 14.9 s\n" ] } ], "source": [ "%time ge, gs = eigh(H, k=1)" ] }, { "cell_type": "markdown", "id": "902666c9-3a4e-44b1-ab04-170c04ed7840", "metadata": {}, "source": [ "Giving us energy:" ] }, { "cell_type": "code", "execution_count": 6, "id": "d2de2242-5737-40aa-9cfd-91dde0924289", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-11.661573929790698" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ge[0]" ] }, { "cell_type": "markdown", "id": "d49f10d7-3e8c-4d2f-817a-9e3a22e7aebe", "metadata": {}, "source": [ "Now let's compute the magnetization at each site.\n", "First we construct the Z spin operators:" ] }, { "cell_type": "code", "execution_count": 7, "id": "a62c50da-6c34-4721-8704-ceb8109f2ab9", "metadata": {}, "outputs": [], "source": [ "Sz = spin_operator('Z', stype='coo')\n", "Sz_ij = [[ikron(Sz, dims, [(i, j)])\n", " for j in range(m)]\n", " for i in range(n)]" ] }, { "cell_type": "markdown", "id": "4d20c530-aa33-4556-acc5-ae471e0b80d9", "metadata": {}, "source": [ "Then we compute the expectation of each at each site:" ] }, { "cell_type": "code", "execution_count": 16, "id": "1dc1e5fb-24c4-42c6-a96c-d8b080483f41", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_2030217/891509411.py:11: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.\n", " plt.colorbar()\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAGiCAYAAAAiDFaYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzp0lEQVR4nO3dfbRddXno++8OkvgWqLcWwwWpjKOXKtUDAhE8DuDcc0EGXu2wOT3gSxFvKgpBzaUgBl8enrYSm2MjvsQbBYfH0he5LVq1JYq1vtASyAWlSoLUWhQ5kFCg0qgYYO91/5hzhZW111p7r7nWYs219/fjmCZzzt9vzt9eia4nz+83nznVaDSQJEmqmyXjHoAkSVInBimSJKmWDFIkSVItGaRIkqRaMkiRJEm1ZJAiSZJqySBFkiTVkkGKJEmqJYMUSZJUSwYpkiSplp40SOfMXAdcBnwoItb2aHcSsBE4ErgH2BARmwe5tyRJWtgqZ1Iy8zjgHOA7c7Q7HLgWuB44miKo+XBmrqp6b0mStPBVyqRk5tOBPwXeBLx7juZvAe5qybTcnpnHAhcC11S5vyRJWviqTvdsAv4mIv42M+cKUk4Arms79mVgdWbuHxGPtnfIzGXAsrbDS4EHIsLXNkuSxiozn0zxvTQMj0TEL4Z0rQWl7yAlM88EXgwcN88uK4Bdbcd2lfd+JnBvhz7rgGg/+PVbn8JXf+O/z3+wAuDMDdeOewgT65o3nzLuIUysVR//yriHMNE+847Txz2EifbNz180NaprZ+aTn/7Uhx/+6c+fMqxL7szMww1UZusrSMnMZwMfAk7t88Nsz35MdTnetJ5ioW3TcuDuPu4nSdKoLP3pz5/CW193DcuWzpoM6MueR/bnI3+6agVFVsYgpU2/mZRjgIOAWzKzeWw/4MTMPB9YFhHTbX12UmRTWh0EPAY80OkmEbEH2NPcb7mXJEm1sGzpowMHKeqt3yDlq8AL2459Cvge8IcdAhSArcAr246dCtzcaT2KJEmTYLoxw3RjZuBrqLu+gpSI2A3c1nosM39GsaD1tnJ/PXBIRJxVNtkMnJ+ZG4ErKBbSrgZeM+DYJUkamxkazHRdtTD/a6i7UVScPRg4rLkTEXcCpwMnA7cC7wHeFhE+fixJmlgzQ/qPuhuo4ixARJzctn92hzbfoHgiSJIkaV4GDlIkSVqMphsNphuDTdcM2n+h8wWDkiRV0FyTMuhW2paZOzJzzTh/prqpdSal/MNag8GUJGlhWxkR/z7uQdRNrYOUiNgEbMrMA4CHxj0eSZKaZmgw7dM9I1XrIEWSpLryEeTRcxpFkiTVkpkUSZIq8Ome0TNIkSSpgplyG/QaVWTmecBFFAVUtwNrI+L6Hu1Ponhx75HAPcCGiNjccv7rwEkdul4bEa8o21wKRNv5XRHR/n6+oXG6R5KkCZKZZwCXA+8DjgauB7Zk5mFd2h8OXFu2Oxq4DPhwZq5qafabFAFPc/t1YBr4i7bLbW9r1/4+v6EykyJJUgXTQ3i6p2L/C4BPRsSV5f7azHw5cC6wrkP7twB3RcTacv/2zDwWuBC4BiAiHmztkJlnAj9ndpDyWETsrDLoKgxSJEmqYLpRbINeo7Q8M1tP7YmIPe3tM3MpcAzw/rZT1wEv7XKbE8rzrb4MrM7M/SPi0Q59VgOfiYiftR1/XmbeA+wBbgIuiYh/6XLfgdV6uicz12TmDmDbuMciSVKrmSFtpbsp6oE1t04ZEYBnAvsBu9qO7wK6rQ1Z0aX9k8rr7SMzV1JM91zZduom4Czg5cCbyuvekJm/3OW+A6t1JsVibpKkReJQYHfL/qwsSpv2HM5Uh2Nzte90HIosym0RsU+CICK2tOx+NzO3Aj8A3kCxKHfoah2kSJJUVzNMMb33u776NUq751kW/36KBa3tWZODmJ0tadrZpf1jwAOtBzPzqcCZwHvnGkhE/Cwzvws8b+5hV1Pr6R5JkupqpjGcrR8R8QhwC3BK26lTgBu6dNvaof2pwM0d1qP8N2AZ8CdzjSUzlwHPB+6dq21VZlIkSZosG4GrMvNmigDkHOAwYDNAZq4HDomIs8r2m4HzM3MjcAXFQtrVwGs6XHs18FcR8UD7icz8APBF4C6KTMy7gQOATw/vR9uXmRRJkiqYLqd7Bt36FRFXA2sppmRuBU4ETo+IH5VNDqYIWprt7wROB04u278HeFtEXNN63cz834CXAZ/scutDgT8H7gA+CzwCHN9y36EzkyJJUgXTTLFkwDUpVde0RMTHgI91OXd2h2PfAF48xzX/CboPKCLO7G+UgzOTIkmSaslMiiRJFcw0pphpDPh0z4D9F7paBymZuQZYgxkfSVLNjHO6Z7Go9Zd/RGyKiBcAK8c9FkmSRmhbZu4o/3GuUq0zKZIk1dU0S1gy4L/1px/vv3KexdwWFYMUSZIqaAxhTUrDNSk9GaRIklTBNFNMuSZlpGq9JkWSJC1eZlIkSapgurGEqcaAa1IG7L/QGaRIklTBDFPMDDghMeN0T0+GcJIkqZbMpEiSVIELZ0ev1kGKFWclSXXlmpTRq/WnY8VZSdIiYcXZDmqdSZEkqa6KhbMDvmDw8f5WnO3AIEWSpApmWNJa1r7yNdSdn44kSaqlvjIpmXkucC7wnPLQduD3ImJLl/YnA1/rcOr5EfG9fu4tSVKdTDeWgAtnR6rfT+du4J3AseX2d8DnM/PIOfodARzcsn2/z/tKklQrMywZyqbu+sqkRMQX2w69q8yuHE+RVenmvoj4SZ9jkySptqYbUzDgW4ynfQtyT5UXzmbmfsBvAU8Dts7R/NuZ+WRgB/AHEdFpCqj12suAZS2HllcdpyRJmkx9BymZ+UKKoOTJwE+BV0fEji7N7wXOAW6hCDp+G/hqZp4cEd/scZt1QPQ7NkmSnijFkz0DrklxuqenKpmUO4CjgF8CVgGfzsyTOgUqEXFH2b5pa2Y+G7gQ6BWkrAc2tuwvp1gPI0lSLcwweMVZ16T01neQEhGPAP9c7t6cmccBbwfePM9L3Ai8fo577AH2NPczs99hSpI0SbZl5gywKSI2jXswdTGMYm5T7Lt+ZC5HU0wDSZI0sYY83WPF2Q76rZNyGbAF+DHFFMyZwMnAaeX59cAhEXFWub8W+CHFkz9LKTIoq8pNkqSJNTOEp3tmfLqnp34zKc8CrqKodfIQ8B3gtIj4Snn+YOCwlvZLgQ8AhwAPUwQrr4iIawcZtCRJWvj6rZOyeo7zZ7ftbwA29D8sSZLqrXg54KALZ82k9OILBiVJqmC6sYRGozHQNWYsi9+Tn44kSaqlWmdSMnMNsAaDKUlSzRRTNQMunHW6p6daBynls+KbMvMAioW6kiTVgtM9o1frIEWSpLqaZgkNBgxSnCjoyU9HkiTVkkGKJEkVzDSmhrKVtmXmjnItpkpO90iSVMHMEKZ7GpbF78lMiiRJqiUzKZIkVTDTWMLUgE/3NHy6pyeDFEmSKphmiqkB65w0rJPSU62DFIu5SZK0eNU6SLGYmySprsY53ZOZ5wEXAQcD24G1EXF9j/YnARuBI4F7gA0Rsbnl/NnApzp0fUpE/KLqfQdlhkKSpAqmmRrK1q/MPAO4HHgfcDRwPbAlMw/r0v5w4Nqy3dHAZcCHM3NVW9N/pwg+9m5tAUpf9x2GWmdSJEnSLBcAn4yIK8v9tZn5cuBcYF2H9m8B7oqIteX+7Zl5LHAhcE1Lu0ZE7BzifQdmkCJJUgUzjSmmBnw6p/F4Mbflmdl6ak9E7Glvn5lLgWOA97edug54aZfbnFCeb/VlYHVm7h8Rj5bHnp6ZPwL2A24F3hMR3x7gvgNzukeSpAqmG0uGspXuplh72dy6ZSaeSRFE7Go7vgtY0aXPii7tn1ReD+B7wNnAq4DXAL8A/iEznzfAfQdmJkWSpAoaTA3hEeK9DzEfCuxuOTErizLr9u0Xmn1srvZ7j0fEjcCNzZOZ+Q/At4C3Am8b4L4DMUiRJGn8ds+zLP79wDSzsxcHMTvL0bSzS/vHgAc6dYiImcz8/4BmJqXKfQfmdI8kSRUMebpnXiLiEeAW4JS2U6cAN3TptrVD+1OBm1vWo+wjM6eAo4B7B7jvwMykSJJUwUxjChoDTvdU678RuCozb6YIQM4BDgM2A2TmeuCQiDirbL8ZOD8zNwJXUCykXU2x9oSyT1BM93wfOIBiiucoioKq87rvKNQ6k5KZazJzB7Bt3GORJKkOIuJqYC3wXoqncE4ETo+IH5VNDqYIHprt7wROB04u278HeFtEtD5+/EvAJ4DbKZ7YOQQ4MSK2tVxnrvsO3VRjwGp5T4Rmxdmv3/oUpmd8z0G/ztxw7biHMLGueXN7ZlPzterjXxn3ECbaZ95x+riHMNG++fmLRvZl0fxOevj/3A77zwx2sUeX8JS/PhLgwHmuSVlUnO6RJKmCMU73LBq1nu6RJEmLl5kUSZIqmBnKv/PNFfRikCJJUgXTw53u2ZaZM8CmiNg04NAWDIMUSZLGb6ULZ2czSJEkqQIXzo6eQYokSRU0GksYtIrHoG9RXuhqHaRk5hqKanf+KUqSamV6CC8YnBr4BYULW62//CNiU0S8AFg57rFIkqQnVq0zKZIk1dVMAxoDrimZqn/R97EySJEkqYKZxhIGfbXMlAtne6r1dI8kSVq8+sqkZOa5wLnAc8pD24Hfi4gtPfqcRPF65yOBe4ANETGy1zpLkvREmGGKQWdrXDjbW7+ZlLuBdwLHltvfAZ/PzCM7Nc7Mw4FrgeuBo4HLgA9n5qrKI5YkqQamG1ND2UrbMnNH+VSrSn1lUiLii22H3lVmV46nyKq0ewtwV0SsLfdvz8xjgQuBa/ocqyRJC5UVZzuovHA2M/cDfgt4GrC1S7MTgOvajn0ZWJ2Z+0fEo12uvQxY1nJoedVxSpI0Ci6cHb2+g5TMfCFFUPJk4KfAqyNiR5fmK4Bdbcd2lfd9JnBvl37rgOh3bJIkPVFmmBq84qxrUnqqkkm5AzgK+CVgFfDpzDypR6DS/kc41eV4q/UUi22blgN3r/r961iy/3TfA17s/uKtp417CBPrj/7Yl5FW9bu/49T6IM78yLXjHsKEu2jcA9AQ9B2kRMQjwD+Xuzdn5nHA24E3d2i+kyKb0uog4DHggR732APsae5nZr/DlCRppBpMMTPgNZaYSelpGMXcpth3/UirrcAr246dCtzcbT2KJEmTYKYxxcygzyC7JqWnfuukXAZsAX5MMQVzJnAycFp5fj1wSEScVXbZDJyfmRuBKygW0q4GXjOMwUuSNC4zjSXMDLooxSClp37rpDwLuIpiXcpXgZcAp0XEV8rzBwOHNRtHxJ3A6RSBzK3Ae4C3RYSPH0uSpJ76rZOyeo7zZ3c49g3gxf0NS5KkenO6Z/R8d48kSRXMMDWUrWTF2Q58C7IkSeNnxdkOah2klBHlGsz4SJJqxume0at1kBIRm4BNmXkA8NC4xyNJUpNByuiZoZAkSbVU60yKJEl1ZSZl9AxSJEmqwCBl9JzukSRJtWQmRZKkChrQWuekEvMovRmkSJJUwTCme6YGnS5a4JzukSSpgiJIGXwrWXG2AzMpkiSNnxVnO6h1kGLFWUlSXTndM3q1DlKsOCtJqiuDlNEzQyFJkmqp1pkUSZLqqtGYojFgJmTQ/gudQYokSRXMMMXMgNewTkpvBimSJE2YzDwPuAg4GNgOrI2I63u0PwnYCBwJ3ANsiIjNLeffBJwF/Hp56BbgkojY1tLmUiDaLr0rIlYM/AN14ZoUSZIqGHKdlHnLzDOAy4H3AUcD1wNbMvOwLu0PB64t2x0NXAZ8ODNXtTQ7Gfhz4D8DJwB3Addl5iFtl9tOERg1txf2/QP0wUyKJEkVDHlNyvLMbD21JyL2dOl2AfDJiLiy3F+bmS8HzgXWdWj/FuCuiFhb7t+emccCFwLXAETE61o7lJmV/wr8F+CPW049FhE75/7JhsMgRZKk8bu7bT+BS9sbZeZS4Bjg/W2nrgNe2uXaJ5TnW30ZWJ2Z+0fEox36PBXYH3iw7fjzMvMeYA9wE8WU0L90ue/Aah2kWMxNklRXQ66Tciiwu+VUtyzKM4H9gF1tx3cB3daGrOjS/knl9e7t0Of9wP8E/rbl2E0U61b+CXgW8G7ghsw8MiIe6HLvgdQ6SLGYmySptoYw3cPj/Xf3WRa//c5THY7N1b7TcTLzHcBrgJMj4hfN4xGxpaXZdzNzK/AD4A0Ui3KHzgyFJEkVjGnh7P3ANLOzJgcxO1vStLNL+8eAfTIgmXkhcAlwakR8p9dAIuJnwHeB581r5BUYpEiSNCEi4hGKx4NPaTt1CnBDl25bO7Q/Fbi5dT1KZl4EvAc4LSJunmssmbkMeD6dp4uGotbTPZIk1VWDwSvGVuy+EbgqM2+mCEDOAQ4DNgNk5nrgkIg4q2y/GTg/MzcCV1AspF1NMaVD2ecdwO8DrwV+mJnNzMtPI+KnZZsPAF+keDz5IIo1KQcAn672Y8zNTIokSRUUFWcH3/oVEVcDa4H3ArcCJwKnR8SPyiYHUwQtzfZ3AqdT1EK5lSJb8raIuKblsucBS4G/pMiMNLcLW9ocSlFL5Q7gs8AjwPEt9x06MymSJE2YiPgY8LEu587ucOwbwIt7XO8587jnmfMf4XAYpEiSVIEvGBw9gxRJkioYRp2UQfsvdLUOUizmJklaJLZl5gywqawRJmoepFjMTZJUV43GEJ7uebz/yj6LuS0KtQ5SJEmqK9ekjJ7TKJIkqZbMpEiSVIGZlNEzSJEkqQKf7hm9voKUzFwH/Cbwa8DDFO8JuDgi7ujR52Tgax1OPT8ivtfP/SVJqoshL5xVB/2uSTkJ2AQcT/GyoicB12Xm0+bR9wiKUr3N7ft93luSJC0ifWVSIuK01v3MfCNwH3AM8M05ut8XET+Zz33KNysuazm0vI9hSpI0cq5JGb1B16QcWP764DzafjsznwzsAP4gIjpNATWtA2LAsUmSNDLFdE//Lwjc9xpGKb1UfgQ5M6coXhf99xFxW4+m91K8RnoVxXqWO4CvZuaJPfqspwiAmtuhVccpSdIE2JaZO8pK6yoNkkn5KPAi4GW9GpWLalsX1m7NzGdTvP654xRRROwB9jT3M3OAYUqSNHyNchv0GiUrznZQKZOSmR8BXgX854i4u8IlbgSeV+XekiTVQbEmZfBN3fX7CPIU8BHg1cDJEXFnxfseTTENJEmS1FG/0z2bgNcCvwHszswV5fGHIuJhgMxcDxwSEWeV+2uBHwLbgaXA6ynWp6wadPCSJI3NkOd7NFu/Qcq55a9fbzv+RuB/lL8/GDis5dxS4APAIRQF4LYDr4iIa/u8tyRJ9TGM6ZrG3v9SB/3WSZnzTyMizm7b3wBs6G9YkiTVmxVnR8+3IEuSpFryBYOSJFUwjKdzzKT0VusgpSxqswYzPpKkumlMFdtA1xjOUBaqWn/5R8SmiHgBsHLcY5EkSU+sWgcpkiTVVXPh7KBbybL4HdR6ukeSpNoabp0Uy+J3YCZFkiTVkpkUSZIqaDCEp3uGNJaFyiBFkqQqLIs/ck73SJKkWjKTIklSBRZzG71aBykWc5Mk1ZbTPSNX6y9/i7lJkuprakibuql1kCJJkhYvgxRJkqpoDGkrWHG2g1qvSZEkqbasODtyZlIkSVItmUmRJKmKxlSxDXSN4QxloTJIkSSpgra3GFe+hrpzukeSJNVSrTMpFnOTJNXWGIu5ZeZ5wEXAwcB2YG1EXN+j/UnARuBI4B5gQ0RsbmuzCvh94D8APwDeFRGfG+S+g6r1l7/F3CRJtdVckzLo1qfMPAO4HHgfcDRwPbAlMw/r0v5w4Nqy3dHAZcCHy6Ck2eYE4GrgKuA/lr/+v5n5kqr3HYZaZ1IkSdIsFwCfjIgry/21mfly4FxgXYf2bwHuioi15f7tmXkscCFwTfMawFciYn25v77MvqwFXlPxvgMzSJEkqYKpRrENeo3S8sxsPbUnIva0t8/MpcAxwPvbTl0HvLTLbU4oz7f6MrA6M/ePiEfLNh/s0GbtAPcdWK2neyRJqq3hVpy9G3ioZeuWmXgmsB+wq+34LmBFlz4rurR/Unm9Xm2a16xy34GZSZEkqYrh1kk5FNjdcmZWFqVrz8JUh2NztW8/Pp9r9nvfgRikSJI0frvnWRb/fmCa2dmLg5id5Wja2aX9Y8ADc7RpXrPKfQfmdI8kSVUMd7pnXiLiEeAW4JS2U6cAN3TptrVD+1OBm8v1KL3a3DDAfQdmJkWSpCrGVydlI3BVZt5MEVycAxwGbAbIzPXAIRFxVtl+M3B+Zm4ErqBYJLuax5/aAfgQ8M3MvBj4PPAbwP8BvGy+9x2FWmdSMnNNZu4Ato17LJIk1UFEXE3x1M17gVuBE4HTI+JHZZODKYKHZvs7gdOBk8v27wHeFhHXtLS5ATgTeCPwHeBs4IyIuKmP+w7dVGMCXhyQmQcAD/3yqutZsv/0uIczcf7iraeNewgT64+u3DTuIUys3/2dNeMewkT7rY98adxDmGjnHfG1AVe0dtf8TvqTX34qjy4Z7Db7zzR4/QM/BzhwnmtSFhWneyRJqsK3II9crad7JEnS4mUmRZKkCoZccVYdmEmRJKmq4T1+vC0zd2Smi7la9JVJycx1wG8CvwY8TPFs9MURcccc/eZ8RbQkSYvYShfOztZvJuUkYBNwPEUBlycB12Xm07p1mM8roiVJktr1lUmJiH2eZc3MNwL3UbwZ8Ztdus3nFdH7yMxlwLKWQ8v7GackSaM21Xj8BTiDXEPdDbom5cDy1wd7tOn2iuhjM3P/Ln3Wse/bIO8eZJCSJA1d8xHkQTd1Vfnpnsycolhn8vcRcVuPpnO9IvreDn3Wl9duWg7c/bl3nsL0tH+g/frQlR8Z9xAm1oXPOX7cQ5hYl//wo+MewkR7+5vfOu4hTLTztox7BBqGQR5B/ijwIvat69/NfF4RvVdE7KHlNdWZWWV8kiSNzjCmapzu6alSkJKZHwFeBZwYEXNNxcznFdGSJE0Wg5SR6/cR5CngI8CrgZPLlxbNZSvwyrZj7a+IliRJ2ke/mZRNwGspXuG8OzObGZKHIuJhqPyKaEmSJopP94xev0/3nEvxRM/XKRa8NrczWtr0/YpoSZImzqDVZvetOmvF2Q76rZMyZ9AYEWd3OPYN4MX93EuSpEXEirMd+IJBSZKqcOHsyBmkSJJUgWtSRq/WQUo5N7cG39YsSdKiU+sv/4jYFBEvAFaOeyySJO3DsvgjV+tMiiRJteWalJEzSJEkqQLXpIxerad7JEnS4mUmRZKkKpzuGTkzKZIkVdEop3wG2Kw425uZFEmSxs+Ksx0YpEiSVIXTPSNX6yDFYm6SpNoySBm5Wn/5W8xNkqTFq9aZFEmS6so6KaNX60yKJElavAxSJElSLTndI0lSFS6cHTmDFEmSKphi8DUlvgO5N6d7JEmqojGkrWDF2Q7MpEiSNH5WnO2g1kGKxdwkSbU1rPUkzvl0Vesvf4u5SZLqatCXC+59yaC6qnWQIkmSFq9aT/dIklRbZkFGziBFkqQKhlIWfygjWbgMUiRJWoAy8xnAh4FXlYe+ALw1In7So88UEMA5wDOAm4A1EbG9PP+/AAmcCjwbuB/4K+A9EfFQy3V+CPxq2+X/MCLe2c/P4JoUSZKqGG6dlFH4M+Ao4LRyOwq4ao4+7wAuAM4HjgN2Al/JzOXl+f+13C4EXgicXV77kx2u9V7g4JbtD/r9AcykSJJUxXADjOWZ2bq/JyL2VL1YZj6fIng4PiJuKo+9CdiamUdExB0d+kwBa4H3RcRny2NvAHYBrwU+HhG3Aatauv0gM98F/ElmPikiHms5tzsidlb9GcBMiiRJdXA38FDLtm7A650APNQMUAAi4sby2i/t0udwYAVwXUufPcA3evQBOBD497YABeDizHwgM2/NzHdl5tJ+f4haZ1Is5iZJqqshL5w9FNjdcqpyFqW0Arivw/H7ynPd+kCROWm1i9nrSwDIzF8G3gN8vO3Uh4BvAf9GUetsPUUQ9DtzDbxVrYOUiNgEbMrMAyiiP0mS6mG40z2751MWPzMvpVjY2stx5a+dRjjV5Xir9vMd+5TfzX8D7KBYTLtXRHywZfc7mflvwF9m5sUR8cAc99+r1kGKJEm1NZ46KR8FPjNHmx8CLwKe1eHcrzA7U9LUXD+yAri35fhB7X3KhbRfAn4KvDoiHp1jTDeWvz4XMEiRJGmhiYj7KR777SkztwIHZubKiNhWHnsJxfqRG7p0u5MiUDkF+HbZZylwEnBxy7UPAL5MMSX1qoj4xTyGfnT56709W7UxSJEkqYI6F3OLiNsz80vAFZn55vLwJ4C/bn2yJzO/B6yLiM9FRCMzLwcuyczvA98HLgF+TvE4czODch3wVOD1wAFl0ALwrxExnZknAMcDX6NYqnEc8EHgCxFxVz8/h0GKJElV1L8s/usoirk1n9b5AkX9k1ZHUGRXmjYATwE+xuPF3E6NiOai3mOAl5S//+e2ax1OMdW0BziDYu3MMuBHwBXltfvSd5CSmScCF5UDPZhiLuqverQ/mSKaavf8iPhev/eXJElzi4gHKbIdvdpMte03gEvLrVP7rzNHAigivkWRSRlYlUzK04B/BD4FXNNHvyOA1pXL/1rh3pIk1UKdp3sWir6DlIjYAmwBaKuON5f7er0vQJKkiVL/6Z6J90SuSfl2Zj6Z4nnqP4iITlNAAGTmMop5rKbl3dpKkrQAbMvMGWBTWSNMPDFByr0Ub1O8hSLw+G3gq5l5ckR8s0ufdcxdrEaSpPEZbiZl5XyKuS02Iw9SykedWl9ktDUzn03xBsVuQcp6YGPL/nKK9xpIklQLU7gmZdTG9QjyjfRYcVy+0Gjvewv6XPsiSZIWgHEFKUfTZ9U5SZJqxYWzI1elTsrTKWrvNx2emUcBD0bEXZm5HjgkIs4q26+lKO6yHVhKkUFZVW6SJE0kH0EevSqZlGPZtzhbc+3Ip4GzKQq8HdZyfinwAeAQ4GGKYOUVEXFthXtLklQPZlJGrkqdlK/TI/iLiLPb9jdQoRSuJEla3Hx3jyRJVQ2aTXG+pyeDFEmSKphqFNtA1xjOUBasWgcpmbkGWAMsGfdYJEnSE6vWX/4RsSkiXgCsHPdYJEnaR2NIW2FbZu4o/3GuUq0zKZIk1dWQp3ssi99BrTMpkiRp8TKTIklSFftO12gEDFIkSarAp3tGz+keSZJUS2ZSJEmqwumekTNIkSSpCoOUkat1kGIxN0lSXbkmZfRq/eVvMTdJkhavWgcpkiTVlhVnR67W0z2SJNXVVKMxhOmevRew4mwHZlIkSVItmUmRJKkKn+4ZOYMUSZIq8Ome0XO6R5Ik1ZKZFEmSqnC6Z+RqHaRYzE2SVFdO94xerb/8LeYmSdLiVetMiiRJteV0z8jVOpMiSVJdNad7Bt1KVpztwEyKJElVDDeTYsXZDsykSJKkWjKTIklSBVMM4ekeH+/pySBFkqQqGo3Bp3sarrztxekeSZJUS2ZSJEmqYgjF3Eb5CHNmPgP4MPCq8tAXgLdGxE969JkCAjgHeAZwE7AmIra3tPk6cFJb16sj4sxB7t1JrTMpmbkmM3cA28Y9FkmS9tEY0jY6fwYcBZxWbkcBV83R5x3ABcD5wHHATuArmbm8rd0VwMEt25uHcO9Zap1JiYhNwKbMPAB4aNzjkSRpEmTm8ymCg+Mj4qby2JuArZl5RETc0aHPFLAWeF9EfLY89gZgF/Ba4OMtzX8eETuHde9uah2kSJJUV1MzxTbQNR7/7fLMbD21JyL2DHDpE4CHmkECQETcmJkPAS8FOgUKhwMrgOta+uzJzG+UfVqDlNdl5uspApgtQEbE7gHu3ZFBiiRJVQxjuubx/ne3nUng0gGuvAK4r8Px+8pz3fpAEXi02gX8asv+nwJ3UkwF/TqwHviPwCkD3LsjgxRJksbvUGB3y37HLEpmXkqxsLWX48pfO4VQU12Ot2o/v0+fiLii5dxtmfl94ObMfHFEfGvAe++j7yAlM08ELgKOoVgs8+qI+Ks5+pwEbASOBO4BNkTE5n7vLUlSXbS9e6fyNUq751kW/6PAZ+Zo80PgRcCzOpz7FWZnSpqaa0xWAPe2HD+oRx+AbwGPAs8rf7+zwr07qpJJeRrwj8CngGvmapyZhwPXUqwEfj3wn4CPZea/RsSc/SVJqqUxFHOLiPuB++dql5lbgQMzc2VEbCuPvQQ4ELihS7fmFM4pwLfLPkspHje+uMftjgT25/HApsq9O+o7SImILRSLZGhb5NPNW4C7ImJtuX97Zh4LXMg8ghxJkupoyJmUoYqI2zPzS8AVmdl8PPgTwF+3Pl2Tmd8D1kXE5yKikZmXA5eUUzjfBy4Bfk7xSDGZ+R+A11EkH+4HXgD8EUVQ8w/93Hs+nog1KSfQslK49GVgdWbuHxGPtnfIzGXAspZD7c9nS5Kk3l5HUVCt+R38BYr6J62OoMhwNG0AngJ8jMeLuZ3a8uTOI8B/Ad4OPB34MfA3FE/3TPd57zk9EUHKCjqvFH4S8Ez2nfdqWkeHhUFv/cDV7Lf0saEPcKG74Ozzxj2EifWVez417iFMrFNe69+7QfzulX8y7iFMuHeM/hbDfbpn6CLiQYplFr3aTLXtNyieKrq0S/sfM7vabKV7z8cTVXG200rhTseb1lNEds3t0BGNS5KkSprTPYNu6u6JyKTsZPZz0QcBjwEPdOpQFrDZ+/jVPNe+SJKkBeSJCFK2Aq9sO3YqcHOn9SiSJE2EMTzds9hUqZPydOC5LYcOz8yjgAcj4q7MXA8cEhFnlec3A+dn5kaKx5BPAFYDrxlo5JIkjdGQn+7ZlpkzwKbyvXWiWiblWOBrLfsby18/DZxNUeDtsObJiLgzM08HPgisoSjm9jZrpEiStNfKeRZzW1Sq1En5Ovu8E2nW+bM7HPsG8OJ+7yVJUm3V/OmehcB390iSVEGdi7ktFE/UI8iSJEl9qXUmJTPXUKxjMZiSJNXLTANmBr2IqZReav3lHxGbIuIFwMpxj0WSpFkaA27qqdaZFEmS6so1KaNX60yKJElavMykSJJUhRVnR85MiiRJFQz5BYPbMnNH+cCISmZSJEkaPyvOdmCQIklSFVacHTmDFEmSKphqNIbwdI9RSi+1DlIs5iZJ0uJV6y9/i7lJkmprZkibuqp1JkWSpLpyumf0ap1JkSRJi5eZFEmSqvDpnpEzSJEkqQorzo6c0z2SJFVgxdnRM5MiSdL4WXG2A4MUSZKqcLpn5AxSJEmqYGqm2Aa6xtRwxrJQ1TpIseKsJEmLV62//K04K0mqrUZjOJu6qnUmRZKk2rJOysjVOpMiSZIWLzMpkiRV4Lt7Rs8gRZKkSobwCLLzPT053SNJ0vhZcbYDMymSJFUxU26DXqNgxdkODFIkSaqgePfOYNM1g65pWehqHaRYzE2SVFtDKYsPYNnZbmr95W8xN0mSFq9aZ1IkSaotMykjZ5AiSVIVw104qw4MUiRJWoAy8xnAh4FXlYe+ALw1In7So88UEMA5wDOAm4A1EbG9PP8c4M4u3f9bRPxF2e6HwK+2nf/DiHhnPz9DpSAlM88DLgIOBrYDayPi+i5tTwa+1uHU8yPie1XuL0nSuA2n4uxwxtLFnwGHAqeV+58ArgJe2aPPO4ALgLOBfwLeDXwlM4+IiN3Ajym++1udU/bb0nb8vcAVLfs/7fcH6DtIycwzgMuB84B/AN4MbMnMF0TEXT26HgG0PgP+r/3eW5Kk2hjamhQAlmdm65k9EbGn6mUz8/kUwcnxEXFTeexNwNYy4LijQ58pYC3wvoj4bHnsDcAu4LXAxyNiGtjZ1u/VwNUR0R6E7I6InQygSiblAuCTEXFlub82M18OnAus69Hvvl4pJkmSFrG72/YTuHSA650APNQMUAAi4sbMfAh4KTArSAEOB1YA17X02ZOZ3yj7fLy9Q2YeAxxFUS6k3cWZ+R6K7MtfAP89Ih7p54foK0jJzKXAMcD7205dR/ED9PLtzHwysAP4g4joNAXUvM8yYFnLoeX9jFOSpJEbbiblUGB3y5nKWZTSCuC+DsfvK8916wNF5qTVLmavL2laDdweETe0Hf8Q8C3g3yjKiKynCIJ+p/ew99VvJuWZwH50/gG6/dD3UsxX3UIRePw28NXMPDkivtmlzzqKhTuSJNXTcIOU3fMpi5+ZlzL39+Nxs67+uKkuxzuPqkefzHwKxTTQ77efi4gPtux+JzP/DfjLzLw4Ih6Y4/57VX26Z14/AEA579WaVtqamc8GLgS6BSnrgY0t+8uZnQqTJGmx+SjwmTna/BB4EfCsDud+hdmJhqbm+pEVFAmGpoO69PmvwFOBP55jPAA3lr8+FxhZkHI/MM3srEm3H6CbG4HXdztZLhbam+pqW0wkSdL4jaFOSkTcT/Fd3FNmbgUOzMyVEbGtPPYS4ECgfWqm6U6KQOUU4Ntln6XAScDFHdqvBr4QEfN5EObo8td7e7Zq01eQEhGPZOYtFD/A51pOnQJ8vo9LHU2fA5UkqU7q/AhyRNyemV8CrsjMN5eHPwH8deuTPZn5PWBdRHwuIhqZeTlwSWZ+H/g+cAnwc4rHmWnp91zgROD09ntn5gnA8RTlRx6imH76IEVA0+sp4FmqTPdsBK7KzJuBrRTrTQ4DNpeDWw8cEhFnlftrKVJP24GlFBmUVeUmSdJkGu6alFF4HUUxt+bTOl8Azm9rcwRFdqVpA/AU4GM8Xszt1LJGSqv/C/ifLddutQc4g2LtzDLgRxT1Ujb0+wP0HaRExNWZ+csURVoOBm4DTo+IH5VNDqYIWpqWAh8ADgEepghWXhER1/Z7b0mSND8R8SA9llaUbaba9hsUjz5fOke/SyiyLJ3OfYsikzKwSgtnI+JjFFFWp3Nnt+1voEL0JElSrc00fHfPiPnuHkmSqqj/dM/EWzLuAUiSJHVS60xKZq6hKLVrMCVJqpkhZFIety0zZ4BNEbFpaFedcLUOUso/qE2ZeQDFY0ySJNXDcKd7Vs6n4uxiY4ZCkiTVUq0zKZIk1ZZP94ycQYokSVU0Zny6Z8Sc7pEkSbVkJkWSpCqskzJyBimSJFXhmpSRM0iRJKkKMykj55oUSZJUS7XOpFhxVpJUWw2KbMpA19j7EmIrznZQ6yDFirOSpNpqNIYQpOz9nRVnOzBDIUmSaqnWmRRJkmprZqZ4wmega0zN3WYRM0iRJKmK4U73qAOneyRJUi2ZSZEkqQozKSNnkCJJUhUzjSGsSRnOUBYqp3skSVIt1TqTYjE3SVJdNRozNAac7mk0fLqnl1p/+UfEpoh4AbBy3GORJGkfjcbjUz5Vt8eDnG2ZuaP8x7lKtc6kSJJUW1acHblaZ1IkSdLiZSZFkqQqrDg7cgYpkiRVYZ2UkXO6R5Ik1ZKZFEmSKmjMzNAYcLqn4XRPTwYpkiRV4XTPyNU6SLGYmyRJi1etg5SI2ARsyswDgIfGPR5Jkvby3T0jZ4ZCkqQqGg1ozAy4WXG2l1pnUiRJWiSsONuBQYokSRU0ZhpDeLrHlbO9GKRIklRFc8pmoGsMZygLVaUgJTPPAy4CDga2A2sj4voe7U8CNgJHAvcAGyJic5V7S5JUB2ZSRq/vhbOZeQZwOfA+4GjgemBLZh7Wpf3hwLVlu6OBy4APZ+aqimOWJEmLQJVMygXAJyPiynJ/bWa+HDgXWNeh/VuAuyJibbl/e2YeC1wIXNPpBpm5DFjWcmg5wPSjzk5Vsd9+RupV7XnEv3NV+fduMNP+3RtIWbpid0SM7C/ifksZeLpmv6VDGcqC1df/CjJzKXAM8P62U9cBL+3S7YTyfKsvA6szc/+IeLRDn3VAtB+87bOv6Ge4Kv2nl3X6iDUfGz51xriHMLH8ezeYf7z6N8Y9hEn3EPArwP0juPYjwM6VF/zaiiFdb2d5TbXpN1R/JrAfsKvt+C6g2x/Wii7tn1Re794OfdZTrGFpWg7cDRwK7O5vyHttA1ZW7Dto/3Hee9I/u3H3H/Tzm+SffdD+i/3vnp/d+P/ujeSLPyJ+US5lGFYe5JGI+MWQrrWgVM0ntie4pjocm6t9p+MARMQeYE9zPzObv91d9TnyzJwZ5Bn0QfqP+d7N307kZzfu/oN+fpP8sw/af7H/3fOzq8XfvZEpgwoDixHrd+Hs/cA0s7MmBzE7W9K0s0v7x4AH+rz/IDaNsf847z0M4x7/uPuP896T3n9Q4x7/OP93P6hJ/uyG0V8LwFSjzzc4ZuZNwC0RcV7LsR3A5yNi1sLZzPxD4JUR8YKWY/8PcFREnDDPezbf3XOgFfn642c3GD+/6vzsqvOzG4yf38JRZbpnI3BVZt4MbAXOAQ4DNgNk5nrgkIg4q2y/GTg/MzcCV1AspF0NvKaPe+4BkpYpIM2bn91g/Pyq87Orzs9uMH5+C0TfmRTYW8ztHRTF3G4D/u+I+GZ57n8Az4mIk1vanwR8kMeLuf2hxdwkSVIvlYIUSZKkUeu74qwkSdITwSBFkiTVkkGKJEmqJYMUSZJUS7V/g1X5JNFFFE8SbQfWRsT14x1V/WXmiRSf2zEUn92rI+KvxjqoCZGZ64DfBH4NeBi4Abg4Iu4Y68AmRGaeS/HC0eeUh7YDvxcRW8Y2qAlV/l28DPhQy0ta1UFmXsrsd77tiohhvV9HY1DrTEpmngFcDrwPOBq4HtiSmYeNc1wT4mnAPwLnj3sgE+gkimqXxwOnUATz12Xm08Y6qslxN/BO4Nhy+zvg85l55FhHNWEy8ziKOlTfGfdYJsh2in+UNbcXjnc4GlTdMykXAJ+MiCvL/bWZ+XKKf6XNqm6rx5X/at0CT8x7LBaSiDitdT8z3wjcR5GV+uZYBjVBIuKLbYfeVWZXjqf4EtEcMvPpwJ8CbwLePebhTJLHImLnuAeh4altkJKZSym+FN7fduo64KVP/Ii0iB1Y/vrgWEcxgTJzP+C3KDJ7W8c8nEmyCfibiPjbzDRImb/nZeY9FJVmbwIuiYh/GfOYNIDaBinAM4H9mP3iwl3MfmGhNBKZOUXxKoi/j4jbxj2eSZGZL6QISp4M/JRiTdSO8Y5qMmTmmcCLgePGPZYJcxNwFvBPwLMoMlA3ZOaREfFEvsxWQ1TnIKWpvSTuVIdj0qh8FHgR8LJxD2TC3AEcBfwSsAr4dGaeZKDSW2Y+G/gQcGpE/GLc45kkbQuzv5uZW4EfAG+g+IeGJlCdg5T7gWlmZ00OYnZ2RRq6zPwI8CrgxIi4e9zjmSQR8Qjwz+XuzeUi0LcDbx7fqCbCMRT/H3dLy1qy/YATM/N8YFlETI9rcJMkIn6Wmd8Fnjfusai62j7dU/6f3C0UT1e0OoXikVBpJDJzKjM/SvEY8v8eEXeOe0wLwBSwbNyDmABfpXgi5aiW7WaKRbRHGaDMX2YuA54P3Dvusai6OmdSoEjRXZWZN1PMb58DHAb4BuU5lE8HPLfl0OGZeRTwYETcNZ5RTYxNwGuB3wB2Z2Yzm/dQRDw8vmFNhsy8jOLJsh8Dy4EzgZOB03p0ExARuyneLL9XZv4MeMA1Ub1l5geALwJ3UWSj3g0cAHx6nOPSYGqbSQGIiKuBtcB7gVuBE4HTI+JHYxzWpDgW+Ha5QRHwfRv4vbGNaHKcS/FEz9cp/hXW3M4Y45gmybOAqyjWpXwVeAlwWkR8Zayj0kJ3KPDnFH/vPgs8Ahzv98Vkm2o0XIMqSZLqp9aZFEmStHgZpEiSpFoySJEkSbVkkCJJkmrJIEWSJNWSQYokSaolgxRJklRLBimSJKmWDFIkSVItGaRIkqRaMkiRJEm19P8D+gcR6ZvOqe8AAAAASUVORK5CYII=", "image/svg+xml": [ "" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "m_ij = [[expec(Sz_ij[i][j], gs)\n", " for j in range(m)]\n", " for i in range(n)]\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "with plt.style.context(NEUTRAL_STYLE):\n", " plt.grid(False)\n", " plt.pcolor(m_ij)\n", " plt.colorbar()" ] }, { "cell_type": "markdown", "id": "65a2a031-56b7-4e9f-952c-33b504b3fd2f", "metadata": {}, "source": [ "Which looks pretty much as expected.\n", "\n", "Alternatively to using global operators, we could also use ``partial_trace`` to look at the state of a few qubits. For example, let's find the correlations between the spin we added the small field to, and every other spin.\n", "\n", "Find the reduced density matrices first:" ] }, { "cell_type": "code", "execution_count": 17, "id": "6eb61d9f-857a-4bf8-8460-d26b7f2c23f3", "metadata": {}, "outputs": [], "source": [ "target = (1, 2)\n", "\n", "rho_ab_ij = [[partial_trace(gs, dims=dims, keep=[target, (i, j)])\n", " for j in range(m)]\n", " for i in range(n)]" ] }, { "cell_type": "markdown", "id": "673c1f72-8ee6-4354-b95f-76a8e24fb181", "metadata": {}, "source": [ "Since one density matrix is just the spin itself, let's purify it when we come across it, meaning we'll find its total entanglement with it's environment:" ] }, { "cell_type": "code", "execution_count": 18, "id": "eb475288-3746-4e0a-8a84-0abfcfcfd125", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_2030217/2652293535.py:9: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.\n", " plt.colorbar()\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhQAAAGiCAYAAAC/AV8QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtWUlEQVR4nO3df7BkZXng8W/P6AwGBtxd1jDFjwqllhHKZBCcSNYaprYCsljK4lQWcQ0Zl1oUQRxdXEOMeXyMcRLKnYSYm6WirrJUagsTI5YlU4BEAoGBEQIqjLDEVQglPxTKcTQwyL29f5xzobvn3jt9zuk7ffre78c6NZzf73Qh/fTzvu/zdrrdLpIkSU2sGHcDJEnS5DOgkCRJjRlQSJKkxgwoJElSYwYUkiSpMQMKSZLUmAGFJElqzIBCkiQ1ZkAhSZIaM6CQJEmNvajJzZl5KfAJ4PKI2LLAdacA24DjgR8Al0XEFU3eLUmS2qN2hiIzXwecD3xrP9cdC1wL3AKcQBGA/Flmbqr7bkmS1C61MhSZeQjwV8B/BX5vP5e/G3i4J4Pxncw8CbgE+GKd90uSpHap2+UxBXw1Ir6WmfsLKE4Grh84dh1wXma+OCJ+PnhDZq4GVg8cXgU8GREujypJGqvMPIjie2kUno2IZ0b0rLGpHFBk5tuA1wKvG/KWI4DHB449Xr77cODROe65FIjBgzu27uLU/M3hGys1tHLNmnE3YWJN79kz7iZoGbth5q87i/XszDzokF94+umf/stLRvXIxzLz2EkPKioFFJl5NHA5cFrFv/hgVqEzz/FZWykGcc5aAzxS4X2SJC2WVT/9l5fw3v/8RVav2ifJXsneZ1/Mp/5q0xEU2Y7lE1AAJwIvA+7KzNljK4ENmXkRsDoipgfueYwiS9HrZcBzwJNzvSQi9gJ7Z/d73iVJUiusXvXzxgHFUlI1oLgReM3Asc8B9wN/PEcwAbADePPAsdOAO+caPyFJ0iSY7s4w3Z1p/IylolJAERF7gHt7j2XmzygGS95b7m8FjoyIc8tLrgAuysxtwKcpBmmeB5zTsO2SJI3NDF1m5u25H/4ZS8ViVMpcCxwzuxMR3wPOADYC9wAfAS6OCKeMSpIm1syI/rdUNKqUCRARGwf2N89xzd9TzAyRJElLUOOAQpKk5Wi622W626zLoun9beLiYJIk1TA7hqLpVtqZmbsy88Jx/p2aaHWGovxgL8TAR5K0tK2PiJ+MuxFNtDqgiIgpYCozDwV2j7s9kiTNmqHLtLM8ntfqgEKSpLZy2mg/uxIkSVJjZigkSarBWR79DCgkSaphptyaPmOpsMtDkiQ1ZoZCkqQapkcwy6Pp/W1iQCFJUg3T3WJr+oylotUBhYWtJEltNeIxFDszcwaYKmswTZxWBxQWtpIkLRNWypQkaTmaocM0ncbPWCoMKCRJqmGmW2xNn7FUODZBkiQ1ZoZCkqQapumwomGXRdMukzYxoJAkqQYDin52eUiSpMbMUEiSVMNMt8NMt+Esj4b3t0mrAwoLW0mS2souj36t/qKOiKmIOA5YP+62SJK0iHZm5q7yh/REanWGQpKktppmBSsa/i6ffuF+K2VKkrQcdUcwhqLrGApJkpa3aTp0HEPxvFaPoZAkSZPBDIUkSTVMd1fQ6TYcQ9Hw/jYxoJAkqYYZOsw0TPRXXW00MzcAHwROBNYCZ0XENQtc/3ngt+c4tSsiji+v2Qx8bo5rXhIRzwzbtqUTGkmStPQdDHwTuGjI699HEXjMbkcDTwF/PXDdTwauW1slmAAzFJIk1TKOQZkRsR3YDpCZw1y/G9g9u5+Z/xH4V+ybkehGxGOVGjOg1QGFlTIlSW014jEUawYChL0RsbfRw+d2HvC1iHho4PghmfkQsBK4B/hIRNxd5cGt/qK2UqYkaZl4hCKTMLtdOuoXZOZa4D8Anxk4dT+wGXgLcA7wDHBrZr6yyvNbnaGQJKmtikGZDRcHe+H+o4A9PacWIzuxGfgxcE3vwYi4Hbh9dj8zbwX+EXgvcPGwDzegkCSphhlW9JbOrv2M0p7FLL2dmR3gvwBXRcSzC10bETOZ+Q2gUoai1V0ekiRpJE4BXgF8dn8XlsHHOuDRKi+olKHIzAuAC4BfKg/dB3ysHHU61/Ubga/PcerVEXF/lXdLktQm090VcIALW2XmIRSBwaxjM3Md8FREPJyZW4EjI+LcgVvPA+6IiHvneGZQdHk8CBxK0c2xjmJSxNCqfhKPAL8DnFRufwd8OTOP3899r6J/fuuDFd8rSVKrzLBiJFtFJwF3lxvAtvKfP1burwWO6b0hMw8DNjF/duKlwF8C3wGuB44ENkTEzioNq5ShiIivDBz6cJm1eD1FtmI+T0TEj6u8S5KkNpvudqDhaqHTFe+PiJtg/pGgEbF5jmO7gV9Y4J73A++v1JA51B6UmZkrgd+kqNq1Yz+X352ZBwG7gI9HxFzdIL3PXg2s7jm0pm47JUnS4qscUGTmaygCiIOAn1LUEd81z+WPAucDd1EECL8F3JiZGyPi5gVecykQVdsmSdKBUszwaDiGYgnNjaiToXiAYrDGSyn6ZK7MzFPmCioi4oHy+lk7MvNo4BJgoYBiK0W/0Kw1FOM3JElqhRmaV8psurhYm1QOKMr5q/9U7t6Zma+jWHzkXUM+4nbgHft5x156inoMU69ckqQJtjMzZ4CpiJgad2PqGEVhqw794x325wQqzm2VJKltRtzlsX4xC1sdCFXrUHyCYpWzf6bohngbsBE4vTzfN/81M7cA36eYAbKKIjOxqdwkSZpYMyOY5THT8P42qZqh+EXgKop5rruBbwGnR8QN5fnB+a+rgE9SzGl9miKweFNEXNuk0ZIkqV2q1qE4bz/nNw/sXwZcVr1ZkiS1W7GwV9NBmcs3QyFJkijKZne73UbPmGk4S6RNls7fRJIkjU2rMxSZeSHF4iQGPpKkVim6KxoOyrTL48Ao5+JOZeahFINAJUlqBbs8+rU6oJAkqa2mWUGXhgHFCwl4C1tJkqTGlldhK0mSVLCwVT8DCkmSapgZQZdHdwnNOVg6fxNJkjQ2ZigkSaphpruCTsNZHl1neUiStLxN06HTsI5E1zoUB4aFrSRJmgytDigsbCVJaiu7PPq1OqCQJKmt7PLot3RCI0mSJtfOzNxVdvVPJDMUkiTVMNPt0GnYZdF9obCVlTIlSVqOprsr6DQbQuEYCkmSlrsunRGMgWg6CqM9lk5oJEmSxsYMhSRJNUx3V9BwKQ/orlgyX8RL5e8hSdIBNYrVRhvf3yKtDiislClJ0mRodUBhpUxJUltNj+S37tL5vdzqgEKSpLayy6Pf0gmNJEmaXFbKlCRpOZoZbZeHlTIlSVqOpu3y6GOXhyRJaswMhSRJNTgos58BhSRJNXS7K+g2rJTZdLXSNml1QGFhK0lSW02PYHGwqkuDZeYG4IPAicBa4KyIuGaB6zcCX5/j1Ksj4v6e6zYBfwC8HPgu8OGI+FKVtrX6izoipiLiOGD9uNsiSVILHAx8E7io4n2voghAZrcHZ09k5snA1cBVwK+Wf34hM3+tygtanaGQJKmtZrrQbTgGolOxyyQitgPbATKzyq1PRMSP5zm3BbghIraW+1sz85Ty+DnDvsCAQpKkGma6K+g2HETReSEgWTMQIOyNiL2NHt7v7sw8CNgFfDwiertBTgb+ZOD66ygCiqG1ustDkqRl4hGKNatmt0tH9NxHgfOBTcBbgQeAG8uxGLOOAB4fuO/x8vjQKmUoMvMC4ALgl8pD9wEfK1Mw891zCrANOB74AXBZRFxR5b2SJLXNDB0aTvLoHZR5FLCn59RIshMR8QBFEDFrR2YeDVwC3NxzfPCv0pnj2IKqdnk8AvwO8E/l/m8DX87MEyLivsGLM/NY4Frg08A7gH8H/EVm/jAivljx3ZIktcZ0tzOCaaPPBxR7DmDp7dspvpNnPca+2YiXsW/WYkGVAoqI+MrAoQ+XWYvXU2QrBr0beDgitpT738nMkygiIwMKSZIOvBMoukJm7QBOpX8cxWnAbVUeWntQZmauBH6TYgrLjnkuOxm4fuDYdcB5mfniiPj5PM9eDazuObSmbjslSVoMIx6UOZTMPAR4Rc+hYzNzHfBURDycmVuBIyPi3PL6LcD3KX70r6LITGwqt1mXAzdn5oeALwNnAr8BvKFK2yoHFJn5GooA4iDgpxRFNXbNc/l8Az1eBBxOf4TU61IgqrZNkqQDZYYRdHlUL4x1Ev2FqraVf14JbKaoMXFMz/lVwCeBI4GnKQKLN0XEtbMXRMRtmfk24OMUxa2+C5wdEXdUaVidDMUDwDrgpRQRzpWZecoCQcVcAz3mOt5rKy98SFBkKB6p3FIBsPKww8bdhInVOeLfjrsJE2vlY04ia2J69+5xN0EtFBE3wfxRSERsHti/DLhsiOf+DfA3TdpWOaCIiGd5YVDmnZn5OuB9wLvmuHy+gR7PAU8u8I699IxwrVi8Q5KkRdelw0zDZ6xoWLq7TUZR2KpD/3iHXjuANw8cOw24c77xE5IkTYKZboeZpvNGl+tqo5n5CYqSn/9M0Q3xNmAjcHp5vm8wCHAFcFFmbqOYOnoycB4VSnlKktRGM90VzDQdRLGEAoqqnZy/SLFoyAPAjcCvAadHxA3l+b7BIBHxPeAMiqDjHuAjwMXWoJAkqc/OzNxVrrI9karWoThvP+c3z3Hs74HXVmuWJEntNuIuj/UHsLDVonBxMEmSapgZwaDMBSZsTBzndUmSpMZanaEo+5IuxMBHktQyzvLo1+qAIiKmgKnMPJRiOVdJklrBgKKfv/wlSVJjrc5QSJLUVmYo+hlQSJJUgwFFP7s8JElSYwYUkiTV0GW2FkX9rSfBsbwqZUqSpMIoujw6L9xvpUxJkpajEQcUE88uD0mS1FirMxRWypQktZUZin6tDiislClJaisDin7+8pckSY21OkMhSVJbdbsdug0zDE3vbxMDCkmSaihqSTSzdOpk2uUhSZJGwIBCkqQaikGZzbeSlTIlSVqORjyGYuIrZZqhkCRJjbU6Q2FhK0lSW1mHol+rAwoLW0mSWmsEXR4YUEiStLyNIkPR9P42sStBkiQ1ZoZCkqQaujSvdLmEEhQGFJIk1TGKSplN728TuzwkSRo/C1tJkrQcWdiqnwGFJEk1OMujX6sDCgtbSZI0GVr9RR0RUxFxHLB+3G2RJKlXtzuabalodYZCkqS2GvEYionX6gyFJEmaDGYoJEmqYRwZiszcAHwQOBFYC5wVEdcscP1bgQuAdcBq4D7goxFxXc81m4HPzXH7SyLimWHbZoZCkqQailkezbeKDga+CVw05PUbgBuAMyiCkK8DX8nMEwau+wlFgPL8ViWYgIoZisy8FHgr8MvA08BtwIci4oEF7tlI8RcY9OqIuL/K+yVJaotRDKrsuX9NZvae2hsRewevj4jtwHaAgevnFBFbBg79bmaeCbwZuLu3KRHx2LDtnkvVLo9TgCngG+W9fwhcn5nHRcTP9nPvqygioFk/rPhuSZKWqkcG9hP46KhfkpkrgDXAUwOnDsnMh4CVwD3ARyLibiqoFFBExOkDDXsn8ARFGuXm/dz+RET8eJj3ZOZqir6eWWsqNFOSpEU34jEURwF7ek7tk50Ykf9G0W3yhZ5j9wObgW8DhwLvA27NzF+NiAeHfXDTQZmHlX8ORjpzuTszDwJ2AR+PiLm6QWZdCkTDtkmStGiKLo/KYyAGnvF8RLFnsUtvZ+Y5FFmPMyPiidnjEXE7cHvPdbcC/wi8F7h42OfXHpSZmR1gG/APEXHvApc+CpwPbKIYf/EAcGM5UnU+WymCldntqLrtlCRpucvMs4HPAv8pIr620LURMUMxtOGVVd7RJEPx58CvAG9Y6KJywGbvoM0dmXk0cAnzdJOUA1GeT/cMM/BEkqQDqVtuTZ+x2MrMxP8CzomIrw5xfYdimum3q7ynVkCRmZ8C3gJsiIjBgSTDuB14R513S5LUBsUYiqZdHlAlrMjMQ4BX9Bw6NjPXAU9FxMOZuRU4MiLOLa8/B/jfFOMibs/MI8r7no6I3eU1QfG9/CDFGIqLKQKKSkupV5022gE+BZwFbIyI71W5v8cJFF0hkiRpeCfRX4phW/nnlRQDK9cCx/ScfxfFd/1UuTFwPcBLgb8EjgB2U0wn3RARO6s0rGqGYgp4O3AmsKcn0tkdEU8DzBEdbQG+T1GdaxVFZmJTuUmSNJnG0OcRETcB86ZFImLzwP7GIZ75fuD91Vqyr6oBxQXlnzcNHH8n8Pnynwejo1XAJ4EjKYph3Qe8KSKurfhuSZLaYwRdHlTs8mizqnUo9vvJzREdXQZcVq1ZkiS124grZU481/KQJGn8dmbmrsysNBCyTVxtVJKkGkY3ywOA9Ytd2GqxtTqgKCO1CzGTIklqm26n2Bo9YzRNaYNWf1FHxFREHAesH3dbJEnS/FqdoZAkqa0clNnPgEKSpDompfb2AdLqLg9JkjQZzFBIklRDlxHM8hhRW9rAgEKSpDrs8uhjl4ckSWrMgEKSpBpmC1s13UpWylxMFraSJLXWaLs8Jr5SZqu/qC1sJUlqr86ItqWh1QGFJEmaDK3u8pAkqbWc5dHHgEKSpDoMKPrY5SFJkhozQyFJUh0uX97HgEKSpBpcbbSfXR6SJKmxVmcoLGwlSWqt0Q7K3JmZM8BUREw1fOpYtDqgKD/Uqcw8FNg97vZIkvS80Y6hsFKmJElSqzMUkiS1VadbbE2fsVQYUEiSVIeFrfoYUEiSVId1KPo4hkKSJDVmhkKSpDrs8uhjQCFJUh0GFH1aHVBY2EqSpMnQ6i/qiJiKiOOA9eNuiyRJfboj2go7M3NX+UN6IrU6QyFJUmtZKbNPqzMUkiRpMpihkCSpBitl9jOgkCSpriUUEDRVKaDIzEuBtwK/DDwN3AZ8KCIe2M99pwDbgOOBHwCXRcQVtVosSZJap2qG4hRgCvhGee8fAtdn5nER8bO5bsjMY4FrgU8D7wD+HfAXmfnDiPhi7ZZLkrTMZOYG4IPAicBa4KyIuGY/9+z3R31mbgL+AHg58F3gwxHxpSptqxRQRMTpAw14J/AExV/s5nluezfwcERsKfe/k5knAZcAcwYUmbkaWN1zaE2VdkqStNg6XWg4x6POGIqDgW8Cn2Oe79Bew/yoz8yTgauBjwBfAs4CvpCZb4iIO4ZtWNMxFIeVfz61wDUnA9cPHLsOOC8zXxwRP5/jnkuBaNg2SZIWT7dD45CiYkAREduB7QCZOcwtw/yo3wLcEBFby/2tZVZjC3DOsG2rHVBkZocihfIPEXHvApceATw+cOzx8t2HA4/Occ/W8tmz1gCPrPiFl9B9kSNgKjtm7bhbMLG+et3V427CxDrjjWePuwkTbcV35/qtpSVszUCAsDci9o7gucP8qD8Z+JM5rtlS5UVNMhR/DvwK8IYhrh2MAjrzHAeg/BCf/yCHjMIkSTpwRvH79oVnPDJwJoGPjuANw/yon++aI6q8qFZAkZmfAt4CbIiIwQ9h0GNzNOplwHPAk3XeL0nS2I02oDgK2NNzZhTZiX3fUpjrR/1c11T6G1adNtoBPkUxYGNjRHxviNt2AG8eOHYacOc84yckSVpu9ixS6e1hftTPd81g1mJBVTMUU8DbgTOBPZk524DdEfE0QGZuBY6MiHPLc1cAF2XmNopRpicD51FhoIckSW0zplkeVQ3zo34HcCr94yhOo6g1NbSqAcUF5Z83DRx/J/D58p/XAsfMnoiI72XmGRQNvZBiDuzF1qCQJE200XZ5DCUzDwFe0XPo2MxcBzwVEQ/X/FF/OXBzZn4I+DJF0uA3GG6M5POq1qHYbzAWEZvnOPb3wGurvEuSJO3jJODrPfuzMyKvBDZT40d9RNyWmW8DPk5R3Oq7wNlValCAa3lIklTPGDIUEXETC/S01P1RHxF/A/xNtdb0M6CQJKmGCRlDccC0OqDIzAspUjQrxt0WSZI0v1Z/UUfEVEQcB6wfd1skSerT7YxmK+zMzF3lD+mJ1OoMhSRJrTXaMRTrF6kOxQFjQCFJUg2OoejX6i4PSZI0GcxQSJJUxximjbaZAYUkSXWMoMtjKQUUdnlIkqTGzFBIklSHXR59Wh1QWNhKktRaBhR9Wv1FbWErSZImQ6sDCkmS2qrTHc1WslKmJElqbOIrZZqhkCRJjZmhkCSpDgdl9jGgkCSphg7N1+JoXBirRQwoJEmqY1TZhSUSVTiGQpIkNdbqDIWFrSRJrWWGok+rv6gtbCVJaqsR16GYeK0OKCRJ0mQwoJAkqY7uiLaClTIlSVqOOt3mwx967rdSpiRJkhkKSZLqWEIDKkfBgEKSpDoMKPrY5SFJkhprdYbCwlaSpLYa8aDMidfqL2oLW0mSWmu000YnXqszFJIktdYSCgZGodUZCkmSNBkMKCRJqmHEa3lYKVOSpGVptF0eE18ps3JAkZkbgA8CJwJrgbMi4poFrt8IfH2OU6+OiPurvl+SJLVPnQzFwcA3gc8BX6xw36uA3ujrhzXeLUlSKzhttF/lgCIitgPbATKzyq1PRMSPq75PkqRWcpZHnwM5huLuzDwI2AV8PCLm6gYBIDNXA6t7Dq1Z7MZJkqT6DkRA8ShwPnAXRZDwW8CNmbkxIm6e555LgTgAbZMkqZ4xZigy8z0U4xnXAvcBWyLilnmu/Tzw23Oc2hURx5fXbKYYyjDoJRHxzDBtWvSAIiIeAB7oObQjM48GLgHmCyi2Att69tcAjyxOCyVJqq7DeMZQZObZwJ8C7wFuBd4FbM/M4yLi4TlueR/wOz37L6IYC/nXA9f9hGK84/OGDSZmHzoOtwPvmO9kROwF9s7uVxyrIUnSUvYB4LMR8Zlyf0tmvhG4gCLD3ycidgO7Z/cz8z8C/4p9MxLdiHisbqPGFVCcQNEVIknSZBptl8eagR/Pe8sf130ycxVF2YY/Gjh1PfDrQ77rPOBrEfHQwPFDMvMhYCVwD/CRiLh7yGfWqkNxCPCKnkPHZuY64KmIeDgztwJHRsS55fVbgO9T9PGsoshMbCo3SZIm0oinjQ526yfw0TluOZziC//xgeOPA0fs732ZuRb4D8DbB07dD2wGvg0cStFNcmtm/mpEPLi/50K9DMVJ9Beqmh3rcGXZmLXAMT3nVwGfBI4EnqYILN4UEdfWeLckSe0w2gzFUcCenv19shP7eXtnyBZtBn4MXNN7MCJupxiOAEBm3gr8I/Be4OIhnlurDsVNLBCURcTmgf3LgMuqvkeSpGVkz5Clt38ETLNvNuJl7Ju16JOZHeC/AFdFxLMLXRsRM5n5DeCVQ7QJcHEwSZLq6zbcKioDgbuAUwdOnQrctp/bT6EYsvDZ/b2nDD7WUWG8o4uDSZJUw8BqofWeUe+2bcBVmXknsIOi1tMxwBUAg2MZe5wH3BER9w4+MDODosvjQYoxFBdTBBRDr37a6oCiXMb1QsykSJIEQERcnZn/Bvh9inGL9wJn9MzaGBzLSGYeRjEZ4n3zPPalwF9SdKXsBu4GNkTEzmHb1el221+MPDMPBXbfcfn3mH62/e1tm87Lj9n/RZrTtdddPe4mTKwz3nj2uJsw0brfnas+kYZ13U+vXLR1t2a/k77yozU81232mhd1urz58D0Ahy275cslSdJYuzxaya4ESZLUmBkKSZLqqDlTY6kyQyFJUg2zXR5Nt9LOzNxVTkaYSGYoJEkav/UOypQkaTmyy6OPAYUkSXUYUPRpdUBhYStJUls5bbRfq7+oI2IqIo4D1o+7LZIkaX6tzlBIktRadnn0MaCQJKmGTrc7gi6PpRORtLrLQ5IkTQYzFJIk1WGXRx8zFJIk1WClzH5mKCRJGj8rZUqStCzZ5dGn1QGFha0kSW1lYat+rf6itrCVJEmTodUZCkmSWssujz4GFJIk1WCXRz8DCkmS6jBD0afVYygkSdJkMKCQJKmGDiMobPXC4yxsJUnSstTtNu/y6D7/gIkvbGWGQpIkNWaGQpKkOkYwy2MpDepsdUBhpUxJUmuNYpbHEgooWv1FbaVMSZImQ6szFJIktVVnptgaPWM0TWkFAwpJkuqwy6NPq7s8JEnSZKicocjMDcAHgROBtcBZEXHNfu45BdgGHA/8ALgsIq6o3FpJklpiJGt5LPMMxcHAN4GLhrk4M48FrgVuAU4APgH8WWZuqvFuSZLaodsdzVZYfpUyI2I7sB0gM4e55d3AwxGxpdz/TmaeBFwCfLHq+yVJaoMRZygmvlLmgRiUeTJw/cCx64DzMvPFEfHzwRsyczWwuufQmkVsnyRJauhABBRHAI8PHHu8fPfhwKNz3HMpEPscfW4anms4R2cZ6jy5e9xNmFinXHD+uJswsQ558uFxN2GidZ97btxN0P44y6PPgZrlMfiRdeY5PmsrcFjPdtQitUuSpFoarzQ6itLdLXIgMhSPUWQper0MeA54cq4bImIvsHd2f8ixGpIkaUwORECxA3jzwLHTgDvnGj8hSdJEGO3y5ZVk5nsoSjisBe4DtkTELfNcuxH4+hynXh0R9/dctwn4A+DlwHeBD0fEl4ZtU506FIcAr+g5dGxmrgOeioiHM3MrcGREnFuevwK4KDO3AZ+mGKR5HnBO1XdLktQW46pDkZlnA38KvAe4FXgXsD0zj4uIhQYvvQronUnyw55nngxcDXwE+BJwFvCFzHxDRNwxTLvqZChOoj/S2Vb+eSWwmSJaOmb2ZER8LzPPAP6EYuXQHwAXR4RTRiVJqu4DwGcj4jPl/pbMfCNwAcWkhvk8ERE/nufcFuCGiNha7m8ti1JuYcgEQJ06FDexwHomEbF5jmN/D7y26rskSWqt0c7yWDMwXnBvOZ6wT2auoqhU/UcDp64Hfn0/b7s7Mw8CdgEfj4je5MDJFD/8e11HEVAMxbU8JEmqYcSzPB4Bdvds82UaDgdWMnc5hsEJELMeBc4HNgFvBR4AbiyX0pg1X4mH+Z65D1cblSRp/I4C9vTs75OdGDBXOYY58yUR8QBFEDFrR2YeTVGx+uY6z5xLqwOKsqb5hZhJkSS1zUwXGtdafP77es+Qpbd/BEwzdzmGwQzDQm4H3tGzP1+Jh6Gf2eov6oiYiojjgPXjboskSfvoNtwqiohngbuAUwdOnQrcVuFRJ9BfqXrHHM88rcozW52hkCSprca4fPk24KrMvJMiEDifYnblFQCD5RsycwvwfYp6FasoMhObym3W5cDNmfkh4MvAmcBvAG8YtlGtzlBIkqR+EXE1xeyL3wfuATYAZ0TEQ+UlfeUbKIKITwLfAm6hCBLeFBF/2/PM24C3Ae8sr9sMnD1sDQqATrdmla4DKTMPBXbf8T/+ielnXRysqhWH/5txN2Fi/fR1x+z/Is3pkG+4OFgTMz+ac2UCDem6Z/5q3vIGTc1+J928czXT081es3Jllw3r9wIc5vLlkiQtQ2Ps8mgluzwkSVJjZigkSapjtJUyJ54ZCkmSauh0uyPZSjszc1dZf2kitTpDYWErSdIysX7SB2W2+ovawlaSpNaaGdG2RLQ6QyFJUlsVXRbNn7FUtDpDIUmSJoMZCkmS6nCWRx8DCkmS6uh2RxBQLJ2IwoBCkqQarJTZzzEUkiSpMTMUkiTVYZdHHzMUkiTV0JkZzVayUuZislKmJGmZsFLmYrJSpiSptbrd0WxLRKszFJIktZZ1KPq0OkMhSZImgxkKSZJqcC2PfgYUkiTVMoJpo0uoz8MuD0mS1JgZCkmS6pgpt6bPWCIMKCRJqqFYy6NZl8VSWsuj1QGFha0kSa01ktLbAB0oKmXOAFMRMdW0aePQ6oCi/FCnMvNQYPe42yNJ0iKZ+EqZrQ4oJElqrdFmKCaeAYUkSXU4KLOPYxMkSVJjtTIUmfke4IPAWuA+YEtE3DLPtRuBr89x6tURcX+d90uSNG6jqZQ5mra0QeUMRWaeDfwp8IfACcAtwPbMPGY/t76KIgCZ3R6s+m5JklrD1Ub71MlQfAD4bER8ptzfkplvBC4ALl3gvici4sc13idJklquUkCRmauAE4E/Gjh1PfDr+7n97sw8CNgFfDwi5uoGmX3PamB1z6E1VdopSdKiG9ksj6WhapfH4cBK4PGB448DR8xzz6PA+cAm4K3AA8CNmblhgfdcSlF3YnZ7pGI7JUlaXKPt8tiZmbvKgo4Tqe600cGYqjPHMQAi4gGKIGLWjsw8GrgEuHme528FtvXsr8GgQpK0dC27wlY/AqbZNxvxMvbNWizkduAd852MiL3A3tn9zKzwaEmSDgDrUPSp1OUREc8CdwGnDpw6FbitwqNOoOgKkSRpIhXTRptvS0WdLo9twFWZeSewg2J8xDHAFQCZuRU4MiLOLfe3AN+nqFexiiIzsancJEmaTA7K7FO5DkVEXA1sAX4fuAfYAJwREQ+Vl6ylCDBmrQI+CXyLombFG4A3RcTf1m61JElqlVqDMiPiL4C/mOfc5oH9y4DL6rxHkqTWmuk6hqKHi4NJklTHGLs8Ki6B8VaK4pPrKGo83Qd8NCKu67lmM/C5OW5/SUQ8M0ybXBxMkqQJUmMJjA3ADcAZFMUpvw58JTNPGLjuJ/QvkbF22GACWp6hKAt8XIiBjySpdUaQoXjBmoESCXvLEgpzqbQERkRsGTj0u5l5JvBm4O6e492IeKxO46HlAUVETAFTmXkoRcVMSZLaYbRdHoPFGxP46ODlDZfAmH3GCoqCkU8NnDokMx+iqIh9D/CRiLibIfnLX5Kk8TsKOKxn2zrPdXWWwBj034CDgS/0HLsf2Ay8BTgHeAa4NTNfOeQz252hkCSptUY7y2NPxdLbQy+B0Sszz6HIfJwZEU/MHo+I2ymqWM9edyvwj8B7gYuHaZAZCkmS6ujOjGarpvYSGOVgzs8C/ykivrbQtRExA3wDGDpDYUAhSdKEqLsERpmZ+Dzw9oj46v7ek5kdimmmQy+TYZeHJEl1jK8ORdUlMM4B/jfwPuD2zJzNbjwdEbvLa4Kiy+NB4FCKbo51FDMth2KGQpKkOma6o9kqqrEExrsoEghTFBmH2e3ynmteCvwl8B2KGSNHAhsiYuew7TJDIUlSHWOslFlxCYyNQzzv/cD767WmYIZCkiQ11uoMhZUyJUmt1aXIUjR6RmckTWmDVn9RR8RURBwHrB93WyRJ6tPtjmYr7MzMXeUP6YnU6gyFJEnLxPqKha1ax4BCkqQ6ZmZqzdLof8bS6fIwoJAkqY7+LouazxhNU9qg1WMoJEnSZDBDIUlSHWYo+hhQSJJUR81Kl/3PGE1T2sAuD0mS1FirMxQWtpIktVW3O0O3YZdHdwkVtmp1QBERU8BUZh4K7B53eyRJel53BF0eS2gMhb/8JUmqw0qZfVqdoZAkaZmwUqYkScuSlTL7GFBIklSHdSj6OIZCkiQ1ZoZCkqQaujMzdBt2eXTt8pAkaZmzy6NPqwMKC1tJkjQZWh1QWNhKktRaruXRp9UBhSRJrdXtQrdhRNBdOgn4pfM3kSRpclkpU5Kk5ag70x3BLI/n77dSpiRJy1J3ZgRdHqNpShvUCigy8z3AB4G1wH3Aloi4ZYHrTwG2AccDPwAui4gr6rxbkqQ2GHGGYuJVHkORmWcDfwr8IXACcAuwPTOPmef6Y4Fry+tOAD4B/FlmbqrZZkmS1DJ1MhQfAD4bEZ8p97dk5huBC4BL57j+3cDDEbGl3P9OZp4EXAJ8ca4XZOZqYHXPoTUAK1d1cBxpdStevHQqsR1oK1csnV8PB9pK/71rpLPK/9Y1UZYb2BMRi/Z/4pWraNxlsXLVSJrSCpUCisxcBZwI/NHAqeuBX5/ntpPL872uA87LzBdHxM/nuOdSIAYPnvTel1dprjQC/zLuBkyuX/nX427BhPPza2g38G+BHy3Cs58FHlv/gV8+YkTPe6x85kSrmqE4HFgJPD5w/HFgvg/2iHmuf1H5vEfnuGcrxZiLWWuAR4CjgD3Vmvy8ncD6mvc2vX+c7570z27c9zf9/Cb57970/uX+756f3fj/3VuUL+mIeKbszh9VfuHZiHhmRM8am7qzPAaTPJ05ju3v+rmOAxARe4G9s/uZOfuPe+pOq8nMmSZTcprcP+Z3z/7jRH52476/6ec3yX/3pvcv93/3/Oxa8e/eoikDgIkPAkapaifdj4Bp9s1GvIx9sxCzHpvn+ueAJyu+v4mpMd4/znePwrjbP+77x/nuSb+/qXG3f5z/v29qkj+7UdyvA6zTrbhSWmbeAdwVEe/pObYL+HJE7DMoMzP/GHhzRBzXc+x/Ausi4uQh3zm7lsdhk17440Dzs2vGz68+P7v6/Oya8fMbjzpdHtuAqzLzTmAHcD5wDHAFQGZuBY6MiHPL668ALsrMbcCnKQZpngecU+Gde4GkpxtEQ/Oza8bPrz4/u/r87Jrx8xuDyhkKeL6w1X+nKGx1L/D+iLi5PPd54JciYmPP9acAf8ILha3+2MJWkiQtHbUCCkmSpF5WTpEkSY0ZUEiSpMYMKCRJUmMGFJIkqbG6lTIPmKpLpauQmRsoPrcTKT67syLimrE2akJk5qXAW4FfBp4GbgM+FBEPjLVhEyIzL6BYLPCXykP3AR+LiO1ja9SEKv9d/ARwec8Ci5pDZn6UfdeAejwiRrXehvaj1RmKqkulq8/BwDeBi8bdkAl0CkWVvtcDp1IE3tdn5sFjbdXkeAT4HeCkcvs74MuZefxYWzVhMvN1FHV+vjXutkyQ+yh+QM1urxlvc5aXtmcoqi6VrlL5a3A7HJi69ktJRJzeu5+Z7wSeoMj23DyWRk2QiPjKwKEPl1mL11P8B1/7kZmHAH8F/Ffg98bcnEnyXEQ8Nu5GLFetDShqLpUuLYbDyj+fGmsrJlBmrgR+kyJjtmPMzZkkU8BXI+JrmWlAMbxXZuYPKCpk3gH8bkT8vzG3adlobUBBvaXSpZHKzA5Fufl/iIh7x92eSZGZr6EIIA4CfkoxhmfXeFs1GTLzbcBrgdeNuy0T5g7gXOD/Ar9Ikdm5LTOPj4gDuRDlstXmgGJW1aXSpVH6c+BXgDeMuyET5gFgHfBSYBNwZWaeYlCxsMw8GrgcOK1cHltDGhj0++3M3AF8F/htih8FWmRtDijqLJUujUxmfgp4C7AhIh4Zd3smSUQ8C/xTuXtnOcDwfcC7xteqiXAixX/j7uoZ+7QS2JCZFwGrI2J6XI2bJBHxs8z8NvDKcbdluWjtLI/yP0h3UYyy73UqxTQ+aVFkZicz/5xi6ui/j4jvjbtNS0AHWD3uRkyAGylmJqzr2e6kGKC5zmBieJm5Gng18Oi427JctDlDAftZKl3zK0eJv6Ln0LGZuQ54KiIeHk+rJsYU8HbgTGBPZs5myXZHxNPja9ZkyMxPUMww+mdgDfA2YCNw+gK3CYiIPRQrOD8vM38GPOkYnoVl5ieBrwAPU2R5fg84FLhynO1aTlqboQCIiKuBLcDvA/cAG4AzIuKhMTZrUpwE3F1uUARndwMfG1uLJscFFDM7bqL4dTO7nT3GNk2SXwSuohhHcSPwa8DpEXHDWFulpe4o4P9Q/Hv3t8CzwOv9vjhwXL5ckiQ11uoMhSRJmgwGFJIkqTEDCkmS1JgBhSRJasyAQpIkNWZAIUmSGjOgkCRJjRlQSJKkxgwoJElSYwYUkiSpMQMKSZLU2P8HZBCIMppiq58AAAAASUVORK5CYII=", "image/svg+xml": [ "" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mi_ij = [[mutinf(rho_ab_ij[i][j] if (i, j) != target else\n", " purify(rho_ab_ij[i][j]))\n", " for j in range(m)]\n", " for i in range(n)]\n", "\n", "with plt.style.context(NEUTRAL_STYLE):\n", " plt.grid(False)\n", " plt.pcolor(mi_ij)\n", " plt.colorbar()" ] }, { "cell_type": "markdown", "id": "fd37a17f-5a17-46ed-9d0e-11005bb95a4f", "metadata": {}, "source": [ "One could also compute: ``concurrence``, ``logneg``, ``quantum_discord``, ``correlation`` ...\n", "\n", "For example we could set up the y-correlation function:" ] }, { "cell_type": "code", "execution_count": 19, "id": "5dbe9530-cdc5-426d-a590-d47b53fb1f23", "metadata": {}, "outputs": [], "source": [ "Sy = spin_operator('y')\n", "z_corr = correlation(None, Sy, Sy, 0, 1, dims=[2, 2], precomp_func=True)" ] }, { "cell_type": "markdown", "id": "b2570cd6-1c67-4a0a-9955-359116665457", "metadata": {}, "source": [ "And compute the correlations:" ] }, { "cell_type": "code", "execution_count": 27, "id": "8e70c7a0-798f-40d8-bce0-294a3c3d6f49", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_2030217/2620079929.py:9: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.\n", " plt.pcolor(cy_ij)\n", "/tmp/ipykernel_2030217/2620079929.py:10: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.\n", " plt.colorbar()\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAGiCAYAAAAiDFaYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1wElEQVR4nO3df5QlZXno+2+DzGh0ICaKQ0Ai9+giQDQgMAIrC8g5AVnkShaSXFAMwUNEETQTL4ijwYfHo4yOZlSkOaPoVQ/HHE3EH6ggGOMPIgNzUQjh5zUGRA7MECBOxqgDTO/7R9WGzZ69d/eu2ptd3f39ZFV6qup9q97eM9JPP+9bT021Wi0kSZKaZodJD0CSJKkXgxRJktRIBimSJKmRDFIkSVIjGaRIkqRGMkiRJEmNZJAiSZIaySBFkiQ1kkGKJElqJIMUSZLUSE+r0zkzVwEXAB+OiJUD2h0BrAX2A+4D1kTEujr3liRJC1vlTEpmHgycDtw8S7u9gCuAa4ADKIKaCzPzhKr3liRJC1+lTEpmPgv4DPA64C9naf4G4J6OTMvtmXkQcDZwWZX7S5Kkha/qdM808LWI+LvMnC1IORS4uuvYVcBpmblTRDza3SEzlwJLuw4vAR6KCF/bLEmaqMx8OsXPpVF4JCJ+OaJrLShDBymZeRLwUuDgOXZZDmzqOrapvPdzgPt79FkFRPfB3znxy3zlX0y+DOvsz7x20kOYt6567ZpJD2Heevkn3zrpIcxrHzj5k5Mewrz2iv/j5qlxXTszn/6sX/nFL37282eM6pIbM3MvA5XtDRWkZObzgQ8DRw/5YXZnP6b6HG9bTbHQtm0ZcO8Q95MkaVyW/Oznz+BNJ1/G0iXbTQYMZesjO/GRz5ywHLgpM2eA6YiYHskoF4BhMykHArsC38/M9rEdgcMz8yxgaURs6+qzkSKb0mlX4DHgoV43iYitwNb2fse9JElqhKVLHq0dpHRYERH/PqqLLRTDBinfBF7cdeyTwB3A+3oEKADrgVd0HTsauKHXehRJkuaDba0ZtrVmal9D/Q0VpETEFuCWzmOZ+R8UC1pvKfdXA7tHxCllk3XAWZm5FriEYiHtacCrao5dkqSJmaHFTN9VC3O/hvobR8XZ3YA92zsRcRdwLHAkcBNwHvDmiHAFrCRp3poZ0f+pv1oVZwEi4siu/VN7tPkOxRNBkiRJc1I7SJEkaTHa1mqxrVVvuqZu/4XOIEWSpApckzJ+jQ5SMvNM4Ex8W7MkSYtOo3/4R8R0ROwLrJj0WCRJ6jRDi201t45MyobMvK385VylRmdSJElqqhFP91jMrYdGZ1IkSdLiZSZFkqQKfLpn/AxSJEmqYKbc6l5D/TndI0mSGslMiiRJFbSf0Kl7DfVnkCJJUgXbWsVW9xrqr9FBisXcJElN5ZqU8Wv0D3+LuUmSFgmLufXQ6EyKJElNNcMU25iqfY2Sxdx6MEiRJKmCmVax1b2G+mv0dI8kSVq8zKRIklTBNqbYoeZ0T93pooXOIEWSpAoMUsbPIEWSpHkmM98InAPsBtwKrIyIawa0PwJYC+wH3AesiYh1HedPBT7Zo+szIuKXVe9bl2tSJEmqYKY1NZJtWJl5IvAh4D3AAcA1wJWZuWef9nsBV5TtDgAuAC7MzBO6mv47RfDx+NYVoAx131FodCbFYm6SpKYa8XTPsszsPLU1Irb26fYW4BMR8fFyf2Vmvhw4A1jVo/0bgHsiYmW5f3tmHgScDVzW0a4VERsHDHfY+9bW6CAlIqaB6czcGdg86fFIkjQm93btJ3B+d6PMXAIcCLy369TVwGF9rn1oeb7TVcBpmblTRDxaHntWZv4Y2BG4CTgvIm6scd/azFBIklTBNnYYyVbaA9ilY1vd57bPoQgiNnUd3wQs79NneZ/2TyuvB3AHcCpwHPAq4JfA9zLzRTXuW1ujMymSJDVVq+Kaku5rlLYMWXG2uwzcVI9js7V//HhEXAdc1z6Zmd8DfgC8CXhzjfvWYiZFkqQKtpVl8etuQ3oQ2Mb22Ytd2T7L0baxT/vHgId6dYiIGeD/BdqZlCr3rc0gRZKkeSIiHgG+DxzVdeoo4No+3db3aH80cEPHepQnycwpYH/g/hr3rc3pHkmSKtjW2oGpVr3f9bdV678WuDQzb6AIQE4H9gTWAWTmamD3iDilbL8OOCsz1wKXUCykPY1i7Qlln6CY7vkhsDPFFM/+FE/Yzum+42AmRZKkCmaYYoYdam7Dr2mJiM8BK4F3UjyFczhwbET8uGyyG0Xw0G5/F3AscGTZ/jzgzRHR+fjxrwIfA26neGJnd+DwiNgwxH1HzkyKJEnzTERcDFzc59ypPY59B3jpgOv9BfAXde47DgYpkiRVsI0ppnx3z1g1Okix4qwkqakmuCZl0Wj0pxMR0xGxL7Bi0mORJGmMNmTmbeUv5yo1OpMiSVJTFQtn603XdPRfMWQxt0XBIEWSpApmnlzWvvI11J+fjiRJaqShMimZeQbFK5lfUB66FXhXRFzZp/2RwLd6nNonIu4Y5t6SJDXJttYO4MLZsRr207kXeBtwULn9PfDlzNxvln57UxSXaW8/HPK+kiQ1Sv1Cbjs43TOLoTIpEfGVrkPvKLMrh1BkVfp5ICJ+OuTYJElqrG2tKaj5FuRtNfsvdJUXzmbmjsAfA8+kqOE/yI2Z+XTgNuDdEdFrCqjz2kuBpR2HllUdpyRJmp+GDlIy88UUQcnTgZ8Bx0fEbX2a30/xAqLvUwQdfwJ8MzOPjIjvDrjNKiCGHZskSU+V4smemmtSnO4ZqEom5U6KNyP+KnAC8OnMPKJXoBIRd5bt29Zn5vOBs4FBQcpqirctti2jWA8jSVIjzFC/4mzHmpQNmTkDTEfEdN2xLRRDBykR8Qjwz+XuDZl5MPDnwOvneInrgNfMco+twNb2fmYOO0xJkuYTi7n1MIpiblM8ef3IbA6gmAaSJGnecrpn/Iatk3IBcCXwE4opmJOAI4FjyvOrgd0j4pRyfyVwN8WTP0soMignlJskSfPWzAie7pnx6Z6Bhs2kPA+4lKLWyWbgZuCYiPhGeX43YM+O9kuADwC7A7+gCFb+ICKuqDNoSZK08A1bJ+W0Wc6f2rW/Blgz/LAkSWq24uWAdRfOmkkZxBcMSpJUwbbWDrRarVrXmLEs/kB+OpIkqZEanUnJzDOBMzGYkiQ1TDFVU3PhrNM9AzX6h39ETEfEvsCKSY9FkqRO21o7jGQrbcjM28pfzlVqdCZFkqSm2sYOtKi5JuWJXIHF3HpodCZFkiQtXmZSJEmqwGJu42eQIklSBTMjmO5pOaExkJ+OJElqJDMpkiRVMNPagamaxdxaFnMbyCBFkqQKtjHFVM06Jy3rpAzU6CDFYm6SJC1ejf7hbzE3SVJTzbR2GMlWsphbD43OpEiS1FQjnu6xmFsPjc6kSJKkxctMiiRJFcy0ppiq+XROy2JuAxmkSJJUwbbWDkzVewLZR5BnYZAiSVIFLaZG8Ahx3VUtC5tBiiRJ80xmvhE4B9gNuBVYGRHXDGh/BLAW2A+4D1gTEes6zr8OOAX47fLQ94G3R8SGjjbnA9F16U0Rsbz2N9SHeSZJkirY1tphJNuwMvNE4EPAe4ADgGuAKzNzzz7t9wKuKNsdAFwAXJiZJ3Q0OxL4X8DvAYcC9wBXZ+buXZe7lSIwam8vHvobGIKZFEmSKhjFW5Ar9n8L8ImI+Hi5vzIzXw6cAazq0f4NwD0RsbLcvz0zDwLOBi4DiIiTOzuUmZU/Av4L8D86Tj0WERurDLqKRgcpVpyVJC0SyzKzc39rRGztbpSZS4ADgfd2nboaOKzPtQ8tz3e6CjgtM3eKiEd79PkVYCfg4a7jL8rM+4CtwPUUU0L/0ue+tTX6h78VZyVJTbWNHUayle4FNndsvTIiAM8BdgQ2dR3fBPRbG7K8T/unldfr5b3A/wb+ruPY9RTrVl4OvK687rWZ+et9rlFbozMpkiQ11Yine/YAtnSc2S6L0t2za3+qx7HZ2vc6Tma+FXgVcGRE/LJ9PCKu7Gj2T5m5HvgR8KcUi3JHziBFkqTJ2zLHsvgPAtvYPmuyK9tnS9o29mn/GPBQ58HMPBt4O/D7EXHzoIFExH9k5j8BL5rDuCtp9HSPJElNNcMOI9mGERGPUDwefFTXqaOAa/t0W9+j/dHADZ3rUTLzHOA84JiIuGG2sWTmUmAf4P65jX54ZlIkSapg2+Se7lkLXJqZN1AEIKcDewLrADJzNbB7RJxStl8HnJWZa4FLKBbSnkYxpUPZ563AfwNeDdydme3My88i4mdlmw8AX6F4PHlX4C+BnYFPV/km5sJMiiRJ80hEfA5YCbwTuAk4HDg2In5cNtmNImhpt78LOJaiFspNFNmSN0fEZR2XfSOwBPg8RWakvZ3d0WYPiloqdwJfAB4BDum478iZSZEkqYIJ1kkhIi4GLu5z7tQex74DvHTA9V4wh3ueNPcRjoZBiiRJFbRaO9Cq+YLBum9RXugaHaRYzE2S1FTbRvCCQV8vOFijf/hbzE2StEhsyMzbyl/OVWp0JkWSpKaaaUGr5pqUqSemi1bMsU7KomKQIklSBTOtHWjVXJQyVXfh7QLX6OkeSZK0eA2VScnMMyheBf2C8tCtwLu66vl39zmCovDMfsB9wJqIWFdptJIkNcQMUwNfljMXLpwdbNhMyr3A24CDyu3vgS9n5n69GmfmXsAVwDXAAcAFwIWZeULlEUuS1ADbWlMj2dTfUJmUiPhK16F3lNmVQyiyKt3eANwTESvL/dsz8yCKCnaX9WgvSZIE1Fg4m5k7An8MPJPi3QG9HApc3XXsKuC0zNyp88VGXddeCiztOLSs6jglSRoHF86O39BBSma+mCIoeTrwM+D4iLitT/PlbP/q6E3lfZ9D/zcnrgJi2LFJkvRUmWGqfsVZ16QMVCWTciewP/CrwAnApzPziAGBSvdf4VSf451WUyy2bVsG3PuOv3kNj035Fzqsb/3X9096CPPWn5z8pkkPYd761mf8d1fH7/0/50x6CPPaK86b9AiGtiEzZ4DpiJie9GCaYuggJSIeAf653L0hMw8G/hx4fY/mGymyKZ12BR4DHhpwj63A1vZ+Zg47TEmSxqrFFDM1r7HDE5kUi7n1MIpiblM8ef1Ip/XAK7qOHQ3c0G89iiRJ88FMa4qZus8guyZloGHrpFwAXAn8hGIK5iTgSOCY8vxqYPeIOKXssg44KzPXApdQLKQ9DXjVKAYvSdKkzLR2YKbuohSDlIGGrZPyPOBSinUp3wReBhwTEd8oz+8G7NluHBF3AcdSBDI3AecBb44IHz+WJEkDDVsn5bRZzp/a49h3gJcONyxJkprN6Z7x8wWDkiRVMDOChbP4CPJAvmBQkiQ1UqMzKZl5JnAmBlOSpIZxumf8Gv3DPyKmI2JfYMWkxyJJUqciSKm/lTZk5m3lL+cqNTqTIknSImExtx4MUiRJqsDpnvEzSJEkqQKDlPFr9JoUSZK0eJlJkSSpghZFrZQ6zKMMZpAiSVIFo5jumao7XbTAGaRIklSBQcr4uSZFkiQ1UqMzKVaclSQ11YgzKRsycwaYjojpelddOBodpJR/UdOZuTOwedLjkSSpbcRBisXcejBDIUmSGqnRmRRJkpqq1ZqiVTOTUrf/QmeQIklSBTNMMVPzGtZJGczpHkmS1EhmUiRJqmCSdVIy843AOcBuwK3Ayoi4ZkD7I4C1wH7AfcCaiFjX1eYE4L8B/wn4EfCOiPhinfvWZSZFkqQKijUp9bdhZeaJwIeA9wAHANcAV2bmnn3a7wVcUbY7ALgAuLAMStptDgU+B1wK/E759W8y82VV7zsKZlIkSZpf3gJ8IiI+Xu6vzMyXA2cAq3q0fwNwT0SsLPdvz8yDgLOBy9rXAL4REavL/dVl9mUl8KqK962t0UGKxdwkSU014umeZZnZeWprRGztbp+ZS4ADgfd2nboaOKzPbQ4tz3e6CjgtM3eKiEfLNh/s0WZljfvW1ugf/hExHRH7AismPRZJkp5kFFM9T0z33EtRtLS99ctMPAfYEdjUdXwTsLxPn+V92j+tvN6gNu1rVrlvbY3OpEiS1FSjyKR09N8D2NJxarssSpfuO0/1ODZb++7jc7nmsPetxSBFkqTJ2zLHsvgPAtvYPnuxK9tnOdo29mn/GPDQLG3a16xy39oaPd0jSVJTtSgqxtbahrxnRDwCfB84quvUUcC1fbqt79H+aOCGcj3KoDbX1rhvbWZSJEmqYBQVZyv2Xwtcmpk3UAQXpwN7AusAMnM1sHtEnFK2XweclZlrgUsoFsmexhNP7QB8GPhuZp4LfBn4Q+D3gd+d633HwUyKJEnzSER8juKpm3cCNwGHA8dGxI/LJrtRBA/t9ncBxwJHlu3PA94cEZd1tLkWOAl4LXAzcCpwYkRcP8R9R85MiiRJFUzyBYMRcTFwcZ9zp/Y49h3gpbNc8/PA56vedxwMUiRJqmDET/eoh0YHKRZzkyRp8Wr0D3+LuUmSmqr2kz2tJ033bMjM28pfzlVqdCZFkqSmGvGalBVzrJOyqDQ6kyJJkhYvMymSJFUwyad7FguDFEmSKvDpnvEbKkjJzFXAK4HfAn5BUQr33Ii4c0CfI4Fv9Ti1T0TcMcz9JUlqiq6Fr5Wvof6GXZNyBDANHEJRr/9pwNWZ+cw59N2bogpee/vhkPeWJEmLyFCZlIg4pnM/M18LPAAcCHx3lu4PRMRP53KfzFwKLO04tGyIYUqSNHauSRm/umtSdim/PjyHtjdm5tOB24B3R0SvKaC2VUDUHJskSWNTTPdM1byGUcoglR9Bzswpijci/kNE3DKg6f0Ub0o8gWI9y53ANzPz8AF9VlMEQO1tj6rjlCRpHrCYWw91MikXAS/hya9x3k65qLZzYe36zHw+cDZ9pogiYiuwtb2fmTWGKUnS6LXKre41ShZz66FSJiUzPwIcB/xeRNxb4RLXAS+qcm9JkpqgWJNSf1N/wz6CPAV8BDgeODIi7qp43wMopoEkSZJ6Gna6Zxp4NfCHwJbMXF4e3xwRvwDIzNXA7hFxSrm/ErgbuBVYAryGYn3KCXUHL0nSxIx4vkfbGzZIOaP8+u2u468FPlX+eTdgz45zS4APALtTFIC7FfiDiLhiyHtLktQco5iuaT3+/9TDsHVSZv3biIhTu/bXAGuGG5YkSc1mxdnx8y3IkiSpkXzBoCRJFYzi6RwzKYM1Okgpi9qciRkfSVLTtKaKrdY1Hv/ThsycAaYjYrreRReORgcp5V/UdGbuDGye9HgkSRoTi7n10OggRZKkpnLh7PgZpEiSVIV1UsbOtR6SJKmRzKRIklRBixE83TOisSxUBimSJFXhdM/YOd0jSZIayUyKJEkVWMxt/BodpFjMTZLUWE73jF2jf/hHxHRE7AusmPRYJEl6sqkRbUBRcfa28pdzlRqdSZEkaZGw4mwPBimSJFXhdM/YGaRIklSFQcrYNXpNiiRJWrzMpEiSVEVrqthqXWM0Q+klM58NXAgcVx66HHhTRPx0QJ8pIIDTgWcD1wNnRsSt5flfAxI4Gng+8CDwJeC8iNjccZ27gd/suvz7IuJtw3wPZlIkSaqg/RbkutsY/TWwP3BMue0PXDpLn7cCbwHOAg4GNgLfyMxl5fnfKLezgRcDp5bX/kSPa70T2K1je/ew34CZFEmSFpjM3IcieDgkIq4vj70OWJ+Ze0fEnT36TAErgfdExBfKY38KbAJeDXw0Im4BTujo9qPMfAfwPzPzaRHxWMe5LRGxsc730eggxWJukqTGGu3C2WWZ2Xlma0RsrXHlQ4HN7QAFICKuy8zNwGHAdkEKsBewHLi6o8/WzPxO2eejfe61C/DvXQEKwLmZeR7wE+BvgfdHxCPDfBON/uFvMTdJUmO116TU3Qr3Aps7tlU1R7cceKDH8QfKc/36QJE56bSpX5/M/HXgPLYPYD4MnAT8HnARRYbm4tkG3a3RmRRJkhaJPYAtHfs9syiZeT7FwtZBDi6/9srzTPU53qn7fM8+mbkz8DXgNorFtI+LiA927N6cmf8GfD4zz42Ih2a5/+MMUiRJqmCqVWx1r1HaMseKsxcBn52lzd3AS4Dn9Tj3XLbPlLS1148sB+7vOL5rd59yIe3XgZ8Bx0fEo7OM6bry6wsBgxRJksZqAsXcIuJBisd+B8rM9cAumbkiIjaUx15GsX7k2j7d7qIIVI4Cbiz7LAGOAM7tuPbOwFUU2Z7jIuKXcxj6AeXX+we26mKQIklSFQ2ukxIRt2fm14FLMvP15eGPAV/tfLInM+8AVkXEFyOilZkfAt6emT8Efgi8Hfg5xePM7QzK1cCvAK8Bdi6DFoB/jYhtmXkocAjwLYr1NQcDHwQuj4h7hvk+DFIkSVqYTqYo5tZ+WudyivonnfamyK60rQGeQbHItV3M7eiIaK+XORB4Wfnnf+661l4UU01bgRMp1s4sBX4MXFJeeygGKZIkVdHwd/dExMMU2Y5Bbaa69lvA+eXWq/23KRbSDrrmDygyKbUZpEiSVEXDg5SFoNFBisXcJElavBr9w99ibpKkxmqNaCtsyMzbyl/OVWp0JkWSpMYa7dM9K+ZYJ2VRaXQmRZIkLV5mUiRJqmDEFWfVg0GKJElVGWSM1VBBSmauAl4J/BbwC4rSuud2Vq/r0+8IYC2wH3AfsCYi1lUasSRJWhSGXZNyBDBNUaTlKIog5+rMfGa/Dpm5F3AFcA1F7f4LgAsz84RKI5YkSYvCUJmUiDimcz8zXws8QFEm97t9ur0BuCciVpb7t2fmQcDZwGW9OmTmUopSum3LhhmnJEnjNtWapfTqHK+h/uo+3dOu9//wgDaH8sR7A9quAg7KzJ369FlF8VKi9nZvnUFKkjRy7UeQ627qq/LC2cycolhn8g8RccuApsuBTV3HNpX3fg69X9u8urx22zLg3s+d/GF2WvJo1SEvWqce/4ZJD2He+tGf9YujNRv/3dVz+WXvn/QQ5rm/mPQAhrUhM2eA6YiYnvRgmqLO0z0XAS8BfncObbsTWlN9jgMQEVsp3qIIQGZWGZ8kSeMziqkai7kNVClIycyPAMcBh0fEbFMxGymyKZ12BR4DHqpyf0mSJm60QYp6GPYR5CngI8DxwJERcdccuq0HXtF17Gjghohw7kaSJPU0bCZlGng18IfAlsxsZ0g2R8QvADJzNbB7RJxSnlsHnJWZa4FLKBbSnga8qu7gJUmaFJ/uGb9hn+45g+KJnm9TLHhtbyd2tNkN2LO9U2ZbjgWOBG4CzgPeHBE9Hz+WJGleGO1bkNXDsHVSZg0aI+LUHse+A7x0mHtJkqTFzXf3SJJUhQtnx84gRZKkClyTMn6NDlIy80zgTOpXxpUkqcks5tZDo4OU8i9qOjN3piiPL0lSM7SmqJ1LsZjbQI0OUiRJaizXpIydQYokSRW4JmX8XOshSZIayUyKJElVON0zdgYpkiRVMYLpHoOUwZzukSRJjWQmRZKkKpzuGbtGBykWc5MkNZZBytg1+od/RExHxL7AikmPRZKkMdqQmbeVv5yr1OhMiiRJTTXiOilWnO2h0ZkUSZK0eBmkSJKkRnK6R5KkKlw4O3YGKZIkVTBF/Xfv1C4GN0BmPhu4EDiuPHQ58KaI+OmAPlNAAKcDzwauB86MiFs72nwbOKKr6+ci4qQ69+7F6R5JkqpojWgbn78G9geOKbf9gUtn6fNW4C3AWcDBwEbgG5m5rKvdJcBuHdvrR3Dv7ZhJkSRpgcnMfSiCg0Mi4vry2OuA9Zm5d0Tc2aPPFLASeE9EfKE89qfAJuDVwEc7mv88IjaO6t79NDpIsZibJKmxRpUFKeZ8lmVm59GtEbG1xlUPBTa3gwSAiLguMzcDhwG9AoW9gOXA1R19tmbmd8o+nUHKyZn5GooA5kogI2JLjXv31OggJSKmgenM3BnYPOnxSJLUNpI6KU/88d6uUwmcX+PSy4EHehx/oDzXrw8UgUenTcBvdux/BriLYirot4HVwO8AR9W4d0+NDlIkSVok9gC2dOz3zKJk5vkUC1sHObj82ivXM9XneKfu80/qExGXdJy7JTN/CNyQmS+NiB/UvPeTGKRIklTFaBe9bpljxdmLgM/O0uZu4CXA83qcey7bZ0ra2mtMlgP3dxzfdUAfgB8AjwIvKv+8scK9ezJIkSSpghFP98xJRDwIPDhbu8xcD+ySmSsiYkN57GXALsC1fbq1p3COAm4s+yyheNz43AG32w/YiScCmyr37skgRZKkBSYibs/MrwOXZGb78eCPAV/tfLomM+8AVkXEFyOilZkfAt5eTuH8EHg78HOKR4rJzP8EnAxcQREs7Qv8FUVQ871h7j0XBimSJFXR/GqxJ1MUVGs/rXM5Rf2TTntTZDja1gDPAC7miWJuR3c8ufMI8F+APweeBfwE+BrF0z3bhrz3rAxSJEmqouFBSkQ8DLxmljZTXfstiqeKzu/T/idsX2220r3nwvojkiSpkRqdSbGYmySpqSaxcHaxafQP/4iYjoh9gRWTHoskSU8y2nf3bMjM28pfzlVqdCZFkqTGGu2alBVzrJOyqDQ6kyJJkhYvMymSJFXgmpTxM0iRJKmKhj+CvBAMHaRk5uHAOcCBwG7A8RHxpQHtjwS+1ePUPhFxx7D3lyRJi0OVTMozgX8EPglcNkS/vYHORUH/WuHekiQ1gtM94zd0kBIRVwJXAmTmMF0fiIifDns/SZIayemesXsq16TcmJlPB24D3h0RvaaAAMjMpcDSjkPLxj04SZLULE/FI8j3A6cDJwCvBO4EvlmubelnFbC5Y7t33IOUJGkoFnMbu7FnUsrXMne+mnl9Zj4fOBv4bp9uq4G1HfvLMFCRJDXIFCNdk2Ixtx4m9QjydQx4O2JEbAW2tveHXPsiSZIWgEkFKQdQTANJkjQ/uXB27KrUSXkW8MKOQ3tl5v7AwxFxT2auBnaPiFPK9iuBu4FbgSUUGZQTyk2SpHnJR5DHr0om5SCeXJytvXbk08CpFAXe9uw4vwT4ALA78AuKYOUPIuKKCveWJKkZzKSMXZU6Kd9mQPAXEad27a8B1gx7H0mStLj57h5Jkqqqm01xvmcggxRJkiqYahVbrWuMZigLVqODlLKozZk8NUXnJElSgzT6h39ETEfEvsCKSY9FkqQnseLs2DU6kyJJUlONeLrHirM9NDqTIkmSFi8zKZIkVfHk6RqNgUGKJEkV+HTP+DndI0mSGslMiiRJVTjdM3YGKZIkVWGQMnaNDlIs5iZJairXpIxfo3/4W8xNkrRIWMyth0ZnUiRJaqzRTvdYzK0HgxRJkiqYarVGMN3jopZBGj3dI0mSFi8zKZIkVeHTPWNnkCJJUgVNf7onM58NXAgcVx66HHhTRPx0QJ8pIIDTgWcD1wNnRsSt5fkXAHf16f5/RcTflu3uBn6z6/z7IuJtw3wPTvdIkrQw/TWwP3BMue0PXDpLn7cCbwHOAg4GNgLfyMxl5fmfALt1bQH8B3Bl17Xe2dXu3cN+A2ZSJEmqosHTPZm5D0VgckhEXF8eex2wPjP3jog7e/SZAlYC74mIL5TH/hTYBLwa+GhEbKMIXDr7HQ98LiJ+1nXJLRGxkRoaHaRYzE2S1FQjnu5Zlpmdp7ZGxNYalz4U2NwOUAAi4rrM3AwcBmwXpAB7AcuBqzv6bM3M75R9PtrdITMPpMjQ9Krvcm5mnkeRfflb4P0R8cgw30Sjf/hbzE2StEjcC2zu2FbVvN5y4IEexx8oz/XrA0XmpNOmAX1OA26PiGu7jn8YOAn4PeAiigzNxYOHvL1GZ1IkSWqs0U737AFs6djvmUXJzPMp1oAMcnD5tdfopvoc79R9vmefzHwGxTTQf+s+FxEf7Ni9OTP/Dfh8Zp4bEQ/Ncv/HGaRIklTBiKd7tsyx4uxFwGdnaXM38BLgeT3OPZftMyVt7fUjy4H7O47v2qfPHwG/AvyPWcYDcF359YWAQYokSWM1gYWzEfEg8OBs7TJzPbBLZq6IiA3lsZcBuwDdUzNtd1EEKkcBN5Z9lgBHAOf2aH8acHlE/Oschn5A+fX+ga26GKRIkrTARMTtmfl14JLMfH15+GPAVzuf7MnMO4BVEfHFiGhl5oeAt2fmD4EfAm8Hfk7xODMd/V4IHA4c233vzDwUOAT4FsX6moOBD1IENPcM830YpEiSVMEUI5juGWc1NziZophb+2mdyynqn3TamyK70rYGeAbFItd2MbejI2JLV7//Cvzvjmt32gqcSLF2ZinwY+CS8tpDMUiRJKmKVqv+dE9rfPNFEfEw8JpZ2kx17beA88ttUL+3U2RZep37AUUmpbZGP4IsSZIWLzMpkiRVMYKne5pasbYpGh2kWHFWktRYo3i654n+GzJzBpiOiOmaV10wGh2klH9R05m5M8UKYUmSFqIVc6yTsqg0OkiRJKmppmaKrdY1RjOUBcsgRZKkKkY73aMeXOshSZIaaehMSmYeDpwDHAjsBhwfEV+apc8RwFpgP+A+YE1ErBt6tJIkNcRI3t1jJmWgKpmUZwL/yPZV63rKzL2AK4BrKGr3XwBcmJknVLi3JEnN0GqNZlNfQ2dSIuJK4EqAzJxLlzcA90TEynL/9sw8CDgbuGzY+0uS1ARmUsbvqVg4eyjb1/a/CjgtM3eKiEe7O2TmUop6/23Lxjg+SZLUQE9FkLIc2NR1bFN57+fQ+7XNqyheTPQkr1n35zzmA1tDu/ny/z7pIcxbL/+N35n0EOatr9/3j5Mewrz2kr86Z9JDmNduef9TcBOLuY3dU/UIcvdf41Sf422rKRbati0D7h31oCRJqmrE0z0Wc+vhqQhSNlJkUzrtCjwGPNSrQ0RspXjVMzDntS+SJGkBeSqClPXAK7qOHQ3c0Gs9iiRJ80KrNYLpHlfODlKlTsqzgBd2HNorM/cHHo6IezJzNbB7RJxSnl8HnJWZa4FLKBbSnga8qtbIJUmaIJ/uGb8qdVIOAm4sNyjWjtwIvKvc3w3Ys904Iu4CjgWOBG4CzgPeHBE+fixJkvqqUifl2wx4J1JEnNrj2HeAlw57L0mSGst394ydLxiUJKkCp3vGzxcMSpKkRmp0JiUzzwTOxGBKktQ0My2YqXsRUymDNPqHf0RMR8S+wIpJj0WSpO20am5P2JCZt5W/nKvU6EyKJElNZcXZ8Wt0JkWSJC1eZlIkSarCirNjZ5AiSVIFPoI8fk73SJKkRjKTIklSFVacHTuDFEmSKphqtUYw3WOUMkijgxSLuUmStHg1+oe/xdwkSY01M6KtYDG3HhqdSZEkqalGPN1jMbceGp1JkSRJi5eZFEmSqvDpnrEzSJEkqQorzo6dQYokSRU0veJsZj4buBA4rjx0OfCmiPjpgD6vBF4PHAj8OnBARNzU1WYp8AHgVcAzgG8Cb4yIe+vcuxfXpEiStDD9NbA/cEy57Q9cOkufZwLfA942oM2HgOOBk4DfBZ4FfDUzd6x57+0YpEiSVEWrNZptDDJzH4rg4M8iYn1ErAdeB/yfmbl3v34RcWlEvAv4uz7X3QU4Dfi/I+LvIuJG4DXAi4Hfr3PvXgxSJEmqYGpmNFtpWWbu3LEtrTm8Q4HNEXF9+0BEXAdsBg6rcd0DgZ2Aqzuuex9wS8d1R3bvRq9JseKsJGmRuLdrP4Hza1xvOfBAj+MPlOfqXPeRiPi3ruObOq47sns3OkiJiGlgOjN3pojAJElqhtE+3bMHsKXjzNZezTPzfCBmuerB7av3ODfV53hd3dcdyb0bHaRIktRYo62TsmWOFWcvAj47S5u7gZcAz+tx7rkUWY+qNgJLMvPZXdmUXYFrO9qM5N4GKZIkzRMR8SDw4GztMnM9sEtmroiIDeWxlwG78EQwUcX3gUeBo4C/Ka+7G/DbwFvLNiO7t0GKJEkVjPjdPSMVEbdn5teBSzLz9eXhjwFfjYg72+0y8w5gVUR8sdz/NWBP4DfKJntnJsDGiNgYEZsz8xPAX2XmQ8DDFDVT/onyiaC53nsuXJAqSVIlo3j8eKwVZ0+mCB6uLrebgT/parM3RYaj7TjgRuBr5f5ny/03dLT5C+BLFJmU7wE/B14REduGvPeszKRIkrQARcTDFDVMBrWZ6tr/FPCpWfr8EnhTuVW+91wYpEiSVMVMudW9hvoySJEkqYLi3T31pmvG+e6ehaDRQYrF3CRJjTWSOilQlA9hQ2bOANNljTDR8CDFYm6SpEVixRzrpCwqjQ5SJElqrNFmUtSDQYokSVW4cHbsXOshSZIaqVImJTPfCJwD7AbcCqyMiGv6tD0S+FaPU/tExB1V7i9J0qSNpuLsaMayUA2dScnME4EPAe8BDgCuAa7MzD1n6bo3RVDT3n447L0lSWqMutVmH686q36qZFLeAnwiIj5e7q/MzJcDZwCrBvR7ICJ+WuF+kiRpERoqSMnMJcCBwHu7Tl0NHDZL9xsz8+nAbcC7I6LXFFD7PkuBpR2Hlg0zTkmSxm5kT/eon2Gne54D7Ahs6jq+CVjep8/9wOnACcArgTuBb2bm4QPus4qiLkp7u3fIcUqSNF6jne7ZkJm3lUVMVar6CHJ37DfV4xgA5WuZO1/NvD4znw+cDXy3z/VXA2s79pdhoCJJWrgs5tbDsEHKg8A2ts+a7Mr22ZVBrmPA2xEjYiuwtb2fmUNcWpKkp4B1UsZuqOmeiHgE+D5wVNepo4Brh7jUARTTQJIkzUvFI8j1N/VXZbpnLXBpZt4ArKdYb7InsA4gM1cDu0fEKeX+SuBuinoqSygyKCeUmyRJ85MLZ8du6DopEfE5YCXwTuAm4HDg2Ij4cdlkN4qgpW0J8AHgZoqaKr8L/EFEfKHyqCVJ0oJXaeFsRFwMXNzn3Kld+2uANVXuI0lSY820XJMyZr5gUJKkKpzuGTtfMChJkhqp0ZmUsqjNmRhMSZIaZwSZFA3U6B/+ETEdEfsCKyY9FkmSnsSKs2PX6EyKJEmLhBVnezBIkSSpCp/uGTuDFEmSqmjN+HTPmDV6TYokSVq8zKRIklSFdVLGziBFkqQqXJMydgYpkiRVYSZl7FyTIkmSGqnRmRQrzkqSGqtFZzG2iteYav9pQ2bOANMRMV3vogtHo4OU8i9qOjN3BjZPejySJD3uyRVjK17j8T9ZzK0HMxSSJKmRGp1JkSSpsWZmiid8al1javY2i5hBiiRJVYx2ukc9GKRIkrQAZeazgQuB48pDlwNvioifDujzSuD1wIHArwMHRMRNHed/DUjgaOD5wIPAl4DzImJzR7u7gd/suvz7IuJtw3wPBimSJFXR/EzKXwN7AMeU+x8DLgVeMaDPM4HvAX8LXNLj/G+U29nAbRSByLry2B91tX1n1zV+NtzwDVIkSapmpjWCNSmjGUq3zNyHIjg5JCKuL4+9DlifmXtHxJ29+kXEpWXbF/Q5fwtwQsehH2XmO4D/mZlPi4jHOs5tiYiNdb4PgxRJkiZvWWZ27m+NiK01rncosLkdoABExHWZuRk4DOgZpFS0C/DvXQEKwLmZeR7wE4rMzPsj4pFhLtzoIMVibpKkpmq1ZmjVnO5pPVHM7d6uUwmcX+PSy4EHehx/oDw3Epn568B5wEe7Tn0Y+AHwb8AKYDWwF/Bnw1y/0UGKxdwkSY3VGsF0zxPd9wC2dJzpmUXJzPOBmOWqB2939SdM9Tk+tPJn89co1qY8KQ0UER/s2L05M/8N+HxmnhsRD831Ho0OUiRJaqzRLpzdMseKsxcBn52lzd3AS4Dn9Tj3XGDTHEfXV2YuA75OsRj2+Ih4dJYu15VfXwgYpEiStNBExIMUj/0OlJnrgV0yc0VEbCiPvYxi/ci1dcZQZlCuosj2HBcRv5xDtwPKr/cPcy+DFEmSqmhwxdmIuD0zvw5ckpmvLw9/DPhq55M9mXkHsCoivlju/xqwJ8UjxQB7lwt6N0bExjKDcjXwK8BrgJ3LoAXgXyNiW2YeChwCfItiqcbBwAeByyPinmG+D4MUSZKqaH6dlJMpirldXe5fDpzV1WZviuxK23HAJzv221NL7YW8BwIvK4/9c9e19qKYatoKnEixdmYp8GOKeilrhv0GDFIkSVqAIuJhimzHoDZTXfufAj41oP23KRbfDrrmDygyKbUZpEiSVEFrZoZWzemeli8YHMggRZKkKpo/3TPvNTpIsZibJGmR2JCZM8B0WSNMNDxIsZibJKmxRvvunhVzrJOyqDQ6SJEkqbFaLWjVfENgy4mCQfx0JElSI5lJkSSpgtZMawRP97hydhCDFEmSqmjNjGC6ZzRDWagqBSmZ+UbgHGA34FZgZURcM6D9EcBaYD/gPmBNRKyrcm9JkprATMr4Db0mJTNPBD4EvIfihUHXAFdm5p592u8FXFG2OwC4ALgwM0+oOGZJkrQIVMmkvAX4RER8vNxfmZkvB84AVvVo/wbgnohYWe7fnpkHAWcDl/W6QWYupaj337asGKwRZxVbH3FWr6odl7i2vCr/3dXjf+/qKUtXbImIsX2QOy6h9nTNjktGMpQFa6j/imTmEoqXC72369TVwGF9uh3KEy83arsKOC0zd4qIR3v0WUXxYqInecWzfjbMcFVa88kTJz2EeevQXmG35mTNJ/ed9BDmteOftWXSQ5jvNgPPBR4cw7UfATaueMtvLR/R9TaW11SXYX/VeQ6wI7Cp6/gmoN9f1vI+7Z9WXu/+Hn1WU6xhaVsG3AvsAVT9X+4GYEXFvnX7T/Le8/2zm3T/up/ffP7e6/Zf7P/2/Owm/29vLD/4I+KX5VKGUeVBHomIX47oWgtK1Xxsd4Jrqsex2dr3Og5ARGyleNUzAJnZ/uOWqhX5MnOmTjW/Ov0nfO/2H+flZzfp/nU/v/n8vdftv9j/7fnZNeLf3tiUQYWBxZgNO+H+ILCN7bMmu7J9tqRtY5/2jwEPDXn/Ouq+C6FO/0neexQmPf5J95/kved7/7omPf5J/u++rvn82Y2ivxaAqdaQb3DMzOuB70fEGzuO3QZ8OSK2m8HPzPcBr4iIfTuO/Xdg/4g4dI73bL+7ZxffbTAcP7t6/Pyq87Orzs+uHj+/haPKdM9a4NLMvAFYD5wO7AmsA8jM1cDuEXFK2X4dcFZmrgUuoVhIexrwqiHuuRVIOqaANGd+dvX4+VXnZ1edn109fn4LxNCZFHi8mNtbKYq53QL8RUR8tzz3KeAFEXFkR/sjgA/yRDG391nMTZIkDVIpSJEkSRo3K1VJkqRGMkiRJEmNZJAiSZIaySBFkiQ1UuPfAFY+SXQOxZNEtwIrI+KayY6q+TLzcIrP7UCKz+74iPjSRAc1T2TmKuCVwG8BvwCuBc6NiDsnOrB5IjPPoHjh6AvKQ7cC74qIKyc2qHmq/Ld4AfDhjpe0qofMPJ/t3/m2KSJG9X4dTUCjMymZeSLwIeA9wAHANcCVmbnnJMc1TzwT+EfgrEkPZB46gqLa5SHAURTB/NWZ+cyJjmr+uBd4G3BQuf098OXM3G+io5pnMvNgijpUN096LPPIrRS/lLW3F092OKqr6ZmUtwCfiIiPl/srM/PlFL+l+X7aAcrfWq+Ep+Y9FgtJRBzTuZ+ZrwUeoMhKfXcig5pHIuIrXYfeUWZXDqH4IaJZZOazgM8ArwP+csLDmU8ei4iNkx6ERqexQUpmLqH4ofDerlNXA4c99SPSIrZL+fXhiY5iHsrMHYE/psjsrZ/wcOaTaeBrEfF3mWmQMncvysz7KCrNXg+8PSL+ZcJjUg2NDVKA5wA7sv2LCzex/QsLpbHIzCmKV0H8Q0TcMunxzBeZ+WKKoOTpwM8o1kTdNtlRzQ+ZeRLwUuDgSY9lnrkeOAX4/4DnUWSgrs3M/SLiqXyZrUaoyUFKW3dJ3Kkex6RxuQh4CfC7kx7IPHMnsD/wq8AJwKcz8wgDlcEy8/nAh4GjI+KXkx7PfNK1MPufMnM98CPgTyl+0dA81OQg5UFgG9tnTXZl++yKNHKZ+RHgOODwiLh30uOZTyLiEeCfy90bykWgfw68fnKjmhcOpPhv3Pc71pLtCByemWcBSyNi26QGN59ExH9k5j8BL5r0WFRdY5/uKf8j932Kpys6HUXxSKg0Fpk5lZkXUTyG/J8j4q5Jj2kBmAKWTnoQ88A3KZ5I2b9ju4FiEe3+Bihzl5lLgX2A+yc9FlXX5EwKFCm6SzPzBor57dOBPQHfoDyL8umAF3Yc2isz9wcejoh7JjOqeWMaeDXwh8CWzGxn8zZHxC8mN6z5ITMvoHiy7CfAMuAk4EjgmAHdBETEFoo3yz8uM/8DeMg1UYNl5geArwD3UGSj/hLYGfj0JMelehqbSQGIiM8BK4F3AjcBhwPHRsSPJzis+eIg4MZygyLguxF418RGNH+cQfFEz7cpfgtrbydOcEzzyfOASynWpXwTeBlwTER8Y6Kj0kK3B/C/KP7dfQF4BDjEnxfz21Sr5RpUSZLUPI3OpEiSpMXLIEWSJDWSQYokSWokgxRJktRIBimSJKmRDFIkSVIjGaRIkqRGMkiRJEmNZJAiSZIaySBFkiQ1kkGKJElqpP8fO1UG+uu4/rIAAAAASUVORK5CYII=", "image/svg+xml": [ "" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "cy_ij = [[z_corr(rho_ab_ij[i][j] if (i, j) != target else\n", " purify(rho_ab_ij[i][j]))\n", " for j in range(m)]\n", " for i in range(n)]\n", "\n", "\n", "with plt.style.context(NEUTRAL_STYLE):\n", " # plt.grid(False)\n", " plt.pcolor(cy_ij)\n", " plt.colorbar()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" }, "vscode": { "interpreter": { "hash": "39c10650315d977fb13868ea1402e99f3e10e9885c2c202e692ae90b8995050d" } } }, "nbformat": 4, "nbformat_minor": 2 }