{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "37a622fc-6c7a-4b52-909e-b6bcc5e0b2e5",
   "metadata": {},
   "source": [
    "# Quantum-Selected Configuration Interaction (QSCI)\n",
    "\n",
    "Quantum-Selected Configuration Interaction (QSCI) [1] is a quantum-classical hybrid algorithm designed for electronic structure calculations, particularly to compoute ground state and excited state energies of molecular systems. This method is one of the quantum subspace diagonalization (QSD). It is well-suited for NISQ-era quantum devices, focusing on using quantum resources efficiently while maintaining high accuracy.\n",
    "\n",
    "![quantum-subspace-diagonalization](images/qsd.png)\n",
    "\n",
    "QSCI selects this vector space from a sampled computational basis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1388e9b0-1b9b-4df4-890c-117622b771c4",
   "metadata": {},
   "source": [
    "## 0. Problem definition\n",
    "\n",
    "Let's define Hamiltonian of $\\mathrm{H_4}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "175153b6-3bb4-4378-a174-cead14589e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install --quiet --break-system-packages openfermionpyscf==0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "93c9a86f-6eab-48fb-9edf-51c7d25e9f1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openfermion\n",
    "import openfermionpyscf\n",
    "from openfermion.transforms import jordan_wigner, get_fermion_operator\n",
    "import cudaq\n",
    "import numpy as np\n",
    "\n",
    "if not hasattr(np, 'string_'):\n",
    "    np.string_ = np.bytes_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "336fd99f-7e4e-4401-a60e-eecb9c869299",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of hydrogen atoms.\n",
    "hydrogen_count = 4\n",
    "\n",
    "# Distance between the atoms in Angstroms.\n",
    "bond_distance = 0.7474\n",
    "\n",
    "# Define a linear chain of Hydrogen atoms\n",
    "geometry = [(\"H\", (0, 0, i * bond_distance)) for i in range(hydrogen_count)]\n",
    "\n",
    "basis = \"sto3g\"\n",
    "multiplicity = 1\n",
    "charge = 0\n",
    "\n",
    "molecule = openfermionpyscf.run_pyscf(\n",
    "    openfermion.MolecularData(geometry, basis, multiplicity, charge), run_fci=True\n",
    ")\n",
    "molecular_hamiltonian = molecule.get_molecular_hamiltonian()\n",
    "fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian)\n",
    "qubit_hamiltonian = jordan_wigner(fermion_hamiltonian)\n",
    "qubit_hamiltonian.compress()\n",
    "\n",
    "hamiltonian = cudaq.SpinOperator(qubit_hamiltonian)\n",
    "\n",
    "num_qubits = hamiltonian.qubit_count\n",
    "electron_count = molecule.n_electrons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "908956ec-a201-416d-bc37-f09c6ba02b85",
   "metadata": {},
   "source": [
    "## 1. Prepare an Approximate Quantum State\n",
    "\n",
    "Start with a variational ansatz (e.g., from VQE) to approximate the ground state wavefunction on a quantum computer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "365q8kclsh",
   "metadata": {},
   "outputs": [],
   "source": [
    "if cudaq.num_available_gpus() > 0 and cudaq.has_target(\"nvidia\"):\n",
    "    cudaq.set_target(\"nvidia\")\n",
    "else:\n",
    "    print(\"CUDA or GPU support is unavailable. Running with CPU simulator. Performance may be significantly reduced.\")\n",
    "    cudaq.set_target(\"qpp-cpu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "04b6e475-72bc-4f4a-bed7-32488d044e0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def kernel(thetas: list[float]):\n",
    "\n",
    "    qubits = cudaq.qvector(num_qubits)\n",
    "\n",
    "    # Hartree-Fock\n",
    "    for i in range(electron_count):\n",
    "        x(qubits[i])\n",
    "    # UCCSD\n",
    "    cudaq.kernels.uccsd(qubits, thetas, electron_count, num_qubits)\n",
    "\n",
    "\n",
    "parameter_count = cudaq.kernels.uccsd_num_parameters(electron_count, num_qubits)\n",
    "\n",
    "exp_vals = []\n",
    "\n",
    "\n",
    "def cost(thetas: list[float]):\n",
    "    exp_val = cudaq.observe(kernel, hamiltonian, thetas).expectation()\n",
    "    exp_vals.append(exp_val)\n",
    "    return exp_val\n",
    "\n",
    "\n",
    "# Initial variational parameters\n",
    "np.random.seed(15)\n",
    "x0 = np.clip(np.random.normal(0, np.pi, parameter_count), -np.pi, np.pi)\n",
    "\n",
    "optimizer = cudaq.optimizers.COBYLA()\n",
    "optimizer.max_iterations = 300\n",
    "optimizer.initial_parameters = x0.tolist()\n",
    "\n",
    "energy, opt_params = optimizer.optimize(parameter_count, cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02dc831a-c41f-4abf-8347-4df7a0d5ca72",
   "metadata": {},
   "source": [
    "## 2 Quantum Sampling to Select Configuration\n",
    "\n",
    "Measure this quantum state multiple times to obtain bitstrings, which correspond to Slater determinants (electronic configurations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0eb10a37-9384-41f5-bb9e-92d91700434a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ 01011010:1 11100100:2 11110000:997 }\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sample_result = cudaq.sample(kernel, opt_params)\n",
    "print(sample_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "836fb668-4bb7-4214-aaaf-72d8f697a384",
   "metadata": {},
   "source": [
    "## 3. Classical Diagonalization on the Selected Subspace\n",
    "\n",
    "Construct a truncated Hamiltonian matrix using only the selected determinants.\n",
    "Diagonalize this matrix classically to obtain improved energy estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "eb220374-44ca-4183-9e1a-f27cc62e8fd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_element(\n",
    "    operator: cudaq.SpinOperator, b1: str | np.ndarray, b2: str | np.ndarray\n",
    ") -> complex:\n",
    "    \"\"\"\n",
    "    Compute an inner product <b1|operator|b2>.\n",
    "    \"\"\"\n",
    "\n",
    "    if isinstance(b1, str):\n",
    "        b1 = np.array([i == \"1\" for i in b1], dtype=bool)\n",
    "    if isinstance(b2, str):\n",
    "        b2 = np.array([i == \"1\" for i in b2], dtype=bool)\n",
    "\n",
    "    num_qubits = hamiltonian.qubit_count\n",
    "    num_terms = hamiltonian.term_count\n",
    "\n",
    "    bsv = np.zeros((num_terms, 2 * num_qubits), dtype=bool)\n",
    "    coeffs = np.zeros(num_terms, dtype=np.complex128)\n",
    "    for i, term in enumerate(operator):\n",
    "        coeffs[i] = term.evaluate_coefficient()\n",
    "        bsf = term.get_binary_symplectic_form()\n",
    "        n = len(bsf) // 2\n",
    "        for j, e in enumerate(bsf):\n",
    "            if j < n:\n",
    "                bsv[i, j] = e\n",
    "            else:\n",
    "                bsv[i, j - n + num_qubits] = e\n",
    "    x = bsv[:, :num_qubits]\n",
    "    z = bsv[:, num_qubits:]\n",
    "\n",
    "    # b1 ⊕ b2 must match x for nonzero contribution\n",
    "    delta = np.bitwise_xor(b1, b2)\n",
    "    match = np.all(x == delta, axis=1)\n",
    "\n",
    "    if not np.any(match):\n",
    "        return 0.0\n",
    "\n",
    "    matched_x = x[match]\n",
    "    matched_z = z[match]\n",
    "    matched_coeffs = coeffs[match]\n",
    "\n",
    "    # i^{x·z}\n",
    "    x_dot_z = np.sum(np.bitwise_and(matched_x, matched_z), axis=1) % 4\n",
    "    i_phases = np.array([1, 1j, -1, -1j], dtype=np.complex128)[x_dot_z]\n",
    "\n",
    "    # (-1)^{b1·z}\n",
    "    b1_dot_z_parity = np.sum(np.bitwise_and(b1, matched_z), axis=1) % 2\n",
    "    sign_phases = np.where(b1_dot_z_parity == 0, 1, -1)\n",
    "\n",
    "    result_terms = matched_coeffs * i_phases * sign_phases\n",
    "    return np.sum(result_terms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5c095c67-1222-4366-bf0b-4a4d847a2841",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "\n",
    "bitstrs = list(k for k, _ in sample_result.items())\n",
    "subspace_dimension = len(bitstrs)\n",
    "\n",
    "values = []\n",
    "row_ids = []\n",
    "column_ids = []\n",
    "\n",
    "for i in range(subspace_dimension):\n",
    "    for j in range(i, subspace_dimension):\n",
    "        elem = pauli_element(hamiltonian, bitstrs[i], bitstrs[j])\n",
    "        if abs(elem) > 1e-6:\n",
    "            values.append(elem)\n",
    "            row_ids.append(i)\n",
    "            column_ids.append(j)\n",
    "            if i != j:\n",
    "                values.append(elem.conjugate())\n",
    "                row_ids.append(j)\n",
    "                column_ids.append(i)\n",
    "\n",
    "values = np.real_if_close(values)\n",
    "\n",
    "truncated_hamiltonian = scipy.sparse.coo_array(\n",
    "    (values, (row_ids, column_ids)), shape=(subspace_dimension, subspace_dimension)\n",
    ")\n",
    "\n",
    "# eigsh requires k < N; fall back to dense solver for small subspaces\n",
    "if subspace_dimension > 1:\n",
    "    eigvals, _ = scipy.sparse.linalg.eigsh(truncated_hamiltonian, 1, which=\"SA\")\n",
    "else:\n",
    "    eigvals = scipy.linalg.eigh(truncated_hamiltonian.toarray(), eigvals_only=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb40b6c3-063a-45df-8867-747fac2bb354",
   "metadata": {},
   "source": [
    "## 5. Compare results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "827cfd92-9208-4057-959f-db70853a32d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.1033232756074405"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## QSCI\n",
    "eigvals[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "df643aee-1967-4d3d-aebd-dda13892c9ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.14355458044603"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## FCI Energy\n",
    "molecule.fci_energy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd08b5ed-dc5a-4f25-86fd-a6c805ee8ee6",
   "metadata": {},
   "source": [
    "## Reference\n",
    "\n",
    "[1] \"Quantum-selected configuration interaction: classical diagonalization of hamiltonians in subspaces selected by quantum computers\" K. Kanno, M. Kohda, R. Imai, S. Koh, K. Mitarai, W. Mizukami, and Y. O. Nakagawa (2023) arXiv: 2302.11320"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
