{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "65adf860-4443-494f-ad90-5d8ea444ac4e",
   "metadata": {},
   "source": [
    "# Scale SQD chemistry workflows with Fulqrum\n",
    "\n",
    "[Fulqrum](https://github.com/Qiskit/fulqrum) is a software package that computes eigenvalues and eigenvectors of Hamiltonians projected into a subspace. It constructs a portable representation of the projected Hamiltonian that you can pass to an eigenvalue solver of your choice, including standard off-the-shelf solvers. These solvers can provide a much faster alternative to the default PySCF solver that the SQD addon uses for diagonalization. For information on how to install and use Fulqrum, see the [Fulqrum documentation](https://qiskit.github.io/fulqrum/).\n",
    "\n",
    "Fulqrum's `addon_sqd_support` module builds a solver that [`diagonalize_fermionic_hamiltonian`](../apidocs/qiskit_addon_sqd.fermion.rst) accepts directly, so you pass the same integrals to both and Fulqrum handles the projection and diagonalization.\n",
    "\n",
    "For more details on the SQD code used in this example, see the [chemistry Hamilton tutorial](https://quantum.cloud.ibm.com/docs/tutorials/sample-based-quantum-diagonalization)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "baf7d074-4443-4695-875e-65f7fb8cef25",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "WARN: Unable to to identify input symmetry using original axes.\n",
      "Different symmetry axes will be used.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -108.835236570774\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pyscf\n",
    "import pyscf.cc\n",
    "import pyscf.mcscf\n",
    "from fulqrum.addon_sqd_support import sci_solver_from_integrals\n",
    "from qiskit_addon_sqd.counts import generate_bit_array_uniform\n",
    "from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian\n",
    "\n",
    "# Build N2 molecule\n",
    "mol = pyscf.gto.Mole()\n",
    "mol.build(\n",
    "    atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n",
    "    basis=\"6-31g\",\n",
    "    symmetry=\"Dooh\",\n",
    ")\n",
    "\n",
    "# Define active space\n",
    "n_frozen = 2\n",
    "active_space = range(n_frozen, mol.nao_nr())\n",
    "\n",
    "# Get molecular integrals\n",
    "scf = pyscf.scf.RHF(mol).run()\n",
    "num_orbitals = len(active_space)\n",
    "n_electrons = int(sum(scf.mo_occ[active_space]))\n",
    "num_elec_a = (n_electrons + mol.spin) // 2\n",
    "num_elec_b = (n_electrons - mol.spin) // 2\n",
    "cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))\n",
    "mo = cas.sort_mo(active_space, base=0)\n",
    "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
    "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)\n",
    "\n",
    "# Compute exact energy\n",
    "exact_energy = cas.run().e_tot\n",
    "\n",
    "# Create a seed to control randomness throughout this workflow\n",
    "rng = np.random.default_rng(24)\n",
    "\n",
    "# Generate random samples\n",
    "bit_array = generate_bit_array_uniform(10_000, num_orbitals * 2, rand_seed=rng)\n",
    "\n",
    "# Build a solver that diagonalizes with Fulqrum. It converts the integrals to a\n",
    "# Fulqrum operator once, then projects onto each subspace the SQD loop requests.\n",
    "sci_solver = sci_solver_from_integrals(hcore, eri, tol=1e-12)\n",
    "\n",
    "# Run SQD\n",
    "result = diagonalize_fermionic_hamiltonian(\n",
    "    hcore,\n",
    "    eri,\n",
    "    bit_array,\n",
    "    samples_per_batch=300,\n",
    "    norb=num_orbitals,\n",
    "    nelec=(num_elec_a, num_elec_b),\n",
    "    num_batches=5,\n",
    "    max_iterations=5,\n",
    "    sci_solver=sci_solver,\n",
    "    symmetrize_spin=True,\n",
    "    seed=rng,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "99709a88-f9ec-4dbc-a561-e682f90aa4c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact energy: -109.04667177808024\n",
      "Estimated energy: -109.03402667559365\n"
     ]
    }
   ],
   "source": [
    "print(f\"Exact energy: {exact_energy}\")\n",
    "print(f\"Estimated energy: {result.energy + nuclear_repulsion_energy}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab599762",
   "metadata": {},
   "source": [
    "## Use a different eigensolver\n",
    "\n",
    "Because Fulqrum exposes the projected Hamiltonian as a standard [`scipy.sparse.linalg.LinearOperator`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html), you can diagonalize it with any eigensolver of your choosing rather than being tied to `scipy.sparse.linalg.eigsh`. This makes it straightforward to substitute a solver that is better suited to your problem size or hardware.\n",
    "\n",
    "Pass the eigensolver you want as the `eigensolver` argument. As an example, this guide demonstrates [PRIMME](https://github.com/primme/primme), a high-performance eigensolver library. Its Python interface, `primme.eigsh`, follows a calling convention similar to `scipy.sparse.linalg.eigsh`, so you can use it as a direct replacement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "48850fb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import primme\n",
    "\n",
    "# Build a solver that diagonalizes with PRIMME instead of the default,\n",
    "# scipy.sparse.linalg.eigsh\n",
    "sci_solver = sci_solver_from_integrals(hcore, eri, eigensolver=primme.eigsh, tol=1e-12)\n",
    "\n",
    "# Run SQD using the PRIMME-based solver\n",
    "result = diagonalize_fermionic_hamiltonian(\n",
    "    hcore,\n",
    "    eri,\n",
    "    bit_array,\n",
    "    samples_per_batch=300,\n",
    "    norb=num_orbitals,\n",
    "    nelec=(num_elec_a, num_elec_b),\n",
    "    num_batches=5,\n",
    "    max_iterations=5,\n",
    "    sci_solver=sci_solver,\n",
    "    symmetrize_spin=True,\n",
    "    seed=rng,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "39b960fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact energy: -109.04667177808024\n",
      "Estimated energy: -109.03937429196168\n"
     ]
    }
   ],
   "source": [
    "print(f\"Exact energy: {exact_energy}\")\n",
    "print(f\"Estimated energy: {result.energy + nuclear_repulsion_energy}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}