{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "db96a82a",
   "metadata": {},
   "source": [
    "Variational Quantum Deflation for excited states \n",
    "============================================================"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "de3c3f91",
   "metadata": {},
   "source": [
    "In this file, we will introduce the methodology for finding electronic excited state energies using the Variational Quantum Deflation (VQD) approach in InQuanto. For the original publication by Higgot, Wang and Brierley (2019) and more technical details, please go [here](https://arxiv.org/abs/1805.08138).\n",
    "\n",
    "In VQE, one finds a minimum of the energy by classically optimizing the function below with respect to the wavefunction parameters, $\\{\\lambda\\}$:\n",
    "\n",
    "\n",
    "$E(\\lambda) = \\langle \\psi(\\lambda) | H | \\psi(\\lambda)\\rangle = \\sum_j c_j \\langle \\psi (\\lambda) | P_j | \\psi (\\lambda) \\rangle$, \n",
    "\n",
    "\n",
    "where $P_j$ are terms in the qubit Hamiltonian, and $c_j$ are classically pre-computed coefficients. In contrast, when executing a VQD simulation, the objective function is modified to include a penalty term,  multiplied by a weight $\\beta$,\n",
    "\n",
    "\n",
    "$F(\\lambda_k) = \\langle \\psi(\\lambda_k) | H | \\psi(\\lambda_k)\\rangle + \\sum_{i=0}^{k-1} \\beta_i |\\langle \\psi (\\lambda_k) | P_j | \\psi (\\lambda_i) \\rangle|^2$. \n",
    "\n",
    "\n",
    "This enforces the requirement that each $|\\psi(\\lambda_k)\\rangle$ be orthogonal to the other $|\\psi(\\lambda_0)\\rangle, \\ldots |\\psi(\\lambda_{k-1})\\rangle$ found by previous optimizations of the objective function, $F(\\lambda_k)$.\n",
    "\n",
    "In order to run a VQD calculation, one needs several things:\n",
    "\n",
    "1. A molecular, electronic Hamiltonian operator, $H$,\n",
    "\n",
    "2. A mapping object for constructing the qubit Hamiltonian,\n",
    "\n",
    "3. An ansatz for computing the ground state energy, $|\\psi(\\lambda)\\rangle$,\n",
    "\n",
    "4. A classical minimizer\n",
    "\n",
    "5. A complete VQE experiment,\n",
    "\n",
    "6. A second ansatz for describing the excited states, $\\{|\\psi(\\lambda_k)\\rangle\\}$,\n",
    "\n",
    "7. A expression for evaluating the weights, $\\{\\beta_i\\}$,\n",
    "\n",
    "8. The number of excited states of interest, $k$.\n",
    "\n",
    "We are concerned with first finding the ground state energy of the second-quantized molecular electronic Hamiltonian using the Variational Quantum Eigensolver (please see [tutorials 1, and 2](https://docs.quantinuum.com/inquanto/tutorials/tutorial_overview.html) for a more in-depth explanation). So let's proceed with points 1-4 in the list above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07648800",
   "metadata": {},
   "source": [
    "First, we need to import a backend, and the appropriate space and state objects. Since we are looking at fermions, these are `inquanto.spaces.FermionSpace` and `inquanto.states.FermionState`. We will use the AerBackend available through the <code>pytket</code> qiskit extension to simulate the quantum hardware."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dd6e897",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pytket.extensions.qiskit import AerBackend\n",
    "from inquanto.spaces import FermionSpace\n",
    "from inquanto.states import FermionState"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d1e1537e",
   "metadata": {},
   "source": [
    "We're going to simulate the dihydrogen molecule in the STO-3G basis. There are 4 spin orbitals -- two of which are occupied --  and the reference state lives in 4-dimensional Fock space. The corresponding objects can be constructed as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6faace1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Coefficient             Term\n",
      "0      0.743018                 \n",
      "1     -1.270293          F0^ F0 \n",
      "2      0.680062  F1^ F0^ F0  F1 \n",
      "3      0.179669  F1^ F0^ F2  F3 \n",
      "4     -1.270293          F1^ F1 \n",
      "5      0.668577  F2^ F0^ F0  F2 \n",
      "6      0.179669  F2^ F0^ F2  F0 \n",
      "7      0.668577  F2^ F1^ F1  F2 \n",
      "8      0.179669  F2^ F1^ F3  F0 \n",
      "9     -0.456807          F2^ F2 \n",
      "10     0.668577  F3^ F0^ F0  F3 \n",
      "11     0.179669  F3^ F0^ F2  F1 \n",
      "12     0.668577  F3^ F1^ F1  F3 \n",
      "13     0.179669  F3^ F1^ F3  F1 \n",
      "14     0.179669  F3^ F2^ F0  F1 \n",
      "15     0.702814  F3^ F2^ F2  F3 \n",
      "16    -0.456807          F3^ F3 \n"
     ]
    }
   ],
   "source": [
    "from inquanto.express import load_h5\n",
    "from inquanto.mappings import QubitMappingJordanWigner\n",
    "h2 = load_h5(\"h2_sto3g.h5\")\n",
    "fermion_hamiltonian = h2[\"hamiltonian_operator\"]\n",
    "qubit_hamiltonian = QubitMappingJordanWigner().operator_map(fermion_hamiltonian)\n",
    "\n",
    "space = FermionSpace(4)\n",
    "state = FermionState([1, 1, 0, 0])\n",
    "\n",
    "print(fermion_hamiltonian.df())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e25461f7",
   "metadata": {},
   "source": [
    "where we have loaded in the molecular Hamiltonian using InQuanto's `express` module, and mapped it to a qubit Hamiltonian using the Jordan-Wigner transformation.\n",
    "\n",
    "Now, for points 3 and 4 in the list, we need to construct an ansatz for our ground state calculation and choose a classical minimizer. For this example we will use the k-UpCCGSD `inquanto.ansatzes.FermionSpaceAnsatzkUpCCGSD` ansatz and the COBYLA minimizer available through the `inquanto.minimizers.MinimizerScipy` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "722520df",
   "metadata": {},
   "outputs": [],
   "source": [
    "from inquanto.ansatzes import FermionSpaceAnsatzkUpCCGSD\n",
    "from inquanto.minimizers import MinimizerScipy\n",
    "\n",
    "ansatz = FermionSpaceAnsatzkUpCCGSD(space, state, k_input=2)\n",
    "minimizer = MinimizerScipy(method=\"COBYLA\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ac258e4b",
   "metadata": {},
   "source": [
    "We're now in a position to address item 5 in the list - running a complete VQE calculation. We know from the first equation in this notebook that the objective function is the expectation value of the qubit Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43fb2048",
   "metadata": {},
   "outputs": [],
   "source": [
    "from inquanto.computables import ExpectationValue\n",
    "from inquanto.algorithms import AlgorithmVQE\n",
    "\n",
    "expectation_value = ExpectationValue(ansatz, qubit_hamiltonian)\n",
    "vqe = AlgorithmVQE(\n",
    "    objective_expression=expectation_value,\n",
    "    minimizer=minimizer,\n",
    "    initial_parameters=ansatz.state_symbols.construct_random(seed=0)\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d434d30e",
   "metadata": {},
   "source": [
    "We have passed in some random $\\{\\lambda\\}$ using the state_symbols attribute of the ansatz object as our starting guess.\n",
    "\n",
    "We now choose our measurement protocol -- in this case, a direct measurement by operator averaging, so we choose `inquanto.protocols.PauliAveraging` -- and the number of shots we wish to simulate in each iteration. Then, we build the algorithm object and execute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e937f844",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# TIMER BLOCK-0 BEGINS AT 2025-06-12 15:11:40.620960\n",
      "# TIMER BLOCK-0 ENDS - DURATION (s): 27.6604800 [0:00:27.660480]\n",
      "VQE Energy:     -1.134404492135646\n",
      "VQE Parameters: {gd0k0: np.float64(1.5163690225114363), gd0k1: np.float64(-1.552957880628078), gs0k0: np.float64(0.35595432868977916), gs0k1: np.float64(1.7188670637293226), gs1k0: np.float64(-1.3603698543071556), gs1k1: np.float64(-0.3503255748500169)}\n"
     ]
    }
   ],
   "source": [
    "from inquanto.protocols import PauliAveraging\n",
    "\n",
    "vqe.build(protocol_objective=PauliAveraging(AerBackend(), shots_per_circuit=10000))\n",
    "\n",
    "vqe.run()\n",
    "\n",
    "# VQE Energy:     -1.1354204303965678\n",
    "# VQE Parameters: [ 1.301 -1.538  0.287  1.699 -1.339 -0.344]\n",
    "print(\"VQE Energy:    \", vqe.final_value)\n",
    "print(\"VQE Parameters:\", vqe.final_parameters)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "08fef881",
   "metadata": {},
   "source": [
    "According to point 6, we now need a second, deflationary ansatz we can use to describe the excited states. To do this, we'll use the same ansatz structure and just make a copy of the first ansatz. We update our symbols from those used in the ground state to some other symbols. Consider this as constructing the symbols in $\\{\\lambda_k\\}$ using those in $\\{\\lambda\\}$ as a reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e75f3b21",
   "metadata": {},
   "outputs": [],
   "source": [
    "ansatz_2 = ansatz.subs(\"{}_2\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1692a8b9",
   "metadata": {},
   "source": [
    "It is almost time to construct, build and execute our VQD algorithm. First, we need to write expressions corresponding to the terms in the functional  \n",
    "\n",
    "$F(\\lambda_k) = \\langle \\psi(\\lambda_k) | H | \\psi(\\lambda_k)\\rangle + \\sum_{i=0}^{k-1} \\beta_i |\\langle \\psi (\\lambda_k) | P_j | \\psi (\\lambda_i) \\rangle|^2$.\n",
    "\n",
    "\n",
    "We will refer to the leading term as `expectation_value`, the weights as `weight_expression`, and the overlap term in the penalty as `overlap_expression`.\n",
    "\n",
    "We will also select the weights as the expectation value of the deflationary ansatz with respect to the sign-flipped qubit Hamiltonian to ensure it is sufficiently large to act as a constraint rather than a weak penalty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4241a9fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from inquanto.computables import OverlapSquared\n",
    "\n",
    "expectation_value = ExpectationValue(ansatz_2, qubit_hamiltonian)\n",
    "\n",
    "weight_expression = ExpectationValue(ansatz_2, -1 * qubit_hamiltonian)\n",
    "\n",
    "overlap_expression = OverlapSquared(ansatz, ansatz_2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7afc14f",
   "metadata": {},
   "source": [
    "As was the case with the previous VQE experiment, we must choose protocols for measuring the overlaps. For this we choose to use the vacuum test, available through the `inquanto.protocols.ComputeUncompute` object.\n",
    "\n",
    "We must also choose the number of excited states we wish to find. For this calculation, we'll choose to find 3 and pass this into the `inquanto.algorithms.AlgorithmVQD` constructor in the `n_vectors` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b29ef55",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# TIMER BLOCK-1 BEGINS AT 2025-06-12 15:12:08.301251\n",
      "# TIMER BLOCK-1 ENDS - DURATION (s): 21.9974220 [0:00:21.997422]\n",
      "# TIMER BLOCK-2 BEGINS AT 2025-06-12 15:12:30.298584\n",
      "# TIMER BLOCK-2 ENDS - DURATION (s): 27.8214365 [0:00:27.821436]\n",
      "# TIMER BLOCK-3 BEGINS AT 2025-06-12 15:12:58.119904\n",
      "# TIMER BLOCK-3 ENDS - DURATION (s): 44.4877955 [0:00:44.487795]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<inquanto.algorithms.vqd._algorithm_vqd.AlgorithmVQD at 0x16a4a6ae0>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from inquanto.algorithms import AlgorithmVQD\n",
    "from inquanto.protocols import ComputeUncompute\n",
    "# instantiate VQD object\n",
    "vqd = AlgorithmVQD(\n",
    "    objective_expression=expectation_value,\n",
    "    overlap_expression=overlap_expression,\n",
    "    weight_expression=weight_expression,\n",
    "    minimizer=MinimizerScipy(method=\"COBYLA\"),\n",
    "    initial_parameters=ansatz_2.state_symbols.construct_random(seed=0),\n",
    "    vqe_value=vqe.final_value,\n",
    "    vqe_parameters=vqe.final_parameters,\n",
    "    n_vectors=3,\n",
    ")\n",
    "# build object\n",
    "backend=AerBackend()\n",
    "vqd.build(\n",
    "    #small number of shots for demonstration purposes leads to large stochastic error\n",
    "    objective_protocol=PauliAveraging(backend, shots_per_circuit=1000),\n",
    "    weight_protocol=PauliAveraging(backend, shots_per_circuit=1000),\n",
    "    overlap_protocol=ComputeUncompute(backend, n_shots=1000),\n",
    ")\n",
    "\n",
    "# execute\n",
    "vqd.run()\n",
    "\n",
    "# print results\n",
    "# VQD excited state energies:    [-1.1354204303965678, -0.4949467734066755, -0.11981345040527547, 0.5478565335949009]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b71d1b55",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VQD excited state energies:    [np.float64(-1.134404492135646), np.float64(-0.4803282017464972), np.float64(-0.14069004669271543), np.float64(0.4880875698728495)]\n"
     ]
    }
   ],
   "source": [
    "print(\"VQD excited state energies:   \", vqd.final_values)"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}