{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Water Molecule with Active Space (CPU vs. GPU)\n",
    "\n",
    "#### A- Classical simulation as a reference: CCSD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install the relevant packages.\n",
    "!pip install pyscf==2.6.2\n",
    "!pip install openfermionpyscf==0.5\n",
    "!pip install matplotlib==3.8.4\n",
    "!pip install openfermion==1.6.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -75.9839755372789\n",
      "Total number of electrons =  10\n",
      "Total number of orbitals =  13\n",
      "ncore occupied orbitals =  3\n",
      "E(CCSD) = -75.98508980454675  E_corr = -0.001114267267875617\n",
      "Total CCSD energy for active space Hamiltonian =  -75.98508980454675 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "from pyscf import gto, scf, mcscf, cc\n",
    "\n",
    "geometry ='O 0.1173 0.0 0.0; H -0.4691 0.7570 0.0; H -0.4691 -0.7570 0.0'\n",
    "mol=gto.M(\n",
    "    atom = geometry,\n",
    "    spin = 0,\n",
    "    charge = 0,\n",
    "    basis = '631g',\n",
    ")\n",
    "\n",
    "myhf = scf.RHF(mol)\n",
    "myhf.max_cycle=100\n",
    "myhf.kernel()\n",
    "nelec = mol.nelectron\n",
    "print('Total number of electrons = ', nelec)\n",
    "norb = myhf.mo_coeff.shape[1]\n",
    "print('Total number of orbitals = ', norb)\n",
    "\n",
    "norb_cas, nele_cas = (4,4)\n",
    "mycasci = mcscf.CASCI(myhf, norb_cas, nele_cas)\n",
    "print('ncore occupied orbitals = ', mycasci.ncore)\n",
    "\n",
    "frozen = []\n",
    "frozen += [y for y in range(0,mycasci.ncore)]\n",
    "frozen += [y for y in range(mycasci.ncore+norb_cas, len(mycasci.mo_coeff))]\n",
    "mycc = cc.CCSD(myhf,frozen=frozen).run()\n",
    "print('Total CCSD energy for active space Hamiltonian = ', mycc.e_tot, '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B- VQE-UCCSD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[\u001b[38;2;255;000;000mwarning\u001b[0m] Target \u001b[38;2;000;000;255mnvidia-fp64\u001b[0m: \u001b[38;2;000;000;255mThis target is deprecating. Please use the 'nvidia' target with option 'fp64'.\u001b[0m\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_23147/4290935201.py:35: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  spin_ham = cudaq.SpinOperator(qubit_hamiltonian)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "UCCSD-VQE energy =   -75.98415928173183\n",
      "Pyscf-CCSD energy =  -75.98508980454675\n",
      "Total number of qubits =  8\n",
      "Total number of parameters =  26\n",
      "Total number of terms in the spin hamiltonian =  105\n",
      "Total elapsed time (s) =  28.929891359000067\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAHHCAYAAACvJxw8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABcd0lEQVR4nO3dd3xUZd4+/utMzaT3RmihhS4CZmkWQIqsAZVnLSBFXB8V0UdYd0FFVH4Kj4jr7n53LY8xYEERBQELSLXQQZCaICWEQBIgZVKnn98fkzmZk0wmIUwymcn1fr3mZeacMyf3DIlcfO7PuY8giqIIIiIiInJJ4e0BEBEREbVmDEtEREREbjAsEREREbnBsERERETkBsMSERERkRsMS0RERERuMCwRERERucGwREREROQGwxIRERGRGwxLRERERG4wLBGR30pLS0NgYCDKysrqPWbKlCnQaDQoLCwEAFRUVGDx4sXo168fAgMDERYWhhEjRuDjjz+Gq7tDCYJQ7+Pxxx9vtvdGRC1H5e0BEBE1lylTpmDjxo1Yt24dpk2bVmd/ZWUl1q9fj3HjxiEqKgoFBQUYNWoUTp06hQceeABPPfUUDAYDvvrqK0ybNg2bNm3Cxx9/DIVC/u/MO++80+X5u3fv3mzvjYhaDsMSEfmttLQ0hISEYNWqVS7DzPr161FRUYEpU6YAAKZPn45Tp05h3bp1SEtLk457+umn8dxzz+HNN9/ETTfdhOeee052nu7du2Pq1KnN+2aIyGs4DUdEfkun0+Hee+/Ftm3bcOXKlTr7V61ahZCQEKSlpWHv3r3YvHkzZsyYIQtKDkuWLEG3bt2wdOlSVFVVtcTwiaiVYFgiIr82ZcoUWCwWfPHFF7LtRUVF2Lx5M+655x7odDps3LgRAFxWoABApVLhoYceQlFREXbv3i3bZzAYcO3atToPk8nUPG+KiFoUwxIR+bWRI0ciISEBq1atkm1fs2YNzGazNAV38uRJAED//v3rPZdjn+NYh/T0dMTExNR5rF271pNvhYi8hD1LROTXlEolHnjgAfz9739HdnY2OnXqBMA+BRcXF4dRo0YBgHTFXEhISL3ncuyrfXXdxIkT8dRTT9U5vm/fvp54C0TkZQxLROT3pkyZgr///e9YtWoVnn/+eeTm5uLnn3/G008/DaVSCUAehMLDw12exxGSYmNjZduTkpIwevTo5nsDRORVnIYjIr83cOBApKSk4LPPPgMAfPbZZxBFUZqCA4BevXoBAI4ePVrveRz7kpOTm3G0RNTaMCwRUZswZcoUHD9+HEePHsWqVavQrVs3DB48WNp/9913AwA++ugjl6+3Wq3S1N2tt97aImMmotaBYYmI2gRHFemll17CkSNHZFUlAPjDH/6AMWPGICMjA998802d17/wwgs4ffo0/vrXv0KlYgcDUVsiiK7W7yci8kPDhg2TLvv//fff0bVrV9n+goICjBw5EpmZmXjooYcwYsQIGI1GrF27Fjt37sTUqVPx0UcfQRAE6TWCINS7gndcXBzuvPPO5n1TRNTsGJaIqM34z3/+g9mzZ+OWW27Bvn37XB5TXl6Ot956C1988QXOnj0Lg8EAAFi4cCFeffXVOsc7B6fabrvtNuzcudMjYyci72FYIiJy49KlSxg6dCgsFgv27NmDDh06eHtIRNTC2LNERORGu3btsGnTJhgMBowfPx7FxcXeHhIRtTBWloiIiIjcYGWJiIiIyA2GJSIiIiI3GJaIiIiI3GBYIiIiInKDy9B6gM1mw+XLlxESEuJ2zRUiIiJqPURRRFlZGRITE6FQ1F8/YljygMuXL6N9+/beHgYRERE1wcWLF5GUlFTvfoYlDwgJCQFg/7BDQ0O9PBoiIiJqjNLSUrRv3176e7w+DEse4Jh6Cw0NZVgiIiLyMQ210LDBm4iIiMgNhiUiIiIiNxiWiIiIiNxgWCIiIiJyg2GJiIiIyA2GJSIiIiI3GJaIiIiI3GBYIiIiInKDYYmIiIjIDYYlIiIiIjcYloiIiIjcYFgiIiIicoNhyQdUmawQRdHbwyAiImqTGJZauROX9ej50iYs+T7T20MhIiJqkxiWWrml1SHp/Z/OeXkkREREbRPDUitXbrR4ewhERERtGsNSK1duYFgiIiLyJoalVo6VJSIiIu9iWGrlWFkiIiLyLoalVq6MlSUiIiKvYlgiIiIicoNhqRWz2bgQJRERkbcxLLVipQaz7DlX8SYiImp5DEut2LVyk+y52cqwRERE1NIYllqxogp5WDJZbV4aCRERUdvFsNSKFZYbZc+NZquXRkJERNR2MSy1YoWsLBEREXkdw1IrVlirZ8loZlgiIiJqaQxLrVhhhXwajpUlIiKilsew1IrVnoZjZYmIiKjlMSy1YrUbvE1WNngTERG1NIalVqz20gGsLBEREbU8lbcHQPXrlRCKALUSvxeUo8pshZE9S0RERC2OlaVW7O0HBmDDU8PRKzEUACtLRERE3sCw5AM0SvsfE6+GIyIiankMSz5Aq7b/MXEFbyIiopbHsOQDWFkiIiLyHoYlH6BVKwGwZ4mIiMgbGJZ8ACtLRERE3sOw5ANqepYYloiIiFoaw5IPqKksscGbiIiopTEs+QBWloiIiLyHYckHaNmzRERE5DUMSz6AV8MRERF5D8OSD+DVcERERN7DsOQDpJ4lCxu8iYiIWppPhKWdO3dCEASXjwMHDgAAsrOzXe7fu3ev23Pn5ORgwoQJCAwMRGxsLJ577jlYLJaWeFuNplVVV5YsrCwRERG1NJW3B9AYQ4cORV5enmzbwoULsW3bNgwaNEi2fevWrejdu7f0PCoqqt7zWq1WTJgwAfHx8di9ezfy8vIwbdo0qNVqvP766559EzdAo3JUlhiWiIiIWppPhCWNRoP4+Hjpudlsxvr16zFnzhwIgiA7NioqSnasOz/88ANOnjyJrVu3Ii4uDjfddBMWL16Mv/3tb3j55Zeh0Wg8+j6aSqtigzcREZG3+MQ0XG0bNmxAYWEhZs6cWWdfWloaYmNjMXz4cGzYsMHtefbs2YO+ffsiLi5O2jZ27FiUlpbixIkTHh93UzkavI1s8CYiImpxPlFZqi09PR1jx45FUlKStC04OBjLly/HsGHDoFAo8NVXX2HSpEn4+uuvkZaW5vI8+fn5sqAEQHqen59f7/c3Go0wGo3S89LS0ht5Ow2qWZSSDd5EREQtzauVpfnz59fbuO14ZGZmyl6Tm5uLzZs3Y9asWbLt0dHRmDt3LlJTUzF48GAsXboUU6dOxbJlyzw+7iVLliAsLEx6tG/f3uPfwxmXDiAiIvIer1aW5s2bhxkzZrg9Jjk5WfY8IyMDUVFR9VaLnKWmpmLLli317o+Pj8f+/ftl2woKCqR99VmwYAHmzp0rPS8tLW3WwMRFKYmIiLzHq2EpJiYGMTExjT5eFEVkZGRIV6w15MiRI0hISKh3/5AhQ/Daa6/hypUriI2NBQBs2bIFoaGh6NWrV72v02q10Gq1jR73jWJliYiIyHt8qmdp+/btOH/+PB599NE6+1auXAmNRoMBAwYAANauXYsPP/wQH3zwgXTMunXrsGDBAmlqb8yYMejVqxcefvhhvPHGG8jPz8eLL76I2bNnt2gYagh7loiIiLzHp8JSeno6hg4dipSUFJf7Fy9ejAsXLkClUiElJQWrV6/G5MmTpf16vR5ZWVnSc6VSiW+++QZPPPEEhgwZgqCgIEyfPh2vvvpqs7+X68HKEhERkfcIoiiK3h6ErystLUVYWBj0ej1CQ0M9fv4rZQbc8to2CAJw7vW76qwtRURERNevsX9/++Q6S22NVmlv8BZFwGJjtiUiImpJDEs+wNGzBPCWJ0RERC2NYckHOHqWAN5Ml4iIqKUxLPkAhUKAWmnvUzJaeEUcERFRS2JY8hHSFXGsLBEREbUohiUfIa3izbBERETUohiWfAQrS0RERN7BsOQjpFW82bNERETUohiWfISjssSb6RIREbUshiUfEVDds2RgZYmIiKhFMSz5iAA1K0tERETewLDkI1hZIiIi8g6GJR+hVVWHJVaWiIiIWhTDko9wTMMZzKwsERERtSSGJR8hTcOxskRERNSiGJZ8hFbFyhIREZE3MCz5CDZ4ExEReQfDko/g0gFERETewbDkIwJUjhvpsrJERETUkhiWfIRzg/fxS3qs2pcDURS9PCoiIiL/p/L2AKhxnJcOeH7dMRzN1aNXYihuah/u3YERERH5OVaWfIRWqixZcbXMCAC4Umrw5pCIiIjaBIYlH+E8DVdhtAAAyqv/S0RERM2HYclHSOssWayoNNmbvBmWiIiImh/Dko9wVJZKq8yw2OyN3WUGhiUiIqLmxrDkIwKqK0tFFSZpGytLREREzY9hyUc4KkslVWZpW5nBXN/hRERE5CEMSz7CEZacl1Yq5zQcERFRs2NY8hGOdZaccRqOiIio+TEs+QhHZcnZ9TR4W20i0n85j5OXSz05LCIiIr/HsOQjHPeGc3Y9laVdZ65h8Tcnsfibk54cFhERkd9jWPIRWhfTcNdTWcqvXu27uNLUwJFERETkjGHJRzgWpXR2PZWl0uqr6Axmq8fGRERE1BYwLPkIQRDqBKbruRqupNIelqoYloiIiK4Lw5IPqd3kbbLaYLQ0LvzopcqSzePjIiIi8mcMSz7E5fIBjawuORazZGWJiIjo+jAs+ZAbWT7AUVkyWWyw2sQGjiYiIiIHhiUfciPLB+idroJjkzcREVHjMSz5EFfTcNdbWQI4FUdERHQ9GJZ8iNbFNFxjK0vON+CtMjEsERERNRbDkg9xvdaS2cWRcjabKK2zBKDRV9ARERERw5JPcW7wjg7WAGjcNFyZ0QLnnu4qE5cPICIiaiyGJR/iHJbiQgMANC4sOVeVAPYsERERXQ+GJR8S4DQNF18dlhrTs+RYvduBYYmIiKjxGJZ8iKyyFFYdlhpRWSqpkt88lw3eREREjcew5EOclw6IC2l8ZUlfaxqO6ywRERE1HsOSD3FUlgQBiAnRAgDKDGaIooiD2UUoqTS5fF3taTiGJSIiosZjWPIhjrAUpFEhVKcCAJQaLDiQXYzJ7+7B/K+OuXxd7cpShcmKT/ddwOmCsuYdMBERkR9QeXsA1HiOdZYCNUqE6dQAAH2lGeeulgNAveGndljanlmAXWcKkdo5Eqv/e0gzjpiIiMj3sbLkQxwreAdrVYgItK+zVFxpQnH1NFtBqcHl6/S1puHOX60AAFwqqWquoRIREfkNnwhLO3fuhCAILh8HDhwAAGRnZ7vcv3fv3nrP+9tvv+HBBx9E+/btodPp0LNnT/zjH/9oqbd13RxLBwRqlYgIsoelkkqz1KtUYbK6bPh2XA0XpLGHrYIyIwCgqMJ1jxMRERHV8IlpuKFDhyIvL0+2beHChdi2bRsGDRok275161b07t1beh4VFVXveQ8dOoTY2Fh88sknaN++PXbv3o3HHnsMSqUSTz31lGffhAc4pt4iAjWICLR/bbLaZBWiglIDgmOCZa9zTMPFhwXg7NUKWKuX8640WVFlskKnqXvPOSIiIrLzibCk0WgQHx8vPTebzVi/fj3mzJkDQRBkx0ZFRcmOdeeRRx6RPU9OTsaePXuwdu3aVhmWbu0eg2dHd8fIlFjo1EpoVAqYLDacv1YhHVOgN6CLU1iqNFmQlW/vZWofGYizVytk5yysMCJJE9gyb4CIiMgH+cQ0XG0bNmxAYWEhZs6cWWdfWloaYmNjMXz4cGzYsOG6z63X6xEZGen2GKPRiNLSUtmjJQSolXhmdDf0TQqDIAhSdUkWlsrkfUuf7s1BcaUZHaMCMaJbTJ1zciqOiIjIPZ8MS+np6Rg7diySkpKkbcHBwVi+fDnWrFmDb7/9FsOHD8ekSZOuKzDt3r0bq1evxmOPPeb2uCVLliAsLEx6tG/fvsnv5UY4mrwrnVbkLig1Sl9Xmax476dzAIDZt3dFsLbudFthOcMSERGRO14NS/Pnz6+3cdvxyMzMlL0mNzcXmzdvxqxZs2Tbo6OjMXfuXKSmpmLw4MFYunQppk6dimXLljVqLMePH8fEiROxaNEijBkzxu2xCxYsgF6vlx4XL168vjfuIeHVlSVnzlfE/fz7VVwrNyIxLAD33NxOdrsUh0JWloiIiNzyas/SvHnzMGPGDLfHJCcny55nZGQgKioKaWlpDZ4/NTUVW7ZsafC4kydPYtSoUXjsscfw4osvNni8VquFVqtt8Ljm5qgsOXMOS47puUGdIqFWKqBzFZbKjXW2ERERUQ2vhqWYmBjExNTto6mPKIrIyMjAtGnToFbXrarUduTIESQkJLg95sSJExg5ciSmT5+O1157rdFjaQ3CXYalmvCTU1QJAOgQaW/gdnXVG3uWiIiI3POJq+Ectm/fjvPnz+PRRx+ts2/lypXQaDQYMGAAAGDt2rX48MMP8cEHH0jHrFu3DgsWLJCm9o4fP46RI0di7NixmDt3LvLz8wEASqXyukKct0Q0MA0nhaWo6rDEaTgiIqLr5lNhKT09HUOHDkVKSorL/YsXL8aFCxegUqmQkpKC1atXY/LkydJ+vV6PrKws6fmXX36Jq1ev4pNPPsEnn3wibe/YsSOys7Ob7X14ivM0nFopwGwVcaXUCFEUIQhCncqSy54lTsMRERG5JYiiKHp7EL6utLQUYWFh0Ov1CA0NbbHvu+bgRTz35VEAQI+4EGRV3xvu14V3IjRAhZSFm2Cxidg9fyQSw3U4c6Uco9/6UXaO/klhWP/U8BYbMxERUWvR2L+/fXLpALJzrizFhmoRVX0LlIJSA/L0BlhsIjRKBeJCAwC47lm6xqUDiIiI3GJY8mGO+8MB9mbvpAgdAOBUXikuVk/BJUXqoFTYVzl37lly3DqFDd5ERETuMSz5MOcG74hANW7tbm9K/+FEAS7U6lcC5GGpS0wQAKDKbMVPp69K4YqIiIjkGJZ8mPM0XHigBmN72++Jt/P0Fel+cM5hSauq+eNuFxEITfXzaR/ux63LduDxjw+hwmhpiaETERH5DIYlHxaqU8NxH+HIQDV6J4aiXbgOBrMNXxy0ryruHJYUCkEKTBGBapgsNmmfKAKbTuRjR9aVlnsDREREPoBhyYcpFYLUexQRpIEgCBjTOw5Azf3inMMSUNPkHa5TY1DHCADA6J6xuH+Q/f52jooUERER2TEs+bjI6iZvx5Tc5IFJUCsFBKgVuLNXHIZ1jZYd7+hbCgvU4P+7pw9ev6cv3pk6ECkJIQCATIYlIiIiGZ9alJLqmjOyK3ZkXsUtnSMBAL0Tw3D4pTHQKBVST5IzR1gK16mREh+KlHj7uhI94u1hiZUlIiIiOYYlH3fPgCTcMyBJti1YW/8fa0iAfV90iPxGwI7QlFNUiQqjBUFuzkFERNSWcBqujfnL2B54dHhnDO0SJdseGaRBTHWA+v1KuTeGRkRE1CoxLLUxI7rF4MU/9oJaWfePvkecYyqutKWHRURE1GoxLJHE0bfEJm8iIqIaDEskYZM3ERFRXQxLJOlYvSZTvt7g5ZEQERG1HgxLJIkKtjd4Xys3enkkRERErQfDEkmig+0LW5YaLLJboRAREbVlDEskCQ1QQ6mw32yuuNLk5dEQERG1DgxLJFEoBOn2KZyKIyIismNYIpmo6rBUWM7KEhEREcCwRLVEVzd5F1awskRERAQwLFEtUcGsLBERETljWCIZR89SYQXDEhEREcCwRLVI03Bs8CYiIgLAsES1sMGbiIhIjmGJZKRVvDkNR0REBIBhiWqpafDmNBwRERHAsES1OKbhilhZIiIiAsCwRLU4puEqTVZUmixeHg0REZH3MSyRTJBGCa3K/mPBJm8iIiKGJapFEASnVbwZloiIiBiWqI5IqW+JTd5EREQMS1RHoEYJAKgwWr08EiIiIu9jWKI6HGGpysSwRERExLBEdQRqVACAKjPDEhEREcMS1aGrrixVsrJERETEsER11UzDcZ0lIiIihiWqg5UlIiKiGgxLVIdOXR2W2LNERETEsER1OabhDKwsERERMSxRXbrqq+E4DUdERMSwRC4EchqOiIhIwrBEdfBqOCIiohpNCksVFRWeHge1IgENXA13/loFV/cmIqI2o0lhKS4uDo888gh++eUXT4+HWgHHNJyrQHT2ajnueHMn5nz2a0sPi4iIyCuaFJY++eQTFBUVYeTIkejevTuWLl2Ky5cve3ps5CXubneSU1hp/29RZYuOiYiIyFuaFJYmTZqEr7/+GpcuXcLjjz+OVatWoWPHjvjjH/+ItWvXwmJhr4svc7copclqAwBYrGKLjomIiMhbbqjBOyYmBnPnzsXRo0fx1ltvYevWrZg8eTISExPx0ksvobKS1QdfVNPg7SIsWexhyWyzteiYiIiIvEV1Iy8uKCjAypUrsWLFCly4cAGTJ0/GrFmzkJubi//93//F3r178cMPP3hqrNRCHGHJZLXBYrVBpazJ1I6wxMoSERG1FU0KS2vXrkVGRgY2b96MXr164cknn8TUqVMRHh4uHTN06FD07NnTU+OkFhRQ3eAN2NdaCnUKS2bHNJyNYYmIiNqGJoWlmTNn4oEHHsCuXbswePBgl8ckJibihRdeuKHBkXdoVQooBMAm2m95EhqglvbV9CxxGo6IiNqGJvUs5eXl4b333qs3KAGATqfDokWLmjwwZzt37oQgCC4fBw4cAABkZ2e73L93795GfY/CwkIkJSVBEASUlJR4ZNy+ShAE6Yq4SpNVqiYBnIYjIqK2p0mVJYvFgtLS0jrbBUGAVquFRqO54YE5Gzp0KPLy8mTbFi5ciG3btmHQoEGy7Vu3bkXv3r2l51FRUY36HrNmzUK/fv1w6dKlGx+wH9BplCg3WjB/7VEcv1SK754egQ5RgVJliQ3eRETUVjSpshQeHo6IiIg6j/DwcOh0OnTs2BGLFi2CzUN/oWo0GsTHx0uPqKgorF+/HjNnzoQgCLJjo6KiZMeq1ep6zlrjnXfeQUlJCf7yl794ZLz+wNHkvfdcEcqNFqzckw0AMFvsFSUre5aIiKiNaFJlacWKFXjhhRcwY8YM3HLLLQCA/fv3Y+XKlXjxxRdx9epVvPnmm9BqtXj++ec9OmAA2LBhAwoLCzFz5sw6+9LS0mAwGNC9e3f89a9/RVpamttznTx5Eq+++ir27duHc+fONer7G41GGI1G6bmrKpuv0zk1eQNARKA9dJqs9uUEzFYRoijWCatERET+pklhaeXKlVi+fDn+9Kc/Sdvuvvtu9O3bF++99x62bduGDh064LXXXmuWsJSeno6xY8ciKSlJ2hYcHIzly5dj2LBhUCgU+Oqrr6TFM+sLTEajEQ8++CCWLVuGDh06NDosLVmyBK+88opH3ktr5ViY0iEySAugpmcJsFeXVEqGJSIi8m9NmobbvXs3BgwYUGf7gAEDsGfPHgDA8OHDkZOT4/Y88+fPr7dx2/HIzMyUvSY3NxebN2/GrFmzZNujo6Mxd+5cpKamYvDgwVi6dCmmTp2KZcuW1fv9FyxYgJ49e2Lq1KmNfevS6/R6vfS4ePHidb3eFwTWCksalf1HxezU2M3lA4iIqC1oUmWpffv2SE9Px9KlS2Xb09PT0b59ewD2q8siIiLcnmfevHmYMWOG22OSk5NlzzMyMhAVFdXg9BoApKamYsuWLfXu3759O44dO4Yvv/wSACCK9r/8o6Oj8cILL9RbPdJqtdBqtQ1+f1+mU8t/NBxLBRidKksMS0RE1BY0KSy9+eab+K//+i98//330vIBBw8eRGZmphQ8Dhw4gPvvv9/teWJiYhATE9Po7yuKIjIyMjBt2rRGNW4fOXIECQkJ9e7/6quvUFVVJT0/cOAAHnnkEfz888/o0qVLo8flj2pXlszVwch5GQGutURERG1Bk8JSWloasrKy8N577yErKwsAMH78eHz99dfo1KkTAOCJJ57w2CAdtm/fjvPnz+PRRx+ts2/lypXQaDTS9ODatWvx4Ycf4oMPPpCOWbduHRYsWCBN7dUORNeuXQMA9OzZU7YaeVtUOyw5gpFzz5KZay0REVEbcN1hyWw2Y9y4cXj33XexZMmS5hhTvdLT0zF06FCkpKS43L948WJcuHABKpUKKSkpWL16NSZPnizt1+v1Urgj9wLUtcOSPRiZZNNwrCwREZH/u+6wpFarcfTo0eYYS4NWrVpV777p06dj+vTpbl8/Y8YMtz1St99+u9S31NbVnYarXoxSNg3Hz4qIiPxfk66Gmzp1KtLT0z09FmpF6k7DVVeWrGzwJiKitqXJtzv58MMPsXXrVgwcOBBBQUGy/W+99ZZHBkfeo9O4vhpONg3HBm8iImoDmhSWjh8/jptvvhkAcPr0adk+rujsH+q7Gs65ssQGbyIiaguaFJZ27Njh6XFQK1P7didWF0sH8P5wRETUFjSpZ8nhzJkz2Lx5s7RWEZuj/Uft252YXS0dwKvhiIioDWhSWCosLMSoUaPQvXt33HXXXcjLywMAzJo1C/PmzfPoAMk7usYGw3lG1eXSAZyGIyKiNqBJYenZZ5+FWq1GTk4OAgMDpe33338/Nm3a5LHBkfd0iQnG3gWj8PSobgBq1lSS3RuODd5ERNQGNKln6YcffsDmzZuRlJQk296tWzdcuHDBIwMj74sLDZB6lxwhifeGIyKitqZJlaWKigpZRcmhqKjI728w29aolfa5OEcVSbYoJXuWiIioDWhSWBoxYgQ++ugj6bkgCLDZbHjjjTdwxx13eGxw5H0qhT0sSUsH8N5wRETUxjRpGu6NN97AqFGjcPDgQZhMJvz1r3/FiRMnUFRUhF27dnl6jORFKqU9T0uLUvJ2J0RE1MY0qbLUp08fnD59GsOHD8fEiRNRUVGBe++9F4cPH0aXLl08PUbyopppOBFWmyhbW4nTcERE1BY0qbIEAGFhYXjhhRc8ORZqhZQKe54220RZvxLAyhIREbUNTQ5LJSUl2L9/P65cuQJbrQrDtGnTbnhg1Do4KktWm002BQewskRERG1Dk8LSxo0bMWXKFJSXlyM0NFR2PzhBEBiW/IjKUVmyirLmbsc2IiIif9eknqV58+bhkUceQXl5OUpKSlBcXCw9ioqKPD1G8iKV09IBtafheG84IiJqC5oUli5duoSnn37a5VpL5F+kBm+bq8pS3Wm4Lw/l4k/v7kFhubFFxkdERNTcmhSWxo4di4MHD3p6LNQKuZuGc7WC9+oDOdifXYTdZwtbZHxERETNrUk9SxMmTMBzzz2HkydPom/fvlCr1bL9aWlpHhkceZ/zNFydBm8XlSWD2b6tymRt/sERERG1gCaFpT//+c8AgFdffbXOPkEQYLXyL0p/oXYsSuliGs5VZclosf/ZV5n5M0BERP6hSWGp9lIB5L+k251YbXWufnO1zpIjUFWyskRERH7iunqW7rrrLuj1eun50qVLUVJSIj0vLCxEr169PDY48j6psuRq6QAXodlYfQwrS0RE5C+uKyxt3rwZRmPNVU6vv/66bKkAi8WCrKwsz42OvE6pcFwNZ4Op1vSqu8pSlcnS/IMjIiJqAdcVlkRRdPuc/I986QD5n7erdZZYWSIiIn/TpKUDqO1wLB1gsYp1roZztc6S1OBtYl8bERH5h+sKS4IgyG5t4thG/suxdIDZaoO59tVwtabhbDZRagKvMnMajoiI/MN1XQ0niiJmzJgBrVYLADAYDHj88ccRFBQEALJ+JvIPsqUDaleWajV4O+/nOktEROQvrissTZ8+XfZ86tSpdY7hTXT9i2PpAKuLdZZq9ywZnfZz6QAiIvIX1xWWMjIymmsc1EqplDUztbUDUO1pOEe/EgAY2OBNRER+gg3e5JbjajgAqKy1HEDtBm/nyhOvhiMiIn/BsERuOa6GA4AKoz0AaZz6mJxxGo6IiPwRwxK55ehZAmqucNNplABchCVzTVjiNBwREfkLhiVyS6EQ4MhLjspSoCMs1Z6G49VwRETkhxiWqEGOJm/H1FpNWKpdWaoJSJVmK1d4JyIiv8CwRA1SV5eWHNNwgRr7RZTu1lkSRXkPExERka9iWKIGOSpLjmk4R89SnXWWzPJwxKk4IiLyBwxL1CDH8gFVtabhzHXWWaoVlm6gydtiteFiUWWTX09EROQpDEvUIMfyARXV6ywFVU/D1W3wloejGwlLL288gRFv7MD+80VNPgcREZEnMCxRg1S1KkstMQ135kq57L9ERETewrBEDXLcTNdRWZKm4dw0eAM3Vlmqqg5etVcNJyIiamkMS9Qgx8KUhuoAEyhNw7mvLN3IKt5V1SGJi1sSEZG3MSxRg5ROq3gDQFA9Dd51Kks3EpaqQxJvm0JERN7GsEQNckzDOYQHaQAA1lrTcMZaVaAbqQpVmRzTcAxLRETkXQxL1CBHg7dDuE4NwMU0nMVz03COoMW1moiIyNsYlqhBakWtylKgPSzVbvD21DpLoihKr72RJnEiIiJPYFiiBtWuLEUE2qfhGqosVTXxSjaT1SYtS8BpOCIi8jaGJWqQqlbPUphjGs4mym6Wa/JQZclgqjmP4350DdmReQUzMvbjWrmxSd+TiIioPgxL1CC109VwKoWAYK1Keu68MKXRIr8dSpWpaTfSdQ5ZjaksWW0iZq44gJ1ZV/HxngtN+p5ERET1YViiBjlPwwVpVbLnFllYsocjRwN4Y6tCtTmHpcY0eO/MuiJ9rVHxR5qIiDyLf7NQg5yn4YK1KtlSAmantZUc03ChjrDUxH4j51W7GzOVt2pfjvS1rdYtWIiIiG6UT4SlnTt3QhAEl48DBw4AALKzs13u37t3b4PnX7FiBfr164eAgADExsZi9uzZzf2WfIrzNFywViWt6A24noZzXC3X1OZsw3VMw2Xml2KHU2WpklfPERGRh6kaPsT7hg4diry8PNm2hQsXYtu2bRg0aJBs+9atW9G7d2/peVRUlNtzv/XWW1i+fDmWLVuG1NRUVFRUIDs722Nj9wdKp6UDggNUshW9nVfxdlSWHFfLNbXB27nXyV11qsJowexPf4VzMYnrMhERkaf5RFjSaDSIj4+XnpvNZqxfvx5z5syBIMgva4+KipId605xcTFefPFFbNy4EaNGjZK29+vXzzMD9xNqpbyyJAgCVAoBFpsIi9NaS46eJcfVck1dwVvWs2S2QhTFOn/OAPD21tM4e7UCcaFaTOibiA93nee95IiIyON8Yhqutg0bNqCwsBAzZ86ssy8tLQ2xsbEYPnw4NmzY4PY8W7Zsgc1mw6VLl9CzZ08kJSXhT3/6Ey5evOj2dUajEaWlpbKHP1PVCkvO25xvnuuoLIXd4DScc8+S1SbWueecw77zRQCA5+/qiXYRuhv6nkRERPXxybCUnp6OsWPHIikpSdoWHByM5cuXY82aNfj2228xfPhwTJo0yW1gOnfuHGw2G15//XW8/fbb+PLLL1FUVIQ777wTJpOp3tctWbIEYWFh0qN9+/YefX+tjUohb/AGalb1vv3NnXjx62MAaipLiWH24HK5pAoVRgue+OQQPt7b+Ev6a1eHXE2tiaKI89cqAAA9E0KhU1cvV8DKEhEReZhXw9L8+fPrbdx2PDIzM2Wvyc3NxebNmzFr1izZ9ujoaMydOxepqakYPHgwli5diqlTp2LZsmX1fn+bzQaz2Yx//vOfGDt2LP7whz/gs88+w++//44dO3bU+7oFCxZAr9dLj4YqUb5OXWvpAABQOm378fRVADUN3v2SwqBVKVBcaca7P57F98fzsfyHLNkCls4aCkeuqkVFFSaUGSwQBKBDZKDT2k4MS0RE5Fle7VmaN28eZsyY4faY5ORk2fOMjAxERUUhLS2twfOnpqZiy5Yt9e5PSEgAAPTq1UvaFhMTg+joaOTk5NT3Mmi1Wmi12ga/v7+QLR0QUD0N51RtKiq3V+Ec03DBWhX6JYXhQHYxVuzKBgCUVJpxsagKHaICZef+302ZSP/5PL6ePQy9EkMBAFXmhlcCd1SVEsN0CFArEcDKEhERNROvhqWYmBjExMQ0+nhRFJGRkYFp06ZBrVY3ePyRI0ekQOTKsGHDAABZWVnSlF5RURGuXbuGjh07Nnpc/s556YAQxzScU2WpwmSFwWyVpuG0KiVu7hiBA9nFKDPW9B/9lltSJyztOVsIk9WGX3OKncJSw9NwjrDUKdp+PkdliT1LRETkaT7Vs7R9+3acP38ejz76aJ19K1euxGeffYbMzExkZmbi9ddfx4cffog5c+ZIx6xbtw4pKSnS8+7du2PixIl45plnsHv3bhw/fhzTp09HSkoK7rjjjhZ5T77AubIUVKvB26GowiRVlrRqBW7uEFHnPEdzS+psK6m0V6UKy2t6xGrfgNdVAMourA5LUUEAAF11WOLVcERE5Gk+sXSAQ3p6OoYOHSoLPM4WL16MCxcuQKVSISUlBatXr8bkyZOl/Xq9HllZWbLXfPTRR3j22WcxYcIEKBQK3Hbbbdi0aVOjKldthexquOppOGOtqbKrZUbp1icapTwsKRUCrDYRR3P1dc5dVFEdlipqboBbu7JUaap725Tsa5UAgM7R1WFJraz3WCIiohvhU2Fp1apV9e6bPn06pk+f7vb1M2bMqNMjFRoaivT0dKSnp3tiiH5JLbsazh5KrpQZZcfk6Q3S11q1AhEaFTpEBiKnqBJ/7JeA9Ucu4/glPaw2UVrU0my1odRgDzfyypI8iLmqFknTcLUqS2zwJiIiT/OpaTjyDqXsdieuK255+irpa031tN0Tt3fBgA7h+Nu4FOjUSlSYrDh3tVw6rqTSLH19rbwmfNUOR7Wn4URRrJmGq1VZYoM3ERF5GsMSNUi+dIDS5TGOypJSIUg9Tg/e0gHrnhyGxHAd+rSzN28fu1QzFVdcWVNNKqwwIaewEqv25aDUYIaz2mHpapkRlSYrFNXLBgA1Dd5mqyi7uS8REdGN8qlpOPIO5wbvkHoqS5dL7JUljdJ1/k6JD8WB7GJkFZShymRFnr5K6lcCgMJyI1795gS2nqq5Ka7jliq1p9bOVFenkiICoVHZv59j6QDAXplS1zMOIiKi68W/UahBTrNwUoP34om90SshFNOG2JdYyK+uLGnVrn+kesSHAABO55dh4frjGLn8R2w7VSDtL6404+Rl+W1jIoPsN+StXVnKyi+TnRMAtCqFNE72LRERkSexskQNcqyfBNRMwz08pBMeHtIJn+6z38bEMQ2nVbkPS5n5ZagwFgOArIoEAJedmsQBe1i6Umas04fkCEspTmFJEASpL4p9S0RE5EkMS9Qg52UCtCp5z1JUdfXncnWDt6aesNQ91h5snK+ac1zRVp+oYPu5a6+7lOmisgQAOo0KFSYrF6YkIiKP4jQcNcjdQo9Rwfbbvjhu+1Y7TDmEBaoRHxpwXd83Msh+bufwY7OJOF1Qt7IEADqN/ceZlSUiIvIkhiVqkPM0XG2OviKHbrHB9R7bvVa4aYijalXpFH4uFlei0mSFRqWQ1lhyCFTbC6XsWSIiIk9iWKIG1b6fm7OoWmHp1u713+uvdiWoIRGB9nMbnMKPYwquW2yw7Co9AAjgwpRERNQM2LNEDbp3QDsU6A0Y0iWqzr7QAPlSAu7CUvc4e1hyLHJprb49SlKEDrnFVXWOd/QsOU/DuboSziHQccsTTsMREZEHsbJEDVIpFZgzqhsGdYqss0+hkN9Qt124rt7zpHaORIBagdu7xyAxvKZ/qb6pO0fV6vglPb48lAuz1YYjF0sAuK5SSTfTZWWJiIg8iGGJPKZ7XP39SgDQPjIQe+aPwr+n3CytvA0A3eJqgk9yTE0fUqjOXrUqM1rwlzW/4alVv2J7pn25gdt7xNY5vyMs8Wa6RETkSQxLdMPuuzkJALDwj70aPDYiSIMAtRLtI2rCUtfqypJKIaB/Uri0PUwnn+LbfMK+iOVdfeOlKT1nNfeH4+1OiIjIc9izRDds8aTeeGZUN7eN4LW1r64sKRUCeiXY7xuXHBOE2FCtdEyvhFA8O7o7OkTp8GPWVXx95DIA4OlR3VyeUwpLrCwREZEHMSzRDQvUqNAh6vp+lBxhKSJQjT7twvCPB25C97gQXC0z4r0fzwGw90M9M9oejG7vHovCChNuah+OlPjQesbhqCxZIYoilm3OQkpCKNL6Jzb1rRERETEskXf0SQyFIABdYuxTcBNvagcA6JkA/OOBm9AtVj7NFhGkwcezUt2eM0BdE5Y2nyjAf3aeBQCGJSIiuiEMS+QVyTHB+P6ZES5X9XYEp+sVKDV4W/Fbbom0XRRFCIJQz6uIiIjcY1gir6lvOq2ppKUDzFZcclq3yWC2SfuIiIiuF8MS+Q1Hg3elyYriypob9uqrzAxLRETUZFw6gPyGIxBVGC3Iyi+VtpcazN4aEhER+QGGJfIbjp6lU3llMDittaSvYlgiIqKmY1giv+G4Gq7cKF9nqZRhiYiIbgDDEvmNQI3rFjxWloiI6EYwLJHf6BIT5HIpAlaWiIjoRvBqOPIbIQFq/PjX27HrzDVcLjHgaG4JvjiYC30Vb39CRERNx7BEfkWrUmJkShwAYOn39rWWeDUcERHdCE7Dkd8K1dn/LcCeJSIiuhEMS+S3wnRqAOxZIiKiG8OwRH4rNMAellhZIiKiG8GwRH5LqiwZ7A3eoigyOBER0XVjWCK/FVprGu6tLadx06s/YM/ZQm8Oi4iIfAzDEvmt2j1LW09dgSgCB7OLvDksIiLyMQxL5LdCA+xXw5UZLagyWfF7QRkA4LK+ypvDIiIiH8OwRH7LMQ0HAAcvFMFiEwEAl0sMKKowYdPxPJittvpeTkREBIBhifyYWqlAoMZ+c13nPqU8fRVe/+4UHv/kV3x7NM9bwyMiIh/BFbzJr4Xp1Kg0WbHbKSxdLjFAgB4AcP5ahbeGRkREPoKVJfJrjrWWjlwskbaVGy04c7UcAHClzOiNYRERkQ9hWCK/FubUtwQAgmD/r7W6f+kqwxIRETWAYYn8WnhgTVgK1CjRPTZEtv9qOcMSERG5x54l8muP3ZoMs9WGKrMVaf3bYdupAmRVLyEAAFdLDV4cHRER+QKGJfJrgzpFImPmLdLzk3l62f6r5UaIogjBMT9HRERUC6fhqE1JCNPJnputvF8cERG5x7BEbUq7cF2dbbwijoiI3GFYojYlISwAAKBUCOgYFQgAOHulHNtOFcBWfYUcERGRM4YlalN6JYYiLlSLkSmxUpVpzmeHMWvlQXxx8KKXR0dERK0RwxK1KSEBauz620i8//BAxIRoAUC6Z9wPJwu8OTQiImqleDUctTkqpf3fCLHVYclhz9lCGC1WnL1Sgde+Owmj2YZxfeLx8JCO0KqU3hgqERG1AgxL1GbF1ApLVWYrFn9zEl8cyIXJagMAHLxQjCCtCg/e0sEbQyQiolaA03DUZjmHJU11temTvTkwWW0Y3TMWd/aKA8Cb7RIRtXUMS9RmRQbVhKVpQzpKX4/rHY/3Hx6EIclRAIBLJVUtPjYiImo9fCIs7dy5E4IguHwcOHAAAJCdne1y/969e92e+8CBAxg1ahTCw8MRERGBsWPH4rfffmuJt0Ve1iEyUPr6yTu6IjZEi5vah2P5n/pDoRCQWH213KVihiUiorbMJ8LS0KFDkZeXJ3s8+uij6Ny5MwYNGiQ7duvWrbLjBg4cWO95y8vLMW7cOHTo0AH79u3DL7/8gpCQEIwdOxZmM1d19nedo4Ow8pFb8N3TIxAZpMHu+SPx1RNDEaS1t/I5lha47Kay9NvFEsxdfQR5egYqIiJ/5RMN3hqNBvHx8dJzs9mM9evXY86cOXXu6RUVFSU71p3MzEwUFRXh1VdfRfv27QEAixYtQr9+/XDhwgV07drVc2+CWqXbusdIXzuuknNIDLcvYHmlzAijxeryirh3dp7FphP56BIbjNl38OeFiMgf+URlqbYNGzagsLAQM2fOrLMvLS0NsbGxGD58ODZs2OD2PD169EBUVBTS09NhMplQVVWF9PR09OzZE506dar3dUajEaWlpbIH+Z/IIA20KvuvSIHe9S1RTheUAXBffSIiIt/mk2EpPT0dY8eORVJSkrQtODgYy5cvx5o1a/Dtt99i+PDhmDRpktvAFBISgp07d+KTTz6BTqdDcHAwNm3ahO+//x4qVf1FtyVLliAsLEx6OKpS5F8EQZCm4lw1eRvMVmQX2q+Uy9cbGjzf5hP52Huu0LODJCKiZufVsDR//vx6G7cdj8zMTNlrcnNzsXnzZsyaNUu2PTo6GnPnzkVqaioGDx6MpUuXYurUqVi2bFm937+qqgqzZs3CsGHDsHfvXuzatQt9+vTBhAkTUFVVf6VgwYIF0Ov10uPiRd4mw18luulbOne1Ao7byeU1EJYKy4144pND+O+PD3l8jERE1Ly82rM0b948zJgxw+0xycnJsucZGRmIiopCWlpag+dPTU3Fli1b6t2/atUqZGdnY8+ePVAoFNK2iIgIrF+/Hg888IDL12m1Wmi1Wpf7yL84+pZchaXfr5RJXzfU4H213AibCOirzPX2PxERUevk1bAUExODmJiYhg+sJooiMjIyMG3aNKjV6gaPP3LkCBISEurdX1lZCYVCIWsSdzy32WyNHhf5L6myVB2GsvLLoFYKSI4JlvqVAKC40gyD2YoAtesQpK+subqy3GCBNphhiYjIV/jE1XAO27dvx/nz5/Hoo4/W2bdy5UpoNBoMGDAAALB27Vp8+OGH+OCDD6Rj1q1bhwULFkhTe3feeSeee+45zJ49G3PmzIHNZsPSpUuhUqlwxx13tMybolZNWmupxIC95wrx0P/thU0ERqXE4lqFSXZsnt6AztFBqDJZoVIKUDtdXVdqsEhflxstiApmZZKIyFf4VFhKT0/H0KFDkZKS4nL/4sWLceHCBahUKqSkpGD16tWYPHmytF+v1yMrK0t6npKSgo0bN+KVV17BkCFDoFAoMGDAAGzatMltRYrajqTqsJSVX4q5q49IPUrbMq/UOTZPX4Ws/DI8t+Y3JMcEYf1Tw6V9pVVOlSWjpc5riYio9RJEURS9PQhfV1pairCwMOj1eoSGhnp7OORBefoqDFmyXXreOToIb/2pPx76v32oMlsBAD3iQpBVUIbhXaPxy5lr0rHHXxmL4OoFLjN2nccrG08CAFY/9gekVt9KhYiIvKexf3/75NIBRC0lIUyHFTMH494B7XBzh3D8v4cGYECHCLz38EAIAtAxKhB92oUBgCwoAcDFokrp69Iq+TQcERH5Dp+ahiPyhtt7xOL2HrGybbd2j8G2ubchWKvCx3svSNsVAtAuQoeLRVXIKapEzwT7v1T0nIYjIvJZDEtETZQcEwzAXn1yuKVzJKKDtbhYVCWvLBlqwlKZgWGJiMiXMCwR3aCEsADp6zG94nGt3H5rlBzZNBwrS0REvoo9S0Q3KDa0ZhmAMb3j0CEyEECtsGSQr7NERES+g5UlohuUEh+Ku/rGIz5Uh6SIQNdhiQ3eREQ+i2GJ6AYpFQL+M2Wg9Lx9dVjKLaqCzSZCoRDY4E1E5MM4DUfkYQlhAVApBJisNhSU2W+wy2k4IiLfxbBE5GEqpQLtIuxXyOUUVsJmE2XVJFaWiIh8C8MSUTNw7lsqM1rgvE5+GcMSEZFPYVgiagZJEfawdLG4SrZsAACUG8yuXkJERK0UwxJRM0iqnoa7VFwla+4GOA1HRORrGJaImkG78OqwVFIpNXdrlPZfNzZ4ExH5FoYlomaQWB2WLpcYpDWWEsLtK31XmKyw2cR6X+tpV8oMuFpmbLHvR0TkbxiWiJqB42q4PH0V9FUmAECi0z3kKkwtU10yWqwY//bPGP+Pn2Gx2lrkexIR+RuGJaJmEBeihVIhwGwVceZKOQAgJkQLtVIA0HJ9S5dLDCisMOFaubFO7xQRETUOwxJRM1ApFYgPtU+7ncorAwCE6dQI1toXzW+pvqXLJVXS1yUMS0RETcKwRNRMHE3ep/JKAQChOhWCA+xhqaXWWrpUXBOWWFkiImoahiWiZuLoWyqssPcsxYcGIEhTt7J0/JJemqrztEslDEtERDeKYYmomTgqSw6jesYhpLqy5OhZKq4w4d53duP+9/Y0yxVyztNwtRfHJCKixmFYImomiU5h6eYO4UgM19X0LFWHpROXS2Gy2FBYYUJeqcHjY7isZ2WJiOhGMSwRNRPHNBwA3NU3AQAQHKAGAJRU2qfmMvNLpWMuXKvw+Bice5ZKKhmWiIiagmGJqJk4T8M5wlLnKPs94/617QwO5xQjM79MOuZ8oT0sncorxexVvyK3uPKGvr/NJuKyvqZaxcoSEVHTqLw9ACJ/1SUmCI/dmoyoII00Jffft3XB3nNF2J9dhD9/dAhRQRrp+AuF9nD09y2n8cPJAsQEa7Ho7l6oMFml6bvrUVhhgslSsxClvsosNXzX7qciIqL6sbJE1EwEQcDzd/XEf9/WRdoWpFUhY+ZgxIRoca3ciKwCp8rStQrYbCL2nS8CABzOKcYn+3LQZ9FmbDqef93f37m5GwCulRvxx3/+jLv/9YssRBERkXsMS0QtLEirwoTqaTlnFworcDKvVJouO3G5FCt3ZwMAtp4quO7vc6lWWDp5uRTFlWYUVZiQp6+q51VERFQbwxKRF6TdlCh9HR1sn4q7UFiJ3WevSdsttppbpWQ59TbV52JRpSwgOSpLEYH2pvIrTjfTzS1mWCIiaiyGJSIvGNA+HO0j7X1Dd/SIhUohwGixYd3hywAAQZAf//uVMljdrMN0uaQK497+CWn/+gUGsxUAcPaqPWj1Tgyrc/wlhiUiokZjWCLyAkEQMGtYZwDA+L7x6BBpv0rOcWuUu/rIp+kMZhtyiuq/Ou7vW06jwmRFYYUJB7LtPU97z9n/OzIlts7xuSUMS0REjcWwROQlM4Z1xqlXx2FkShzaV4clAEiK0GHKHzpIzx1XwmXl2/uZRFHEp/suYNK/d2HvuUL8XlCGr37NlY7/+fdryNNX4fy1CigE4M5ecXW+940uS0BE1JZw6QAiL9JplACAWzpH4sfTVzGsaxSW3tsPMSFadIwKREiACt1jQ7D28CXM++I3VJisiArSSPebm/PZYUQHa2ET7b1P18pN+On0VaTEhwAA+iaFy1YSd7hUXIVDF4ogisCgTpHSdoPZigC1sgXeORGR72BYImoFnry9Cyb0TUDHqEAI1Q1L2+fdDlEUkbErGzh8CRUmey9SYYUJCgGICtbiapkRV8uMiA7WIGPGLbj7//2CzPwyfH3E3vs0JDkKSoWA0AAVSp1u3nu6oAwP/d8+AMD+50cjLFCNtb/m4i9rfsPEm9rh9Xv6SkGOiKit4zQcUSsgCAI6RQdJQQkAlAoBKqUCPaqrRIC9ArVi5mBseGo40qcPgkohQKNU4L2HB6JvUhj6tAsFAPx0+ioAYGiXKABAWPUVcQ7FlWYYLTYYLTb8mlMMURTxr+1nYBOBdYcvYWr6PrcN5c42Hc/D14cv3dD7JyJqzVhZImrleiaESl8vuruX7Oq2758ZAZVSgc7RQQCAsb3icfySvUk8RKvCoE4RAIAwnRoXUQWVQoAgAGZrTRA6kF0ErUqB89cqEFRdTTp0oRj7zhdiaJdot2PL1xvw5Ke/wiYCfdqFomtsiNvjiYh8EcMSUSsXE6LFu1NvhkalqLMMQLc4eTh54vYuGNgpAoXlJqTEhyBQY/8VD9PZK0vxYQFQCILsyrqD2cW4UP180oB2sFhFrD54ERt/u9xgWFp/5BIcBajvjuXj6VEMS0TkfzgNR+QDxvVJwMiUule11aZSKjC0SzTu7p8oC1LhOvvCl+3CdUiKsDd8O2b8jlwswebq26k8lNpBWjDzu2P5stuimCw2aQ0nABBFUXYV3nfH8pr47oiIWjdWlojagNDqylK7CB2U1SlpaJconMorQ1H1lXWje8aid2IYrDYRMSH25vEvDl7Ebd1jsO7wJfzfz+cQolVhw5zh2HqyAAeyi3G6oBwapQI2UURmfhnOXi1Hl5hgmK02LP7mJEqrzFg8qQ9CAtT1jo2IqLVjWCJqAxxLCdzUPhzxoQFYf+QyHv5DR3z16yVsOVmAQI0Sr0zsA8DeWD6xfyI++OU8Xvz6uOw8ZQYLHvq/vThdUC5tG90rFhVGK348fRWvf3sKC+5KwX92nMXa6qbv3OIqvD9tECKDNC30bomIPEsQRbFxl7xQvUpLSxEWFga9Xo/Q0NCGX0DUwkRRxLlrFegcFQSFQoAoihAEAd8fy8P/rD6CVyf2xv2DaxbCrDJZ8cbmTHxx4CIMFhtu7hCO23vE4s0fsuD4P8aEvglIDA/A9KGdcOZKOR5ZcQDOF9ApFQIC1UqUGS1QKwWM7R2Pl9N6IzpY28Lv3rUjF0vw4Pt7MW9Mdzw6ItnbwyEiL2js398MSx7AsES+zGYToVAILvcZLVZYbaLUKL74m5NI/+U8xvSKw7tTB8pe99vFEixcfxyZ+WWIC9Xib+NS0DEyCPPXHsWJy/Yr9KKDtZiS2gHj+8YjJb5pvyuiKGL/+SL0bx9+Qwtozv/qKD4/cBHd44Lxw7O3Nfk8ROS7GJZaEMMStRVWmz2oDOoUAbWy8deHHMvVY96aI9L0nUohYOUjt2BYV/nVdgazFZtP5KNTVBCCtEoczdVjVEqcbJ2olbuzsWjDCfwhORIfz0q9rnE4iKKIIUu2I7/UAAD47aUxddaiIiL/19i/v3k1HBE1mlIhYEiXqOsOKH2TwrDhqeFYem9fpHaOhMUm4olPDuHbo3nQV5kB2Ctc//P5ETzz+RFM/PcujH7rJ8z94jc8/skhOP5NZ7WJ+L+fzwGw3yj4b18dxdHcEjTm33z/3nEGz64+gjKDGZn5ZVJQAoBfc4qlr0/llWLCP3/GmoMXr+s9EpH/YmXJA1hZImo8g9mKB/9vLw7nlAAAQgJUWDa5H3afLcRHey5ArRSgVSmlZQtMVhtm39EFgRoVTBYb/rHtdwRqlKg01SxjcFv3GDw3tgeO5urRt10YdBoldmZdwYAO4RjYMRIbf7uMOZ8dBgCMTInFwI4RWLY5S3r9vTe3g1IQcFuPGPxnx1mczCtFSIAKv/x1JEICVNJ04/lrFUgMD4BWpYS+yozQABXMVhFfH76EWzpHIkirwttbT+NauREdo4LwlzE9oFHx36RErRWn4VoQwxLR9SmuMOH/7TiDbacKkF1YKdu3/L/6Y+JNibCJwIrd5/H6d5l1Xj/7ji7olRCGLw5exN5zhTA6rQdVW/+kMJy7VoEyp3vjKQTAJgJ924Xh2CV9va8d1DECmfll+NOg9kiK0OHVb06ie1wwUjtH4eO9F/DHfgmIDtZixe5shOnUiAvVyq4UfCWtN6YP7XQdnwwRtSSGpRbEsETUNCaLDS9+fQxfHMxFcnQQ5o7pjj/2S5T2W6w2TM/Yj1N5ZbipfTh2n70GrUqJzf9zK+LDAgAAxy/p8eePDiK/1IB+SeE4dbkUZpsNN3eIwG8XS2CpvkSvT7tQPDo8GfPXHoXBbEOAWoEPpg3G1HT7DYU1SgVMVnvomtA3Ad82cZHN6GAt7uwVh8/25yAmRIufnrujzk2JTRYblnx/CiculWLhH3uhb1JYPWcjoubEsNSCGJaIbsylkirEhwZAWc9VeQ5lBjOsNhHhgfI1m6pMVlSZrYgM0kBfZYbFakNUsBa5xZX4NacE+iozxvaOQ2xIAMoMZuQWVyE8UI340ACkvr4NV8qMeGNyP7QL16Gg1IBJN7XD9Iz9OJJTgjt7xUlrRo3rHY+SKhPOXq3AmF5x+HRfDgAgrX8iSg1mZOaV4YPpg9A9LgQjl+9EbnEVRqXEYmjXaJgsNiTHBCFIo8K/d5zBnnOFAAC1UsCzd3bHzKGdcbG4Eu3CdQjS2q8+NFttMFtt0tWItV0pNeCl9SdwMq8Uf7//JgzsGHFDfw5EbQ3DUgtiWCLyXcdy9bhYXInxfeIhCDVhzVpdkVIqBHy8JxtZBWV4/q6eCNSopIby5T+cxonLerx9/wCEBaplyzCsP3IJz3x+pN7vG6RR4uaOEfj592sAaqYGAzVKDOsaDYVgb2IvM5hxW/cYDOsajbjQAFhtIkxWG05eLsVXv+ZK04uBGiXm3tkd/duHo31EIBSKmhBpNNtgstpQbrBgy6kCXCsz4v7B7XGxqBLZhZX4Y78E9E4MgyBAthyDKIr4Nce+UvugjhHoGhsMAMgpqkRplQVatQIJYQH46fQ1nLtajqFdo7Du8CUcOF+Mp0Z2xaiesag0WREdrMWlkiqcv1qBmzqEw2C2IqeoEh0jAxESoIYIERqlQvb5u+JYH4zIUxiWWhDDEhG5svdcIXZkXkFuSRXUCgHHLulRZrDgjh6xmDWiM7rFBuOrXy9hyXenUFhhgkalkN2PrzH6tgtDSIAKu88WemTM0cEaAAJMFit0GiUKSo3SvpAAFXRqJa6UGes/QT3nvFZuv62OWinAbK37145CAAI1KgRplQjXaSBChMlir6qVGc24WmaEwWxDiFaF2FAt1EoFFIIApUKAQiFAKdiDrWOb0WJD9rUKCAIQExKAuFAtYkO0iAzSQoQIm02ETQRUSgEapQJWm4jc4ipoVQokhutgstpQabTAYhORGK5DSIAKl4qrsOlEPgI1KozrHQ+rzQYRgFalQFGFGSJEBFVXAW2iCJVCQJfYYIQHalBmMKPMYEGZwQyD2YaIIA1igrWICFSj0mxFmcECg9mKyEANyoxmXC4xoGNUIKw2EReLKhEfpkO4Tg2z1YYAjRLa6nCpEOz3hFQrBVSarLCJIrQqBbQqJRSCAKtNhFUUYbWJsDn+67RNo1IgTKeGwWyFKNr/DMxWGyw2EYIACAAUggBBqPmvAAEKRfV/Bft9JgVBgKdirCMQC6i5h6UAAfFhAR6/YIJhqQUxLBHRjag0WZCvN6BjVBAO5xTj+CU9RAC9EkIRFazFd8fykJVfhqvlRmiq/2KMDNIi7aZEDO8aDYvNhhW7srH3XCFOF5QjT18FEYBOrUSgRgmtSgmNSgGNUoG+SWEI1qqw+sBFJIQFoF9SGL47nl9vSNOqFOiXFIbfcvXSMRqlAlHBGpQbLSgzWBATokWfxFD8cuYakqODMaJbNFbuya4TiuJCtVL4ig3RXnfoorZt+7zbkBwT7NFz+lVY2rlzJ+644w6X+/bv34/BgwcDsJdoly9fjvfffx8XLlxAdHQ0nnzySbzwwgv1nruoqAhz5szBxo0boVAocN999+Ef//gHgoMb/wfCsERErYnVJlb/i79x/9Y3WqywWEWYrTbkFldBIQjQqAToqyzoEhOE8EANjBYrsq9VotRgRt92YQhQKyGKIkqrLAjSKqFSKmTTkOVGi71yoVTgVH4p4kMDkBAWgAuFlQgOUCE6WAuD2SpdyWgwW1FpsqLCaEFxpQlKQYBKqUCFyWKvJoUEIECjQGmVBVfLjFLFxGYT5V9XV0yUCgGdooKgEAQUlBlwpdSAK6VGFFWaoBAEqBT2som1+n2LgL2iZLEhv9SAAJUSQVolBNjvb2iwWBGkUeGOlFgUV5qw71wRQnUqKAQBBrMVEUEaKAVBWtJCEOzvKSu/DFVmG0ICVAgNUCFYq0KAWomiChOulhtRUmlGoEYpbS+uNCFApURieACyCyuhVAjoGBWIfL0BFSYr1AoBVWYrzFYbbKJ9fTKLzf4edBollIK9qmYwWyECUErVN0ClUEDhogKnrzJDVz39WmmyQqtSQKkQIEKEKKL6IUKEvWJmc7FNFO3/vdHqkuj0heNrR0z55ukR6BwddIPfQc6vwpLJZEJRUZFs28KFC7Ft2zacPXtW+h/C008/jR9++AFvvPEG+vbti6KiIhQVFeHOO++s99zjx49HXl4e3nvvPZjNZsycORODBw/GqlWrGj0+hiUiIiLf41dhqTaz2Yx27dphzpw5WLhwIQDg1KlT6NevH44fP44ePXo06jynTp1Cr169cODAAQwaNAgAsGnTJtx1113Izc1FYmJiA2ewY1giIiLyPX59u5MNGzagsLAQM2fOlLZt3LgRycnJ+Oabb9C5c2d06tQJjz76aJ2KlLM9e/YgPDxcCkoAMHr0aCgUCuzbt69Z3wMRERH5Bp8MS+np6Rg7diySkpKkbefOncOFCxewZs0afPTRR1ixYgUOHTqEyZMn13ue/Px8xMbGyrapVCpERkYiPz+/3tcZjUaUlpbKHkREROSfvBqW5s+fb7/c0M0jM1N+q4Pc3Fxs3rwZs2bNkm232WwwGo346KOPMGLECNx+++1IT0/Hjh07kJWVBU9asmQJwsLCpEf79u09en4iIiJqPVwvC9tC5s2bhxkzZrg9Jjk5WfY8IyMDUVFRSEtLk21PSEiASqVC9+7dpW09e/YEAOTk5LjsY4qPj8eVK1dk2ywWC4qKihAfH1/vmBYsWIC5c+dKz0tLSxmYiIiI/JRXw1JMTAxiYmIafbwoisjIyMC0adOgVqtl+4YNGwaLxYKzZ8+iS5cuAIDTp08DADp27OjyfEOGDEFJSQkOHTqEgQMHAgC2b98Om82G1NTUeseh1Wqh1WobPW4iIiLyXT51Ndy2bdswevRonDp1CikpKbJ9NpsNgwcPRnBwMN5++23YbDbMnj0boaGh+OGHHwDY12SaNm0atm3bhnbt2gGwLx1QUFCAd999V1o6YNCgQVw6gIiIyM/55dVw6enpGDp0aJ2gBAAKhQIbN25EdHQ0br31VkyYMAE9e/bE559/Lh1TWVmJrKwsmM1madunn36KlJQUjBo1CnfddReGDx+O999/v0XeDxEREbV+PlVZaq1YWSIiIvI9fllZIiIiImppDEtEREREbjAsEREREbnBsERERETkBsMSERERkRteXZTSXzguKOQ94oiIiHyH4+/thhYGYFjygLKyMgDgLU+IiIh8UFlZGcLCwurdz3WWPMBms+Hy5csICQmBIAgeO6/jnnMXL17k+k0N4Gd1ffh5NR4/q8bjZ9V4/Kwarzk/K1EUUVZWhsTERCgU9XcmsbLkAQqFAklJSc12/tDQUP4yNRI/q+vDz6vx+Fk1Hj+rxuNn1XjN9Vm5qyg5sMGbiIiIyA2GJSIiIiI3GJZaMa1Wi0WLFkGr1Xp7KK0eP6vrw8+r8fhZNR4/q8bjZ9V4reGzYoM3ERERkRusLBERERG5wbBERERE5AbDEhEREZEbDEtEREREbjAstWL//ve/0alTJwQEBCA1NRX79+/39pC87uWXX4YgCLJHSkqKtN9gMGD27NmIiopCcHAw7rvvPhQUFHhxxC3np59+wt13343ExEQIgoCvv/5atl8URbz00ktISEiATqfD6NGj8fvvv8uOKSoqwpQpUxAaGorw8HDMmjUL5eXlLfguWkZDn9WMGTPq/JyNGzdOdkxb+ayWLFmCwYMHIyQkBLGxsZg0aRKysrJkxzTm9y4nJwcTJkxAYGAgYmNj8dxzz8FisbTkW2l2jfmsbr/99jo/W48//rjsmLbwWb3zzjvo16+ftNDkkCFD8P3330v7W9vPFMNSK7V69WrMnTsXixYtwq+//or+/ftj7NixuHLlireH5nW9e/dGXl6e9Pjll1+kfc8++yw2btyINWvW4Mcff8Tly5dx7733enG0LaeiogL9+/fHv//9b5f733jjDfzzn//Eu+++i3379iEoKAhjx46FwWCQjpkyZQpOnDiBLVu24JtvvsFPP/2Exx57rKXeQotp6LMCgHHjxsl+zj777DPZ/rbyWf3444+YPXs29u7diy1btsBsNmPMmDGoqKiQjmno985qtWLChAkwmUzYvXs3Vq5ciRUrVuCll17yxltqNo35rADgz3/+s+xn64033pD2tZXPKikpCUuXLsWhQ4dw8OBBjBw5EhMnTsSJEycAtMKfKZFapVtuuUWcPXu29NxqtYqJiYnikiVLvDgq71u0aJHYv39/l/tKSkpEtVotrlmzRtp26tQpEYC4Z8+eFhph6wBAXLdunfTcZrOJ8fHx4rJly6RtJSUlolarFT/77DNRFEXx5MmTIgDxwIED0jHff/+9KAiCeOnSpRYbe0ur/VmJoihOnz5dnDhxYr2vaauflSiK4pUrV0QA4o8//iiKYuN+77777jtRoVCI+fn50jHvvPOOGBoaKhqNxpZ9Ay2o9mcliqJ42223ic8880y9r2mrn5UoimJERIT4wQcftMqfKVaWWiGTyYRDhw5h9OjR0jaFQoHRo0djz549XhxZ6/D7778jMTERycnJmDJlCnJycgAAhw4dgtlsln1uKSkp6NChQ5v/3M6fP4/8/HzZZxMWFobU1FTps9mzZw/Cw8MxaNAg6ZjRo0dDoVBg3759LT5mb9u5cydiY2PRo0cPPPHEEygsLJT2teXPSq/XAwAiIyMBNO73bs+ePejbty/i4uKkY8aOHYvS0lKpkuCPan9WDp9++imio6PRp08fLFiwAJWVldK+tvhZWa1WfP7556ioqMCQIUNa5c8Ub6TbCl27dg1Wq1X2QwAAcXFxyMzM9NKoWofU1FSsWLECPXr0QF5eHl555RWMGDECx48fR35+PjQaDcLDw2WviYuLQ35+vncG3Eo43r+rnynHvvz8fMTGxsr2q1QqREZGtrnPb9y4cbj33nvRuXNnnD17Fs8//zzGjx+PPXv2QKlUttnPymaz4X/+538wbNgw9OnTBwAa9XuXn5/v8mfPsc8fufqsAOChhx5Cx44dkZiYiKNHj+Jvf/sbsrKysHbtWgBt67M6duwYhgwZAoPBgODgYKxbtw69evXCkSNHWt3PFMMS+ZTx48dLX/fr1w+pqano2LEjvvjiC+h0Oi+OjPzJAw88IH3dt29f9OvXD126dMHOnTsxatQoL47Mu2bPno3jx4/L+gTJtfo+K+e+tr59+yIhIQGjRo3C2bNn0aVLl5Yeplf16NEDR44cgV6vx5dffonp06fjxx9/9PawXOI0XCsUHR0NpVJZp/O/oKAA8fHxXhpV6xQeHo7u3bvjzJkziI+Ph8lkQklJiewYfm6Q3r+7n6n4+Pg6FxBYLBYUFRW1+c8vOTkZ0dHROHPmDIC2+Vk99dRT+Oabb7Bjxw4kJSVJ2xvzexcfH+/yZ8+xz9/U91m5kpqaCgCyn6228llpNBp07doVAwcOxJIlS9C/f3/84x//aJU/UwxLrZBGo8HAgQOxbds2aZvNZsO2bdswZMgQL46s9SkvL8fZs2eRkJCAgQMHQq1Wyz63rKws5OTktPnPrXPnzoiPj5d9NqWlpdi3b5/02QwZMgQlJSU4dOiQdMz27dths9mk/6G3Vbm5uSgsLERCQgKAtvVZiaKIp556CuvWrcP27dvRuXNn2f7G/N4NGTIEx44dkwXMLVu2IDQ0FL169WqZN9ICGvqsXDly5AgAyH622sJn5YrNZoPRaGydP1Mebxknj/j8889FrVYrrlixQjx58qT42GOPieHh4bLO/7Zo3rx54s6dO8Xz58+Lu3btEkePHi1GR0eLV65cEUVRFB9//HGxQ4cO4vbt28WDBw+KQ4YMEYcMGeLlUbeMsrIy8fDhw+Lhw4dFAOJbb70lHj58WLxw4YIoiqK4dOlSMTw8XFy/fr149OhRceLEiWLnzp3Fqqoq6Rzjxo0TBwwYIO7bt0/85ZdfxG7duokPPvigt95Ss3H3WZWVlYl/+ctfxD179ojnz58Xt27dKt58881it27dRIPBIJ2jrXxWTzzxhBgWFibu3LlTzMvLkx6VlZXSMQ393lksFrFPnz7imDFjxCNHjoibNm0SY2JixAULFnjjLTWbhj6rM2fOiK+++qp48OBB8fz58+L69evF5ORk8dZbb5XO0VY+q/nz54s//vijeP78efHo0aPi/PnzRUEQxB9++EEUxdb3M8Ww1Ir961//Ejt06CBqNBrxlltuEffu3evtIXnd/fffLyYkJIgajUZs166deP/994tnzpyR9ldVVYlPPvmkGBERIQYGBor33HOPmJeX58URt5wdO3aIAOo8pk+fLoqiffmAhQsXinFxcaJWqxVHjRolZmVlyc5RWFgoPvjgg2JwcLAYGhoqzpw5UywrK/PCu2le7j6ryspKccyYMWJMTIyoVqvFjh07in/+85/r/EOlrXxWrj4nAGJGRoZ0TGN+77Kzs8Xx48eLOp1OjI6OFufNmyeazeYWfjfNq6HPKicnR7z11lvFyMhIUavVil27dhWfe+45Ua/Xy87TFj6rRx55ROzYsaOo0WjEmJgYcdSoUVJQEsXW9zMliKIoer5eRUREROQf2LNERERE5AbDEhEREZEbDEtEREREbjAsEREREbnBsERERETkBsMSERERkRsMS0RERERuMCwREXmAIAj4+uuvvT0MImoGDEtE5PNmzJgBQRDqPMaNG+ftoRGRH1B5ewBERJ4wbtw4ZGRkyLZptVovjYaI/AkrS0TkF7RaLeLj42WPiIgIAPYpsnfeeQfjx4+HTqdDcnIyvvzyS9nrjx07hpEjR0Kn0yEqKgqPPfYYysvLZcd8+OGH6N27N7RaLRISEvDUU0/J9l+7dg333HMPAgMD0a1bN2zYsEHaV1xcjClTpiAmJgY6nQ7dunWrE+6IqHViWCKiNmHhwoW477778Ntvv2HKlCl44IEHcOrUKQBARUUFxo4di4iICBw4cABr1qzB1q1bZWHonXfewezZs/HYY4/h2LFj2LBhA7p27Sr7Hq+88gr+9Kc/4ejRo7jrrrswZcoUFBUVSd//5MmT+P7773Hq1Cm88847iI6ObrkPgIiarlluz0tE1IKmT58uKpVKMSgoSPZ47bXXRFG03w3+8ccfl70mNTVVfOKJJ0RRFMX3339fjIiIEMvLy6X93377rahQKMT8/HxRFEUxMTFRfOGFF+odAwDxxRdflJ6Xl5eLAMTvv/9eFEVRvPvuu8WZM2d65g0TUYtizxIR+YU77rgD77zzjmxbZGSk9PWQIUNk+4YMGYIjR44AAE6dOoX+/fsjKChI2j9s2DDYbDZkZWVBEARcvnwZo0aNcjuGfv36SV8HBQUhNDQUV65cAQA88cQTuO+++/Drr79izJgxmDRpEoYOHdqk90pELYthiYj8QlBQUJ1pMU/R6XSNOk6tVsueC4IAm80GABg/fjwuXLiA7777Dlu2bMGoUaMwe/ZsvPnmmx4fLxF5FnuWiKhN2Lt3b53nPXv2BAD07NkTv/32GyoqKqT9u3btgkKhQI8ePRASEoJOnTph27ZtNzSGmJgYTJ8+HZ988gnefvttvP/++zd0PiJqGawsEZFfMBqNyM/Pl21TqVRSE/WaNWswaNAgDB8+HJ9++in279+P9PR0AMCUKVOwaNEiTJ8+HS+//DKuXr2KOXPm4OGHH0ZcXBwA4OWXX8bjjz+O2NhYjB8/HmVlZdi1axfmzJnTqPG99NJLGDhwIHr37g2j0YhvvvlGCmtE1LoxLBGRX9i0aRMSEhJk23r06IHMzEwA9ivVPv/8czz55JNISEjAZ599hl69egEAAgMDsXnzZjzzzDMYPHgwAgMDcd999+Gtt96SzjV9+nQYDAb8/e9/x1/+8hdER0dj8uTJjR6fRqPBggULkJ2dDZ1OhxEjRuDzzz/3wDsnouYmiKIoensQRETNSRAErFu3DpMmTfL2UIjIB7FniYiIiMgNhiUiIiIiN9izRER+j90GRHQjWFkiIiIicoNhiYiIiMgNhiUiIiIiNxiWiIiIiNxgWCIiIiJyg2GJiIiIyA2GJSIiIiI3GJaIiIiI3GBYIiIiInLj/wf+ryIaxpISlQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import openfermion\n",
    "import openfermionpyscf\n",
    "from openfermion.transforms import jordan_wigner, get_fermion_operator\n",
    "\n",
    "import timeit\n",
    "\n",
    "\n",
    "import cudaq\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize\n",
    "import numpy as np\n",
    "\n",
    "# GPU\n",
    "cudaq.set_target(\"nvidia-fp64\")\n",
    "# CPU\n",
    "#cudaq.set_target(\"qpp-cpu\")\n",
    "\n",
    "# 1- Classical pre-processing:\n",
    "\n",
    "geometry = [('O', (0.1173,0.0,0.0)), ('H', (-0.4691,0.7570,0.0)), ('H', (-0.4691,-0.7570,0.0))]\n",
    "basis = '631g'\n",
    "multiplicity = 1\n",
    "charge = 0\n",
    "ncore = 3\n",
    "norb_cas, nele_cas = (4,4)\n",
    "\n",
    "molecule = openfermionpyscf.run_pyscf(openfermion.MolecularData(geometry, basis, multiplicity,charge))\n",
    "\n",
    "molecular_hamiltonian = molecule.get_molecular_hamiltonian(\n",
    "        occupied_indices=range(ncore), active_indices=range(ncore, ncore + norb_cas))\n",
    "\n",
    "fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian)\n",
    "qubit_hamiltonian = jordan_wigner(fermion_hamiltonian)\n",
    "\n",
    "spin_ham = cudaq.SpinOperator(qubit_hamiltonian)\n",
    "\n",
    "# 2- Quantum computing using UCCSD ansatz\n",
    "\n",
    "electron_count = nele_cas\n",
    "qubit_count = 2*norb_cas\n",
    "\n",
    "@cudaq.kernel\n",
    "def kernel(qubit_num:int, electron_num:int, thetas: list[float]):\n",
    "        qubits = cudaq.qvector(qubit_num)\n",
    "\n",
    "        for i in range(electron_num):\n",
    "                x(qubits[i])\n",
    "\n",
    "        cudaq.kernels.uccsd(qubits, thetas, electron_num, qubit_num)\n",
    "\n",
    "parameter_count = cudaq.kernels.uccsd_num_parameters(electron_count,qubit_count)\n",
    "\n",
    "# Define a function to minimize\n",
    "def cost(theta):\n",
    "\n",
    "        exp_val = cudaq.observe(kernel, spin_ham, qubit_count, electron_count, theta).expectation()\n",
    "\n",
    "        return exp_val\n",
    "\n",
    "\n",
    "exp_vals = []\n",
    "\n",
    "def callback(xk):\n",
    "        exp_vals.append(cost(xk))\n",
    "\n",
    "# Initial variational parameters.\n",
    "np.random.seed(42)\n",
    "x0 = np.random.normal(0, 1, parameter_count)\n",
    "\n",
    "# Use the scipy optimizer to minimize the function of interest\n",
    "start_time = timeit.default_timer()\n",
    "result = minimize(cost, x0, method='COBYLA', callback=callback, options={'maxiter': 300})\n",
    "end_time = timeit.default_timer()\n",
    "\n",
    "print('UCCSD-VQE energy =  ', result.fun)\n",
    "print('Pyscf-CCSD energy = ', mycc.e_tot)\n",
    "print('Total number of qubits = ', qubit_count)\n",
    "print('Total number of parameters = ', parameter_count)\n",
    "print('Total number of terms in the spin hamiltonian = ',spin_ham.get_term_count())\n",
    "print('Total elapsed time (s) = ', end_time-start_time)\n",
    "\n",
    "plt.plot(exp_vals)\n",
    "plt.xlabel('Epochs')\n",
    "plt.ylabel('Energy')\n",
    "plt.title('VQE')\n",
    "plt.show()"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}