vqe_dcs_qiskit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "adfa7d5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Installation of the requirements\n",
    "#!python -m pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "74dcd404",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SparsePauliOp(['IIII', 'IIIZ', 'IIZZ', 'IZII', 'ZZZI', 'IIZI', 'IZIZ', 'ZXIX', 'IXZX', 'ZXZX', 'IXIX', 'ZZZZ', 'IZZZ', 'ZZIZ', 'ZIZI'],\n",
      "              coeffs=[-0.81054798+0.j,  0.17218393+0.j, -0.22575349+0.j,  0.17218393+0.j,\n",
      " -0.22575349+0.j,  0.12091263+0.j,  0.16892754+0.j,  0.0452328 +0.j,\n",
      " -0.0452328 +0.j, -0.0452328 +0.j,  0.0452328 +0.j,  0.16614543+0.j,\n",
      "  0.16614543+0.j,  0.17464343+0.j,  0.12091263+0.j])\n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "(C) Renata Wong 2023\n",
    "\n",
    "Qiskit code for testing fidelity of derandomised classical shadow on the ground state energy of molecules.\n",
    "\n",
    "Procedure:\n",
    "1. Derandomize the molecule-in-question's Hamiltonian.\n",
    "2. Choose a variational ansatz with initial parameters selected at random.\n",
    "3. Apply the derandomized Hamiltonian as basis change operators to the ansatz.\n",
    "4. Measure the ansatz in the Pauli Z basis and store the results as a shadow.\n",
    "5. Obtain the expectation value of the molecular Hamiltonian from the shadow.\n",
    "6. Optimize for minimum Hamiltonian expectation value. \n",
    "7. Feed the calculated angles/parameters back to the ansatz.\n",
    "8. Repeat steps 3-7 till the optimization is completed. \n",
    "9. Output the minimized expectation value of the molecular Hamiltonian and the mean-square-root-error. \n",
    "\n",
    "Note: Below we perform calculations on the molecular Hamiltonian of H_2.\n",
    "To perform calculations on other molecules, you will need to specify their geometry, charge and spin \n",
    "to replace the values in the driver. \n",
    "\n",
    "Note: predicting_quantum_properties module comes from https://github.com/hsinyuan-huang/predicting-quantum-properties\n",
    "'''\n",
    "\n",
    "from qiskit import QuantumCircuit, execute\n",
    "from qiskit_aer import QasmSimulator\n",
    "\n",
    "from qiskit_nature.units import DistanceUnit\n",
    "from qiskit_nature.second_q.drivers import PySCFDriver\n",
    "from qiskit_nature.second_q.mappers import BravyiKitaevMapper\n",
    "\n",
    "from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver\n",
    "from qiskit_nature.second_q.algorithms import GroundStateEigensolver\n",
    "\n",
    "from qiskit.circuit.library import EfficientSU2\n",
    "from qiskit.algorithms.optimizers import SLSQP, COBYLA, SPSA\n",
    "\n",
    "import time\n",
    "import numpy as np\n",
    "from functools import partial\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import Counter\n",
    "\n",
    "from predicting_quantum_properties.prediction_shadow import estimate_exp\n",
    "from modified_derandomization import modified_derandomized_classical_shadow\n",
    "\n",
    "# handle deprecation issues\n",
    "import qiskit_nature\n",
    "qiskit_nature.settings.use_pauli_sum_op = False\n",
    "import h5py\n",
    "H5PY_DEFAULT_READONLY=1\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# classically obtained ground state energy\n",
    "EXPECTED_EIGENVALUE = -1.86\n",
    "\n",
    "# specifying the geometry of the molecule in question\n",
    "driver = PySCFDriver(\n",
    "    atom=\"H 0 0 0; H 0 0 0.735\",\n",
    "    basis=\"sto3g\",\n",
    "    charge=0,\n",
    "    spin=0,\n",
    "    unit=DistanceUnit.ANGSTROM,\n",
    ")\n",
    "\n",
    "\n",
    "problem = driver.run()\n",
    "hamiltonian = problem.hamiltonian\n",
    "\n",
    "# electronic Hamiltonian of the system\n",
    "second_q_op = hamiltonian.second_q_op()\n",
    "\n",
    "# The Bravyi-Kitaev repserentation of the Fermionic Hamiltonian\n",
    "mapper = BravyiKitaevMapper()\n",
    "bkencoded_hamiltonian = mapper.map(second_q_op)\n",
    "\n",
    "print(bkencoded_hamiltonian)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0b205cb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_hamiltonian(hamiltonian, derandomize = False):\n",
    "\n",
    "    hamiltonian_observables = []\n",
    "    hamiltonian_coefficients = []\n",
    "\n",
    "    for observable in hamiltonian.paulis:\n",
    "        op_list = []\n",
    "        for op_index, pauli_op in enumerate(observable):\n",
    "            pauli_op = str(pauli_op)\n",
    "            if pauli_op == 'I' or pauli_op == 'X' or pauli_op == 'Y' or pauli_op == 'Z':\n",
    "                op_list.append((pauli_op, op_index))\n",
    "        hamiltonian_observables.append(op_list)\n",
    "    \n",
    "    hamiltonian_coefficients = hamiltonian.coeffs.real\n",
    "    \n",
    "    system_size = len(hamiltonian_observables[0])\n",
    "\n",
    "    # removing all occurrences of Pauli-I, for all-Pauli-I there is an empty list left\n",
    "    # these observables are needed for estimate_exp()\n",
    "    observables_xyze = []\n",
    "\n",
    "    for observable in hamiltonian_observables:\n",
    "        XYZE = []\n",
    "        for pauli in observable:\n",
    "            if pauli[0] != 'I':\n",
    "                XYZE.append(pauli)\n",
    "        observables_xyze.append(XYZE)   \n",
    "    \n",
    "    # derandomisation procedure requires that coefficients are non-negative\n",
    "    if derandomize == True:\n",
    "        absolute_coefficients = [abs(coeffcient) for coeffcient in hamiltonian_coefficients]\n",
    "\n",
    "        # removing the empty list as well\n",
    "        # these observables are needed for derandomisation procedure\n",
    "        observables_xyz = []\n",
    "        for idx, observable in enumerate(observables_xyze):\n",
    "            if observable:\n",
    "                observables_xyz.append(observable)\n",
    "            else:\n",
    "                absolute_coefficients.pop(idx)\n",
    "        \n",
    "        return system_size, observables_xyze, observables_xyz, hamiltonian_coefficients, absolute_coefficients\n",
    "    \n",
    "    return system_size, observables_xyze, hamiltonian_coefficients\n",
    "\n",
    "\n",
    "\n",
    "# process the Hamiltonian to obtain properly formatted data\n",
    "hamiltonian_data = process_hamiltonian(bkencoded_hamiltonian, derandomize = True)\n",
    "system_size, observables_xyze, observables_xyz, hamiltonian_coefficients, absolute_coefficients = hamiltonian_data\n",
    "\n",
    "#print('HAMILTONIAN\\n', observables_xyze)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "265d3af8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAADWCAYAAABIdCwuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvXElEQVR4nO3deVxUZf//8dfMALILiEqCC8hSopDiiguYlubXr+JtWmreaX7TwF+5ttyhlWHaz6y4v7/b7G4x7rtU7iRLM80WAffKShM3VFxSMfcFxQWG3x8TKAIzA8zMmXP6PB+PeSRnzlzzhs65PnOuuc45urKysjKEEEIIDdMrHUAIIYSwNyl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNM9F6QDObt86uHxKmff2aQJR99XttUrlrk9m4TiyXdeObNfqJ8XOgsun4MIxpVPUnlpzC8dQ6/ah1txCeTKMKYQQQvOk2AkhhNA8GcYUqnf1Oly+DjrA1wPcXZVOJIRwNlLshOqUlcHhM7B5Pxw8Beeu3HpOBwT6QGQQdI+AZv6KxRRCOBEpdjYwbWEie45swWBwRa83EOQfysg+qSTEDlM6mllqzP37RcjcCofOVP98GXD6sumxaT9EB8OwzuDn6dCYmqDG7QPUm1vYlxQ7GxnVdyaj+s6gtLSEFZv/wdwlIwkPbk9wYLjS0cxSU+6tByDrRygxWv+aXcehYBWM7g5tgu2XTavUtH3cTq25hf3IBBUbMxhceLDLE5QaSzh4YrvScazm7Lk35UPm97UrdOWKb8L7uZAnU9brzNm3j5qoNbewPSl2Nnaz5AarNi8EICQwUuE01nPm3IfPmI7ozEkfZXrUxFgGH22Cs0W2zfZn4czbhzlqzS1sT4qdjSz57lWSZvox8AUPPlw7g6nD3iesWQwAx88cICU9jpslNwD4JOd1Mta+qGTcCuZyz1k8kq27V1Ws+1JGEtv2fe3QfCWlsGSL6bu4+rpeYvq+r8wWjdmIM2WpjmzXQiucutgZjUbmz59PREQE7u7uxMbGkpubS1RUFOPHj1c6XiUj+6TyedoFsl4+Q+e7B7DjQHbFc8GB4fRoN5TMdXMpPHeInO2ZjOyTqmDaW8zlTh6cTsbamRRfL2LDzuV4uTekY9QDDs338xE4dcl27e3/HQpO2669uvj9InzyAzz/CUxZAi8thzW/QtE1ZXNVR7ZroRVOXezGjRtHWloaEyZMYM2aNQwfPpwRI0ZQUFBAXFyc0vGq5ePpz9Rh7/P93i/ZnLeiYvnwxGfYumcVcxaPIHlQOm4uDRRMWVV1uf29mzCkxyQWrHiaJd/N5slBbzk816Z8dbRprT0n4PXVsGU/XLtpWnaxGNbuNC0/c1m5bObIdi3UzmmL3dKlS8nIyGDlypVMnz6d3r17k5qaSrdu3SgpKaFDhw5KR6yRr2cAQ3tOZdFXL2A0mmZUuBhcaRfWi6Li87QN7aFwwupVl7tfpzEcO51PUven8fUMcGieq9fhyFnbt7u3UJnhw6JrsGg9lBqrH5a9dA0+WO+8Q5uyXQs1c9piN2fOHPr3709CQkKl5eHh4bi6uhITYxp/P3z4MAkJCURGRtKuXTs2bNigRNwqhvScxLlLhXzz078BOHxyF7sOb6J9eF9Wf/+ewulqdmdugGaNwhWZsn3svH3avXqj8onojrL1INwsrfn7x7IyKLyg/DCrObJdC7XSlZU53+fIY8eO0bx5cz744AMef/zxSs+NGDGCvXv38ssvvwDQr18/Bg8eTEpKCps3b2bYsGEcOnQINzc3s++h0+msyjL/yWxiWyfW6fcoZzQamfZOAsmD0gkJjGTSgnjmjf8Wf5+mZl+342AO09/pXaf3tEXucvMyxzCgy/9Y9cm9Ppnv1KbXWO4fv6jSMnMzLs2ZvLjyz5++2ptje3Lq1lgdDU3NJjiqFzp9zZ8xy8qMbPviNTZ/Yt/vvmS7Vm67FrZlbQlzyiO7Y8dMJ0QFBQVVWl5cXExubm7FEOaZM2fYuHEj48aNAyA+Pp5mzZqRnZ2NM/liy0IiguOIDInD092HMf3SeHvlZKVjOT2d3qDKtmticGlgttABUFaG3mD+g5qzkO1aqIlTHtkdOHCAiIgI3nrrLSZPnlyxfNasWbz88sssWLCAlJQUfv75Zx566CEKCgoq1hk+fDh9+/a12WzNbZnK3T/LLwQ6PlK31yqVuz6Z77TjKHxo5ah0+RHfnUdwNZn2IDR38Fc1mVvh+4OWT6MY0RW6tLZvFtmua8eW27VQhlNeLiwsLIyYmBjmzJlDQEAAwcHBZGVlsXr1agCnnYkpbCvYThdxNujhrob2aduc7hGm7+3MaeAC7Vs6Jo8QfyZOOYyp1+tZtmwZ0dHRJCcnM3bsWAIDA5k4cSIGg6FickqLFi34/fffuX79esVrDx06RMuW0ltoQSNv8HW3fbvNA8DF8aOYNG8EPWu4iEf5N8jDOoObU34EFULdnLLYAURGRpKdnc2VK1c4evQoaWlp7Ny5kzZt2uDh4QFAYGAg3bt354MPPgBg8+bNHD9+nN695YtkLdDpoKsdJsvZo01rDekIA+8Fzzu+lmvkDY/3go6hisQSQvNU9Rly27ZtdO3atdKyd955hzFjxpCeno6bmxtLly61OBNTqEd8BKzbY7psmC34uEMHBQ/89TroGw2Jd8P0TNOyp++H0Mam4i6EsA/VFLuioiLy8/NJSUmptDwsLIz169crlMrkg9V/Y9fhTUS36k5I4ygys+cyeei7xLZO4JOc19m8awVN/VvyzMMZ3Cy5zrPv9iW4UTjPj/zYaXMDbNi5nIUrJrFkxm8UXy9SJLefJwyMhc9/tk17zjJMePswalgT5XKYU9P20S60J/My/8qpC0dxMbiR+mgmbi7uim7XNWWNbhXP1Ld7cejkTt6Zsp3gwHAKzx1iXuZf0aEjsGEIz434CIPewIxFAykqvkD6xI0Ozy/sz2mHMe/k7e1NaWkpTz31lNJRKjl0Mo8r1y7xZsp6Ll09y7UbVxiW8AyxrRM4X3SK7QezSZ+4kdC7YtiU9zkeDbxJHZWpdGyzuctt+DWLxn7NARTN3SsK7mlmfp3Jiy3PxIyPgJjmtsulZea2j4MntuPi4sabKevp12ks3/28WNHtw1xWg96FWWM+p2e7hyrW93b3Y/bYVbyZsp6ggFB+2Gua+Db78VU1vYXQANUUO2eVd2gjHSNNF5HtEHE/+tvO38r/bRuxYYl/PNeXPUe2KBGxWuZyA3y/ZzUdIvqi0ym/iej1MLYn3H1X3dvoHAYPdbRdJq0zt30ENgzGaDSNKxcVX8DXs5EiGcuZy6rT6aqc5O7j6Y+Xh2k6rovBFb1OgdlKwuGcYEBH3S5fPceqLe/w6Ya3KCq+QELscPy8TeNSV65dwNPdFwAv94YUXbugYNLKzOUG+Oanf/HcIx/x1Y+LzLTiOG4u8EQirNsNX+00XV/SGu6ukNTBdN6afCdmPXPbh69XINdvFvP46/dg0Bn4f0//4LRZzTlz8QQ/5X/DqD4zHJBSKE2KXT35eAbwWL9XiI8exNbdqzh98dYZr17uDTn9xxmwV69dwtvdT6GUVZnL/cuBdbRp2Q1XF+ea6GPQw/1tIbYFbNgHPxSY7lFXHe8GplmXPSJN3/uJ2jG3ffyU/zUNvRqz6Jk9rP81i2W58xl9v3L3sTOXtSY3Sq7z+n8eY+qw9zAYpBv8M1B+jErl2ob2YGeBaYLMjoM5FcM7AJHNO/FrQS4AP+//lntadq22DSWYy334ZB5bdq3kb+/158jvu/jwK+f65NvEF4Z2grShMLkfDL1tePLhLqaro7w8xDTFXwpd3ZjbPsrKyiruFNDQK5Ar1y4qkrGcuaw1Sc8az6D4ibRs2sbe8YSTkGJXT6FBbXExuDJtYSIuBlfc3bwqnvP3bkK7sF5MXtCDgye2Ex+dpFzQO5jLPaTH07z+5DrmPvEVLZtGM7b/bAWT1szNBVoFQs+oW8u6hSt30riWmNs+OkY+wNFTe5i2MJF/rX2RQd1SzLRkf+ayAqR9NJyf9n/NvMzH2Jy3gt2Ht7AxbznLN6QzbWEiG3d+plBy4Uhy/G4D4wbMrfj3+l+zyMx+jeDACGJbJ/BI7+d4pPdzFc8XXy/itaWPEtW8kxJRKzGXu1z5NGxnyi0cw9z28eJfsyqtq/T2YS7rzNGfVFl/5eyqd8mdsWggAb71mAUlnJpTXgjamcgFc2tHqQvmlp92UNdbAClBycyyXdeOXAha/eTIzgIfBU/4rc97K5Vbyb+XsJ5s1+p4X2E7UuwsiLpP6QR1o9bcwjHUun2oNbdQnkxQEUIIoXlS7IQQQmieFDshhBCaJ8VOCCGE5kmxE0IIoXlS7IQQQmieFDshhBCaJ8VOCCGE5kmxE0IIoXlS7IQQQmieFDshhBCaJ8VOCCGE5kmxE0IIoXly1wML9q2Dy6eUeW+fJnW/yrtSueuTWQhL1Lhdq7UP0RopdhZcPqXcTS7rQ625hTBHjdu1GjNrkQxjCiGE0DwpdkIIITRPhjGFcKBrN2FvIfx29tayjA3QzA9aBEJEUzDIR1AhbE6KnQ1MW5jIniNbMBhc0esNBPmHMrJPKgmxw5SOZpZac6vRhavwdR5sOwQ3Sio/t/2o6QHQ0APiI6D3PeAme2edqHG7VmNmtZHdyUZG9Z3JqL4zKC0tYcXmfzB3yUjCg9sTHBiudDSz1JpbTX4ogOXbTEd1llwshjW/moriyG4Q2tj++bRIjdu1GjOriQyY2JjB4MKDXZ6g1FjCwRPblY5jNbXmdnard8CSLdYVutudvgz/+BZ2/mafXH8Watyu1ZhZDaTY2djNkhus2rwQgJDASIXTWE+tuZ3Z+n2mocu6KjVCxkYoUOgcLS1Q43atxsxqIMXORpZ89ypJM/0Y+IIHH66dwdRh7xPWLAaA42cOkJIex82SGwB8kvM6GWtfVDJuBXO55yweydbdqyrWfSkjiW37vlYqqqr8fhFW/mx+nfRRpoc5pUZYsrXq93zCPDXuj7Iv2pdTFzuj0cj8+fOJiIjA3d2d2NhYcnNziYqKYvz48UrHq2Rkn1Q+T7tA1stn6Hz3AHYcyK54LjgwnB7thpK5bi6F5w6Rsz2TkX1SFUx7i7ncyYPTyVg7k+LrRWzYuRwv94Z0jHpAwbTqsfwnKDHapq0zl2HdHtu0VR83S03fJX7xi2l4tuAUlJUpnap6atwfZV+0L6cuduPGjSMtLY0JEyawZs0ahg8fzogRIygoKCAuLk7peNXy8fRn6rD3+X7vl2zOW1GxfHjiM2zds4o5i0eQPCgdN5cGCqasqrrc/t5NGNJjEgtWPM2S72bz5KC3FE6pDr9fgn2Ftm1z837TUZ5Sth+BFz+FjzfDd7tNw7P/+w3MX2Mqxs5Kjfuj7Iv24bTFbunSpWRkZLBy5UqmT59O7969SU1NpVu3bpSUlNChQwelI9bI1zOAoT2nsuirFzAaTT2Ui8GVdmG9KCo+T9vQHgonrF51uft1GsOx0/kkdX8aX88AhROqw0+HbN/mpWLIP2n7dq2x65jpu8PiaibZnLgA/+8buHzN4bGspsb9UfZF23PaYjdnzhz69+9PQkJCpeXh4eG4uroSE2May37xxReJjIxEr9eTlZWlRNRqDek5iXOXCvnmp38DcPjkLnYd3kT78L6s/v49hdPV7M7cAM0ahcv051o4ctbyOs7UrjllZbDiF9CZef5iMWzMd2isWlPj/ij7om055Xl2x44dIy8vjylTplR57ujRo0RHR9OggWnYoX///owZM4bHH3/c0TErvJGcU2WZl7svy185B5i+e/z78id5asgCQgIjmbQgnvjowfj7NHVw0sos5RZ1U3jBTu2et0+75hw5C6cuWV5v8354MMb+eayhxv1R9kX7c8oju2PHTJcIDwoKqrS8uLiY3NzcSkOY8fHxhIWF1fo9dDqdVY/c3Jx6/S4AX2xZSERwHJEhcXi6+zCmXxpvr5xs8XW5uTlW57RH7rqoT+b6PMop8d53Pk6fq1wdymdd3vmw9vlyK778yuG/y/0DHrLq//vla6DT6+2eR6n9UY37Yn1zq+VhLac8sgsMDAQgPz+fAQMGVCyfN28ehYWFTjs5pSaDu0+s9HP3tkl0b5ukTJg6ePaRDKUjqIrxjyntamnXnBvFF61a7+b1K847NfMOat4fZV+sO6csdmFhYcTExDBnzhwCAgIIDg4mKyuL1atXA9ik2JVZuWNuy1TuXlQJCYmULaxbB6JU7vpkro/Ji03/tfb/qz2lr4XDZ279XJ6tynqjzD9/p3GjBrHyDcf+fiWl8NJncOV6zevogJ5tvFjggL+9GrdrtfYhWuOUw5h6vZ5ly5YRHR1NcnIyY8eOJTAwkIkTJ2IwGCompwjhjJrbaaKcvdo1x8UAve+p+XkdoNdDwt0OiyREnTjlkR1AZGQk2dnZlZaNHj2aNm3a4OHhoVAqISyLaQEbbDw70dUAdzezbZvWuq+N6a4NG/NNxe324wSDHsb2hLv8lMkmhLWc8siuJtu2basyhDlz5kxCQkLYsmULEyZMICQkhIMHDyqUUAgIbwJNfW3bZlwr8HSzbZvW0uvgoU4wtT90um0u2IAYeDEJokOUySVEbTjtkd2dioqKyM/PJyUlpdLytLQ00tLSFEpl8sHqv7Hr8CaiW3UnpHEUmdlzmTz0XaJbxTP17V4cOrmTd6ZsJzgwnOLrRTz7bl+CG4Xz/MiPnTJ3bOsEBs9sSHiz9gC89NhyfD0DmLFoIEXFF0ifuFHR3M5Op4PBHeDdHNu018AF+rWzTVv10aKR6bZDPxSYfn7ACTJVx9x2/Y/Pn+JQ4U7uahTGlIfew6A3KLpd16bvuHjlDC9+OAiDwRUv94bMePQ/GI2lTtOfODvVHNl5e3tTWlrKU089pXSUSg6dzOPKtUu8mbKeS1fPcu3GFYYlPENs6wQMehdmjfmcnu1uTd/2aOBN6qhMBRObmMsNEBrUjjeSc3gjOafiag2zH19lrklxmzbB0Ln2Z8RUKykO/L1s05bWmduu9/32IyUlN3gjOYeWTaP5/o8LKyu1Xde27/D28OetlI28mZxLZHAcW3evcpr+RA1UU+ycVd6hjXSMNF2QtUPE/ej1horndDqd4ieO18RcboCjp/Yw5e2evL/6eaeY4ahGwzpDZFDNz09ebHkm5n1toGtr2+bSMnPbdeHZAkLvMk1ua93sXnYd2axIxnK17TsMegN6vanLLi0rJTgwwnFhNUCKXT1dvnqOf339EtMWJrLku1e5fFUdVzywlDvjuf28mbyeoqvn2bL7C4VSqpurAZ5IrNsRnkFvGgr973tNw6LCOua265DGUfxakAvA9gPruFJ8QaGUJnXpO/Ye/YGUv3dk+4F13BUQ6oCU2qGa7+yclY9nAI/1e4X46EFs3b2K0xcVOqGmlizlLh+6jG+bxIHjvxAfPUiJmKrnajB9zxXbAlb8bN2ltyKawl86ygzHujC3XYcH30uroLZMf6c3rYLa4qfwqEtd+o67W3Tm7UnbWJb7Bl/9uIihvapeUlFUT47s6qltaA92FqwHYMfBHIzGUoUTWcdc7uIbVyj94+ddhzfRrJGMo9VXdDD8bSBM7AM9IqBlI/BwMxVD7wam4c4+beD5gTCxrxS6urK0P46+/0XmP5mNr2cjutzzX0pErFDbvuPmbVfQ8XL3xc1VTsGqDTmyq6fQoLa4GFyZtjCRNi274e7mRanx1m2l0z4aTt7hjRw/s5+HE58lvu1gBdPeYi738dP7eWPZ43i4eRMUEMpfH5ilcFpt0OkgIsj0EPZhbrs2Go0888/70OsNtA/vwz0tujhtVqjadwT43sW7Xz6DXqfHxyOA50Z8pGB69ZFiZwPjBsyt+Pf6X7PIzH6N4MAIYlsnMHP0J5XWLb5exGtLHyWqeSdHx6zCXO6Fk3+usv6MRQMJ8L3LkRGFqDVz23V1dxdQcruuTd8B8GZybqWfnak/cXa6MplqZ5aS17XzC4GOj9TttUrlrk/m+iif1Vjd3QKEbSn5t1bjdq3WPkRr5MjOAp8m6nxvpXIr+fcS2qfG7VqtfYjWSLGzIOo+pRPUjVpzC2GOGrdrNWbWIpmNKYQQQvOk2AkhhNA8KXZCCCE0T4qdEEIIzZNiJ4QQQvOk2AkhhNA8KXZCCCE0T4qdEEIIzZNiJ4QQQvOk2AkhhNA8KXZCCCE0T4qdEEIIzZNiJ4QQQvPkrgcW7FsHl08p894+Tep+xXSlctcnsxBapMY+RI2ZLZFiZ8HlU8rdeLE+1JpbCK1R476oxsyWyDCmEEIIzZNiJ4QQQvOk2AkhrFJqhBPnb/18tgjKypTLI0RtyHd2Qoga3SyFX47A1gNw9CyUGG89l7YCvBpAVBD0iITQxqDTKZdVCHOk2NnAtIWJ7DmyBYPBFb3eQJB/KCP7pJIQO0zpaGapNbdwjN3H4T/fw8Ximte5ch1+PmJ6RAbBI10gwNtxGbVCrfuimnLLMKaNjOo7ky9eLWL5y2fp12kMc5eM5PiZA0rHskituYX9GMtg+TZ4N8d8obtT/kn4v1/CLo3N4nMUte6Lasktxc7GDAYXHuzyBKXGEg6e2K50HKupNbewrbIyyPoB1u+r2+uvl8AH62HXcdvm+jNR677o7Lml2NnYzZIbrNq8EICQwEiF01hPrbmFbf14CDZb+FCePsr0qImxDD7aBBeu2jbbn4Va90Vnz+3Uxc5oNDJ//nwiIiJwd3cnNjaW3NxcoqKiGD9+vNLxKlny3askzfRj4AsefLh2BlOHvU9YsxgAjp85QEp6HDdLbgDwSc7rZKx9Ucm4FczlnrN4JFt3r6pY96WMJLbt+1qpqMLOLhXDZ9ts09a1m/DJ97Zpyxau3YTLxaYZpc5K+hD7cupiN27cONLS0pgwYQJr1qxh+PDhjBgxgoKCAuLi4pSOV8nIPql8nnaBrJfP0PnuAew4kF3xXHBgOD3aDSVz3VwKzx0iZ3smI/ukKpj2FnO5kwenk7F2JsXXi9iwczle7g3pGPWAgmmFPW3aD8U3bdfe7hNw/Lzl9expzwlY8C08/wnMXA6pWfDZT7X7LtJRpA+xL6ctdkuXLiUjI4OVK1cyffp0evfuTWpqKt26daOkpIQOHTooHbFaPp7+TB32Pt/v/ZLNeSsqlg9PfIate1YxZ/EIkgel4+bSQMGUVVWX29+7CUN6TGLBiqdZ8t1snhz0lsIphb2UGmHLftu3u8kObVordy/8MxsO3HaNx2s3TcvfXGM6T9AZSR9iH05b7ObMmUP//v1JSEiotDw8PBxXV1diYmI4f/48AwcOJDIyktjYWB544AEOHFB+FpCvZwBDe05l0VcvYDSaxk1cDK60C+tFUfF52ob2UDhh9arL3a/TGI6dziep+9P4egYonFDYy+8X4dI127e7/6Tt27TGsXOmIzio/sT3S8Xw8WbHZqoN6UNszymL3bFjx8jLy2PYsKrnahw9epTo6GgaNGiATqdj8uTJ5Ofns2PHDgYOHMjYsWMVSFzVkJ6TOHepkG9++jcAh0/uYtfhTbQP78vq799TOF3N7swN0KxROMGB4QqmEvb22zn7tHv6MhTfsE/b5mzMB3Pnt5cBh04rP8xqjvQhtuWUJ5UfO2Y6UScoKKjS8uLiYnJzc3nwwQcB8PPzo2/fvhXPx8fHM2/ePKveQ2flpR7mP5lNbOtEs+u8kZxTZZmXuy/LXzH1IEajkb8vf5KnhiwgJDCSSQviiY8ejL9PU7Pt5ubm0GlEb6ty2iN3XdQnc31M+tj08d3a/6+iss5JM+j2UFqlZeZmXJp7fvLiyj8Hh97D+RN765Gu9sa8VUDDxqEW1xs4YhLb1/6vXbOosQ+xJrM1ueuitpnLrLxmnVMe2QUGBgKQn59fafm8efMoLCyscXJKeno6SUlJ9o5Xa19sWUhEcByRIXF4uvswpl8ab6+crHQsIW5jvw8JOju2XRO93mDVejor11Oa9CH1pyuztiw6kNFopH379hQWFjJ//nyCg4PJyspi9erVHD16lK1bt9KlS5dKr5k1axZr1qxh3bp1eHp62izLtkzl7uvkFwIdH6nba5XKXZ/M9VF+NGHpaERUb1M+LPvRunXL/8Z3HsHVZPZQ8HavW666+nA9/HrM8oWq/09fCDd/cFRvauxD1JjZEqc8stPr9Sxbtozo6GiSk5MZO3YsgYGBTJw4EYPBQExMTKX1Z8+ezapVq/jqq69sWuiE+LMIsdO8AT9Pxxc6gO6R5gudTgdNfKF1E8dlEspyyu/sACIjI8nOzq60bPTo0bRp0wYPD4+KZbNmzWL16tV88803+Pn5OTilENrQzB/cXU1T821JqWIS0dR0J4aN+VWf0+nARQ+PxstdGv5MnLbYVWfbtm107dq14uddu3bx8ssv07p1axITEyuWb9++3fHhhFAxVwN0CoMNdbwmZk26KTQBT6eDoR2hsQ9k76l86bK774KB90KwvzLZhDJUU+yKiorIz88nJSWlYll0dLTVM3GEEOb1ioTN+213Sa3mAcoOE+p0kHA39IyEqUtNy15KAn8v5TIJ5aim2Hl7e1NaWqp0jGp9sPpv7Dq8iehW3QlpHEVm9lwmD32XJv4tmJf5V3ToCGwYwnMjPsKgNzBj0UCKii+QPnGjU+aObZ3AT/nfkLluLsYyIxP++w2aN47i2Xf7EtwonOdHfqxobmEfjX2hfzv4ckf92zLoYURX5xgm1N82M8HZCl1N+2B0q3imvt2LQyd38s6U7RXnqH2S8zqbd62gqX9Lnnk4g5sl1x2+X9Y2M8CGnctZuGISS2b8RvH1IkX6EqecoKImh07mceXaJd5MWc+lq2e5duMKwxKeIbZ1At7ufsweu4o3U9YTFBDKD3tXAzD78VUWWrU/c7mv3yzmy63/5LXx3/BGcg6RIXF4NPAmdVSm0rGFnd3XxjTMZ87kxZZnYg6JM30PKGpmbh806F2YNeZzerZ7qGL980Wn2H4wm/SJGwm9K4ZNeZ87fL+sbeZyG37NorFfcwDF+hIpdvWUd2gjHSNNFzbtEHF/pfN7fDz98fJoCJgu9aPXOc85PeZy7z6yBZ1OzwvvP8hrS0dTfOOKUjGFgxn08HgvaNOsbq/XYSp0PZzvDi9Ox9w+qNPpqpwwnv/bNmLDEv9Yvy97jmxxWNZytc0M8P2e1XSI6ItOp2y5kWJXT5evnuNfX7/EtIWJLPnuVS5frXrlgDMXT/BT/jcVG4kzMJf7/OXfOXe5kDn/s4bolvF8ueWfCiYVjubmAv+TAEkdwKUWn88CfeCp+03fkwnLrOk7bnfl2gU83X0B8HJvSNG1Cw5IWVltMwN889O/6NPhUQekM08139k5Kx/PAB7r9wrx0YPYunsVpy9WPhPzRsl1Xv/PY0wd9h4Gg/P8uc3l9nJvSNtWPTDoDdwbfh/LcucrmFQoQa+HxHsgtoVp0sqWg1BUw4WiQwKgRwR0aGUqlMI6lvqOO3m5N+T0H2d6X712CW93PwekrKy2mX85sI42Lbvh6uLmoIQ1kyO7emob2oOdBesB2HEwB6Ox8iSa9KzxDIqfSMumbZSIVyNzuaOad+LoqT0AHDyxnaAAy9cYFNrk7wX/dS+88hdI/W94rAcM6wTDO8OE3qaro0x/ELqGS6GrLUt9x50im3fi14JcAH7e/y33tOxqdn17qG3mwyfz2LJrJX97rz9Hft/Fh1/NcETMasnmWU+hQW1xMbgybWEibVp2w93Ni1JjCQC7D29hY95yfj9/hOUb0hnSYxI92g1ROLGJudx+3o2JCUtg6tu9aODqyd9GLVE4rVCaXmeardnYV+kk2mFuHwRI+2g4eYc3cvzMfh5OfJb4toNpF9aLyQt60MSvBX/pOdnpMw/p8TRDejwNwOQFPRjbf7bDM5eTYmcD4wbMrfj3+l+zyMx+jeDACGJbJ7By9uUq689YNJAAXwtT3hzAXO6hvaYwtNeUiueLrxfx2tJHiWreSYmoQmiSuX1w5uhPqqz/SO/neKT3cxU/K7Ff1jZzufJTrZTqS5zyQtDORK0XRJULQQtxi5Lbhxr7EDVmtkSO7CzwUfAKEPV5b6VyK/n3EsIZqbEPUWNmS6TYWRB1n9IJ6katuYXQGjXui2rMbInMxhRCCKF5UuyEEEJonhQ7IYQQmifFTgghhOZJsRNCCKF5UuyEEEJonhQ7IYQQmifFTgghhOZJsRNCCKF5UuyEEEJonhQ7IYQQmifFTgghhOZJsRNCCKF5ctcDC/atg8unlHlvnyZ1v/q4Urnrk1kILVJjH6LGzJZIsbPg8inlbmJYH2rNLYTWqHFfVGNmS2QYUwghhOZJsRNCCKF5MowphNCkU5dg13H47eytZf/4Fpr5QatAaBsCbtID/mnI/2obmLYwkT1HtmAwuKLXGwjyD2Vkn1QSYocpHc0steYWwpzDZ2DNDth3supzB343PdbvAw836BYOD7QFd1fH57ydWvdFNeWWYmcjo/rOZFTfGZSWlrBi8z+Yu2Qk4cHtCQ4MVzqaWWrNLcSdSo3w5Q7I3g1lVqxffAPW7YZfjsCobhDe1O4RzVLrvqiW3PKdnY0ZDC482OUJSo0lHDyxXek4VlNrbiHAVOg+3GAqXtYUutudvwILv4Odv9klWq2pdV909txS7GzsZskNVm1eCEBIYKTCaayn1txCAGT9CHn1mCpfWgYZG+HoWcvr2pta90Vnzy3FzkaWfPcqSTP9GPiCBx+uncHUYe8T1iwGgONnDpCSHsfNkhsAfJLzOhlrX1QybgVzuecsHsnW3asq1n0pI4lt+75WKqoQ1dp9HLYcML9O+ijTw5xSIyzZAjdLbZetNqQPsS+nLnZGo5H58+cTERGBu7s7sbGx5ObmEhUVxfjx45WOV8nIPql8nnaBrJfP0PnuAew4kF3xXHBgOD3aDSVz3VwKzx0iZ3smI/ukKpj2FnO5kwenk7F2JsXXi9iwczle7g3pGPWAgmmrKr4BG/bd+vmHArhRolwe4VjGMli+zXbtnbwIm/Jt115tSB9iX05d7MaNG0daWhoTJkxgzZo1DB8+nBEjRlBQUEBcXJzS8arl4+nP1GHv8/3eL9mct6Ji+fDEZ9i6ZxVzFo8geVA6bi4NFExZVXW5/b2bMKTHJBaseJol383myUFvKZyysp8Pw4vL4dPbOrslW+Dlz2D/74rFEg60rxDOFNm2zU37TUVUKdKH2IfTFrulS5eSkZHBypUrmT59Or179yY1NZVu3bpRUlJChw4dlI5YI1/PAIb2nMqir17AaDQC4GJwpV1YL4qKz9M2tIfCCatXXe5+ncZw7HQ+Sd2fxtczQOGEt+wthI82VT/kVHwT/pkNJ847PpdwrJ8O277N05crn5unBOlDbM9pi92cOXPo378/CQkJlZaHh4fj6upKTIxpTDgpKYmYmBjat29P586d+fbbb5WIW8WQnpM4d6mQb376NwCHT+5i1+FNtA/vy+rv31M4Xc3uzA3QrFG4000jXrOj5ufKyqC0FL7b7bg8Qhn2mlDy2zn7tFsb0ofYllOeZ3fs2DHy8vKYMmVKleeOHj1KdHQ0DRqYDuEzMjLw8/MD4JdffiExMZFz585hMBgclveN5Jwqy7zcfVn+immPMRqN/H35kzw1ZAEhgZFMWhBPfPRg/H2UPbHHUm5ndfoyHLHQyZUB24/Aw13kKhlaVWqE05fs03bhBfu0WxPpQ+zPKY/sjh0zzSEOCgqqtLy4uJjc3NxKQ5jlhQ7g4sWL6HQ6ysosD7jrdDqrHrm5OfX+fb7YspCI4DgiQ+LwdPdhTL803l452eLrcnNzrM5pj9x1UZ/M1j5i47pblaW0DAKaBNs9jzyUebh7elU5p6581uWdD2ufL/feBx/aLKca+xCl+o+6ZLaWU37mDQwMBCA/P58BAwZULJ83bx6FhYVVJqdMnDiRNWvWcPHiRT799FNcXJzr1xrcfWKln7u3TaJ72yRlwtTBs49kKB2hkuKiM1atZzSWcP3qBfuGEYox/jENv6ysrFadnjVKblyzaXv1JX1I/enKrDkMcjCj0Uj79u0pLCxk/vz5BAcHk5WVxerVqzl69Chbt26lS5cuVV6Xm5vLlClTWL9+Pd7e3jbJsi1Tufs6+YVAx0fq9lqlctcns7XKymD+GtMElJo2Xh0Q0xzG9rJvFqGsV1eahrUtKT96m7zYunb/0hF6RdU91+3U2IeoMbMlTjmMqdfrWbZsGdHR0SQnJzN27FgCAwOZOHEiBoOhYnLKnRISEtDr9WzatMnBiYUj6XTQv535QqfTQZ9oR6YSSmhup8l99mpXKMe5xvtuExkZSXZ2dqVlo0ePpk2bNnh4eABQVFTE2bNnadmyJWCaoHLw4EHuueceh+cVjtWuuWnyybIfTOdElY9ilZWBqwEe6wktGimbUdhf+5bw8xHbtunvBS1l29Ecpy121dm2bRtdu3at+PnKlSs8/PDDFBUV4eLigru7Ox9//DEtWrRQMKVwlG7hEB0MWw+azovS6aB1E+gUBp5uSqcTjtAmGPw94fxV27XZPQL0TjnmJepDNcWuqKiI/Px8UlJSKpY1bdqUrVu3KpjK5IPVf2PX4U1Et+pOSOMoMrPnMnnou7QKiubFDwdhMLji5d6QGY/+B6OxlGff7Utwo3CeH/mxU+b2cm/Iwj9mep06f4QhPSfxl56TmbFoIEXFF0ifuFHR3Lfz9TDdj0z8ORn0kBRnuuOBLQR6Q08bfVdnjZr2wehW8Ux9uxeHTu7knSnbCQ4Mp6T0ZpVlxdeLHN6f1CYzwOCZDQlv1h6Alx5bjq9ngCJ9iWo+v3h7e1NaWspTTz2ldJRKDp3M48q1S7yZsp5LV89y7cYVhiU8Q2zrBLw9/HkrZSNvJucSGRzH1t2r8GjgTeqoTKVjm80dHnwvbyTn8EZyDqF3xdDlnoEAzH58lYVWhXC82BYQ18r8OpMXW56cotfBiG7QwEGHAOb2QYPehVljPqdnu4cq1q9umaP7k9pmBggNalfRn5RfQUWJvkQ1xc5Z5R3aSMdI04VNO0Tcj15/62R2g96A/o/xkNKyUoIDIxTJWB1zucsV37jC+csnFb/ygRCWPNIVIoMsr1cTHTCym2kY3FHM7YM6na7KCePVLXO02mYGOHpqD1Pe7sn7q5+36hxoe1HNMKazunz1HKu2vMOnG96iqPgCCbHD8fO+tcfsPfoD//tZCm4u7gzrNU3BpJVZyg3w4941dIzqr1BCIaznaoAnEuGzbbDZwu1+7uTjDiO6mr7/cyRr9kFnU5fMGc/tx8fDn79/+iRbdn9BfPQgB6WtTIpdPfl4BvBYv1eIjx7E1t2rOH2x8skpd7fozNuTtrEs9w2++nERQ3tVvQSaEizlBtiU9xnDE59VIJ0QtedqgOFdTMOaq3dYvqScqwE6hcJ/3QteCtxAwJp90NnUJXP50GV82yQOHP9FsWInw5j11Da0BzsL1gOw42AORuOty/CX32gRTNeLc3P1cHi+mpjLDVBSepOjp/bQulmsEvGEqLOou2BKf5jWH/q1gzbNoIkvNPKGEH/oHAbDO8OsIabiqEShA8v7oDOqbebiG1co/WOdXYc30axRa7tnrIkc2dVTaFBbXAyuTFuYSJuW3XB386LUaLp76MET23n3y2fQ6/T4eATw3IiPFE57i7ncAL8cWMe9re9TMKEQ9dO8kenhrCztg2kfDSfv8EaOn9nPw4nPEt92cLXLnDlzE/+WvLHscTzcvAkKCOWvD8xyaN7bSbGzgXED5lb8e/2vWWRmv0ZwYASxrRN4Mzm30rrF14t4bemjRDXv5OiYVZjL3SmqH52i+lVaf8aigQT43uXomEJolrl9cOboT6qsf+cyJfqT2mZeOPnnKsuU6Euc8tqYzkSt14jT8rUxhVATNfYhasxsiRzZWeCj4OSo+ry3UrmV/HsJ4YzU2IeoMbMlcmQnhBBC82Q2phBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELz/j+SiJIgyOB6dAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 568.197x264.88 with 1 Axes>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "VARIATIONAL ANSATZ\n",
    "Note that for molecules other than H_2 you may need to specify a different number of reps.\n",
    "'''\n",
    "\n",
    "reps = 1   \n",
    "ansatz = EfficientSU2(system_size, su2_gates=['rx', 'ry'], entanglement='circular', reps=reps, skip_final_rotation_layer=False)  \n",
    "    \n",
    "ansatz.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d5f6e337",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "COST FUNCTION\n",
    "'''\n",
    "\n",
    "\n",
    "def basis_change_circuit(pauli_op):\n",
    "\n",
    "    #   Generating circuit with just the basis change operators\n",
    "    #\n",
    "    #       pauli_op: n-qubit Pauli operator\n",
    "\n",
    "    basis_change = QuantumCircuit(ansatz.num_qubits, ansatz.num_qubits)\n",
    "    \n",
    "    for idx, op in enumerate(pauli_op):\n",
    "        if op == 'X':\n",
    "            basis_change.h(idx)\n",
    "        elif op == 'Y':\n",
    "            basis_change.h(idx)\n",
    "            basis_change.p(-np.pi/2, idx)\n",
    "\n",
    "    return basis_change\n",
    "\n",
    "\n",
    "\n",
    "def min_cost():\n",
    "    \n",
    "    return min(cost_history)\n",
    "\n",
    "\n",
    "def log_cost(cost):\n",
    "\n",
    "    global cost_history\n",
    "    cost_history.append(cost)\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "def objective_function(operators, params):\n",
    "    \n",
    "    backend = QasmSimulator(method='statevector', shots=1)\n",
    "    \n",
    "    pauli_op_dict = Counter(tuple(x) for x in operators)    \n",
    "    \n",
    "    shadow = []\n",
    "    for pauli_op in pauli_op_dict:\n",
    "        \n",
    "        qc = ansatz.bind_parameters(params)\n",
    "        qc = qc.compose(basis_change_circuit(pauli_op))                                   \n",
    "        qc.measure(reversed(range(system_size)), range(system_size))\n",
    "        result = execute(qc, backend, shots=pauli_op_dict[pauli_op]).result()\n",
    "        counts = result.get_counts()   \n",
    "\n",
    "        # generate a separate circuit for each basis change Pauli operator\n",
    "        output_str = list(list(counts.keys())[list(counts.values()).index(1)])\n",
    "        output = [int(i) for i in output_str]\n",
    "        eigenvals = [x+1 if x == 0 else x-2 for x in output]\n",
    "        snapshot = [(op, eigenval) for op, eigenval in zip(pauli_op, eigenvals)] \n",
    "                \n",
    "        shadow.append(snapshot)\n",
    "    \n",
    "    \n",
    "    cost = 0.0\n",
    "    for term, weight in zip(observables_xyze, hamiltonian_coefficients):\n",
    "        sum_product, match_count = estimate_exp(shadow, term)\n",
    "        if match_count != 0:\n",
    "            exp_val = sum_product / match_count\n",
    "            cost += weight * exp_val\n",
    "            \n",
    "    log_cost(cost)\n",
    "    \n",
    "    return cost\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b963d524",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NUMBER OF OPERATORS | DERANDOMISED OPERATORS | AVERAGE RMSE ERROR\n",
      "\n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "RUNNING EXPERIMENTS\n",
    "''' \n",
    "\n",
    "start_time = time.time()\n",
    "rmse_errors = []\n",
    "\n",
    "print('NUMBER OF OPERATORS | DERANDOMISED OPERATORS | AVERAGE RMSE ERROR\\n')\n",
    "measurement_range = [50, 100, 200, 600, 1000, 1600]\n",
    "\n",
    "for num_operators in measurement_range:   \n",
    "    \n",
    "    derandomized_hamiltonian = modified_derandomized_classical_shadow(observables_xyz, num_operators, \n",
    "                                                                      system_size, weight=absolute_coefficients)\n",
    "\n",
    "    tuples = (tuple(pauli) for pauli in derandomized_hamiltonian)\n",
    "    counts = Counter(tuples)\n",
    "\n",
    "    optimizer = SPSA(maxiter=2000) \n",
    "    cost_function = partial(objective_function, derandomized_hamiltonian)\n",
    "    \n",
    "    expectation_values = []\n",
    "    num_experiments = 5\n",
    "\n",
    "\n",
    "    for iteration in range(num_experiments):\n",
    "        \n",
    "        cost_history = []\n",
    "        \n",
    "        params = np.random.rand(ansatz.num_parameters)\n",
    "        result = optimizer.minimize(fun=cost_function, x0=params)\n",
    "        minimal_cost = min_cost()\n",
    "        expectation_values.append(minimal_cost)\n",
    "                \n",
    "        print(\"EXPERIMENT {}: GROUND STATE ENERGY FOUND = {}\".format(iteration, minimal_cost))\n",
    "\n",
    "\n",
    "    rmse_derandomised_cs = np.sqrt(np.sum([(expectation_values[i] - EXPECTED_EIGENVALUE)**2 \n",
    "                                           for i in range(num_experiments)])/num_experiments)\n",
    "\n",
    "    rmse_errors.append(rmse_derandomised_cs)\n",
    "    \n",
    "    print('{} | {} | {}'.format(num_operators, counts, rmse_derandomised_cs))\n",
    "\n",
    "\n",
    "elapsed_time = time.time() - start_time\n",
    "print(\"Execution time = \", time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\n",
    "\n",
    "          \n",
    "points = measurement_range\n",
    "num_points = len(measurement_range)\n",
    "\n",
    "plt.plot([i for i in points], [rmse_errors[i] for i in range(num_points)], 'r', label='SPSA(maxiter=2000)')\n",
    "plt.xlabel('Number of measurements')\n",
    "plt.ylabel('Average RMSE error')\n",
    "plt.legend(loc=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8841c2fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Electronic ground state energy:\n",
      " === GROUND STATE ENERGY ===\n",
      " \n",
      "* Electronic ground state energy (Hartree): -1.857275030202\n",
      "  - computed part:      -1.857275030202\n",
      "~ Nuclear repulsion energy (Hartree): 0.719968994449\n",
      "> Total ground state energy (Hartree): -1.137306035753\n",
      " \n",
      "=== MEASURED OBSERVABLES ===\n",
      " \n",
      "  0:  # Particles: 2.000 S: 0.000 S^2: 0.000 M: 0.000\n",
      " \n",
      "=== DIPOLE MOMENTS ===\n",
      " \n",
      "~ Nuclear dipole moment (a.u.): [0.0  0.0  1.3889487]\n",
      " \n",
      "  0: \n",
      "  * Electronic dipole moment (a.u.): [0.0  0.0  1.388948701555]\n",
      "    - computed part:      [0.0  0.0  1.388948701555]\n",
      "  > Dipole moment (a.u.): [0.0  0.0  -0.000000001555]  Total: 0.000000001555\n",
      "                 (debye): [0.0  0.0  -0.000000003953]  Total: 0.000000003953\n",
      " \n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "Above we have assumed a particular ground state energy for the molecule of interest.\n",
    "Below we corroborate this assumption using a classical minimum eigensolver on our Hamiltonian.\n",
    "'''\n",
    "\n",
    "converter = BravyiKitaevMapper()\n",
    "numpy_solver = NumPyMinimumEigensolver()   \n",
    "\n",
    "calc = GroundStateEigensolver(converter, numpy_solver)\n",
    "res = calc.solve(problem)\n",
    "print('Electronic ground state energy:\\n', res) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00222479",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}