1001 lines
48 KiB
Plaintext
1001 lines
48 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "090b6884",
|
|
"metadata": {},
|
|
"source": [
|
|
"{/* cspell:ignore hcore ccsd Motta */}\n",
|
|
"\n",
|
|
"# SQD for energy estimation of a chemistry Hamiltonian"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ee7de1bc",
|
|
"metadata": {},
|
|
"source": [
|
|
"In this lesson, we will apply SQD to estimate the ground state energy of a molecule.\n",
|
|
"\n",
|
|
"In particular, we will discuss the following topics using the $4$-step Qiskit pattern approach:\n",
|
|
"\n",
|
|
"1. Step 1: Map problem to quantum circuits and operators\n",
|
|
" - Setup the molecular Hamiltonian for $N_2$.\n",
|
|
" - Explain the chemistry-inspired and hardware-friendly local unitary cluster Jastrow (LUCJ) [\\[1\\]](#references)\n",
|
|
"2. Step 2: Optimize for target hardware\n",
|
|
" - Optimize gate counts and layout of the ansatz for hardware execution\n",
|
|
"3. Step 3: Execute on target hardware\n",
|
|
" - Run the optimized circuit on a real QPU to generate samples of the subspace.\n",
|
|
"4. Step 4: Post-process results\n",
|
|
" - Introduce the self-consistent configuration recovery loop [\\[2\\]](#references)\n",
|
|
" - Post-process the full set of bitstring samples, using prior knowledge of particle number and the average orbital occupancy calculated on the most recent iteration.\n",
|
|
" - Probabilistically create batches of subsamples from recovered bitstrings.\n",
|
|
" - Project and diagonalize the molecular Hamiltonian over each sampled subspace.\n",
|
|
" - Save the minimum ground state energy found across all batches and update the avg orbital occupancy.\n",
|
|
"\n",
|
|
"We will use several software packages throughout the lesson.\n",
|
|
"- `PySCF` to define the molecule and setup the Hamiltonian.\n",
|
|
"- `ffsim` package to construct the LUCJ ansatz.\n",
|
|
"- `Qiskit` for transpiling the ansatz for hardware execution.\n",
|
|
"- `Qiskit IBM Runtime` to execute the circuit on a QPU and collect samples.\n",
|
|
"- `Qiskit addon SQD` configuration recovery and ground state energy estimation using subspace projection and matrix diagonalization."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "719a9c0e-8c00-4fab-ba23-f2c8a7ebb573",
|
|
"metadata": {},
|
|
"source": [
|
|
"## 1. Map problem to quantum circuits and operators"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d02f97af",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Molecular Hamiltonian"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5dac27ce-56b1-4f7f-ab83-ac153c004e82",
|
|
"metadata": {},
|
|
"source": [
|
|
"A molecular Hamiltonian takes the generic form:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n",
|
|
"+\n",
|
|
"\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n",
|
|
"\\frac{(pr|qs)}{2} \\,\n",
|
|
"\\hat{a}^\\dagger_{p\\sigma}\n",
|
|
"\\hat{a}^\\dagger_{q\\tau}\n",
|
|
"\\hat{a}_{s\\tau}\n",
|
|
"\\hat{a}_{r\\sigma}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"$\\hat{a}^\\dagger_{p\\sigma}$/$\\hat{a}_{p\\sigma}$ are the fermionic creation/annihilation operators associated to the $p$-th basis set element and the spin $\\sigma$. $h_{pr}$ and $(pr|qs)$ are the one- and two-body electronic integrals. Using pySCF, we will define the molecule and compute the one- and two-body integrals of the Hamiltonian for basis set `6-31g`."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "cc987c08-7261-4c4c-a06b-609d7003efe9",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"converged SCF energy = -108.835236570774\n",
|
|
"CASCI E = -109.046671778080 E(CI) = -32.8155692383188 S^2 = 0.0000000\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"import warnings\n",
|
|
"import pyscf\n",
|
|
"import pyscf.cc\n",
|
|
"import pyscf.mcscf\n",
|
|
"\n",
|
|
"warnings.filterwarnings(\"ignore\")\n",
|
|
"\n",
|
|
"# Specify molecule properties\n",
|
|
"open_shell = False\n",
|
|
"spin_sq = 0\n",
|
|
"\n",
|
|
"# Build N2 molecule\n",
|
|
"mol = pyscf.gto.Mole()\n",
|
|
"mol.build(\n",
|
|
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]], # Two N atoms 1 angstrom apart\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) # hcore: one-body integrals\n",
|
|
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals) # eri: two-body integrals\n",
|
|
"\n",
|
|
"# Compute exact energy for comparison\n",
|
|
"exact_energy = cas.run().e_tot"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "3aa4e0f0",
|
|
"metadata": {},
|
|
"source": [
|
|
"In this lesson, we will use Jordan-Wigner (JW) transformation to map a fermionic wavefunction to a qubit wavefunction so that it can be prepared using a quantum circuit. The JW transformation maps the Fock space of fermions in M spatial orbitals onto the Hilbert space of 2M qubits, i.e., a spatial orbital is split into two _spin orbitals_, one associated with a spin up ($\\alpha$) electron and another with spin down ($\\beta$). A spin orbital can be occupied or unoccupied. Usually, when we refer to number of orbitals, we will be using number of _spatial_ orbitals. The number of spin orbitals will be double. In quantum circuits, we will represent each spin orbital with one qubit. Thus, a set of qubits will represent spin-up or $\\alpha$-orbitals, and another set will represent spin-up or $\\beta$-orbitals. For example, $N_2$ molecule for `6-31g` basis set has $16$ spatial orbitals (i.e., $16$ $\\alpha$ + $16$ $\\beta$ = $32$ spin orbitals). Thus, we will need a $32$-qubit quantum circuit (we may need extra ancilla qubits as discussed later). The qubits are measured in computational basis to generate bitstrings, which represent electronic configurations or (Slater) determinants. Throughout this lesson, we will use the terms bitstrings, configurations, and determinants interchangeably. The bitstrings tell us electron occupancy in spin orbitals: a $1$ in a bit position means the corresponding spin orbital is occupied, while a $0$ means the spin orbital is empty. As electronic structure problems are particle preserving, only a fixed number of spin orbitals must be occupied. The $N_2$ molecule has $5$ spin-up ($\\alpha$) and $5$ spin-down ($\\beta$) electrons. Thus, any bitstring representing the $\\alpha$ and $\\beta$ orbitals must have five $1\\text{s}$ each for $N_2$ molecule."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "f6a89c89-85e1-4c0d-abeb-8ab8b154a7ba",
|
|
"metadata": {},
|
|
"source": [
|
|
"### 1.1 Quantum circuit for sample generation: The LUCJ ansatz"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "a464c865-1528-45c2-8ede-325583b15976",
|
|
"metadata": {},
|
|
"source": [
|
|
"In this lesson, we will use the local unitary coupled cluster Jastrow (LUCJ) [\\\\[1\\\\]](#references) ansatz for quantum state preparation and subsequent sampling. First, we will explain different building blocks of the full UCJ ansatz and the approximations made in the local version of it. Next, by using ffsim package, we will construct the LUCJ ansatz and optimize it using Qiskit transpiler for hardware execution.\n",
|
|
"\n",
|
|
"The UCJ ansatz has the following form (for a product of $L$ layers or repetitions of the UCJ operator.)\n",
|
|
"\n",
|
|
"$$\n",
|
|
"|\\psi\\rangle = \\prod_{\\mu=1}^{L}{(e^{K^{\\mu}} \\times {e^{iJ^{\\mu}}} \\times {e^{-K^{\\mu}}})} |\\Phi_{0}\\rangle\n",
|
|
"$$\n",
|
|
"\n",
|
|
"where, $\\vert \\Phi_{0} \\rangle$ is a reference state, typically taken as the Hartree-Fock (HF) state. As the Hartree-Fock state is defined as having the lowest numbered orbitals occupied, the HF state preparation will involve applying X gates to set qubits corresponding to occupied orbitals to one. For example, the HF state preparation block for 4 spatial orbitals and 2 up- and 2 down-spin may look like the following:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "0f0bb614",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "4a32d88b",
|
|
"metadata": {},
|
|
"source": [
|
|
"A single repetition of the UCJ operator ${(e^{K^{(\\mu)}} \\times {e^{iJ^{(\\mu)}}} \\times {e^{-K^{(\\mu)}}})}$ consists of a diagonal Coulomb evolution ($e^{iJ^{(\\mu)}}$) sandwiched by orbital rotations ($e^{K^{(\\mu)}}$ and $e^{-K^{(\\mu)}}$)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "963e8386-39b6-40b7-9740-fffcf1573fe6",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "271b1825-bc78-4957-bbf8-21a714e45ed3",
|
|
"metadata": {},
|
|
"source": [
|
|
"Orbital rotation blocks work on a single spin species ($\\alpha$ (up-spin)/$\\beta$ (down-spin)). For each electron species, orbital rotation consists of a layer of single-qubit $R_{z}$ gates followed by a sequence of 2-qubit Given's rotation gates ($XX + YY$ gates).\n",
|
|
"\n",
|
|
"The 2-qubit gates act on adjacent spin-orbitals (nearest neighbor qubits), and therefore, are implementable on IBM QPUs without the need for SWAP gates."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "2e8ac1d2-8f04-4591-921b-8ba0174e4ad0",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ebc9b48c-e77d-4af5-b8cf-460daeeadcdb",
|
|
"metadata": {},
|
|
"source": [
|
|
"The $e^{iJ^{(\\mu)}}$, also known as the diagonal Coulomb operator, consists of three blocks. Two of them work on same spin sectors ($e^{iJ_{\\alpha \\alpha}^{(\\mu)}}$ and $e^{iJ_{\\beta \\beta}^{(\\mu)}}$), and one works between two spin sectors ($e^{iJ_{\\alpha \\beta}^{(\\mu)}}$).\n",
|
|
"\n",
|
|
"All the blocks in $e^{iJ^{(\\mu)}}$ consists of number-number gates $U_{nn}(\\phi)$ [\\[1\\]](#references). A $U_{nn}(\\phi)$ gate can be further broken down into a $R_{ZZ}(\\frac{\\phi}{2})$ gate followed by two single-qubit $Rz(-\\frac{\\phi}{2})$ gates acting on two separate qubits.\n",
|
|
"\n",
|
|
"Same-spin components ($J_{\\alpha \\alpha}$ and $J_{\\beta \\beta}$) have $U_{nn}$ gates between all possible pairs of qubits. However, as superconducting QPUs have restrictive connectivity, qubits must be swapped to realize gates between non-adjacent qubits.\n",
|
|
"\n",
|
|
"For example, consider the following $e^{iJ_{\\alpha \\alpha}^{(\\mu)}}$ (or $e^{iJ_{\\beta \\beta}^{(\\mu)}}$) block for $N = 4$ spatial orbitals. For a linear qubit connectivity, the last three gates are not directly implementable as they work between non-adjacent qubits (e.g., Q0 and Q2 are not directly connected). Therefore, we need SWAP gates to make them adjacent (following figure shows an example with $3$ SWAP gates)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d3e24a20-1c86-4aea-8300-13bd2ead4e8f",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ef3f1d36-e8ab-46f0-bddb-fa8fb884e531",
|
|
"metadata": {},
|
|
"source": [
|
|
"Next, the $J_{\\alpha \\beta}$ implements gates between same indexed orbitals from different spin sectors (e.g., between $0\\alpha$ and $0\\beta$). Similarly, if the qubits are not physically adjacent on a QPU, these gates will also require SWAPs."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9afe0036-318c-43b9-9b2c-39a808bda82c",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "6219aa8e-8a63-4bfc-b52b-64507292471d",
|
|
"metadata": {},
|
|
"source": [
|
|
"From the above discussion, the UCJ ansatz faces some hurdles for HW execution as it needs SWAP gates due to non-adjacent qubit interactions. The local variant of the UCJ ansatz, LUCJ, addresses this challenge by removing some $U_{nn}$ from the diagonal Coulomb operator.\n",
|
|
"\n",
|
|
"In the same electron species blocks, $J_{\\alpha \\alpha}$ and $J_{\\beta \\beta}$), we only keep the $U_{nn}$ gates compatible with nearest-neighbor connectivity and remove gates between non-adjacent qubits in the LUCJ version. Following figure shows the LUCJ block after removal of non-adjacent gates."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "24f54a55-4a09-4508-997a-96306835c7e6",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "165d0a19-0366-4af8-8ea0-20df351f49f5",
|
|
"metadata": {},
|
|
"source": [
|
|
"Next, the LUCJ version of the $J_{\\alpha \\beta}$ block that works between different electron species can take different shape based on the device topology.\n",
|
|
"\n",
|
|
"Here, also, the LUCJ version gets rid of non-compatible gates. The figure below shows variants of the $J_{\\alpha \\beta}$ block for different qubit topology including grid, hexagonal, heavy-hex, and linear.\n",
|
|
"\n",
|
|
"- **Grid**: we can have $U_{nn}$ gates between all $\\alpha$ and $\\beta$ orbitals without any SWAPs, and therefore, do not need to remove any $U_{nn}$ gates.\n",
|
|
"- **Hexagonal**: Every other orbital (0th, 2nd, 4th, etc. indexed orbitals) becomes nearest neighbors when $\\alpha$ and $\\beta$ are laid out in two adjacent linear chains.\n",
|
|
"- **Linear**: Only one $\\alpha$ and one $\\beta$ orbital are connected, which means the $J_{\\alpha \\beta}$ block will have only one gate.\n",
|
|
"- **Heavy-hex**: The $\\alpha$-$\\beta$ interactions are kept between every $4$-th indexed (0th, 4th, 8th, etc.) spin orbitals and are need _ancilla_ mediated, i.e., we need ancilla qubits between the linear chains representing $\\alpha$ and $\\beta$ orbitals. This arrangement needs a limited number of SWAPs."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "8ec1e433-bc00-42bc-bdbb-288c56c32f9d",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "3e070615-2b4a-4866-8c6c-d2d038447f45",
|
|
"metadata": {},
|
|
"source": [
|
|
"While removing gates from the UCJ ansatz to construct the LUCJ version makes it more HW compatible, the ansatz loses some expressivity. Therefore, more repetitions ($L$) of the modified UCJ operator may be needed when using the LUCJ ansatz."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "04367dac",
|
|
"metadata": {},
|
|
"source": [
|
|
"### 1.2 LUCJ ansatz initialization"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "f5eae55f",
|
|
"metadata": {},
|
|
"source": [
|
|
"The LUCJ is a parameterized ansatz, and we need to initialize the parameters before hardware execution. One way to initialize ansatz is by using `t1` and `t2` amplitudes from classical coupled cluster singles and doubles (CCSD) method, where `t1` amplitudes are the coefficient of single excitation operators and `t2` amplitudes are for double excitation operators.\n",
|
|
"\n",
|
|
"Note that while initializing the LUCJ ansatz with `t1` and `t2` amplitudes generate decent results, the ansatz parameters may need further optimization."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "1835e74e-3354-425f-8596-c574f03e7a6e",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"E(CCSD) = -109.0398256929733 E_corr = -0.20458912219883\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# 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",
|
|
")\n",
|
|
"ccsd.run()\n",
|
|
"\n",
|
|
"t1 = ccsd.t1\n",
|
|
"t2 = ccsd.t2"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9f0842fa",
|
|
"metadata": {},
|
|
"source": [
|
|
"### 1.3 Constructing the LUCJ ansatz using `ffsim`"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e5f7f7e8-30e3-40e8-9b85-bf057b35b766",
|
|
"metadata": {},
|
|
"source": [
|
|
"We will use the [ffsim](https://github.com/qiskit-community/ffsim/tree/main) package to create and initialize the ansatz with `t1` and `t2` amplitudes computed above. Since our molecule has a closed-shell Hartree-Fock state, we will use the spin-balanced variant of the UCJ ansatz, [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced).\n",
|
|
"\n",
|
|
"As IBM hardware has a heavy-hex topology, we will adopt the _zig-zag_ pattern used in [\\[1\\]](#references) and explained above for qubit interactions. In this pattern, orbitals (qubits) with the same spin are connected with a line topology (red and blue circles). Due to the heavy-hex topology, orbitals for different spins have connections between every 4th orbital (0th, 4th, 8th, etc.) (purple circles).\n",
|
|
"\n",
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"id": "6799421a-4425-404a-9994-47215957054d",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import ffsim\n",
|
|
"from qiskit import QuantumCircuit, QuantumRegister\n",
|
|
"\n",
|
|
"n_reps = 2\n",
|
|
"alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]\n",
|
|
"alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]\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",
|
|
")\n",
|
|
"\n",
|
|
"nelec = (num_elec_a, num_elec_b)\n",
|
|
"\n",
|
|
"# create an empty quantum circuit\n",
|
|
"qubits = QuantumRegister(2 * num_orbitals, 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(num_orbitals, 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",
|
|
"# circuit.decompose().draw(\"mpl\", scale=0.5, fold=-1)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "3b2de316",
|
|
"metadata": {},
|
|
"source": [
|
|
"The LUCJ ansatz with repeated layers can be optimized by merging some adjacent blocks. Consider a case for `n_reps=2`. The two orbital rotation blocks in the middle can be merged into a single orbital rotation block. The `ffsim` package has a pass manager named `ffsim.qiskit.PRE_INIT` to optimize the circuit by merging such adjacent blocks."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7cb99cd9",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e9ad460d",
|
|
"metadata": {},
|
|
"source": [
|
|
"## 2. Optimize for target hardware"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "cecf2994",
|
|
"metadata": {},
|
|
"source": [
|
|
"First, we fetch a backend of our choice. We will optimize our circuit for the backend, and then execute the optimized circuit on the same backend to generate samples for the subspace."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"id": "4c8dbbef-378e-4ae3-9290-bde58b72024d",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
|
|
"\n",
|
|
"service = QiskitRuntimeService()\n",
|
|
"backend = service.backend(\"ibm_kyiv\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "1b8ba388-6904-4725-ad97-2a31c0adaa07",
|
|
"metadata": {},
|
|
"source": [
|
|
"Next, 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 adheres to the zig-zag pattern (two linear chains with ancilla qubit in-between them) described above. Laying out qubits in this pattern leads to an efficient hardware-compatible circuit with less gates.\n",
|
|
"- Generate a staged pass manager using the [generate_preset_pass_manager](https://docs.quantum.ibm.com/api/qiskit/transpiler_preset#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": "code",
|
|
"execution_count": 5,
|
|
"id": "ecd2675c-a302-43fb-80f4-d658d56360d5",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})\n",
|
|
"Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
|
|
"\n",
|
|
"spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]\n",
|
|
"spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]\n",
|
|
"\n",
|
|
"initial_layout = spin_a_layout + spin_b_layout\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": "72a780f9",
|
|
"metadata": {},
|
|
"source": [
|
|
"## 3. Execute on target hardware"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "4cd9164c-697a-4aa6-b71f-86720e3d5b66",
|
|
"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](https://docs.quantum.ibm.com/guides/execution-modes) and execute our circuit."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"id": "9f542448-5767-4e12-85c8-dd8584545dbb",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
|
|
"\n",
|
|
"sampler = Sampler(mode=backend)\n",
|
|
"sampler.options.dynamical_decoupling.enable = True\n",
|
|
"\n",
|
|
"job = sampler.run([isa_circuit], shots=10_000) # Takes approximately 5sec of QPU time"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "8c5eb3e9-2a5a-423b-ac18-7bba7f69f18f",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Run cell after IQX job completion\n",
|
|
"primitive_result = job.result()\n",
|
|
"pub_result = primitive_result[0]\n",
|
|
"counts = pub_result.data.meas.get_counts()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e3ebf77a-8efd-43d6-bd19-e26423921c6f",
|
|
"metadata": {},
|
|
"source": [
|
|
"## 4. Post-process results"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "1e7f0ecd",
|
|
"metadata": {},
|
|
"source": [
|
|
"The post-processing part of the SQD workflow can be summarized using the following diagram.\n",
|
|
"\n",
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "8c51a527",
|
|
"metadata": {},
|
|
"source": [
|
|
"Sampling the LUCJ ansatz in the computational basis generates a pool of noisy configurations $\\tilde{\\mathcal{\\chi}}$, which are used in the post-processing routine. It involves a method called (details discussed later) _configuration recovery_ to probabilistically correct configurations with incorrect electron numbers. Configurations only with correct electron numbers $\\tilde{\\mathcal{\\chi}}_{R}$ are then subsampled and distributed into multiple batches based on the frequency of appearance of each unique configuration. Each batch of samples defines a subspace ($\\mathcal{S^{(k)}}$). Next, the molecular Hamiltonian, $H$, is projected onto subspaces:\n",
|
|
"$$\n",
|
|
"H_{\\mathcal{S}^{(k)}} = P_{\\mathcal{S}^{(k)}} H _{\\mathcal{S}^{(k)}} \\text{ with } P_{\\mathcal{S}^{(k)}} = \\sum_{x \\in \\mathcal{S}^{(k)}} \\vert x \\rangle \\langle x \\vert\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Each projected Hamiltonian $H_{\\mathcal{S}^{(k)}}$ is then fed into an Eigensolver, where it is diagonalized to compute eigenvalues and eigenvectors to reconstruct an eigenstate. In this lesson, we project and diagonalize the Hamiltonian using the `qiskit-addon-sqd` package which uses the Davidson's method from PySCF for diagonalization.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"H_{\\mathcal{S}^{(k)}} \\vert \\psi^{(k)} \\rangle = E^{(k)} \\vert \\psi^{(k)} \\rangle\n",
|
|
"$$\n",
|
|
"\n",
|
|
"We then collect the lowest eigenvalue (energy) from the batches, and also compute average orbital occupancy, $\\text{n}$. The average occupancy information is used in the configuration recovery step to probabilistically correct noise configurations.\n",
|
|
"\n",
|
|
"Next, we explain the self-consistent configuration recovery loop in detail and show concrete code examples to implement the above-mentioned steps to estimate the ground state energy of $N_2$ Hamiltonian."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "39997b5d-8bd3-4bc8-9ade-11fa5cfb34f9",
|
|
"metadata": {},
|
|
"source": [
|
|
"### 4.1 Configuration recovery: overview"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "b05881ac-80ce-47a6-ac28-c1237f85bc0a",
|
|
"metadata": {},
|
|
"source": [
|
|
"Each bit in a bitstring (Slater determinant) represents a spin orbital. The right half of a bitstring represents spin-up orbitals, and the left half represents spin-down orbitals. A `1` means the orbital is occupied by an electron, and a `0` means the orbital is empty. We know the correct number of particles (both up-spin electron and down-spin electron) a priori. Suppose we have a determinant $x$ with $N_x$ electrons (i.e., there are $N_x$ numbers of $1$s in the bitstring) in it. The correct number of particles is $N$. If $N_x \\neq N$, then we know that the bitstring is corrupted by noise. The self-consistent configuration routine attempts to correct the bitstring by probabilistically flipping $|N_x - N|$ bits by leveraging average orbital occupancy information. The average orbital occupancy ($n$) tells us how likely an orbital be occupied by an electron. If $N_x < N$, we have fewer electrons and need to flip some $0$s to $1$s and vice versa.\n",
|
|
"\n",
|
|
"The probability of flipping can be $|x[i] - avg\\_occupancy[i]|$ for `i`-th spin orbital. In [\\[2\\]](#references), the authors used a weighted probability of flipping using modified ReLU function.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{align}\n",
|
|
" w(y) = \\begin{cases}\n",
|
|
"\n",
|
|
" \\delta \\frac{y}{h} & \\text{if } y \\leq h\\\\ \\nonumber\n",
|
|
"\n",
|
|
" \\delta + (1 - \\delta) \\frac{y - h}{1 - h} & \\text{if } y > h\n",
|
|
"\n",
|
|
"\\end{cases}\n",
|
|
"\\end{align}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Here $h$ defines the location of \"corner\" of the ReLU function, and the parameter $\\delta$ defines the value of the ReLU function at the corner. For $\\delta = 0$, $w$ becomes true ReLU function, and for $\\delta >0$, it becomes _modified_ ReLU. In the paper, the authors used $\\delta = 0.01$ and $h =$ number of alpha (or beta) particles/number of alpha (or beta) spin orbitals $= N/M$ (filling factor).\n",
|
|
"\n",
|
|
"The average orbital occupancy ($n$) is not known a priori. The first iteration of the ground state estimation starts with configurations with only correct particle numbers in both spin species. After the first iteration, we have an estimate of the ground state, and using the estimate, we can construct the first guess of $n$. This guess of $n$ is used to recover configurations, run the next iteration of ground state estimation, and self-consistently refine the guess of $n$. The process repeats until the a stopping criterion is met.\n",
|
|
"\n",
|
|
"Consider the following example for $N = 2$ and $x = |1000\\rangle$ ($N_x = 1$). We need to flip one of the 0s to 1 to correct it for particle numbers, and the choices are `1100`, `1010`, and `1001`. Based on the probability of flipping, one of the choices will be selected as _recovered configuration_ (or the bitstring with correct number of particles).\n",
|
|
"\n",
|
|
"Suppose in the first iteration we run two batches, and the estimated ground states from them are:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{align}\\nonumber\n",
|
|
" \\text{Batch0: } \\vert \\psi \\rangle &= 0.8 \\times \\vert 1001 \\rangle + 0.6 \\times \\vert 0101 \\rangle \\\\ \\nonumber\n",
|
|
" \\text{Batch1: } \\vert \\psi \\rangle &= \\frac{1}{\\sqrt{3}} \\left( \\vert 1001 \\rangle + \\vert 0101 \\rangle + \\vert 0110 \\rangle \\right) \\nonumber\n",
|
|
"\\end{align}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Using the computational basis states and their amplitudes, we can compute probability of electron occupancies (in short _occupancies_) per spin-orbital (qubit) (note that probability = |amplitude|$^2$). Below we tabulate qubit-wise occupancies for each bitstring appearing in the estimated ground state and compute total orbital occupancy for a batch. Note that, as per Qiskit ordering convention, the right most bit represents qubit-0 (Q0), and the left most bit represents Q3.\n",
|
|
"\n",
|
|
"Occupancy (Batch0):\n",
|
|
"| | Q3 | Q2 | Q1 | Q0 |\n",
|
|
"|:---------:|:--------:|:--------:|:--------:|:--------:|\n",
|
|
"| 1001 | 0.64 | 0.0 | 0.0 | 0.64 |\n",
|
|
"| 0110 | 0.0 | 0.36 | 0.36 | 0.0 |\n",
|
|
"| **n** _(Batch0)_ | **0.64** | **0.36** | **0.36** | **0.64** |\n",
|
|
"\n",
|
|
"Occupancy (Batch1)\n",
|
|
"| | Q3 | Q2 | Q1 | Q0 |\n",
|
|
"|:---------:|:--------:|:--------:|:--------:|:--------:|\n",
|
|
"| 1001 | 0.33 | 0.00 | 0.00 | 0.33 |\n",
|
|
"| 0101 | 0.0 | 0.33 | 0.00 | 0.33 |\n",
|
|
"| 0110 | 0.0 | 0.33 | 0.33 | 0.00 |\n",
|
|
"| **n** _(Batch1)_ | **0.33** | **0.66** | **0.33** | **0.66** |\n",
|
|
"\n",
|
|
"Occupancy (average across batches)\n",
|
|
"| | Q3 | Q2 | Q1 | Q0 |\n",
|
|
"|:---------:|:--------:|:--------:|:--------:|:--------:|\n",
|
|
"| **n** _(Batch0)_ | 0.64 | 0.36 | 0.36 | 0.64 |\n",
|
|
"| **n** _(Batch1)_ | 0.33 | 0.66 | 0.33 | 0.66 |\n",
|
|
"| **n** _(average)_ | **0.49** | **0.51** | **0.35** | **0.65** |\n",
|
|
"\n",
|
|
"Using the average orbital occupancy computed above, we can find the probabilities of flip for different orbitals in the configuration $x = \\vert 1000 \\rangle$. As the orbital represented by Q3 is already occupied and need not to be flipped, we set its p(flip) to $0$. For the remaining orbitals, which are unoccupied, the probability of flip is $\\vert x[i] - \\text{n}[i] \\vert$ each. Along with p(flip), we also compute the probability weight associated with flipping using the modified ReLU function described above.\n",
|
|
"\n",
|
|
"Probability of flip ($x = \\vert 1000 \\rangle$, $\\delta = 0.01$, $h = N/M = 2/4 = 0.50$)\n",
|
|
"| | Q3 | Q2 | Q1 | Q0 |\n",
|
|
"|:---------:|:--------:|:--------:|:--------:|:--------:|\n",
|
|
"| p(flip) ($\\vert x[i] - \\text{n}[i] \\vert$) | 0 | 0.51 | 0.35 | 0.65 |\n",
|
|
"| w(p(flip)) | 0 | 0.03 | 0.007 | 0.31 |\n",
|
|
"\n",
|
|
"Finally, using weighted probabilities above, we can flip one of the unoccupied Q2, Q1, and Q0 orbitals. Based on the values above, Q0 will be flipped most likely, and a possible recovered configuration can be $\\vert \\text{1001} \\rangle$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "c1940235",
|
|
"metadata": {},
|
|
"source": [
|
|
""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5c614290",
|
|
"metadata": {},
|
|
"source": [
|
|
"The complete self-consistent configuration recovery process can be summarized as follows:\n",
|
|
"\n",
|
|
"**First iteration:** Suppose the bitstrings (configurations or Slater determinants) generated by the quantum computer form a set $\\widetilde{\\chi}$, which includes both configurations with correct ($\\widetilde{\\chi}_{correct}$) and incorrect ($\\widetilde{\\chi}_{incorrect}$) number of particles in each spin sector.\n",
|
|
"1. Configurations from ($\\widetilde{\\chi}_{correct}$) are randomly sampled to create batches $(\\mathcal{S}^{(1)}, \\cdots, \\mathcal{S}^{(K)})$ of vectors for subspace projection. The number of batches and samples in each batch are user defined parameters. The larger the number of samples in each batch, the larger the subspace dimension and more computationally demanding the diagonalization becomes. On the other hand, too small number of samples may miss the ground state support vectors and lead to incorrect estimation.\n",
|
|
"2. Run the eigenstate solver (i.e. projection onto subspace and diagonalization) on the batches and obtain approximate eigenstates. $|\\psi^{(1)}\\rangle, \\cdots, |\\psi^{(K)}\\rangle$.\n",
|
|
"3. From the approximate eigenstates construct the first guess for $n$.\n",
|
|
"\n",
|
|
"**Subsequent iterations:**\n",
|
|
"1. Using $n$ correct the configurations with wrong particle number in $\\widetilde{\\chi}_{incorrect}$. Suppose we name them $\\widetilde{\\chi}_{correct\\_new}$. Then, $\\widetilde{\\chi}_{recovered} (\\widetilde{\\chi}_{R}) = \\widetilde{\\chi}_{correct} \\cup \\widetilde{\\chi}_{correct\\_new}$ forms the new set of configurations with correct particle numbers.\n",
|
|
"2. $\\widetilde{\\chi}_{R}$ is sampled to create batches $\\mathcal{S}^{(1)}, \\cdots, \\mathcal{S}^{(K)}$.\n",
|
|
"3. Eigenstate solver runs with new batches and generates new estimates of ground states $|\\psi^{(1)}\\rangle, \\cdots, |\\psi^{(K)}\\rangle$.\n",
|
|
"4. From the approximate eigenstates construct refined guess for $n$.\n",
|
|
"5. If the stopping criterion is not met, go back to step `2.1`."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5eea548a",
|
|
"metadata": {},
|
|
"source": [
|
|
"### 4.2 Ground state estimation"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "31dfe00e",
|
|
"metadata": {},
|
|
"source": [
|
|
"First, we will transform the counts into a bitstring matrix and probability array for post-processing.\n",
|
|
"\n",
|
|
"Each row in the matrix represents one unique bitstring. Since qubits are indexed from the right of a bitstring in Qiskit, column ``0`` represents qubit ``N-1``, and column ``N-1`` represents qubit ``0``, where ``N`` is the number of qubits.\n",
|
|
"\n",
|
|
"The alpha orbitals are represented in the column index range ``(N, N/2]`` (right half), and the beta orbitals are represented in the column range ``(N/2, 0]`` (left half)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"id": "71550274",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from qiskit_addon_sqd.counts import counts_to_arrays\n",
|
|
"\n",
|
|
"# Convert counts into bitstring and probability arrays\n",
|
|
"bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "acbdbdc5",
|
|
"metadata": {},
|
|
"source": [
|
|
"There are a few user-controlled options which are important for this technique:\n",
|
|
"\n",
|
|
"- ``iterations``: Number of self-consistent configuration recovery iterations\n",
|
|
"- ``n_batches``: Number of batches of configurations used by the different calls to the eigenstate solver\n",
|
|
"- ``samples_per_batch``: Number of unique configurations to include in each batch\n",
|
|
"- ``max_davidson_cycles``: Maximum number of Davidson cycles run by each eigensolver"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"id": "a2d22e0b-8f51-42a0-858f-ad0297cb0bae",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting configuration recovery iteration 0\n",
|
|
" Batch 0 subspace dimension: 21609\n",
|
|
" Batch 1 subspace dimension: 21609\n",
|
|
" Batch 2 subspace dimension: 21609\n",
|
|
" Batch 3 subspace dimension: 21609\n",
|
|
" Batch 4 subspace dimension: 21609\n",
|
|
"Starting configuration recovery iteration 1\n",
|
|
" Batch 0 subspace dimension: 609961\n",
|
|
" Batch 1 subspace dimension: 616225\n",
|
|
" Batch 2 subspace dimension: 627264\n",
|
|
" Batch 3 subspace dimension: 633616\n",
|
|
" Batch 4 subspace dimension: 624100\n",
|
|
"Starting configuration recovery iteration 2\n",
|
|
" Batch 0 subspace dimension: 564001\n",
|
|
" Batch 1 subspace dimension: 605284\n",
|
|
" Batch 2 subspace dimension: 582169\n",
|
|
" Batch 3 subspace dimension: 559504\n",
|
|
" Batch 4 subspace dimension: 591361\n",
|
|
"Starting configuration recovery iteration 3\n",
|
|
" Batch 0 subspace dimension: 550564\n",
|
|
" Batch 1 subspace dimension: 549081\n",
|
|
" Batch 2 subspace dimension: 531441\n",
|
|
" Batch 3 subspace dimension: 527076\n",
|
|
" Batch 4 subspace dimension: 531441\n",
|
|
"Starting configuration recovery iteration 4\n",
|
|
" Batch 0 subspace dimension: 544644\n",
|
|
" Batch 1 subspace dimension: 580644\n",
|
|
" Batch 2 subspace dimension: 527076\n",
|
|
" Batch 3 subspace dimension: 531441\n",
|
|
" Batch 4 subspace dimension: 537289\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"from qiskit_addon_sqd.configuration_recovery import recover_configurations\n",
|
|
"from qiskit_addon_sqd.fermion import (\n",
|
|
" bitstring_matrix_to_ci_strs,\n",
|
|
" solve_fermion,\n",
|
|
")\n",
|
|
"from qiskit_addon_sqd.subsampling import postselect_and_subsample\n",
|
|
"\n",
|
|
"rng = np.random.default_rng(24)\n",
|
|
"# SQD options\n",
|
|
"iterations = 5\n",
|
|
"\n",
|
|
"# Eigenstate solver options\n",
|
|
"n_batches = 5\n",
|
|
"samples_per_batch = 500\n",
|
|
"max_davidson_cycles = 300\n",
|
|
"\n",
|
|
"# Self-consistent configuration recovery loop\n",
|
|
"e_hist = np.zeros((iterations, n_batches)) # energy history\n",
|
|
"s_hist = np.zeros((iterations, n_batches)) # spin history\n",
|
|
"occupancy_hist = []\n",
|
|
"avg_occupancy = None\n",
|
|
"for i in range(iterations):\n",
|
|
" print(f\"Starting configuration recovery iteration {i}\")\n",
|
|
" # On the first iteration, we have no orbital occupancy information from the\n",
|
|
" # solver, so we begin with the full set of noisy configurations.\n",
|
|
" if avg_occupancy is None:\n",
|
|
" bs_mat_tmp = bitstring_matrix_full\n",
|
|
" probs_arr_tmp = probs_arr_full\n",
|
|
"\n",
|
|
" # If we have average orbital occupancy information, we use it to refine the full set of noisy configurations\n",
|
|
" else:\n",
|
|
" bs_mat_tmp, probs_arr_tmp = recover_configurations(\n",
|
|
" bitstring_matrix_full,\n",
|
|
" probs_arr_full,\n",
|
|
" avg_occupancy,\n",
|
|
" num_elec_a,\n",
|
|
" num_elec_b,\n",
|
|
" rand_seed=rng,\n",
|
|
" )\n",
|
|
"\n",
|
|
" # Create batches of subsamples. We post-select here to remove configurations\n",
|
|
" # with incorrect hamming weight during iteration 0, since no config recovery was performed.\n",
|
|
" batches = postselect_and_subsample(\n",
|
|
" bs_mat_tmp,\n",
|
|
" probs_arr_tmp,\n",
|
|
" hamming_right=num_elec_a,\n",
|
|
" hamming_left=num_elec_b,\n",
|
|
" samples_per_batch=samples_per_batch,\n",
|
|
" num_batches=n_batches,\n",
|
|
" rand_seed=rng,\n",
|
|
" )\n",
|
|
"\n",
|
|
" # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.\n",
|
|
" e_tmp = np.zeros(n_batches)\n",
|
|
" s_tmp = np.zeros(n_batches)\n",
|
|
" occs_tmp = []\n",
|
|
" coeffs = []\n",
|
|
" for j in range(n_batches):\n",
|
|
" strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])\n",
|
|
" print(f\" Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}\")\n",
|
|
" energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(\n",
|
|
" batches[j],\n",
|
|
" hcore,\n",
|
|
" eri,\n",
|
|
" open_shell=open_shell,\n",
|
|
" spin_sq=spin_sq,\n",
|
|
" max_davidson=max_davidson_cycles,\n",
|
|
" )\n",
|
|
" energy_sci += nuclear_repulsion_energy\n",
|
|
" e_tmp[j] = energy_sci\n",
|
|
" s_tmp[j] = spin\n",
|
|
" occs_tmp.append(avg_occs)\n",
|
|
" coeffs.append(coeffs_sci)\n",
|
|
"\n",
|
|
" # Combine batch results\n",
|
|
" avg_occupancy = tuple(np.mean(occs_tmp, axis=0))\n",
|
|
"\n",
|
|
" # Track optimization history\n",
|
|
" e_hist[i, :] = e_tmp\n",
|
|
" s_hist[i, :] = s_tmp\n",
|
|
" occupancy_hist.append(avg_occupancy)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e4f6ec2c-032d-4e18-8ea4-cf867cba5054",
|
|
"metadata": {},
|
|
"source": [
|
|
"### 4.3 Discussion of results"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "592eabb7-18f3-4622-8710-bfc43ad6cdec",
|
|
"metadata": {},
|
|
"source": [
|
|
"The first plot shows that after a few iterations we estimate the ground state energy within ~24 mH (chemical accuracy is typically accepted to be 1 kcal/mol $\\approx$ 1.6 mH). 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.\n",
|
|
"\n",
|
|
"Although the estimated ground state energy is decent, it is not within the chemical accuracy limit ($\\pm \\approx 1.6$ mH). This gap can be attributed to the small subspace dimension we used above for projection and diagonalization. As we used `samples_per_batch=500`, the subspace is spanned by max $500$ vectors, which is missing vectors from ground state support. Increasing the `samples_per_batch` parameter should improve the accuracy at the expense of more classical compute resources and runtime."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"id": "61959e69-a182-4636-abcb-a32349fc9076",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Data for energies plot\n",
|
|
"x1 = range(iterations)\n",
|
|
"min_e = [np.min(e) for e in e_hist]\n",
|
|
"e_diff = [abs(e - exact_energy) for e in min_e]\n",
|
|
"yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]\n",
|
|
"\n",
|
|
"# Chemical accuracy (+/- 1 milli-Hartree)\n",
|
|
"chem_accuracy = 0.001\n",
|
|
"\n",
|
|
"# Data for avg spatial orbital occupancy\n",
|
|
"y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]\n",
|
|
"x2 = range(len(y2))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"id": "8cd90034-6ef3-41bd-a847-c115cade82f7",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Exact energy: -109.04667 Ha\n",
|
|
"SQD energy: -109.02234 Ha\n",
|
|
"Absolute error: 0.02434 Ha\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<Image src=\"/learning/images/courses/quantum-diagonalization-algorithms/qda-4-sqd-implementation/extracted-outputs/8cd90034-6ef3-41bd-a847-c115cade82f7-1.avif\" alt=\"Output of the previous code cell\" />"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\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-6)\n",
|
|
"axs[0].axhline(\n",
|
|
" y=chem_accuracy, color=\"#BF5700\", linestyle=\"--\", 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",
|
|
"print(f\"Exact energy: {exact_energy:.5f} Ha\")\n",
|
|
"print(f\"SQD energy: {min_e[-1]:.5f} Ha\")\n",
|
|
"print(f\"Absolute error: {e_diff[-1]:.5f} Ha\")\n",
|
|
"plt.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e1530472",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Exercise for the reader\n",
|
|
"Progressively increase the `samples_per_batch` parameter (e.g., from $1000$ to $10000$ at a step of $1000$; permitted my your computer's memory) and compare the estimated ground state energies."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d2b2241f",
|
|
"metadata": {},
|
|
"source": [
|
|
"## References\n",
|
|
"\n",
|
|
"\\[1] M. Motta et al., “Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure” (2023). [Chem. Sci., 2023, 14, 11213](https://pubs.rsc.org/en/content/articlehtml/2023/sc/d3sc02516k).\n",
|
|
"\n",
|
|
"\\[2] J. Robledo-Moreno et al., \"Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer\" (2024). [arXiv:quant-ph/2405.05068](https://arxiv.org/abs/2405.05068)."
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"description": "Sample-based quantum diagonalization (SQD) is implemented in the context of solving the ground state of a nitrogen molecule.",
|
|
"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"
|
|
},
|
|
"title": "SQD Implementation"
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|