diff --git a/docs/tutorials/sample-based-quantum-diagonalization.ipynb b/docs/tutorials/sample-based-quantum-diagonalization.ipynb
index e0a29f29850..93207a8bece 100644
--- a/docs/tutorials/sample-based-quantum-diagonalization.ipynb
+++ b/docs/tutorials/sample-based-quantum-diagonalization.ipynb
@@ -10,10 +10,10 @@
"description: Use the sample-based quantum diagonalization algorithm to simulate a nitrogen molecule using noisy quantum hardware.\n",
"---\n",
"\n",
- "\n",
"{/* cspell:ignore fontdict fontsize milli LUCJ CCSD ccsd hcore pvdz */}\n",
"\n",
"# Sample-based quantum diagonalization of a chemistry Hamiltonian\n",
+ "\n",
"*Usage estimate: under one minute on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)*"
]
},
@@ -22,31 +22,34 @@
"id": "9701917b",
"metadata": {},
"source": [
- "## Background\n",
+ "## Learning outcomes\n",
"\n",
- "In this tutorial, we show how to post-process noisy quantum samples to find an approximation to the ground state of the nitrogen molecule $\\text{N}_2$ at equilibrium bond length, using the [sample-based quantum diagonalization (SQD) algorithm](https://arxiv.org/abs/2405.05068), for samples taken from a 59-qubit quantum circuit (52 system qubits + 7 ancilla qubits). A Qiskit-based implementation is available in the [SQD Qiskit addons](https://github.com/Qiskit/qiskit-addon-sqd), more details can be found in the corresponding [docs](/docs/guides/qiskit-addons-sqd) with a [simple example](/docs/guides/qiskit-addons-sqd-get-started) to get started.\n",
- "SQD is a technique for finding eigenvalues and eigenvectors of quantum operators, such as a quantum system Hamiltonian, using quantum and distributed classical computing together. Classical distributed computing is used to process samples obtained from a quantum processor, and to project and diagonalize a target Hamiltonian in a subspace spanned by them. This allows SQD to be robust to samples corrupted by quantum noise and deal with large Hamiltonians, such as chemistry Hamiltonians with millions of interaction terms, beyond the reach of any exact diagonalization methods. An SQD-based workflow has the following steps:\n",
+ "After going through this tutorial, users should understand:\n",
+ "- How to use the [SQD Qiskit addon](https://github.com/Qiskit/qiskit-addon-sqd) to approximate the ground state energy of a molecular system using bitstrings sampled from a quantum processing unit (QPU).\n",
+ "- How to use [ffsim](https://github.com/qiskit-community/ffsim) to construct a local unitary cluster Jastrow (LUCJ) circuit for quantum chemistry simulation.\n",
"\n",
- "1. Choose a circuit ansatz and apply it on a quantum computer to a reference state (in this case, the [Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) state).\n",
- "2. Sample bitstrings from the resulting quantum state.\n",
- "3. Run the *self-consistent configuration recovery* procedure on the bitstrings to obtain the approximation to the ground state.\n",
+ "## Prerequisites\n",
"\n",
- "SQD is known to work well when the target eigenstate is sparse: the wave function is supported in a set of basis states $\\mathcal{S} = \\{|x\\rangle \\}$ whose size does not increase exponentially with the size of the problem.\n",
+ "We recommend that users familiarize themselves with the following topics before going through this tutorial:\n",
+ "- Quantum chemistry and second quantization\n",
+ "- Using the Sampler primitive to sample from quantum circuits\n",
"\n",
- "### Quantum chemistry\n",
+ "## Background\n",
+ "In this tutorial, we show how to post-process noisy quantum samples to approximate the ground state of the nitrogen molecule $\\text{N}_2$ at equilibrium bond length, by using the [SQD Qiskit addon](https://github.com/Qiskit/qiskit-addon-sqd) to implement the [sample-based quantum diagonalization (SQD) algorithm](https://arxiv.org/abs/2405.05068). More details on the software can be found in the corresponding [documentation](/docs/guides/qiskit-addons-sqd), including a [simple example](/docs/guides/qiskit-addons-sqd-get-started) to get started.\n",
"\n",
- "The properties of molecules are largely determined by the structure of the electrons within them. As fermionic particles, electrons can be described using a mathematical formalism called second quantization. The idea is that there are a number of *orbitals*, each of which can be either empty or occupied by a fermion. A system of $N$ orbitals is described by a set of fermionic annihilation operators $\\{\\hat{a}_p\\}_{p=1}^N$ that satisfy the fermionic anticommutation relations,\n",
+ "This tutorial is recommended for users who are familiar with quantum chemistry: specifically, familiarity with finding ground state energies of a molecule. For a detailed walkthrough on the workflow, refer to the [quantum diagonalization algorithm course](https://quantum.cloud.ibm.com/learning/courses/quantum-diagonalization-algorithms).\n",
"\n",
- "$$\n",
- "\\begin{align*}\n",
- "\\hat{a}_p \\hat{a}_q + \\hat{a}_q \\hat{a}_p &= 0, \\\\\n",
- "\\hat{a}_p \\hat{a}_q^\\dagger + \\hat{a}_q^\\dagger \\hat{a}_p &= \\delta_{pq}.\n",
- "\\end{align*}\n",
- "$$\n",
"\n",
- "The adjoint $\\hat{a}_p^\\dagger$ is called a creation operator.\n",
+ "SQD is a technique for finding eigenvalues and eigenvectors of quantum operators, such as a quantum system Hamiltonian, by using quantum and distributed classical computing together. Classical distributed computing is used to process samples obtained from a quantum processor, and to project and diagonalize a target Hamiltonian in a subspace that they span. An SQD-based workflow has the following steps:\n",
"\n",
- "So far, our exposition has not accounted for spin, which is a fundamental property of fermions. When accounting for spin, the orbitals come in pairs called *spatial orbitals*. Each spatial orbital is composed of two *spin orbitals*, one that is labeled \"spin-$\\alpha$\" and one that is labeled \"spin-$\\beta$\". We then write $\\hat{a}_{p\\sigma}$ for the annihilation operator associated with the spin-orbital with spin $\\sigma$ ($\\sigma \\in \\{\\alpha, \\beta\\}$) in spatial orbital $p$. If we take $N$ to be the number of spatial orbitals, then there are a total of $2N$ spin-orbitals. The Hilbert space of this system is spanned by $2^{2N}$ orthonormal basis vectors labeled with two-part bitstrings $\\lvert z \\rangle = \\lvert z_\\beta z_\\alpha \\rangle = \\lvert z_{\\beta, N} \\cdots z_{\\beta, 1} z_{\\alpha, N} \\cdots z_{\\alpha, 1} \\rangle$.\n",
+ "1. Choose a circuit ansatz and apply it on a quantum computer to a reference state (in this case, the [Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) state).\n",
+ "2. Sample bitstrings from the resulting quantum state.\n",
+ "3. Run the *self-consistent configuration recovery* procedure on the bitstrings to obtain the ground state approximation.\n",
+ "\n",
+ "SQD is known to work well when the target eigenstate is sparse: the wave function is supported in a set of basis states $\\mathcal{S} = \\{|x\\rangle \\}$ whose size does not increase exponentially with the size of the problem.\n",
+ "\n",
+ "\n",
+ "### Quantum chemistry\n",
"\n",
"The Hamiltonian of a molecular system can be written as\n",
"\n",
@@ -63,72 +66,24 @@
"\n",
"where the $h_{pr}$ and $h_{prqs}$ are complex numbers called molecular integrals that can be calculated from the specification of the molecule using a computer program. In this tutorial, we compute the integrals using the [PySCF](https://pyscf.org/) software package.\n",
"\n",
- "For details about how the molecular Hamiltonian is derived, consult a textbook on quantum chemistry (for example, *Modern Quantum Chemistry* by Szabo and Ostlund). For a high-level explanation of how quantum chemistry problems are mapped onto quantum computers, check out the lecture [*Mapping Problems to Qubits*](https://youtube.com/watch?v=TyFU6r8uEsE&t=900) from Qiskit Global Summer School 2024.\n",
+ "For details about how the molecular Hamiltonian is derived, consult a textbook on quantum chemistry (for example, *Modern Quantum Chemistry* by Szabo and Ostlund). For a high-level explanation of how quantum chemistry problems are mapped onto quantum computers, check out the lecture [*Mapping Problems to Qubits*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900) from Qiskit Global Summer School 2024.\n",
"\n",
"### Local unitary cluster Jastrow (LUCJ) ansatz\n",
"\n",
- "SQD requires a quantum circuit ansatz to draw samples from. In this tutorial, we'll use the [local unitary cluster Jastrow (LUCJ)](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) ansatz due to its combination of physical motivation and hardware-friendliness.\n",
- "\n",
- "The LUCJ ansatz is a specialized form of the general unitary cluster Jastrow (UCJ) ansatz, which has the form\n",
- "\n",
- "$$\n",
- " \\lvert \\Psi \\rangle = \\prod_{\\mu=1}^L e^{\\hat{K}_\\mu} e^{i \\hat{J}_\\mu} e^{-\\hat{K}_\\mu} | \\Phi_0 \\rangle\n",
- "$$\n",
- "\n",
- "where $\\lvert \\Phi_0 \\rangle$ is a reference state, often taken to be the Hartree-Fock state, and the $\\hat{K}_\\mu$ and $\\hat{J}_\\mu$ have the form\n",
- "\n",
- "$$\n",
- "\\hat{K}_\\mu = \\sum_{pq, \\sigma} K_{pq}^\\mu \\, \\hat{a}^\\dagger_{p \\sigma} \\hat{a}^{\\phantom{\\dagger}}_{q \\sigma}\n",
- "\\;,\\;\n",
- "\\hat{J}_\\mu = \\sum_{pq, \\sigma\\tau} J_{pq,\\sigma\\tau}^\\mu \\, \\hat{n}_{p \\sigma} \\hat{n}_{q \\tau}\n",
- "\\;,\n",
- "$$\n",
- "\n",
- "where we have defined the number operator\n",
- "\n",
- "$$\n",
- "\\hat{n}_{p \\sigma} = \\hat{a}^\\dagger_{p \\sigma} \\hat{a}^{\\phantom{\\dagger}}_{p \\sigma}.\n",
- "$$\n",
- "\n",
- "The operator $e^{\\hat{K}_\\mu}$ is an orbital rotation, which can be implemented using known algorithms in linear depth and using linear connectivity.\n",
- "Implementing the $e^{i \\mathcal{J}_k}$ term of the UCJ ansatz requires either all-to-all connectivity or the use of a fermionic swap network, making it challenging for noisy pre-fault-tolerant quantum processors that have limited connectivity. The idea of the *local* UCJ ansatz is to impose sparsity constraints on the $\\mathbf{J}^{\\alpha\\alpha}$ and $\\mathbf{J}^{\\alpha\\beta}$ matrices which allow them to be implemented in constant depth on qubit topologies with limited connectivity. The constraints are specified by a list of indices indicating which matrix entries in the upper triangle are allowed to be nonzero (since the matrices are symmetric, only the upper triangle needs to be specified). These indices can be interpreted as pairs of orbitals that are allowed to interact.\n",
- "\n",
- "As an example, consider a square lattice qubit topology. We can place the $\\alpha$ and $\\beta$ orbitals in parallel lines on the lattice, with connections between these lines forming \"rungs\" of a ladder shape, like this:\n",
- "\n",
- "\n",
- "\n",
- "With this setup, orbitals with the same spin are connected with a line topology, while orbitals with different spins are connected when they share the same spatial orbital. This yields the following index constraints on the $\\mathbf{J}$ matrices:\n",
- "\n",
- "$$\n",
- "\\begin{align*}\n",
- "\\mathbf{J}^{\\alpha\\alpha} &: \\set{(p, p+1) \\; , \\; p = 0, \\ldots, N-2} \\\\\n",
- "\\mathbf{J}^{\\alpha\\beta} &: \\set{(p, p) \\;, \\; p = 0, \\ldots, N-1}\n",
- "\\end{align*}\n",
- "$$\n",
- "\n",
- "In other words, if the $\\mathbf{J}$ matrices are nonzero only at the specified indices in the upper triangle, then the $e^{i \\mathcal{J}_k}$ term can be implemented on a square topology without using any swap gates, in constant depth. Of course, imposing such constraints on the ansatz makes it less expressive, so more ansatz repetitions may be required.\n",
- "\n",
- "The IBM® hardware has a heavy-hex lattice qubit topology, in which case we can adopt a \"zig-zag\" pattern, depicted below. In this pattern, orbitals with the same spin are mapped to qubits with a line topology (red and blue circles), and a connection between orbitals of different spin is present at every 4th spatial orbital, with the connection being facilitated by an ancilla qubit (purple circles). In this case, the index constraints are\n",
+ "SQD requires a quantum circuit ansatz to draw samples from. In this tutorial, we'll use the [local unitary cluster Jastrow (LUCJ)](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) ansatz due to its combination of physical motivation and hardware-friendliness. We'll use [ffsim](https://qiskit-community.github.io/ffsim/) to construct the ansatz circuit.\n",
"\n",
- "\n",
- "$$\n",
- "\\begin{align*}\n",
- "\\mathbf{J}^{\\alpha\\alpha} &: \\set{(p, p+1) \\; , \\; p = 0, \\ldots, N-2} \\\\\n",
- "\\mathbf{J}^{\\alpha\\beta} &: \\set{(p, p) \\;, \\; p = 0, 4, 8, \\ldots (p \\leq N-1)}\n",
- "\\end{align*}\n",
- "$$\n",
+ "The LUCJ ansatz adapts to QPUs with restricted qubit connectivity. The spin orbitals are mapped to qubits such that the ansatz does not require routing with SWAP gates. IBM® hardware has a heavy-hex lattice qubit topology, in which case we can adopt a \"zig-zag\" pattern, depicted below. In this pattern, orbitals with the same spin are mapped to qubits with a line topology (red and blue circles), and a connection between orbitals of different spin is present at every 4th spatial orbital, with the connection being facilitated by an ancilla qubit (purple circles).\n",
"\n",
"\n",
"\n",
- "\n",
"### Self-consistent configuration recovery\n",
"\n",
"The self-consistent configuration recovery procedure is designed to extract as much signal as possible from noisy quantum samples. Because the molecular Hamiltonian conserves particle number and spin Z, it makes sense to choose a circuit ansatz that also conserves these symmetries. When applied to the Hartree-Fock state, the resulting state has a fixed particle number and spin Z in the noiseless setting. Therefore, the spin-$\\alpha$ and spin-$\\beta$ halves of any bitstring sampled from this state should have the same [Hamming weight](https://en.wikipedia.org/wiki/Hamming_weight) as in the Hartree-Fock state. Due to the presence of noise in current quantum processors, some measured bitstrings will violate this property. A simple form of postselection would discard these bitstrings, but this is wasteful because those bitstrings might still contain some signal. The self-consistent recovery procedure attempts to recover some of that signal in post-processing. The procedure is iterative and requires as input an estimate of the average occupancies of each orbital in the ground state, which is first computed from the raw samples. The procedure is run in a loop, and each iteration has the following steps:\n",
"\n",
- "1. For each bitstring that violates the specified symmetries, flip its bits with a probabilistic procedure designed to bring the bitstring closer to the current estimate of the average orbital occupancies, to obtain a new bitstring.\n",
- "2. Collect all of the old and new bitstrings that satisfy the symmetries, and subsample subsets of a fixed size, chosen in advance.\n",
- "3. For each subset of bitstrings, project the Hamiltonian into the subspace spanned by the corresponding basis vectors (see the [previous section](#quantum-chemistry) for a description of these basis vectors), and compute a ground state estimate of the projected Hamiltonian on a classical computer.\n",
- "4. Update the estimate of the average orbital occupancies with the ground state estimate with the lowest energy.\n",
+ "1. For each bitstring that violates the specified symmetries, flip its bits with a probabilistic procedure designed to bring the bitstring closer to the current estimate of the average orbital occupancies, to obtain a new bitstring.\n",
+ "2. Collect all of the old and new bitstrings that satisfy the symmetries, and subsample subsets of a fixed size, chosen in advance.\n",
+ "3. For each subset of bitstrings, project the Hamiltonian into the subspace spanned by the corresponding basis vectors (see the [previous section](#quantum-chemistry) for a description of these basis vectors), and compute a ground state estimate of the projected Hamiltonian on a classical computer.\n",
+ "4. Update the estimate of the average orbital occupancies with the ground state estimate with the lowest energy.\n",
"\n",
"### SQD workflow diagram\n",
"\n",
@@ -146,10 +101,10 @@
"\n",
"Before starting this tutorial, ensure that you have the following installed:\n",
"\n",
- "- Qiskit SDK v1.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
- "- Qiskit Runtime v0.22 or later (`pip install qiskit-ibm-runtime`)\n",
- "- SQD Qiskit addon v0.11 or later (`pip install qiskit-addon-sqd`)\n",
- "- ffsim v0.0.58 or later (`pip install ffsim`)"
+ "* Qiskit SDK v1.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
+ "* Qiskit Runtime v0.22 or later (`pip install qiskit-ibm-runtime`)\n",
+ "* SQD Qiskit addon v0.11 or later (`pip install qiskit-addon-sqd`)\n",
+ "* ffsim v0.0.72 or later (`pip install ffsim`)"
]
},
{
@@ -176,19 +131,31 @@
"import pyscf.cc\n",
"import pyscf.mcscf\n",
"from qiskit import QuantumCircuit, QuantumRegister\n",
+ "from qiskit.primitives import StatevectorSampler\n",
+ "from qiskit.providers.fake_provider import GenericBackendV2\n",
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"from qiskit_ibm_runtime import SamplerV2 as Sampler"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "4bc6ee26-4371-4cd2-80a7-60752bf8775d",
+ "metadata": {},
+ "source": [
+ "## Small-scale simulator example\n",
+ "\n",
+ "In this tutorial, we will find an approximation to the ground state of a nitrogen molecule at equilibrium. We first use a small STO-6G basis set so that we can simulate the experiment and make sure it's working."
+ ]
+ },
{
"cell_type": "markdown",
"id": "afeb054c",
"metadata": {},
"source": [
- "## Step 1: Map classical inputs to a quantum problem\n",
+ "### Step 1: Map classical inputs to a quantum problem\n",
"\n",
- "In this tutorial, we will find an approximation to the ground state of the molecule at equilibrium in the cc-pVDZ basis set. First, we specify the molecule and its properties."
+ "First, we specify the molecule and its properties."
]
},
{
@@ -201,7 +168,10 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "converged SCF energy = -108.929838385609\n"
+ "converged SCF energy = -108.464957764796\n",
+ "CASCI E = -108.595987350986 E(CI) = -32.4115475088426 S^2 = 0.0000000\n",
+ "norb = 8\n",
+ "nelec = (5, 5)\n"
]
}
],
@@ -214,7 +184,7 @@
"mol = pyscf.gto.Mole()\n",
"mol.build(\n",
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n",
- " basis=\"cc-pvdz\",\n",
+ " basis=\"sto-6g\",\n",
" symmetry=\"Dooh\",\n",
")\n",
"\n",
@@ -228,13 +198,17 @@
"n_electrons = int(sum(scf.mo_occ[active_space]))\n",
"n_alpha = (n_electrons + mol.spin) // 2\n",
"n_beta = (n_electrons - mol.spin) // 2\n",
- "cas = pyscf.mcscf.CASCI(scf, norb, (n_alpha, n_beta))\n",
+ "nelec = (n_alpha, n_beta)\n",
+ "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\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), norb)\n",
"\n",
- "# Store reference energy from SCI calculation performed separately\n",
- "reference_energy = -109.22802921665716"
+ "# Compute exact energy using FCI\n",
+ "reference_energy = cas.run().e_tot\n",
+ "\n",
+ "print(f\"norb = {norb}\")\n",
+ "print(f\"nelec = {nelec}\")"
]
},
{
@@ -255,7 +229,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "E(CCSD) = -109.2177884185545 E_corr = -0.2879500329450042\n"
+ "E(CCSD) = -108.5933309085008 E_corr = -0.1283731437052352\n"
]
}
],
@@ -301,8 +275,6 @@
" options=dict(maxiter=1000),\n",
")\n",
"\n",
- "nelec = (n_alpha, n_beta)\n",
- "\n",
"# create an empty quantum circuit\n",
"qubits = QuantumRegister(2 * norb, name=\"q\")\n",
"circuit = QuantumCircuit(qubits)\n",
@@ -320,117 +292,441 @@
"id": "db11bf6d",
"metadata": {},
"source": [
- "## Step 2: Optimize problem for quantum hardware execution"
+ "### Step 2: Optimize for quantum hardware execution\n",
+ "\n",
+ "Next, we optimize the circuit for a target hardware. For now, we'll create a generic backend with a specified number of qubits and a gate set that the LUCJ ansatz naturally decomposes to."
]
},
{
- "cell_type": "markdown",
- "id": "0760b3f3",
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "53a039d8",
"metadata": {},
+ "outputs": [],
"source": [
- "Next, we optimize the circuit for a target hardware."
+ "backend = GenericBackendV2(\n",
+ " 2 * norb, basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\"]\n",
+ ")"
]
},
{
"cell_type": "code",
- "execution_count": 5,
- "id": "53a039d8",
+ "execution_count": 6,
+ "id": "7d554aa5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
- "Using backend ibm_fez\n"
+ "Gate counts: OrderedDict({'xx_plus_yy': 86, 'cp': 16, 'p': 16, 'measure': 16, 'x': 10, 'barrier': 1})\n"
]
}
],
"source": [
- "service = QiskitRuntimeService()\n",
- "backend = service.least_busy(\n",
- " operational=True, simulator=False, min_num_qubits=133\n",
+ "pass_manager = generate_preset_pass_manager(\n",
+ " optimization_level=3, backend=backend\n",
")\n",
- "\n",
- "print(f\"Using backend {backend.name}\")"
+ "# Set the pre-initialization stage of the pass manager with passes suggested by ffsim\n",
+ "pass_manager.pre_init = ffsim.qiskit.PRE_INIT\n",
+ "isa_circuit = pass_manager.run(circuit)\n",
+ "print(f\"Gate counts: {isa_circuit.count_ops()}\")"
]
},
{
"cell_type": "markdown",
- "id": "057ebbf6",
+ "id": "0cc1edef",
"metadata": {},
"source": [
- "We recommend the following steps to optimize the ansatz and make it hardware-compatible.\n",
+ "### Step 3: Execute using Qiskit primitives\n",
"\n",
- "- Select physical qubits (`initial_layout`) from the target hardware that adheres to the \"zig-zag\" pattern described previously. Laying out qubits in this pattern leads to an efficient hardware-compatible circuit with less gates. Here, we include some code to automate the selection of the \"zig-zag\" pattern that uses a scoring heuristic to minimize the errors associated with the selected layout.\n",
- "- Generate a staged pass manager using the [generate_preset_pass_manager](/docs/api/qiskit/qiskit.transpiler.generate_preset_pass_manager) function from qiskit with your choice of `backend` and `initial_layout`.\n",
- "- Set the `pre_init` stage of your staged pass manager to `ffsim.qiskit.PRE_INIT`. `ffsim.qiskit.PRE_INIT` includes qiskit transpiler passes that decompose gates into orbital rotations and then merges the orbital rotations, resulting in fewer gates in the final circuit.\n",
- "- Run the pass manager on your circuit."
+ "After optimizing the circuit for hardware execution, we are ready to run it on the target hardware and collect samples for ground state energy estimation. As we only have one circuit, we will use Qiskit Runtime's [Job execution mode](/docs/guides/execution-modes) and execute our circuit."
]
},
{
- "cell_type": "markdown",
- "id": "59eaa216",
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "93c1cef3-298e-4deb-8512-769fe94cd5a5",
"metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/kjs/projects/documentation/.venv/lib/python3.12/site-packages/qiskit/circuit/quantumcircuit.py:4625: UserWarning: Trying to add QuantumRegister to a QuantumCircuit having a layout\n",
+ " circ.add_register(qreg)\n"
+ ]
+ }
+ ],
"source": [
- "\n",
- " \n",
- " \n",
- " \n",
- ""
+ "rng = np.random.default_rng()\n",
+ "sampler = StatevectorSampler(seed=rng)\n",
+ "job = sampler.run([isa_circuit], shots=100_000)"
]
},
{
"cell_type": "code",
- "execution_count": 6,
- "id": "d9f56529",
- "metadata": {
- "tags": [
- "id-automated-layout-code"
- ]
- },
+ "execution_count": 8,
+ "id": "332ecab3-77e6-473f-b0e7-af30f983393a",
+ "metadata": {},
"outputs": [],
"source": [
- "from typing import Sequence\n",
- "\n",
- "import rustworkx\n",
- "from qiskit.providers import BackendV2\n",
- "from rustworkx import NoEdgeBetweenNodes, PyGraph\n",
+ "primitive_result = job.result()\n",
+ "pub_result = primitive_result[0]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6df05b6e",
+ "metadata": {},
+ "source": [
+ "### Step 4: Post-process and return result in desired classical format\n",
"\n",
- "IBM_TWO_Q_GATES = {\"cx\", \"ecr\", \"cz\"}\n",
+ "A useful metric to judge the quality of the QPU output is the number of valid configurations returned. A valid configuration has the correct particle number and spin Z, which means that the right half of the bitstring has Hamming weight equal to the number of spin-up electrons, and the left half has Hamming weight equal to the number of spin-down electrons. The following cell computes the fraction of sampled configurations that are valid."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "718f8517",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Fraction of sampled configurations that are valid: 1.0\n"
+ ]
+ }
+ ],
+ "source": [
+ "def is_valid_bitstring(\n",
+ " bitstring: str, norb: int, nelec: tuple[int, int]\n",
+ ") -> bool:\n",
+ " n_alpha, n_beta = nelec\n",
+ " return (\n",
+ " len(bitstring) == 2 * norb\n",
+ " and bitstring[norb:].count(\"1\") == n_alpha\n",
+ " and bitstring[:norb].count(\"1\") == n_beta\n",
+ " )\n",
"\n",
"\n",
- "def create_linear_chains(num_orbitals: int) -> PyGraph:\n",
- " \"\"\"In zig-zag layout, there are two linear chains (with connecting qubits between\n",
- " the chains). This function creates those two linear chains: a rustworkx PyGraph\n",
- " with two disconnected linear chains. Each chain contains `num_orbitals` number\n",
- " of nodes, that is, in the final graph there are `2 * num_orbitals` number of nodes.\n",
+ "bit_array = pub_result.data.meas\n",
+ "num_valid = sum(\n",
+ " is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
+ ")\n",
+ "valid_fraction = num_valid / bit_array.num_shots\n",
+ "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7b126f3a",
+ "metadata": {},
+ "source": [
+ "All of the bitstrings are valid because we are sampling the circuit on a noiseless simulator. When running on a noisy QPU, the fraction will be less than one, but hopefully it will be larger than the fraction you would expect if the bitstrings were sampled uniformly at random, which is computed in the following cell."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "6b3e4bca",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625\n"
+ ]
+ }
+ ],
+ "source": [
+ "expected_fraction_random = (\n",
+ " math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
+ ")\n",
+ "print(\n",
+ " f\"Expected fraction of valid configurations from uniformly random bitstrings: {expected_fraction_random}\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "eb704101-0fe8-4d12-b572-b1d844e35a90",
+ "metadata": {},
+ "source": [
+ "Now, we estimate the ground state energy of the Hamiltonian using the `diagonalize_fermionic_hamiltonian` function. This function performs the self-consistent configuration recovery procedure to iteratively refine the noisy quantum samples to improve the energy estimate. We pass a callback function so that we can save the intermediate results for later analysis. See the [API documentation](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) for explanations of the arguments to `diagonalize_fermionic_hamiltonian`.\n",
"\n",
- " Args:\n",
- " num_orbitals (int): Number orbitals or nodes in each linear chain. They are\n",
- " also known as alpha-alpha interaction qubits.\n",
+ "Here, we use the `initial_occupancies` argument to `diagonalize_fermionic_hamiltonian` to specify the Hartree-Fock configuration as the initial guess for the orbital occupancies in the ground state. This approach is sensible for systems where the ground state has significant support on the Hartree-Fock configuration, but it might not be appropriate in other situations, though more advanced computational methods might yield better initial guesses in those cases. Specifying `initial_occupancies` also allows configuration recovery to run even if no valid configurations were sampled, as may be the case when sampling a large circuit on a noisy QPU. Without this argument, the configuration recovery would fail and raise an error if no valid configurations were provided."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "2f32a352",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Iteration 1\n",
+ "\tSubsample 0\n",
+ "\t\tEnergy: -108.59570174319286\n",
+ "\t\tSubspace dimension: 1849\n",
+ "\tSubsample 1\n",
+ "\t\tEnergy: -108.59570174319286\n",
+ "\t\tSubspace dimension: 1849\n",
+ "\tSubsample 2\n",
+ "\t\tEnergy: -108.59570174319286\n",
+ "\t\tSubspace dimension: 1849\n",
+ "Iteration 2\n",
+ "\tSubsample 0\n",
+ "\t\tEnergy: -108.59570174319286\n",
+ "\t\tSubspace dimension: 1849\n",
+ "\tSubsample 1\n",
+ "\t\tEnergy: -108.59570174319286\n",
+ "\t\tSubspace dimension: 1849\n",
+ "\tSubsample 2\n",
+ "\t\tEnergy: -108.59570174319286\n",
+ "\t\tSubspace dimension: 1849\n",
+ "Final energy: -108.59570174319286\n",
+ "Final energy error: 0.0002856077932023027\n"
+ ]
+ }
+ ],
+ "source": [
+ "from functools import partial\n",
"\n",
- " Returns:\n",
- " A rustworkx.PyGraph with two disconnected linear chains each with `num_orbitals`\n",
- " number of nodes.\n",
- " \"\"\"\n",
- " G = rustworkx.PyGraph()\n",
+ "from qiskit_addon_sqd.fermion import (\n",
+ " SCIResult,\n",
+ " diagonalize_fermionic_hamiltonian,\n",
+ " solve_sci_batch,\n",
+ ")\n",
"\n",
- " for n in range(num_orbitals):\n",
- " G.add_node(n)\n",
+ "# SQD options\n",
+ "energy_tol = 1e-3\n",
+ "occupancies_tol = 1e-3\n",
+ "max_iterations = 5\n",
"\n",
- " for n in range(num_orbitals - 1):\n",
- " G.add_edge(n, n + 1, None)\n",
+ "# Eigenstate solver options\n",
+ "num_batches = 3\n",
+ "samples_per_batch = 300\n",
+ "symmetrize_spin = True\n",
+ "carryover_threshold = 1e-4\n",
+ "max_cycle = 200\n",
"\n",
- " for n in range(num_orbitals, 2 * num_orbitals):\n",
- " G.add_node(n)\n",
+ "# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies\n",
+ "initial_occupancies = (\n",
+ " np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
+ " np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
+ ")\n",
"\n",
- " for n in range(num_orbitals, 2 * num_orbitals - 1):\n",
- " G.add_edge(n, n + 1, None)\n",
+ "# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
+ "# you can omit this step, in which case you would not specify the sci_solver argument\n",
+ "# in the call to diagonalize_fermionic_hamiltonian below.\n",
+ "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n",
"\n",
- " return G\n",
+ "# List to capture intermediate results\n",
+ "result_history = []\n",
"\n",
"\n",
- "def create_lucj_zigzag_layout(\n",
+ "def callback(results: list[SCIResult]):\n",
+ " result_history.append(results)\n",
+ " iteration = len(result_history)\n",
+ " print(f\"Iteration {iteration}\")\n",
+ " for i, result in enumerate(results):\n",
+ " print(f\"\\tSubsample {i}\")\n",
+ " print(f\"\\t\\tEnergy: {result.energy + nuclear_repulsion_energy}\")\n",
+ " print(\n",
+ " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
+ " )\n",
+ "\n",
+ "\n",
+ "result = diagonalize_fermionic_hamiltonian(\n",
+ " hcore,\n",
+ " eri,\n",
+ " bit_array,\n",
+ " samples_per_batch=samples_per_batch,\n",
+ " norb=norb,\n",
+ " nelec=nelec,\n",
+ " num_batches=num_batches,\n",
+ " energy_tol=energy_tol,\n",
+ " occupancies_tol=occupancies_tol,\n",
+ " max_iterations=max_iterations,\n",
+ " sci_solver=sci_solver,\n",
+ " symmetrize_spin=symmetrize_spin,\n",
+ " initial_occupancies=initial_occupancies,\n",
+ " carryover_threshold=carryover_threshold,\n",
+ " callback=callback,\n",
+ " seed=rng,\n",
+ ")\n",
+ "\n",
+ "final_energy = result.energy + nuclear_repulsion_energy\n",
+ "energy_error = final_energy - reference_energy\n",
+ "print(f\"Final energy: {final_energy}\")\n",
+ "print(f\"Final energy error: {energy_error}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9d78906b-4759-4506-9c69-85d4e67766b3",
+ "metadata": {},
+ "source": [
+ "#### Visualize the results\n",
+ "\n",
+ "The first plot shows that in this simulation we are already within `1 mH` of the exact answer after the first iteration (chemical accuracy is typically accepted to be `1 kcal/mol` $\\approx$ `1.6 mH`). This is a small system, though, and because the samples are noiseless, configuration recovery is not needed. On a larger system run on a noisy QPU, multiple configuration recovery iterations might be needed, and the final accuracy might be worse. Generally, the energy can be improved by allowing more configuration recovery iterations or by increasing the number of samples per batch.\n",
+ "\n",
+ "The second plot shows the average occupancy of each spatial orbital after the final iteration. We can see that both the spin-up and spin-down electrons occupy the first five orbitals with high probability in our solutions."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "caffd888-e89c-4aa9-8bae-4d1bb723b35e",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# Data for energies plot\n",
+ "x1 = range(len(result_history))\n",
+ "min_e = [\n",
+ " min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
+ " for result in result_history\n",
+ "]\n",
+ "e_diff = [abs(e - reference_energy) for e in min_e]\n",
+ "yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n",
+ "\n",
+ "# Chemical accuracy (+/- 1 milli-Hartree)\n",
+ "chem_accuracy = 0.001\n",
+ "\n",
+ "# Data for avg spatial orbital occupancy\n",
+ "y2 = np.sum(result.orbital_occupancies, axis=0)\n",
+ "x2 = range(len(y2))\n",
+ "\n",
+ "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
+ "\n",
+ "# Plot energies\n",
+ "axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n",
+ "axs[0].set_xticks(x1)\n",
+ "axs[0].set_xticklabels(x1)\n",
+ "axs[0].set_yticks(yt1)\n",
+ "axs[0].set_yticklabels(yt1)\n",
+ "axs[0].set_yscale(\"log\")\n",
+ "axs[0].set_ylim(1e-4)\n",
+ "axs[0].axhline(\n",
+ " y=chem_accuracy,\n",
+ " color=\"#BF5700\",\n",
+ " linestyle=\"--\",\n",
+ " label=\"chemical accuracy\",\n",
+ ")\n",
+ "axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
+ "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
+ "axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n",
+ "axs[0].legend()\n",
+ "\n",
+ "# Plot orbital occupancy\n",
+ "axs[1].bar(x2, y2, width=0.8)\n",
+ "axs[1].set_xticks(x2)\n",
+ "axs[1].set_xticklabels(x2)\n",
+ "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
+ "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
+ "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
+ "\n",
+ "plt.tight_layout()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ce0eecb3-8a23-4118-aa1e-a28afcec6334",
+ "metadata": {},
+ "source": [
+ "## Large-scale hardware example\n",
+ "\n",
+ "Now, we run a larger example on real quantum hardware. Here, we'll derive an active space for the nitrogen molecule from the cc-pVDZ basis set.\n",
+ "\n",
+ "When executing the LUCJ circuit on hardware, you should carefully map the qubits onto the hardware to minimize the overhead of qubit routing with SWAP gates. We recommend the following steps to optimize the ansatz and make it hardware-compatible.\n",
+ "\n",
+ "* Select physical qubits (`initial_layout`) from the target hardware that adhere to the \"zig-zag\" pattern described previously. Laying out qubits in this pattern leads to an efficient hardware-compatible circuit with fewer gates. Here, we include some code to automate the selection of the \"zig-zag\" pattern that uses a scoring heuristic to minimize the errors associated with the selected layout.\n",
+ "* Generate a staged pass manager using the [generate\\_preset\\_pass\\_manager](/docs/api/qiskit/qiskit.transpiler.generate_preset_pass_manager) function from qiskit with your choice of `backend` and `initial_layout`.\n",
+ "* Set the `pre_init` stage of your staged pass manager to `ffsim.qiskit.PRE_INIT`. `ffsim.qiskit.PRE_INIT` includes qiskit transpiler passes that decompose gates into orbital rotations and then merges the orbital rotations, resulting in fewer gates in the final circuit.\n",
+ "* Run the pass manager on your circuit."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "195e8307",
+ "metadata": {},
+ "source": [
+ "\n",
+ " \n",
+ " \n",
+ " \n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "c2a20804",
+ "metadata": {
+ "tags": [
+ "id-automated-layout-code"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "from typing import Sequence\n",
+ "\n",
+ "import rustworkx\n",
+ "from qiskit.providers import BackendV2\n",
+ "from rustworkx import NoEdgeBetweenNodes, PyGraph\n",
+ "\n",
+ "IBM_TWO_Q_GATES = {\"cx\", \"ecr\", \"cz\"}\n",
+ "\n",
+ "\n",
+ "def create_linear_chains(num_orbitals: int) -> PyGraph:\n",
+ " \"\"\"In zig-zag layout, there are two linear chains (with connecting qubits between\n",
+ " the chains). This function creates those two linear chains: a rustworkx PyGraph\n",
+ " with two disconnected linear chains. Each chain contains `num_orbitals` number\n",
+ " of nodes, that is, in the final graph there are `2 * num_orbitals` number of nodes.\n",
+ "\n",
+ " Args:\n",
+ " num_orbitals (int): Number orbitals or nodes in each linear chain. They are\n",
+ " also known as alpha-alpha interaction qubits.\n",
+ "\n",
+ " Returns:\n",
+ " A rustworkx.PyGraph with two disconnected linear chains each with `num_orbitals`\n",
+ " number of nodes.\n",
+ " \"\"\"\n",
+ " G = rustworkx.PyGraph()\n",
+ "\n",
+ " for n in range(num_orbitals):\n",
+ " G.add_node(n)\n",
+ "\n",
+ " for n in range(num_orbitals - 1):\n",
+ " G.add_edge(n, n + 1, None)\n",
+ "\n",
+ " for n in range(num_orbitals, 2 * num_orbitals):\n",
+ " G.add_node(n)\n",
+ "\n",
+ " for n in range(num_orbitals, 2 * num_orbitals - 1):\n",
+ " G.add_edge(n, n + 1, None)\n",
+ "\n",
+ " return G\n",
+ "\n",
+ "\n",
+ "def create_lucj_zigzag_layout(\n",
" num_orbitals: int, backend_coupling_graph: PyGraph\n",
") -> tuple[PyGraph, int]:\n",
" \"\"\"This function creates the complete zigzag graph that 'can be mapped' to an IBM QPU with\n",
@@ -603,242 +899,206 @@
" return layouts[0][:-num_alpha_beta_qubits], num_alpha_beta_qubits"
]
},
- {
- "cell_type": "code",
- "execution_count": 7,
- "id": "7d554aa5",
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Gate counts (w/o pre-init passes): OrderedDict({'rz': 12407, 'sx': 12155, 'cz': 3986, 'x': 581, 'measure': 52, 'barrier': 1})\n",
- "Gate counts (w/ pre-init passes): OrderedDict({'sx': 7059, 'rz': 6980, 'cz': 1876, 'measure': 52, 'x': 35, 'barrier': 1})\n"
- ]
- }
- ],
- "source": [
- "initial_layout, _ = get_zigzag_physical_layout(norb, backend=backend)\n",
- "\n",
- "pass_manager = generate_preset_pass_manager(\n",
- " optimization_level=3, backend=backend, initial_layout=initial_layout\n",
- ")\n",
- "\n",
- "# without PRE_INIT passes\n",
- "isa_circuit = pass_manager.run(circuit)\n",
- "print(f\"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}\")\n",
- "\n",
- "# with PRE_INIT passes\n",
- "# We will use the circuit generated by this pass manager for hardware execution\n",
- "pass_manager.pre_init = ffsim.qiskit.PRE_INIT\n",
- "isa_circuit = pass_manager.run(circuit)\n",
- "print(f\"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}\")"
- ]
- },
{
"cell_type": "markdown",
- "id": "0cc1edef",
+ "id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062",
"metadata": {},
"source": [
- "## Step 3: Execute using Qiskit primitives"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cbf7ef9f",
- "metadata": {},
- "source": [
- "After optimizing the circuit for hardware execution, we are ready to run it on the target hardware and collect samples for ground state energy estimation. As we only have one circuit, we will use Qiskit Runtime's [Job execution mode](/docs/guides/execution-modes) and execute our circuit."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 8,
- "id": "93c1cef3-298e-4deb-8512-769fe94cd5a5",
- "metadata": {},
- "outputs": [],
- "source": [
- "sampler = Sampler(mode=backend)\n",
- "job = sampler.run([isa_circuit], shots=100_000)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 9,
- "id": "332ecab3-77e6-473f-b0e7-af30f983393a",
- "metadata": {},
- "outputs": [],
- "source": [
- "primitive_result = job.result()\n",
- "pub_result = primitive_result[0]"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "6df05b6e",
- "metadata": {},
- "source": [
- "## Step 4: Post-process and return result in desired classical format"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a8736983",
- "metadata": {},
- "source": [
- "A useful metric to judge the quality of the QPU output is the number of valid configurations returned. A valid configuration has the correct particle number and spin Z, which means that the right half of the bitstring has Hamming weight equal to the number of spin-up electrons, and the left half has Hamming weight equal to the number of spin-down electrons. The following cell computes the fraction of sampled configurations that are valid."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 10,
- "id": "718f8517",
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Fraction of sampled configurations that are valid: 0.00065\n"
- ]
- }
- ],
- "source": [
- "def is_valid_bitstring(\n",
- " bitstring: str, norb: int, nelec: tuple[int, int]\n",
- ") -> bool:\n",
- " n_alpha, n_beta = nelec\n",
- " return (\n",
- " len(bitstring) == 2 * norb\n",
- " and bitstring[norb:].count(\"1\") == n_alpha\n",
- " and bitstring[:norb].count(\"1\") == n_beta\n",
- " )\n",
- "\n",
+ "### Steps 1-4\n",
"\n",
- "bit_array = pub_result.data.meas\n",
- "num_valid = sum(\n",
- " is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
- ")\n",
- "valid_fraction = num_valid / bit_array.num_shots\n",
- "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7b126f3a",
- "metadata": {},
- "source": [
- "While this fraction might seem fairly small, it is significantly higher than the fraction of bitstrings you would expect to be valid if the bitstrings were sampled uniformly at random, as shown by the following cell."
+ "Here we put all the steps together into a single workflow at a larger scale, which is then run on real quantum hardware."
]
},
{
"cell_type": "code",
- "execution_count": 11,
- "id": "6b3e4bca",
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n"
- ]
- }
- ],
- "source": [
- "expected_fraction_random = (\n",
- " math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
- ")\n",
- "print(\n",
- " f\"Expected fraction of valid configurations from uniformly random bitstrings: {expected_fraction_random}\"\n",
- ")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "eb704101-0fe8-4d12-b572-b1d844e35a90",
- "metadata": {},
- "source": [
- "Now, we estimate the ground state energy of the Hamiltonian using the `diagonalize_fermionic_hamiltonian` function. This function performs the self-consistent configuration recovery procedure to iteratively refine the noisy quantum samples to improve the energy estimate. We pass a callback function so that we can save the intermediate results for later analysis. See the [API documentation](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) for explanations of the arguments to `diagonalize_fermionic_hamiltonian`.\n",
- "\n",
- "Here, we use the `initial_occupancies` argument to `diagonalize_fermionic_hamiltonian` to specify the Hartree-Fock configuration as the initial guess for the orbital occupancies in the ground state. This approach is sensible for systems where the ground state has significant support on the Hartree-Fock configuration, but it might not be appropriate in other situations, though more advanced computational methods might yield better initial guesses in those cases. Specifying `initial_occupancies` also allows configuration recovery to run even if no valid configurations were sampled, as may be the case when sampling a large circuit on a noisy QPU. Without this argument, the configuration recovery would fail and raise an error if no valid configurations were provided."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 12,
- "id": "2f32a352",
+ "execution_count": 14,
+ "id": "3858949c-a55d-4ff8-a0fc-54fb53e131b5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
+ "converged SCF energy = -108.929838385609\n",
+ "norb = 26\n",
+ "nelec = (5, 5)\n",
+ "E(CCSD) = -109.2177884185543 E_corr = -0.2879500329450042\n",
+ "Using backend ibm_fez\n",
+ "Gate counts: OrderedDict({'sx': 7070, 'rz': 7004, 'cz': 1874, 'x': 52, 'measure': 52, 'barrier': 1})\n",
+ "Fraction of sampled configurations that are valid: 0.00018\n",
+ "Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n",
"Iteration 1\n",
"\tSubsample 0\n",
- "\t\tEnergy: -109.04961576529212\n",
- "\t\tSubspace dimension: 291600\n",
+ "\t\tEnergy: -109.06896153670648\n",
+ "\t\tSubspace dimension: 281961\n",
"\tSubsample 1\n",
- "\t\tEnergy: -109.05419845387246\n",
- "\t\tSubspace dimension: 298116\n",
+ "\t\tEnergy: -109.02861595000536\n",
+ "\t\tSubspace dimension: 263169\n",
"\tSubsample 2\n",
- "\t\tEnergy: -109.03983430633693\n",
- "\t\tSubspace dimension: 291600\n",
+ "\t\tEnergy: -109.03904563251682\n",
+ "\t\tSubspace dimension: 288369\n",
"Iteration 2\n",
"\tSubsample 0\n",
- "\t\tEnergy: -109.13341224824521\n",
- "\t\tSubspace dimension: 432964\n",
+ "\t\tEnergy: -109.13773131837928\n",
+ "\t\tSubspace dimension: 419904\n",
"\tSubsample 1\n",
- "\t\tEnergy: -109.12829311613213\n",
- "\t\tSubspace dimension: 434281\n",
+ "\t\tEnergy: -109.11957805994658\n",
+ "\t\tSubspace dimension: 404496\n",
"\tSubsample 2\n",
- "\t\tEnergy: -109.13563153146234\n",
- "\t\tSubspace dimension: 444889\n",
+ "\t\tEnergy: -109.13928526896287\n",
+ "\t\tSubspace dimension: 405769\n",
"Iteration 3\n",
"\tSubsample 0\n",
- "\t\tEnergy: -109.15805816816365\n",
- "\t\tSubspace dimension: 644809\n",
+ "\t\tEnergy: -109.17341837077211\n",
+ "\t\tSubspace dimension: 603729\n",
"\tSubsample 1\n",
- "\t\tEnergy: -109.17104671394021\n",
- "\t\tSubspace dimension: 588289\n",
+ "\t\tEnergy: -109.16337313138857\n",
+ "\t\tSubspace dimension: 579121\n",
"\tSubsample 2\n",
- "\t\tEnergy: -109.17084043977475\n",
- "\t\tSubspace dimension: 585225\n",
+ "\t\tEnergy: -109.17116470035694\n",
+ "\t\tSubspace dimension: 582169\n",
"Iteration 4\n",
"\tSubsample 0\n",
- "\t\tEnergy: -109.18154982106486\n",
- "\t\tSubspace dimension: 776161\n",
+ "\t\tEnergy: -109.1805148352\n",
+ "\t\tSubspace dimension: 769129\n",
"\tSubsample 1\n",
- "\t\tEnergy: -109.17757205183068\n",
- "\t\tSubspace dimension: 762129\n",
+ "\t\tEnergy: -109.17712191946752\n",
+ "\t\tSubspace dimension: 801025\n",
"\tSubsample 2\n",
- "\t\tEnergy: -109.17890903245448\n",
- "\t\tSubspace dimension: 760384\n",
+ "\t\tEnergy: -109.18252582523913\n",
+ "\t\tSubspace dimension: 776161\n",
"Iteration 5\n",
"\tSubsample 0\n",
- "\t\tEnergy: -109.18924157389883\n",
- "\t\tSubspace dimension: 966289\n",
+ "\t\tEnergy: -109.18759713210412\n",
+ "\t\tSubspace dimension: 976144\n",
"\tSubsample 1\n",
- "\t\tEnergy: -109.18449630323602\n",
- "\t\tSubspace dimension: 952576\n",
+ "\t\tEnergy: -109.18538156651687\n",
+ "\t\tSubspace dimension: 986049\n",
"\tSubsample 2\n",
- "\t\tEnergy: -109.18632692423985\n",
- "\t\tSubspace dimension: 942841\n",
- "Final energy: -109.18924157389883\n",
- "Final energy error: 0.03878764275833646\n"
+ "\t\tEnergy: -109.18622955530022\n",
+ "\t\tSubspace dimension: 1022121\n",
+ "Final energy: -109.18759713210412\n",
+ "Final energy error: 0.040432084553046366\n"
]
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
}
],
"source": [
- "from functools import partial\n",
+ "# ------------------------------ Step 1 ------------------------------\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=\"cc-pvdz\",\n",
+ " symmetry=\"Dooh\",\n",
+ ")\n",
"\n",
- "from qiskit_addon_sqd.fermion import (\n",
- " SCIResult,\n",
- " diagonalize_fermionic_hamiltonian,\n",
- " solve_sci_batch,\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",
+ "norb = len(active_space)\n",
+ "n_electrons = int(sum(scf.mo_occ[active_space]))\n",
+ "n_alpha = (n_electrons + mol.spin) // 2\n",
+ "n_beta = (n_electrons - mol.spin) // 2\n",
+ "nelec = (n_alpha, n_beta)\n",
+ "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\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), norb)\n",
+ "\n",
+ "# Store reference energy from SCI calculation performed separately\n",
+ "reference_energy = -109.22802921665716\n",
+ "\n",
+ "print(f\"norb = {norb}\")\n",
+ "print(f\"nelec = {nelec}\")\n",
+ "\n",
+ "# Get CCSD t2 amplitudes for initializing the ansatz\n",
+ "ccsd = pyscf.cc.CCSD(\n",
+ " scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
+ ").run()\n",
+ "t1 = ccsd.t1\n",
+ "t2 = ccsd.t2\n",
+ "n_reps = 1\n",
+ "alpha_alpha_indices = [(p, p + 1) for p in range(norb - 1)]\n",
+ "alpha_beta_indices = [(p, p) for p in range(0, norb, 4)]\n",
+ "\n",
+ "\n",
+ "ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
+ " t2=t2,\n",
+ " t1=t1,\n",
+ " n_reps=n_reps,\n",
+ " interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),\n",
+ " # Setting optimize=True enables the \"compressed\" factorization\n",
+ " optimize=True,\n",
+ " # Limit the number of optimization iterations to prevent the code cell from running\n",
+ " # too long. Removing this line may improve results.\n",
+ " options=dict(maxiter=1000),\n",
+ ")\n",
+ "\n",
+ "\n",
+ "# create an empty quantum circuit\n",
+ "qubits = QuantumRegister(2 * norb, name=\"q\")\n",
+ "circuit = QuantumCircuit(qubits)\n",
+ "\n",
+ "# prepare Hartree-Fock state as the reference state and append it to the quantum circuit\n",
+ "circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n",
+ "\n",
+ "# apply the UCJ operator to the reference state\n",
+ "circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
+ "circuit.measure_all()\n",
+ "\n",
+ "\n",
+ "# ------------------------------ Step 2 ------------------------------\n",
+ "\n",
+ "service = QiskitRuntimeService()\n",
+ "backend = service.least_busy(\n",
+ " operational=True, simulator=False, min_num_qubits=133\n",
")\n",
"\n",
+ "print(f\"Using backend {backend.name}\")\n",
+ "initial_layout, _ = get_zigzag_physical_layout(norb, backend=backend)\n",
+ "\n",
+ "pass_manager = generate_preset_pass_manager(\n",
+ " optimization_level=3, backend=backend, initial_layout=initial_layout\n",
+ ")\n",
+ "pass_manager.pre_init = ffsim.qiskit.PRE_INIT\n",
+ "isa_circuit = pass_manager.run(circuit)\n",
+ "\n",
+ "print(f\"Gate counts: {isa_circuit.count_ops()}\")\n",
+ "\n",
+ "\n",
+ "# ------------------------------ Step 3 ------------------------------\n",
+ "sampler = Sampler(mode=backend)\n",
+ "job = sampler.run([isa_circuit], shots=100_000)\n",
+ "primitive_result = job.result()\n",
+ "pub_result = primitive_result[0]\n",
+ "\n",
+ "\n",
+ "# ------------------------------ Step 4 ------------------------------\n",
+ "\n",
+ "bit_array = pub_result.data.meas\n",
+ "num_valid = sum(\n",
+ " is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
+ ")\n",
+ "valid_fraction = num_valid / bit_array.num_shots\n",
+ "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")\n",
+ "expected_fraction_random = (\n",
+ " math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
+ ")\n",
+ "print(\n",
+ " f\"Expected fraction of valid configurations from uniformly random bitstrings: {expected_fraction_random}\"\n",
+ ")\n",
"# SQD options\n",
"energy_tol = 1e-3\n",
"occupancies_tol = 1e-3\n",
@@ -866,18 +1126,6 @@
"result_history = []\n",
"\n",
"\n",
- "def callback(results: list[SCIResult]):\n",
- " result_history.append(results)\n",
- " iteration = len(result_history)\n",
- " print(f\"Iteration {iteration}\")\n",
- " for i, result in enumerate(results):\n",
- " print(f\"\\tSubsample {i}\")\n",
- " print(f\"\\t\\tEnergy: {result.energy + nuclear_repulsion_energy}\")\n",
- " print(\n",
- " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
- " )\n",
- "\n",
- "\n",
"result = diagonalize_fermionic_hamiltonian(\n",
" hcore,\n",
" eri,\n",
@@ -894,44 +1142,14 @@
" initial_occupancies=initial_occupancies,\n",
" carryover_threshold=carryover_threshold,\n",
" callback=callback,\n",
- " seed=12345,\n",
+ " seed=rng,\n",
")\n",
"\n",
"final_energy = result.energy + nuclear_repulsion_energy\n",
"energy_error = final_energy - reference_energy\n",
"print(f\"Final energy: {final_energy}\")\n",
- "print(f\"Final energy error: {energy_error}\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9d78906b-4759-4506-9c69-85d4e67766b3",
- "metadata": {},
- "source": [
- "### Visualize the results\n",
- "\n",
- "The first plot shows that after a couple of iterations we estimate the ground state energy within ``~40 mH`` (chemical accuracy is typically accepted to be ``1 kcal/mol`` $\\approx$ ``1.6 mH``). The energy can be improved by allowing more iterations of configuration recovery or increasing the number of samples per batch.\n",
+ "print(f\"Final energy error: {energy_error}\")\n",
"\n",
- "The second plot shows the average occupancy of each spatial orbital after the final iteration. We can see that both the spin-up and spin-down electrons occupy the first five orbitals with high probability in our solutions."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 13,
- "id": "caffd888-e89c-4aa9-8bae-4d1bb723b35e",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
"# Data for energies plot\n",
"x1 = range(len(result_history))\n",
"min_e = [\n",
@@ -983,14 +1201,24 @@
},
{
"cell_type": "markdown",
- "id": "4c3c04c1",
+ "id": "405e89ea-57da-4021-bb18-91e8d583d310",
"metadata": {},
"source": [
- "## Tutorial survey\n",
+ "## Next steps\n",
"\n",
- "Please take this short survey to provide feedback on this tutorial. Your insights will help us improve our content offerings and user experience.\n",
- "\n",
- "[Link to survey](https://your.feedback.ibm.com/jfe/form/SV_8cECHVBBWBjt72S)"
+ "If you found this work interesting, you might be interested in the following material:\n",
+ "- [Sample-based Krylov quantum diagonalization of a fermionic lattice model](/docs/tutorials/sample-based-krylov-quantum-diagonalization) - a related tutorial using time evolution circuits instead of a variational ansatz\n",
+ "- [Scale SQD chemistry workflows with Dice solver](https://qiskit.github.io/qiskit-addon-sqd/how_tos/integrate_dice_solver.html) - a page showing how to use the more efficient Dice software for diagonalization\n",
+ "- [SQD addon API documentation](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) - reference for the `diagonalize_fermionic_hamiltonian` function\n",
+ "- [*Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer*](https://www.science.org/doi/10.1126/sciadv.adu9991) - the paper this tutorial is based on"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a1b8767d",
+ "metadata": {},
+ "source": [
+ "© IBM Corp., 2017-2026"
]
}
],
diff --git a/public/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/baad4e53-5bfd-4cb4-8027-2ec26d27ecdd.avif b/public/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/baad4e53-5bfd-4cb4-8027-2ec26d27ecdd.avif
deleted file mode 100644
index 0a33f70b47d..00000000000
Binary files a/public/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/baad4e53-5bfd-4cb4-8027-2ec26d27ecdd.avif and /dev/null differ
diff --git a/public/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/3858949c-a55d-4ff8-a0fc-54fb53e131b5-1.avif b/public/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/3858949c-a55d-4ff8-a0fc-54fb53e131b5-1.avif
new file mode 100644
index 00000000000..b83e3c68cca
Binary files /dev/null and b/public/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/3858949c-a55d-4ff8-a0fc-54fb53e131b5-1.avif differ
diff --git a/public/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/caffd888-e89c-4aa9-8bae-4d1bb723b35e-0.avif b/public/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/caffd888-e89c-4aa9-8bae-4d1bb723b35e-0.avif
index f9059adfcbb..7302971aa03 100644
Binary files a/public/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/caffd888-e89c-4aa9-8bae-4d1bb723b35e-0.avif and b/public/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/caffd888-e89c-4aa9-8bae-4d1bb723b35e-0.avif differ