{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Necessary dependencies for this Notebook:\n",
"# Needs tequila version>1.7.0\n",
"# e.g. pip install --upgrade git+https://github.com/aspuru-guzik-group/tequila.git@devel\n",
"# PySCF: pip install pyscf\n",
"# Recommended: pip install qulacs\n",
"import tequila as tq\n",
"import numpy\n",
"import time\n",
"\n",
"# convenience to avoid \"tq.gates\" in circuit construction\n",
"from tequila.circuit.gates import X, CNOT, Ry, QubitExcitation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Molecular Quantum Circuit Design: A Graph Based Approach\n",
"\n",
"\n",
"\n",
"- This notebook showcases the results of the article with the same title. \n",
"- You can find the original article [here](https://arxiv.org/abs/2207.12421). \n",
"\n",
"- If you are using this for your research, considering citing dependencies as well (see the \"Computational Details & Code Availability\" section of the article).\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Two-Electron Systems\n",
"\n",
"We will start by creating a function that generates a general version of the circuit displayed in Eq.(11) in the paper\n",
"\n",
"\n",
"\n",
"it generates a wavefunction for a spin-restriced electron-pair in $R_\\text{e}$ spatial orbitals. Spin-restriction in this context means, that electrons are always forced to form spin pairs. They move as quasi-particles sometimes referred to as *hardcore Bosons* - (see for example Section I.A. of [arxiv/2105.03836](https://arxiv.org/pdf/2105.03836.pdf)). \n",
"\n",
"The $N$-qubit wavefunction looks like\n",
"$$\n",
"U_2^N \\lvert 0\\dots0 \\rangle = \\sum_{k=0}^{R_\\text{e}-1} \\cos\\left(\\frac{\\theta_k}{2}\\right)\\prod_{l=0}^{k-1}\\left(\\sin\\left(\\frac{\\theta_l}{2}\\right)\\right)\\lvert2^{2k+1}+2^{2k}\\rangle\\nonumber = c_0 \\lvert 1100\\dots0 \\rangle + c_1 \\lvert 0011\\dots0 \\rangle + \\dots\n",
"$$\n",
"with $N=2R_\\text{e}$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def make_Upair(qubits: list, label=None):\n",
" # need even number of qubits\n",
" assert len(qubits) % 2 == 0\n",
" U = X(qubits[0])\n",
" c = None\n",
" for q in range(len(qubits)//2-1):\n",
" angle = (q,q+1,label)\n",
" U+= Ry(angle,qubits[2*q+2],control=c)\n",
" U+= CNOT(qubits[2*q+2], qubits[2*q])\n",
" c = 2*q+2\n",
" # the CNOT sequence that transforms from Hardcore Boson to Jordan-Wigner\n",
" U += sum([CNOT(qubits[2*q],qubits[2*q+1]) for q in range(len(qubits)//2)], tq.QCircuit())\n",
" return U"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's run this circuit for some angles and see how the wavefunction looks. It should look like described above:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"+0.8776|11000000> +0.4207|00110000> +0.2017|00001100> +0.1102|00000011> \n"
]
}
],
"source": [
"U = make_Upair(qubits=[0,1,2,3,4,5,6,7])\n",
"variables = {k:1.0 for k in U.extract_variables()}\n",
"wfn = tq.simulate(U, variables)\n",
"print(wfn)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's reproduce the left plot in Fig.2 of the paper where this particular circuit is tried out for Hydrogen-Molecules in an 6-31G basis. First we will create the molecule with tequila and then compute two VQEs, one with the standard Hartree-Fock orbitals and one with optimized orbitals for this circuit - note that the optimized orbitals are still in the same basis set, we are just optimizing the orbital representation within the given 6-31G basis. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -1.05802481296927\n",
"converged SCF energy = -1.12530860416883\n",
"converged SCF energy = -1.08435007825667\n",
"converged SCF energy = -1.02926301467656\n",
"converged SCF energy = -0.977519224346808\n",
"converged SCF energy = -0.932367146464438\n",
"converged SCF energy = -0.894183762983294\n",
"converged SCF energy = -0.862514623706141\n",
"converged SCF energy = -0.836584370668952\n",
"converged SCF energy = -0.815591772270656\n",
"runtime: 592.6103935241699s\n"
]
}
],
"source": [
"# loop over bond distances\n",
"# feel free to reduce points for faster runtime\n",
"# 10 points should take about 400s (reference: Dell XPS from 2019 ,intel i5, 4 threads)\n",
"# Time consuming part: orbital-optimization (evaluation of density matrices, takes a few iterations since no clever guess is used here)\n",
"npoints=10\n",
"geom = \"H 0.0 0.0 0.0\\nH 0.0 0.0 {}\"\n",
"start = time.time()\n",
"# collecting data\n",
"spa_hf=[]\n",
"spa_opt=[]\n",
"exact = []\n",
"x = []\n",
"\n",
"guess = None\n",
"for R in numpy.linspace(0.5, 3.0, npoints):\n",
" mol = tq.Molecule(geometry=geom.format(R), basis_set=\"6-31G\")\n",
" # no orbital optimization\n",
" H = mol.make_hamiltonian()\n",
" E = tq.ExpectationValue(H=H, U=U)\n",
" result = tq.minimize(E, silent=True)\n",
" delta=1.0\n",
" energy = result.energy\n",
" # orbital optimization\n",
" while(delta > 1.e-4):\n",
" opt = tq.chemistry.optimize_orbitals(circuit=U, molecule=mol, silent=True, initial_guess=guess)\n",
" guess = opt.mo_coeff\n",
" delta=abs(opt.energy-energy)\n",
" energy = opt.energy\n",
" # save data\n",
" spa_hf.append(result.energy)\n",
" spa_opt.append(opt.energy)\n",
" exact.append(mol.compute_energy(\"fci\"))\n",
" x.append(R)\n",
"end = time.time()\n",
"print(\"runtime: {}s\".format(end-start))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoLElEQVR4nO3de3hU1b3/8fc3kxsQIEC4BMIlgHJERdQIop7W9mhLVQRvFdSCgmC1eK8WtdVWEbXVHhRQRMErh4vWUrSIx/5OrfYAClhOVRDLTQkJ5EK4E5LMrN8fk8QQCASyk72T+byeJ8/MrNmz93dG+cyatdesMeccIiLS9MX5XYCIiDQMBb6ISIxQ4IuIxAgFvohIjFDgi4jEiHi/CziStLQ016NHD7/LEBFpVFauXFngnGtfvT3Qgd+jRw9WrFjhdxkiIo2KmX19uHYN6YiIxAgFvohIjFDgi4jEiEAGvpkNMbMZO3fu9LsUEZEmI5CB75x72zk3rnXr1n6XIiLSZAQy8EVExHsKfBGRGBHoefgitRWJOA4cKKO4uIwDB8Lllwdfr7ivpCRMXJwRCln5ZVyN1yu2q+n60R5/uH2Zmd8vV50453Au+po754hEDv9Xsc3h76vdY6pu15D3Offt86y4rPrcvbyM7vfQ+8aMOZ3u3VM9/W+nwJfjFg5HOHAgXBmsR7s8XAgfLpArLmsb4MXFZZSVRfx+OWrNDOLiom8EFeFf8R7g1+1jCWb9hEbDuOCCno078M1sGHAx0AqY6Zz774Y8fqyIRBz5+XvJzd1DTs5utm7dw549JQeF5qGXxxbYBw6UEQ578y8/Ls5ISgqRnBxPUlJ8+eW3t5OSQrRsmUhaWvNq94Vq3P5I+0pMDFWGVzjsCIcjldejl5FaXT/c42u6XtPjgYN6j37drvgkUvXP7NC2g+/nCPcd2+MOt331N8aGuK+ireJ69Uugxvvqchnd78HHqQ+1DnwzmwVcAuQ5506p0j4YeBoIAS865x6vaR/OuQXAAjNrAzwJKPCPQUWQ5+Tsrgzz3Nzd5OTsJidnT+X1bdv2HrHHWxGwFQFY02VKSvNDwrI2jzvSZfUATk6OJz5ep5JEGsKx9PBfBqYCr1Y0mFkImAZcCGQDy81sIdHwf6za40c75/LKr/+y/HFCdGgkP3/fQQH+baB/G+xbt+45bK86La05nTu3JD09hVNO6UB6ekr57ZaV7S1bJlWGrgJWJDbVOvCdcx+aWY9qzQOAdc65DQBmNhcY6px7jOingYNY9LPK48C7zrlPD3ccMxsHjAPo1q1bbcsLpHA4Ql7e3kN649XDfNu2wwd5+/bNK0P71FM7VIZ31TDv1CmFxMSQD89ORBqbuo7hdwE2V7mdDQw8wva3AhcArc2st3NuevUNnHMzzCwXGJKYmHhmHetrUNu372fZsmyWLt3MkiXZfPLJFvbsKTlku/btm1eGdr9+HQ7qiVe0K8hFxGsNetLWOfcM8EwttnsbeDsrK2ts/Vd1fCIRx+rV+SxdupmlS7NZsmQza9cWAtEx8n79OjJyZD9OPvngnnnHjgpyEfFHXQN/C9C1yu2M8rY6MbMhwJDevXvXdVee2bGjmI8/zmbp0ujfsmXZ7Np1AIB27ZoxaFBXRo48jUGDMjjrrC6kpCT6XLGIyMHqGvjLgRPMLJNo0A8HrqlrUX738CMRx1dfFbJkyebKHvzq1fk4F52SdcopHRgx4hQGDcpg0KCunHBC20b/ZRoRafqOZVrmHOB8IM3MsoGHnHMzzWw88B7RmTmznHNf1LWohu7h7959gI8/3lIZ7suWZVNUVAxAamoyZ5+dwdVXn8ygQV0ZMKALrVolNUhdIiJeMhfgr81lZWU5r3/i0DnHunXby3vv0eGZzz/Pq/wSTN++7TnnnGjPfdCgDPr0SSMuTr13EWk8zGylcy6rensgl1bwsoe/d28Jy5fnVM6cWbYsm4KCfQC0apXEwIFduOyy7zBoUAYDB2aQmppc52OKiARRk+zhL168jnfe+YolSzbzz39uq5zj3qdPO845p2vl2PtJJ6URCulLSCLStMRUD3/Bgi+ZPfszBgzown33ncegQV05++wM2rZt5m2hIiKNSJPs4e/cWUxKSqJ67yISkxpVD7+uWrfWOLyISHWB7ALrR8xFRLwXyMDXj5iLiHgvkIEvIiLeC2Tga0hHRMR7gQx8DemIiHgvkIEvIiLeU+CLiMQIBb6ISIxQ4IuIxIhABr5m6YiIeC+Qga9ZOiIi3gtk4IuIiPcU+CIiMUKBLyISIxos8M3sJDObbmZvmtnNDXVcERGJqlXgm9ksM8szs8+rtQ82s7Vmts7MJhxpH865Nc65nwI/Bs49/pJFROR41LaH/zIwuGqDmYWAacCPgL7ACDPra2anmtk71f46lD/mUuDPwCLPnoGIiNRKrX7xyjn3oZn1qNY8AFjnnNsAYGZzgaHOuceAS2rYz0JgoZn9Gfiv465aRESOWV1+4rALsLnK7WxgYE0bm9n5wOVAEkfo4ZvZOGAcQLdu3epQnoiIVNVgv2nrnPsA+KAW280ws1xgSGJi4pn1XZeISKyoyyydLUDXKrczyttERCSA6hL4y4ETzCzTzBKB4cBCL4rS0goiIt6r7bTMOcBSoI+ZZZvZGOdcGTAeeA9YA8x3zn3hRVFaPE1ExHvmnPO7hhplZWW5FStW+F2GiEijYmYrnXNZ1dsDubSCevgiIt4LZOBrDF9ExHuBDHz18EVEvBfIwFcPX0TEe4EMfBER8V4gA19DOiIi3gtk4GtIR0TEe4EMfBER8Z4CX0QkRgQy8DWGLyLivUAGvsbwRUS8F8jAFxER7ynwRURihAJfRCRGBDLwddJWRMR7gQx8nbQVEfFeIANfRES8p8AXEYkRCnwRkRihwBcRiRENGvhm1sLMVpjZJQ15XBERqWXgm9ksM8szs8+rtQ82s7Vmts7MJtRiV78A5h9PoSIiUjfxtdzuZWAq8GpFg5mFgGnAhUA2sNzMFgIh4LFqjx8NnAasBpLrVrKIiByPWgW+c+5DM+tRrXkAsM45twHAzOYCQ51zjwGHDNmY2flAC6AvsN/MFjnnIofZbhwwDqBbt261fiIiInJkte3hH04XYHOV29nAwJo2ds49AGBm1wMFhwv78u1mADMAsrKyXB3qExGRKuoS+MfFOffy0bYxsyHAkN69e9d/QSIiMaIus3S2AF2r3M4obxMRkQCqS+AvB04ws0wzSwSGAwu9KEpr6YiIeK+20zLnAEuBPmaWbWZjnHNlwHjgPWANMN8594UXRWm1TBER75lzwT0vmpWV5VasWOF3GSIijYqZrXTOZVVvD+TSCurhi4h4L5CBrzF8ERHvBTLw1cMXEfFeIANfPXwREe8FMvBFRMR7gQx8DemIiHgvkIGvIR0REe8FMvBFRMR7CnwRkRgRyMDXGL6IiPcCGfgawxcR8V4gA19ERLynwBcRiREKfBGRGBHIwNdJWxER7wUy8HXSVkTEe4EMfBER8Z4CX0QkRijwRURihAJfRCRGNFjgm9n5ZvaRmU03s/Mb6rgiIhJVq8A3s1lmlmdmn1drH2xma81snZlNOMpuHLAHSAayj69cERE5XvG13O5lYCrwakWDmYWAacCFRAN8uZktBELAY9UePxr4yDn3NzPrCPweuLZupYuIyLGoVeA75z40sx7VmgcA65xzGwDMbC4w1Dn3GHDJEXZXBCTVdKeZjQPGAXTr1q025YmISC3Utod/OF2AzVVuZwMDa9rYzC4HfgikEv20cFjOuRlmlgsMSUxMPLMO9YmISBUNdtLWOfeWc+4m59zVzrkPjrKtvmkrIuKxugT+FqBrldsZ5W11prV0RES8V5fAXw6cYGaZZpYIDAcWelOWiIh4rbbTMucAS4E+ZpZtZmOcc2XAeOA9YA0w3zn3hRdFaUhHRMR7tZ2lM6KG9kXAIk8rIjqkAwzp3bu317sWEYlZgVxaQT18ERHvBTLwRUTEe4EMfM3SERHxXiADX0M6IiLeC2Tgi4iI9wIZ+BrSERHxXiADX0M6IiLeC2Tgi4iI9wIZ+BrSERHxXiADX0M6IiLeC2Tgi4iI9xT4IiIxQoEvIhIjAhn4OmkrIuK9QAa+TtqKiHgvkIEvIiLeU+CLiMQIBb6ISIxokoFfmpNDaU6O32WIiARKgwW+mcWZ2aNmNsXMRtXnsbY9/gTrL7qYgunTiRw4UJ+HEhFpNGoV+GY2y8zyzOzzau2DzWytma0zswlH2c1QIAMoBbKPr9za6TjhF6R897vkT36aDZcMYfdf/1qfhxMRaRRq28N/GRhctcHMQsA04EdAX2CEmfU1s1PN7J1qfx2APsAS59xdwM3ePYVDJXTuTMbTk+k2ayaWkED2zbew+aafUvL11/V5WBGRQKtV4DvnPgS2V2seAKxzzm1wzpUAc4GhzrnPnHOXVPvLI9qrLyp/bLimY5nZODNbYWYr8vPzj/0ZVdHinHPoueCPdLj3XvatWMGGS4aQN3kykX376rRfEZHGqC5j+F2AzVVuZ5e31eQt4IdmNgX4sKaNnHMznHNZzrms9u3b16G8KEtMpN3oG+j57iJaXfQjCqc/z/qLL2HX4sU45+q8fxGRxqLBTto65/Y558Y45251zk070rb1sbRCQocOdH7iCbrPfp1Q69ZsueNOvrlhNAfWrfPsGCIiQVaXwN8CdK1yO6O8LdCan3kmmX94k44P/oriNWvYMOwytj3+BOE9e/wuTUSkXtUl8JcDJ5hZppklAsOBhV4UVd9r6VgoRNtrrqHX4ndJvfxytr/yCusH/4gdCxbgIpF6OaaIiN9qOy1zDrAU6GNm2WY2xjlXBowH3gPWAPOdc194UVRDrZYZ36YN6Q//hh7z55PQpTO5E+7j62uvo3j16no9rvinrKiIPR99RKSkxO9SRBqcBfnEZVZWlluxYkWDHMtFIuz84wLynnqKcFERqcOvpsPttxNKTW2Q40v9Ks3NpfCll9jxxpu4/ftJyMig/e230+rii7C4JvmFc4lhZrbSOZd1SHsQA9/MhgBDevfuPfZf//pXgx47vGsX+VOnUjT7vwi1bEn7O+4g9aorsVCoQesQbxzYuJHCmTPZ+aeFEInQ+pJLaHHeeRTOnMmBL78k6aST6HDXXbQ471zMzO9yRTzRqAK/QkP28KsrXvsV2yZOZN/y5ST37UvHX/2S5qef7kstcuyKV6+mYMYL7H7vPSwxkdQrrqDt6NEkZkRnDrtIhF1//jP5k5+mdMsWmp99Nh3uvptmp57ic+UiddeoAt/PHn5Vzjl2LVpE3hO/pSwvj9aXXUaHu+8iPi3Nt5rkyPatWEHB8zPY+9FHxKWk0GbECNqOGlnjf7NISQk75s6l4NnnCO/YQcsfDabD7beT2KNHwxYu4qFGFfgV/OzhVxXZu5eC6dMpfPkV4pKSaH/brbQZMQJLSPC7NCH6xrz3ww8peH4G+z/9lFDbtrQdOZI214wg1KpVrfYR3rOH7bNmUfjSy7jSUtr8+CrSbr6ZeA++/CfS0BT4HjiwYSPbJk1i79//TtIJven4wC9pcfZAv8uKWS4cZtfixRS+8CIHvvyS+PR02o0eTeqVVxDXrNlx7bMsP5/8Z59lxxtvRr+lff0o2o4eTSglxePqRepPowr8oAzpHI5zjj3/8z9sm/QYpVu20OqiH9Hh3ntJ6NTJ79JiRqSkhJ0LFlA4cyalX39DYmYm7caOpfUlF2OJiZ4co2TTJvKefprd7y4m1KYNaTffTOrwq4nzaP8i9alRBX6FoPXwq4oUF1P44kwKX3gB4uJIu/lm2l4/SoFQjyJ791I0/w22v/QSZXl5JJ98Mu3GjaPlBf9Rb7Oo9n/2GXlPPsW+jz/WVE5pNBT49aQkewt5TzzO7vf/QmL37nR84H5SvvMdv8tqUsI7drD99dkUvfYa4Z07aT5gAO3GjaPFuec0yFRK5xx7//6/5D31lKZySqOgwK9nez76O9sefZSSTZtI+f736XjfBBK7dj36A6VGpdvy2P7yyxTNm4fbt4+U732PduPG+jY99vBTOe+i2amn+lKPSE0aVeAHeQz/SFxJCdtffZX8Z5+DsjLa3XgjbUeNJFRPawI1VSXffEPhizPZ+cc/4sJhWl10Ee3GjiW5z4l+lwZUTOWcR8FzzxEuKtJUTgmcRhX4FRpTD7+q0m3byPvt79j15z8DEJeSQkJ6OvHpnUhI70xCejoJndPL2zqT0LGDpngCxWvXUjjjBXa9+y4WCtH68stpN2Y0id26+V3aYUWncr5E4csv40pKSL3qStrfcoumcorvFPg+2L9qFftWfkppbi6lW3Mpy8mlNDeXcFHRwRuaEd+hAwmdOhHfOf0wbwrphFJTm+x48b5P/0Hh88+z529/I655c1JHDKftqFEkdOjgd2m1UpafT8Fzz1E0/w0sIYF2N1yvqZziKwV+gET276c0dyuluTmU5eZSmpNL6dby2+VvCq7aao7WrFn0TaDyk0L5G0PFm0KnTsQlJfn0jI5dxYnQwuefZ9+KFYRSU2kz8ie0veaaRrtgnaZySlAo8BsR5xzh7dujbwS5OZRt3Vp+veIvh3B+wSGPC6WlRd8IOnUioXP0k0H000In4lJaYokJWEIClpgYvaz4a8Aphi4cZvf7f6FgxvMcWL2G+I4daTf6BlKvuoq45s0brI76tP+zz8l76in2LVtWPpXzNlpdfLGmckqDaVSB31hP2jakSEnJt28EW3O//aSQ++2fq+2PtcfHfxv+Vd8MEhOwhMQq12t4wzjoMd9ej0tMhCrXw7t2UzR7NiWbNpHQvRtpY8fS6tJLm2QP2DnH3v9dEp3KuWaNpnJKg2pUgV8hVnv4XnDOEdm5MzpUlJNLZP8+XEkprrQEV1pafr38duX1UlxJSd2ul5TAEf6fSjrpJNLGjaXlD34QE0tOu0iEXYveJX/yZEqzszWVUxqEAl8ajAuHD/umgHMkdOsWkz1cV1JC0bz5FDz7bHQq5+DBdLjzDhK7d/e7NGmCFPgiAVB1KqeZ0XPhn0jo0sXvsqSJqSnwdRZJpAGFUlJof9ut9FzwR3CO3F89SJA7XdK0KPBFfJDYrRvtf343e5csYecf/uB3ORIjGizwzezfzWy6mb1oZksa6rgiQdVm+HCan3UW2x5/gtJt2/wuR2JAfG02MrNZwCVAnnPulCrtg4GngRDwonPu8Zr24Zz7CPjIzIYBy+tStEhTYHFxpE98hA1Dh7H1wYfImP5cTJ7QbiqccxCJ4MJhKCvDRSLRy3A42lbt0pWVRbcvC0MkHL0Ml+HCEVy4jGannOL5Oly1CnzgZWAq8GpFg5mFgGnAhUA2sNzMFhIN/8eqPX60cy6v/Po1wJg61CzSZCR2706HO+9g22OPs2vhQloPHep3STFr75IlFM56ifD27dGwDpfhysK4SBjKqoX1YQKccNjTerq/9irNzzrL033WKvCdcx+aWY9qzQOAdc65DQBmNhcY6px7jOingUOYWTdgp3Nud03HMrNxwDiAbgFdNEvES22uu45di99j66THaHHOOVp8rYHt/+wz8n7/e/YtXUZ8ejrJffpAKISFQlh8COKi1yvaiA9hceX3heKxUFz5ffEQisMq2+LLH1feVrGvqpcVxwlV21d8PEknnOD5c61tD/9wugCbq9zOBo72A69jgJeOtIFzboaZ5QJDEhMTz6xDfSKNgoVCpD86kY3DLmPrww/T5ZlnNLTTAA6sX0/+5KfZ/f77hNq0oeP995E6fHiT/OZ3hQadpeOce8g5d9QTts65t51z41prHXmJEUk9e9L+tlvZ/f5f2L14sd/lNGmlOTnk3P8AG4Zcyt4lS0i7dTy93n+ftiNHNumwh7r18LcAVX/SKaO8rc6qrKXjxe5EGoW211/Prvf+m60PP0LzgQOJb9vW75KalLLt2yl8/nmK/msOmNH2Jz+h3U3jYup1rksPfzlwgpllmlkiMBxY6EVR6uFLLLL4eNIfnUh4zx62TZzodzlNRnjPXvKnTmP9BRey/bXXaXXpEHotfpeO902IqbCHWga+mc0BlgJ9zCzbzMY458qA8cB7wBpgvnPuCy+KMrMhZjZj586dXuxOpNFIPvFE2t9yM7sWvcuu99/3u5xGLXLgANtfeYX1F15IwdSptDjvPHq+vZDOjz5KQufOfpfnC62lIxIwrrSUjT++mrL8fHq983aj/UEYv7iyMnb+aSH506ZSlpNLi3MG0f7OO2NqhdJGtZaOevgSyywhgc6THiW8YwfbHqv+lRapiXOOXe+/z4ahw8h94AHi27aj26yZdJs1K6bC/kgCGfgaw5dYl1z+uwE7/7SQ3R984Hc5gbd32TI2XT2cLbfeBs7R5Zmn6fHGfFqcc47fpQVKIId09ItXItFfNdt0xRWEd+2m5ztvE2rZ0u+SAmf/Z5+T/5//yd4lS4hPT6f9+J/ReuhQLL4uExAbv0Y1pKMevgjEJSaSPmkSZfn55P32t36XEygHNmwk+/Y72HTVVRSvXk2HCb+g1+J3Sb3iipgP+yPRKyMSYM1OPZV2o2+g8MWZtBw8mJRzz/W7JF+Vbt1KwbRp7Hjrj8QlJZH2s5/R9obrCaWk+F1ao6AhHZGAixw4wMZhlxE5UEzPhW8TSmnhd0kNrqyoiMIZL1A0ezY4R5trRtDupptibh59bWlIR6SRiktKIv3RRynL3Ur+75/yu5wGFdm7l/xnn41+aeqVV2h18cXlX5q6T2F/HDSkI9IIND/jdNqO/AnbX3mVlj8cTIuBA/wuqV5FSkrYMXceBdOnE96+nZYXXkD7228nScut1ImGdEQaicj+/WwYOgyco+efFhDXvLnfJXnOhcPsfPttCp6ZQmlODs0HDqTDXXfS7LTT/C6tUdGQjkgjF9esGekTH6F082byn37a73I85Zxj9//7f2wcNozcCfcRatOGrjNfpNvLLynsPaQhHZFGpMWAAbS5ZgTbX32Nlj/8Ic3POMPvkuqsrKCALXfdzb5PPiGxRw+6TJ5Myx/+QL8JUA8COaRT4XBr6ZSWlpKdnU1xcbFPVTUeycnJZGRkkJCQ4Hcp4qHwnr1svPRSLDGRzAV/JC452e+Sjlvxl1+y+ZZbCG8vouOECaReqXn0XqhpSKfRvbLZ2dm0bNmSHj16qAdwBM45CgsLyc7OJjMz0+9yxEOhlBZ0euRhNo+5kYKpU+nw85/7XdJx2fX+++Tc+wtCrVvTffbrNDv5ZL9LavICOYZ/JMXFxbRr105hfxRmRrt27fRJqIlKOfdcUq+6ksJZL7H/n//0u5xj4pyjYPp0ttx6G0knnkCP+fMU9g0kkIF/tNUyFfa1o9epaetw773Et29P7gMPECkp8bucWokUF5Nz98/Jn/w0rS4dQvdXXyWhQwe/y4oZgQx8zdIRObpQy5akP/wbDvxrHQXPPed3OUdVui2Pr38ykl3vvkv7u++i8xNPEJeU5HdZMSWQgd8YPProo5x88sn069eP/v378/HHH3P++efTp08fTjvtNM4991zWrl1buf2wYcM4++yza9zfpk2bOOWUUw5q+/Wvf82TTz4JwPXXX09mZib9+/enf//+PPPMM/XzxKRRSfnud2k9dCiFM16gePVqv8up0f7PPmfTVVdRsn49GdOmkjZ2rD6B+qDRnbQNgqVLl/LOO+/w6aefkpSUREFBASXlH6lnz55NVlYWM2bM4J577mHhwoXs2LGDlStXkpKSwoYNG+jZs+dxHfd3v/sdV155pZdPRZqAjvdNYM+S/yXngV+SOX8eFrBZWbsWLSLnvvuJT0uj+5w5JPc50e+SYlajDvw77ljMqlVbPd1n//6dmDx58BG3yc3NJS0tjaTyj6NpaWmHbPOd73yHyZMnA/DWW28xZMgQOnbsyNy5c7n//vs9rVliWyg1lfSHHiJ7/K0UvPAC7W+5xe+SAHCRCAVTp1Lw7HM0O/NMMqY8o/VvfNZgQzpm1s3MFpjZLDOb0FDHrQ8/+MEP2Lx5MyeeeCK33HILf/vb3w7Z5u233+bU8p9VmzNnDiNGjGDEiBHMmTOnxv2uX7++csimf//+TJ8+/aD777nnnsr7PvvsM2+flDRqLS+4gFYXXUTBc9Mp/uorv8shsm8fW26/g4Jnn6P1FZfT/aVZCvsAqFUP38xmAZcAec65U6q0DwaeBkLAi865x4+wm1OBN51zr5vZvDrUXOloPfH6kpKSwsqVK/noo4/461//ytVXX83jj0ef+rXXXkuzZs3o0aMHU6ZMYdu2bfzrX//ivPPOw8xISEjg888/P2S8HqBXr16sWrWq8vavf/3rg+7XkI4cScdfPsDeZcvIvf8Besyd49sXmEpzctj8s/EcWLuWjvdNoM3IkRqvD4ja/h/xMjAVeLWiwcxCwDTgQiAbWG5mC4mGf/VfXh4NLAPeNLPRwGt1K9t/oVCI888/n/PPP59TTz2VV155Bfh2DL/ClClTKCoqqvzy065du5gzZw6XXnopN910EwAPP/ww/fr1a/gnIU1KfNu2dPrVL9ly510UvvQSaWPHNngN+/7xD7JvvQ1XXEzX56eT8u//3uA1SM1qNaTjnPsQ2F6teQCwzjm3wTlXAswFhjrnPnPOXVLtLw+4AXjIOfd94GIvn0RDW7t2LVVX8Vy1ahXdu3c/7LZz5sxh8eLFbNq0iU2bNrFy5Urmzp3LwIEDWbVqFatWreLSSy9tqNKliWs5eDAtL7yQgilTObBhQ4Mee8eCBXwzchRxLZrTY95chX0A1WUMvwuwucrt7PK2miwGbjOz6cCmmjYys3FmtsLMVuTn59ehvPqzZ88eRo0aRd++fenXrx+rV68+ZPgFolMtv/7664OmY2ZmZtK6dWs+/vjjBqxYYoWZ0enBXxHXrBm59z+AC4fr/ZguHCbvySfJnXAfzc48k8x580jq1avejyvHrtaLp5lZD+CdijF8M7sSGOycu7H89k+Agc658XUu6gjr4a9Zs4aTTjqproeIGXq9YtPOhQvJufcXdLxvAm1Hjaq344T37CXn5z9nzwcfkDpiOJ3uvz9w00JjUX2sh78F6FrldkZ5m4j4rNWQIaR897vk/edkSr7+ul6OUZKdzdcjRrDno4/o+OCvSH/oIYV9wNUl8JcDJ5hZppklAsOBhV4UpaUVROrGzOj08G+whARyH/glLhLxdP/7li9n05VXUZqXR7cXX6DtNdd4un+pH7UKfDObAywF+phZtpmNcc6VAeOB94A1wHzn3BdeFHW0xdNE5OgSOnak44RfsG/FCormzvVsv0VvvMHXN4wm1KYNmfPm0mLQIM/2LfWrVtMynXMjamhfBCzytKLoft8G3s7Kymr4eWUiTUjryy9n16J3yXvyKVK+810SM440r+LIXFkZ2377W4pefY0W551Hl98/RahVKw+rlfoWyMXT1MMX8YaZkf7Iwxiw9cFfcby/cBfetYvNN/2Uoldfo+2oUXSd/pzCvhEKZOBrDF/EOwmdO9Ph3nvYu2QpO95885gff2DjRjZdPZy9n3xC+sRH6HjfBP0MYSMVyMBvDD18r5dHBti5cycjR46kd+/e9OrVi5EjR1Kb12DSpEl1fj7StKX++Mc0HzCAvCd+S+nW2i84uHfJEjZdPZzwjh10nzWTVC3t0agFMvCD3sOvujzyP//5T/7yl7/QtWt0hurs2bP5v//7P0aNGsU999wDULk88s6dO9lwhG8/jhkzhp49e7Ju3TrWr19PZmYmN95441HrUeDL0VhcHOkTH8GFw+Q+9FCthna2z57NN2PHkdCxIz3emE/zs85qgEqlPjXqz2VbJ03iwJovPd1n0kn/RqejLF9cH8sjr1u3jpUrVzJv3rfryj344IP07t2b9evXs3nzZh588EFatmzJunXr+N73vsezzz7L/fffz/79++nfvz8nn3wys2fPrsOzl6YssVs3Otx5B9smPcbOP/2J1GHDDrudKy1l66RJ7Jgzl5TvfY/Ov/sdoZQWDVus1ItA9vCDPqRTH8sjr169mv79+xMKhSrbQqEQ/fv354svorNdP/nkE6ZMmcLq1atZv349b731Fo8//jjNmjVj1apVCns5qjbXXUezM85g26THKM3LO+T+sqIivrlxLDvmzKXd2LFkTJ2isG9CAtnDr+20zKP1xOtLfS2PfDQDBgyo/LWsESNG8Pe//13LJcsxsbg40h+dyMZhl7H1Nw+TMXVK5dLFB9avZ/PNt1CWm0vnJx6n9dChPlcrXgtk4DcGXi+P3LdvX1atWkUkEiEuLvrBKxKJsGrVKvr27Ut2dvYha4prjXE5HkmZmbS/7VbyfvckuxYtovXFF7Pnww/ZctfdWHIy3V97lWb9+/tdptSDQA7pBF19LI/cu3dvTj/9dCZOnFj52IkTJ3LGGWfQu3dvIDqks3HjRiKRCPPmzeO8884DICEhgdLS0np8xtLUtL3+epL79WPbIxPJnzaNzT+9mYRuXcl8Y77CvgkLZOAHfQy/vpZHnjlzJl999RW9evWiV69efPXVV8ycObPy/rPOOovx48dz0kknkZmZyWWXXQbAuHHj6NevH9dee633T1aaJAuF6PzoRCJ791IwZSotL7iAHq+/TkJ6ut+lST2q9fLIfsjKynIrVqw4qC1Wl/v94IMPePLJJ3nnnXeO6XGx+npJ7exatIiyggLaXHcdFhfI/p8ch5qWR9YYvkgMa3XRRX6XIA1Igd9IVJwgFhE5Xo3yM1yQh6GCRK+TiFQVyMA/0knb5ORkCgsLFWZH4ZyjsLCQ5ORkv0sRkYBodCdtS0tLyc7Opri42KeqGo/k5GQyMjJI0M/OicSUJnPSNiEhofJLTCIiUnuBHNIRERHvKfBFRGKEAl9EJEYE+qStmeUDXx/nw9OAAg/Laez0enxLr8XB9HocrCm8Ht2dc+2rNwY68OvCzFYc7ix1rNLr8S29FgfT63Gwpvx6aEhHRCRGKPBFRGJEUw78GX4XEDB6Pb6l1+Jgej0O1mRfjyY7hi8iIgdryj18ERGpQoEvIhIjmmTgm9lgM1trZuvMbILf9fjFzLqa2V/NbLWZfWFmt/tdUxCYWcjM/mFmx/bzYU2QmaWa2Ztm9qWZrTGzQX7X5Bczu7P838nnZjbHzJrcUrNNLvDNLARMA34E9AVGmFlff6vyTRlwt3OuL3A28LMYfi2quh1Y43cRAfE0sNg592/AacTo62JmXYDbgCzn3ClACBjub1Xea3KBDwwA1jnnNjjnSoC5wFCfa/KFcy7XOfdp+fXdRP8xd/G3Kn+ZWQZwMfCi37X4zcxaA98BZgI450qcczt8Lcpf8UAzM4sHmgM5PtfjuaYY+F2AzVVuZxPjIQdgZj2A04GPfS7Fb5OBe4GIz3UEQSaQD7xUPsT1opm18LsoPzjntgBPAt8AucBO59x/+1uV95pi4Es1ZpYC/AG4wzm3y+96/GJmlwB5zrmVftcSEPHAGcBzzrnTgb1ATJ7zMrM2REcCMoHOQAszu87fqrzXFAN/C9C1yu2M8raYZGYJRMN+tnPuLb/r8dm5wKVmtonoUN/3zex1f0vyVTaQ7Zyr+NT3JtE3gFh0AbDROZfvnCsF3gLO8bkmzzXFwF8OnGBmmWaWSPTEy0Kfa/KFmRnR8dk1zrnf+12P35xz9znnMpxzPYj+f/E/zrkm14urLefcVmCzmfUpb/oPYLWPJfnpG+BsM2te/u/mP2iCJ7Ab3U8cHo1zrszMxgPvET3TPss594XPZfnlXOAnwGdmtqq87X7n3CL/SpKAuRWYXd452gDc4HM9vnDOfWxmbwKfEp3d9g+a4BILWlpBRCRGNMUhHREROQwFvohIjFDgi4jECAW+iEiMUOCLiMQIBb6ISIxQ4IuIxIj/D7itf+tT4dSgAAAAAElFTkSuQmCC\n",
"text/plain": [
"