{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2c935a4a-aa7d-499a-b358-3ef7fac417c2",
   "metadata": {},
   "source": [
    "# Obtaining the deuteron energy descibed on two qubits (Dumitrescu et al., PRL 120, 210501 (2018))\n",
    "\n",
    "Our starting point is the 2‑qubit Hamiltonian\n",
    "\n",
    "$H_2 = 5.906709\\, I + 0.218291\\, Z_0 - 6.125\\, Z_1\n",
    "      - 2.143304\\left(X_0X_1 + Y_0Y_1\\right).$\n",
    "      \n",
    "The trial state is \n",
    "$\n",
    "|\\psi(\\theta)\\rangle =\n",
    "e^{\\,i\\frac{\\theta}{2}\\left(X_0 Y_1 - X_1 Y_0\\right)} |00\\rangle,$\n",
    "\n",
    "1. Define the Hamiltonian and the two qubit circuit\n",
    "2. For a given theta, prepare the ansatz |\\psi(\\theta)\\rangle\n",
    "3. Compute the expectation values for $Z_0$, $Z_1$, $X_0X_1$, $Y_0Y_1$ with a finite number of measurements.  \n",
    "4. Deduce the expectation value of $E(\\theta)=\\langle H_2\\rangle$\n",
    "5. Then reproduce the Fig. 2 of (Dumitrescu et al., PRL 120, 210501 (2018)) where the different Pauli strings are shown as a function of \\(theta\\), show also \\(E(\\theta)\\) as a function of \\(theta\\), and deduce the ground state energy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4766d432-a17d-4e44-8f69-1734e11c0ff6",
   "metadata": {},
   "source": [
    "## 0. Imports\n",
    "\n",
    "- `QuantumCircuit` to build circuits\n",
    "- `StatevectorSampler` for **shot-based** counts (like running the circuit many times)\n",
    "- `Statevector` for exact (ideal) checks (no shot noise)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "09303518",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit\n",
    "from qiskit.quantum_info import Statevector\n",
    "from qiskit.primitives import StatevectorSampler\n",
    "from qiskit.visualization import plot_histogram\n",
    "import numpy as np\n",
    "\n",
    "# for vizualization \n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae79692c-9c5b-4cf1-bcd1-049d2d064d72",
   "metadata": {},
   "source": [
    "## 1. Circuit for the trial state vector\n",
    "\n",
    "The trial state is \n",
    "\\[\n",
    "|\\psi(\\theta)\\rangle =\n",
    "e^{\\,i\\frac{\\theta}{2}\\left(X_0 Y_1 - X_1 Y_0\\right)} |00\\rangle,\n",
    "\\]\n",
    "It could be prepared by \n",
    "1) performing a X gate on the first qubit\n",
    "2) a Y-rotation with angle \\( theta \\) on the first qubit\n",
    "3) a CNOT where 1 control 0\n",
    "\n",
    "**Task: For a fixed $\\theta$, make the circuit and perform measurements by counts** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "47e471f4-5345-405d-91cf-c474a4257ecb",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Here Make the circuit\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a86c0e-08a2-468f-8c8f-7782c1289d2d",
   "metadata": {},
   "source": [
    "## 3. Computing expectation values of observables for one single angle value \n",
    "\n",
    "To measure on one qubit:\n",
    "\n",
    "- **Z**: do nothing, measure in computational basis\n",
    "- **X**: apply **H**, then measure Z\n",
    "- **Y**: apply **Sdg** then **H**, then measure Z\n",
    "\n",
    "For two-qubit Pauli strings (like $X_0X_1$), apply the rule on each qubit.\n",
    "\n",
    "### From counts to expectation values\n",
    "\n",
    "Eigenvalues:\n",
    "- bit `0` → +1\n",
    "- bit `1` → −1\n",
    "\n",
    "Qiskit bitstrings are shown as `'b1b0'` (highest classical bit on the left). \n",
    "Tips: You need to convert counts to expectation values, you can use a similar \n",
    "routine used in the two-qubit tutorial. Examples of routines are also given below.  \n",
    "\n",
    "**Task: For a fixed $\\theta$, make the circuits and computes the different expectation values** \n",
    "\n",
    "In this \"hard\" version, for a fixed $\\theta$, the expectation of $\\langle Z_0 \\rangle$, $\\langle Z_1 \\rangle$, $\\langle X_0 X_1 \\rangle$ and $\\langle Y_0 Y_1 \\rangle$ should be computed.\n",
    "\n",
    "The routines from counts to expectations are given. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "63288326-c79b-4d9c-9e3a-b62473ee9ec2",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Here, run the code and complete it to get all expectation values. \n",
    "\n",
    "sampler = StatevectorSampler()\n",
    "Nshot = 1000\n",
    "theta = np.pi/3\n",
    "\n",
    "# -----------------------\n",
    "# 3) Convert counts -> expectation values for Z0, Z1, X0X1, Y0Y1\n",
    "# -----------------------\n",
    "def bit_to_ev(b: int) -> int:\n",
    "    # Z-eigenvalue mapping\n",
    "    return 1 if b == 0 else -1\n",
    "\n",
    "def exp_Z0(counts: dict, shots: int) -> float:\n",
    "    # Qiskit bitstrings are 'b1b0' (bit 1 left, bit 0 right)\n",
    "    return sum(bit_to_ev(int(s[1])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_Z1(counts: dict, shots: int) -> float:\n",
    "    return sum(bit_to_ev(int(s[0])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_two_qubit_product(counts: dict, shots: int) -> float:\n",
    "    # After basis rotations, Z-measurement corresponds to X or Y measurement.\n",
    "    # The product eigenvalue is (ev on qubit0) * (ev on qubit1).\n",
    "    total = 0.0\n",
    "    for s, c in counts.items():\n",
    "        b1 = int(s[0])\n",
    "        b0 = int(s[1])\n",
    "        total += (bit_to_ev(b0) * bit_to_ev(b1)) * c\n",
    "    return total / shots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f647fec-3bbe-4311-8e78-9d8aeca8e207",
   "metadata": {},
   "source": [
    "## 3. Computing expectation values of observables for a set of angles\n",
    "\n",
    "Tips: \n",
    "\n",
    "### Making loops on $\\theta$\n",
    "\n",
    "- You should define a list of angle values and make a loop on the angles, for instance:  \n",
    " \n",
    "``` \n",
    "n_theta = 41  # number of theta points\n",
    "theta_list = np.linspace(-np.pi, np.pi, n_theta)\n",
    "```\n",
    "\n",
    "### Example of useful routines for the $\\langle Z_0 \\rangle $ expectation value\n",
    "\n",
    "- Prior to the loop, you can define a routine that builds the circuit for a given $\\theta$ value. \n",
    "An example of such a circuit to measure $Z_0$ (or $Z_1$) is \n",
    "\n",
    "```\n",
    "def circuit_z(theta: float) -> QuantumCircuit:\n",
    "    qc = QuantumCircuit(2, 2, name=\"Measure_Z0Z1\")\n",
    "    qc.ry(theta, 1)\n",
    "    qc.x(0)\n",
    "    qc.cx(1, 0)\n",
    "    qc.measure([0, 1], [0, 1])\n",
    "    return qc\n",
    "```\n",
    "\n",
    "The routine can be called in a loop on $\\theta$ value for instance using the syntax\n",
    "  \n",
    "```\n",
    "counts_z  = sampler.run([circuit_z(theta)],  shots=Nshot).result()[0].data.c.get_counts()\n",
    "```\n",
    "\n",
    "### Storing\n",
    "\n",
    "To store results as a function of $\\theta$ this could be done by defining a list: \n",
    "\n",
    "```\n",
    "Z0_list = []\n",
    "```\n",
    "\n",
    "Results could be stored in the $\\theta$ loop using \n",
    "\n",
    "```\n",
    "    Z0_list.append(exp_Z0(counts_z, Nshot))\n",
    "```\n",
    "\n",
    "with `exp_Z0` the routine to transform bits into qubits and compute expectation values (given below).  \n",
    "After the loop, you do \n",
    "```\n",
    "Z0_list = np.array(Z0_list)\n",
    "```\n",
    "that are the $\\langle Z_0 \\rangle$ results.\n",
    "\n",
    "- As an intermediate step, just do the code that shows the expectation value of $\\langle Z_0 \\rangle$ and display this expectation value as a function of $\\theta$ and display it\n",
    "\n",
    "To display the result, use matplotib standard routines. For instance, a minimal script is:\n",
    "\n",
    "```\n",
    "fig, ax = plt.subplots(figsize=(6.0, 4.0))\n",
    "ax.plot(theta_list, Z_0_list, marker = 'o', label=r\"$\\langle Z_0\\rangle$\")\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "- Then, do the same for $\\langle Z_1 \\rangle$, $\\langle X_0X_1\\rangle $ and $\\langle Y_0Y_1\\rangle $.\n",
    "- Note that, while for $\\langle Z_1 \\rangle$, the same circuit can be used,for $\\langle X_0X_1\\rangle $ and $\\langle Y_0Y_1 \\rangle$ additional circuit should be built with proper rotations before measurements.   \n",
    "\n",
    "## Task: \n",
    "In this \"hard\" version, for a fixed $\\theta$, using the tips above,  the expectation of $\\langle Z_0 \\rangle$, $\\langle Z_1 \\rangle$, $\\langle X_0 X_1 \\rangle$ and $\\langle Y_0 Y_1 \\rangle$ should be computed as a function of $\\theta$, results should be stored in lists and displayed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "cbbd0d97-b74c-4fb3-9e3b-3bc89749a707",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Here, run the code and complete it to get all expectation values.   \n",
    "\n",
    "Nshot = 1000\n",
    "n_theta = 41  # number of theta points\n",
    "theta_list = np.linspace(-np.pi, np.pi, n_theta)\n",
    "\n",
    "sampler = StatevectorSampler()\n",
    "\n",
    "# -----------------------\n",
    "# Helpers: counts -> expectation values\n",
    "# -----------------------\n",
    "def bit_to_ev(b: int) -> int:\n",
    "    # Z-eigenvalue mapping\n",
    "    return 1 if b == 0 else -1\n",
    "\n",
    "def exp_Z0(counts: dict, shots: int) -> float:\n",
    "    # Qiskit bitstrings are 'b1b0' (bit 1 left, bit 0 right)\n",
    "    return sum(bit_to_ev(int(s[1])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_Z1(counts: dict, shots: int) -> float:\n",
    "    return sum(bit_to_ev(int(s[0])) * c for s, c in counts.items()) / shots\n",
    "\n",
    "def exp_two_qubit_product(counts: dict, shots: int) -> float:\n",
    "    # After basis rotations, Z-measurement corresponds to X or Y measurement.\n",
    "    total = 0.0\n",
    "    for s, c in counts.items():\n",
    "        b1 = int(s[0])\n",
    "        b0 = int(s[1])\n",
    "        total += (bit_to_ev(b0) * bit_to_ev(b1)) * c\n",
    "    return total / shots\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "970f568b-3ee4-417a-b34d-115b78f04765",
   "metadata": {},
   "source": [
    "## 4. Expectation value of $H_2$ from measurements as a function of $\\theta$\n",
    "\n",
    "## Task \n",
    "From the different expectation values above, compute the expectation value of $\\langle H_2 \\rangle$ \n",
    "as a list and display the result as a function of $\\theta$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e1fe9799-f4b7-4102-a9ea-b71d8f609454",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Here show the result \n"
   ]
  }
 ],
 "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.11.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}