Sample-basierte Quantendiagonalisierung von'm Chemie-Hamiltonian

Geschätzter Rechenaufwand: unter einer Minute auf'm Heron-r2-Prozessor (HINWEIS: Das is nur a Schätzung. Deine Laufzeit kann abweichen.)

Hintergrund

In dem Tutorial zeigen wir, wie ma verrauschte Quantensamples nachbearbeitet, um a Annäherung an'n Grundzustand vom Stickstoffmolekül $\text{N}_2$ bei Gleichgewichtsbindungslänge zu finden — mit'n sample-basierten Quantendiagonalisierungs-(SQD-)Algorithmus, für Samples aus'm 59-Qubit-Quantenschaltkreis (52 System-Qubits + 7 Ancilla-Qubits). A Qiskit-basierte Implementierung gibt's bei de SQD-Qiskit-Addons; mehr Details findst du in de entsprechenden Docs mit'n einfachen Beispiel zum Einstieg.

SQD is a Technik zum Finden von Eigenwerten und Eigenvektoren von Quantenoperatoren, wie zum Beispiel'n Quantensystem-Hamiltonian, indem Quantencomputing und verteiltes klassisches Computing gemeinsam genutzt werden. Verteiltes klassisches Computing verarbeitet Samples vom Quantenprozessor und projiziert und diagonalisiert einen Ziel-Hamiltonian in'n durch sie aufgespannten Unterraum. Dadurch is SQD robust gegenüber durch Quantenrauschen korrumpierten Samples und kann große Hamiltonians behandeln — etwa Chemie-Hamiltonians mit Millionen von Wechselwirkungstermen, die für exakte Diagonalisierungsmethoden unerreichbar san. A SQD-basierter Workflow hat die folgenden Schritte:

1. Wähl an Schaltkreis-Ansatz und wend'n auf'm Quantencomputer auf'n Referenzzustand an (in dem Fall den Hartree-Fock-Zustand).
2. Sampel Bitstrings aus'm resultierenden Quantenzustand.
3. Führ die selbstkonsistente Konfigurations-Recovery-Prozedur auf de Bitstrings aus, um die Annäherung an'n Grundzustand zu erhalten.

SQD funktioniert bekanntlich gut, wenn der Ziel-Eigenzustand spärlich is: Die Wellenfunktion is in'ner Menge von Basiszuständen $\mathcal{S} = \{|x\rangle \}$ konzentriert, deren Größe net exponentiell mit der Problemgröße wächst.

Quantenchemie

Die Eigenschaften von Molekülen werden größtenteils durch die Elektronenstruktur bestimmt. Als fermionische Teilchen lassen sich Elektronen mit'm mathematischen Formalismus der zweiten Quantisierung beschreiben. Die Idee is, dass es a Reihe von Orbitalen gibt, die entweder leer oder von einem Fermion besetzt sein können. Ein System mit $N$ Orbitalen wird durch a Menge fermionischer Vernichtungsoperatoren $\{\hat{a}_p\}_{p=1}^N$ beschrieben, die die fermionischen Antikommutationsrelationen erfüllen:

$$\begin{align*} \hat{a}_p \hat{a}_q + \hat{a}_q \hat{a}_p &= 0, \\ \hat{a}_p \hat{a}_q^\dagger + \hat{a}_q^\dagger \hat{a}_p &= \delta_{pq}. \end{align*}$$

Das Adjungierte $\hat{a}_p^\dagger$ heißt Erzeugungsoperator.

Bis jetzt haben wir den Spin außer Acht gelassen, der a fundamentale Eigenschaft von Fermionen is. Berücksichtigt man den Spin, kommen die Orbitale in Paaren, die Raumorbitale genannt werden. Jedes Raumorbital setzt sich aus zwei Spinorbitalen zusammen — eines mit der Bezeichnung "Spin-$\alpha$" und eines mit "Spin-$\beta$". Wir schreiben dann $\hat{a}_{p\sigma}$ für den Vernichtungsoperator, der dem Spinorbital mit Spin $\sigma$ ($\sigma \in \{\alpha, \beta\}$) im Raumorbital $p$ zugeordnet is. Nehmen wir $N$ als die Anzahl der Raumorbitale, dann gibt es insgesamt $2N$ Spinorbitale. Der Hilbert-Raum dieses Systems wird von $2^{2N}$ orthonormalen Basisvektoren aufgespannt, die mit zweiteiligen 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$ bezeichnet werden.

Den Hamiltonian eines Molekülsystems kann ma schreiben als

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},$$

wobei die $h_{pr}$ und $h_{prqs}$ komplexe Zahlen san, die Molekülintegrale genannt werden und aus der Spezifikation des Moleküls mit'm Computerprogramm berechnet werden können. In dem Tutorial berechnen wir die Integrale mit'm Softwarepaket PySCF.

Für Details darüber, wie der molekulare Hamiltonian hergeleitet wird, schau in a Lehrbuch über Quantenchemie (zum Beispiel Modern Quantum Chemistry von Szabo und Ostlund). Für a Erklärung auf hohem Niveau, wie Quantenchemieprobleme auf Quantencomputer abgebildet werden, schau dir die Vorlesung Mapping Problems to Qubits von der Qiskit Global Summer School 2024 an.

Lokaler unitärer Cluster-Jastrow-(LUCJ-)Ansatz

SQD braucht an Schaltkreis-Ansatz, aus dem Samples gezogen werden. In dem Tutorial verwenden wir den lokalen unitären Cluster-Jastrow-(LUCJ-)Ansatz wegen seiner Kombination aus physikalischer Motivation und Hardware-Freundlichkeit.

Der LUCJ-Ansatz is a spezialisierte Form des allgemeinen unitären Cluster-Jastrow-(UCJ-)Ansatzes, der die Form hat:

$$\lvert \Psi \rangle = \prod_{\mu=1}^L e^{\hat{K}_\mu} e^{i \hat{J}_\mu} e^{-\hat{K}_\mu} | \Phi_0 \rangle$$

wobei $\lvert \Phi_0 \rangle$ an Referenzzustand is, oft der Hartree-Fock-Zustand, und die $\hat{K}_\mu$ und $\hat{J}_\mu$ die Form haben:

$$\hat{K}_\mu = \sum_{pq, \sigma} K_{pq}^\mu \, \hat{a}^\dagger_{p \sigma} \hat{a}^{\phantom{\dagger}}_{q \sigma} \;,\; \hat{J}_\mu = \sum_{pq, \sigma\tau} J_{pq,\sigma\tau}^\mu \, \hat{n}_{p \sigma} \hat{n}_{q \tau} \;,$$

wobei wir den Zahlenoperator definiert haben:

$$\hat{n}_{p \sigma} = \hat{a}^\dagger_{p \sigma} \hat{a}^{\phantom{\dagger}}_{p \sigma}.$$

Der Operator $e^{\hat{K}_\mu}$ is a Orbitalrotation, die mit bekannten Algorithmen in linearer Tiefe und mit linearer Konnektivität implementiert werden kann. Die Implementierung des $e^{i \mathcal{J}_k}$-Terms vom UCJ-Ansatz erfordert entweder All-to-all-Konnektivität oder ein fermionisches Swap-Netzwerk, was's für verrauschte, vorfehlertolerierende Quantenprozessoren mit eingeschränkter Konnektivität schwierig macht. Die Idee vom lokalen UCJ-Ansatz is, Sparsitätsbeschränkungen auf die $\mathbf{J}^{\alpha\alpha}$- und $\mathbf{J}^{\alpha\beta}$-Matrizen aufzuerlegen, die es ermöglichen, sie in konstanter Tiefe auf Qubit-Topologien mit eingeschränkter Konnektivität zu implementieren. Die Einschränkungen werden durch a Liste von Indizes angegeben, die angeben, welche Matrixeinträge im oberen Dreieck ungleich null sein dürfen (da die Matrizen symmetrisch san, muss nur das obere Dreieck angegeben werden). Diese Indizes können als Paare von Orbitalen interpretiert werden, die miteinander wechselwirken dürfen.

Als Beispiel: Stell dir a quadratische Qubit-Topologie vor. Ma kann die $\alpha$- und $\beta$-Orbitale in parallelen Linien auf'm Gitter anordnen, wobei Verbindungen zwischen diesen Linien "Sprossen" einer Leiterform bilden, so:

Mit dieser Anordnung san Orbitale mit gleichem Spin mit einer Linientopologie verbunden, während Orbitale mit unterschiedlichem Spin verbunden san, wenn sie dasselbe Raumorbital teilen. Das ergibt die folgenden Indexbeschränkungen für die $\mathbf{J}$-Matrizen:

$$\begin{align*} \mathbf{J}^{\alpha\alpha} &: \set{(p, p+1) \; , \; p = 0, \ldots, N-2} \\ \mathbf{J}^{\alpha\beta} &: \set{(p, p) \;, \; p = 0, \ldots, N-1} \end{align*}$$

Anders gesagt: Wenn die $\mathbf{J}$-Matrizen nur an den angegebenen Indizes im oberen Dreieck ungleich null san, dann kann der $e^{i \mathcal{J}_k}$-Term auf a quadratischer Topologie ohne Swap-Gates in konstanter Tiefe implementiert werden. Natürlich macht das Auferlegen solcher Einschränkungen den Ansatz weniger ausdrucksstark, sodass mehr Ansatz-Wiederholungen nötig sein können.

Die IBM®-Hardware hat a Heavy-Hex-Gitter-Qubit-Topologie. In dem Fall kann ma a "Zickzack"-Muster verwenden, wie unten dargestellt. In dem Muster werden Orbitale mit gleichem Spin auf Qubits mit einer Linientopologie abgebildet (rote und blaue Kreise), und eine Verbindung zwischen Orbitalen mit unterschiedlichem Spin gibt's bei jedem 4. Raumorbital — vermittelt durch an Ancilla-Qubit (lila Kreise). Die Indexbeschränkungen san in dem Fall:

$$\begin{align*} \mathbf{J}^{\alpha\alpha} &: \set{(p, p+1) \; , \; p = 0, \ldots, N-2} \\ \mathbf{J}^{\alpha\beta} &: \set{(p, p) \;, \; p = 0, 4, 8, \ldots (p \leq N-1)} \end{align*}$$

Selbstkonsistente Konfigurations-Recovery

Die selbstkonsistente Konfigurations-Recovery-Prozedur is dazu gedacht, so viel Signal wie möglich aus verrauschten Quantensamples herauszuholen. Da der molekulare Hamiltonian die Teilchenzahl und den Spin-Z erhält, macht's Sinn, an Schaltkreis-Ansatz zu wählen, der diese Symmetrien ebenfalls erhält. Wenn ma'n auf'n Hartree-Fock-Zustand anwendet, hat der resultierende Zustand im rauschfreien Fall eine feste Teilchenzahl und einen festen Spin-Z. Deshalb sollten die Spin-$\alpha$- und Spin-$\beta$-Hälften jedes aus diesem Zustand gesampelten Bitstrings dasselbe Hamming-Gewicht wie im Hartree-Fock-Zustand haben. Wegen Rauschen in aktuellen Quantenprozessoren verletzen manche gemessenen Bitstrings diese Eigenschaft. Eine einfache Form der Nachselektion würde diese Bitstrings verwerfen, aber das is verschwenderisch, weil sie trotzdem noch a Signal enthalten können. Die selbstkonsistente Recovery-Prozedur versucht, a Teil dieses Signals bei der Nachbearbeitung wiederherzustellen. Die Prozedur is iterativ und braucht als Eingabe a Schätzung der durchschnittlichen Besetzung jedes Orbitals im Grundzustand, die zunächst aus den Roh-Samples berechnet wird. Die Prozedur läuft in aner Schleife, und jede Iteration hat die folgenden Schritte:

1. Für jeden Bitstring, der die angegebenen Symmetrien verletzt: Flippe seine Bits mit aner probabilistischen Prozedur, die darauf ausgelegt is, den Bitstring näher an die aktuelle Schätzung der durchschnittlichen Orbitalbesetzungen zu bringen, um an neuen Bitstring zu erhalten.
2. Sammle alle alten und neuen Bitstrings, die die Symmetrien erfüllen, und subsample Teilmengen fester Größe, die im Voraus gewählt werden.
3. Für jede Teilmenge von Bitstrings: Projiziere den Hamiltonian in den Unterraum, der von den entsprechenden Basisvektoren aufgespannt wird (schau in'n vorigen Abschnitt für a Beschreibung dieser Basisvektoren), und berechne auf'm klassischen Computer a Grundzustandsschätzung des projizierten Hamiltonians.
4. Aktualisiere die Schätzung der durchschnittlichen Orbitalbesetzungen mit der Grundzustandsschätzung mit der niedrigsten Energie.

SQD-Workflow-Diagramm

Der SQD-Workflow is im folgenden Diagramm dargestellt:

Voraussetzungen

Bevor du mit dem Tutorial anfängst, stell sicher, dass du Folgendes installiert hast:

• Qiskit SDK v1.0 oder neuer, mit Visualisierungsunterstützung
• Qiskit Runtime v0.22 oder neuer (pip install qiskit-ibm-runtime)
• SQD-Qiskit-Addon v0.11 oder neuer (pip install qiskit-addon-sqd)
• ffsim v0.0.58 oder neuer (pip install ffsim)

Setup

# Added by doQumentation — installs packages not in the Binder environment
!pip install -q ffsim qiskit-addon-sqd

import pyscf
import pyscf.cc
import pyscf.mcscf
import ffsim
import numpy as np
import matplotlib.pyplot as plt

from qiskit import QuantumCircuit, QuantumRegister
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Schritt 1: Klassische Eingaben auf's Quantenproblem abbilden

In dem Tutorial finden wir a Annäherung an'n Grundzustand des Moleküls im Gleichgewicht im cc-pVDZ-Basissatz. Zuerst legen wir das Molekül und seine Eigenschaften fest.

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="cc-pvdz",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Store reference energy from SCI calculation performed separately
exact_energy = -109.22690201485733

converged SCF energy = -108.929838385609

Bevor wir den LUCJ-Ansatz-Schaltkreis aufbauen, führen wir zunächst in der folgenden Code-Zelle a CCSD-Berechnung durch. Die $t_1$- und $t_2$-Amplituden aus dieser Berechnung werden verwendet, um die Parameter des Ansatzes zu initialisieren.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -109.2177884185543  E_corr = -0.2879500329450045

Jetzt verwenden wir ffsim, um den Ansatz-Schaltkreis zu erstellen. Da unser Molekül an geschlossenschaligen Hartree-Fock-Zustand hat, verwenden wir die spinausgeglichene Variante des UCJ-Ansatzes, UCJOpSpinBalanced. Wir übergeben Wechselwirkungspaare, die für a Heavy-Hex-Gitter-Qubit-Topologie geeignet san (schau in'n Hintergrundabschnitt zum LUCJ-Ansatz). Wir setzen optimize=True in der from_t_amplitudes-Methode, um die "komprimierte" Doppelfaktorisierung der $t_2$-Amplituden zu aktivieren (schau in The local unitary cluster Jastrow (LUCJ) ansatz in ffsims Dokumentation für Details).

n_reps = 1
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
    # Limit the number of optimization iterations to prevent the code cell from running
    # too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

nelec = (num_elec_a, num_elec_b)

# create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

Schritt 2: Problem für die Quantenhardware-Ausführung optimieren

Als nächstes optimieren wir den Schaltkreis für a Zielhardware.

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)

print(f"Using backend {backend.name}")

Using backend ibm_fez

Wir empfehlen die folgenden Schritte, um den Ansatz zu optimieren und hardwarekompatibel zu machen.

• Wähl physikalische Qubits (initial_layout) von der Zielhardware aus, die dem zuvor beschriebenen "Zickzack"-Muster entsprechen. Das Anordnen von Qubits in diesem Muster führt zu'nem effizienten, hardwarekompatiblen Schaltkreis mit weniger Gates. Hier gibt's a bisschen Code, der die Auswahl des "Zickzack"-Musters automatisiert und a Scoring-Heuristik verwendet, um die mit dem gewählten Layout verbundenen Fehler zu minimieren.
• Generiere an gestaffelten Pass-Manager mit der Funktion generate_preset_pass_manager aus Qiskit mit deiner Wahl von backend und initial_layout.
• Setz die pre_init-Stage deines gestaffelten Pass-Managers auf ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT enthält Qiskit-Transpiler-Passes, die Gates in Orbitalrotationen zerlegen und diese dann zusammenführen, was zu weniger Gates im endgültigen Schaltkreis führt.
• Führ den Pass-Manager auf deinem Schaltkreis aus.

Code zur automatischen Auswahl des "Zickzack"-Layouts

from typing import Sequence

import rustworkx
from qiskit.providers import BackendV2
from rustworkx import NoEdgeBetweenNodes, PyGraph

IBM_TWO_Q_GATES = {"cx", "ecr", "cz"}

def create_linear_chains(num_orbitals: int) -> PyGraph:
    """In zig-zag layout, there are two linear chains (with connecting qubits between
    the chains). This function creates those two linear chains: a rustworkx PyGraph
    with two disconnected linear chains. Each chain contains `num_orbitals` number
    of nodes, that is, in the final graph there are `2 * num_orbitals` number of nodes.

    Args:
        num_orbitals (int): Number orbitals or nodes in each linear chain. They are
            also known as alpha-alpha interaction qubits.

    Returns:
        A rustworkx.PyGraph with two disconnected linear chains each with `num_orbitals`
            number of nodes.
    """
    G = rustworkx.PyGraph()

    for n in range(num_orbitals):
        G.add_node(n)

    for n in range(num_orbitals - 1):
        G.add_edge(n, n + 1, None)

    for n in range(num_orbitals, 2 * num_orbitals):
        G.add_node(n)

    for n in range(num_orbitals, 2 * num_orbitals - 1):
        G.add_edge(n, n + 1, None)

    return G

def create_lucj_zigzag_layout(
    num_orbitals: int, backend_coupling_graph: PyGraph
) -> tuple[PyGraph, int]:
    """This function creates the complete zigzag graph that 'can be mapped' to an IBM QPU with
    heavy-hex connectivity (the zigzag must be an isomorphic sub-graph to the QPU/backend
    coupling graph for it to be mapped).
    The zigzag pattern includes both linear chains (alpha-alpha interactions) and connecting
    qubits between the linear chains (alpha-beta interactions).

    Args:
        num_orbitals (int): Number of orbitals, that is, number of nodes in each alpha-alpha linear chain.
        backend_coupling_graph (PyGraph): The coupling graph of the backend on which the LUCJ ansatz
            will be mapped and run. This function takes the coupling graph as a undirected
            `rustworkx.PyGraph` where there is only one 'undirected' edge between two nodes,
            that is, qubits. Usually, the coupling graph of a IBM backend is directed (for example, Eagle devices
            such as ibm_brisbane) or may have two edges between two nodes (for example, Heron `ibm_torino`).
            A user needs to be make such graphs undirected and/or remove duplicate edges to make them
            compatible with this function.

    Returns:
        G_new (PyGraph): The graph with IBM backend compliant zigzag pattern.
        num_alpha_beta_qubits (int): Number of connecting qubits between the linear chains
            in the zigzag pattern. While we want as many connecting (alpha-beta) qubits between
            the linear (alpha-alpha) chains, we cannot accommodate all due to qubit and connectivity
            constraints of backends. This is the maximum number of connecting qubits the zigzag pattern
            can have while being backend compliant (that is, isomorphic to backend coupling graph).
    """
    isomorphic = False
    G = create_linear_chains(num_orbitals=num_orbitals)

    num_iters = num_orbitals
    while not isomorphic:
        G_new = G.copy()
        num_alpha_beta_qubits = 0
        for n in range(num_iters):
            if n % 4 == 0:
                new_node = 2 * num_orbitals + num_alpha_beta_qubits
                G_new.add_node(new_node)
                G_new.add_edge(n, new_node, None)
                G_new.add_edge(new_node, n + num_orbitals, None)
                num_alpha_beta_qubits = num_alpha_beta_qubits + 1
        isomorphic = rustworkx.is_subgraph_isomorphic(
            backend_coupling_graph, G_new
        )
        num_iters -= 1

    return G_new, num_alpha_beta_qubits

def lightweight_layout_error_scoring(
    backend: BackendV2,
    virtual_edges: Sequence[Sequence[int]],
    physical_layouts: Sequence[int],
    two_q_gate_name: str,
) -> list[list[list[int], float]]:
    """Lightweight and heuristic function to score isomorphic layouts. There can be many zigzag patterns,
    each with different set of physical qubits, that can be mapped to a backend. Some of them may
    include less noise qubits and couplings than others. This function computes a simple error score
    for each such layout. It sums up 2Q gate error for all couplings in the zigzag pattern (layout) and
    measurement of errors of physical qubits in the layout to compute the error score.

    Note:
        This lightweight scoring can be refined using concepts such as mapomatic.

    Args:
        backend (BackendV2): A backend.
        virtual_edges (Sequence[Sequence[int]]): Edges in the device compliant zigzag pattern where
            nodes are numbered from 0 to (2 * num_orbitals + num_alpha_beta_qubits).
        physical_layouts (Sequence[int]): All physical layouts of the zigzag pattern that are isomorphic
            to each other and to the larger backend coupling map.
        two_q_gate_name (str): The name of the two-qubit gate of the backend. The name is used for fetching
            two-qubit gate error from backend properties.

    Returns:
        scores (list): A list of lists where each sublist contains two items. First item is the layout, and
            second item is a float representing error score of the layout. The layouts in the `scores` are
            sorted in the ascending order of error score.
    """
    props = backend.properties()
    scores = []
    for layout in physical_layouts:
        total_2q_error = 0
        for edge in virtual_edges:
            physical_edge = (layout[edge[0]], layout[edge[1]])
            try:
                ge = props.gate_error(two_q_gate_name, physical_edge)
            except Exception:
                ge = props.gate_error(two_q_gate_name, physical_edge[::-1])
            total_2q_error += ge
        total_measurement_error = 0
        for qubit in layout:
            meas_error = props.readout_error(qubit)
            total_measurement_error += meas_error
        scores.append([layout, total_2q_error + total_measurement_error])

    return sorted(scores, key=lambda x: x[1])

def _make_backend_cmap_pygraph(backend: BackendV2) -> PyGraph:
    graph = backend.coupling_map.graph
    if not graph.is_symmetric():
        graph.make_symmetric()
    backend_coupling_graph = graph.to_undirected()

    edge_list = backend_coupling_graph.edge_list()
    removed_edge = []
    for edge in edge_list:
        if set(edge) in removed_edge:
            continue
        try:
            backend_coupling_graph.remove_edge(edge[0], edge[1])
            removed_edge.append(set(edge))
        except NoEdgeBetweenNodes:
            pass

    return backend_coupling_graph

def get_zigzag_physical_layout(
    num_orbitals: int, backend: BackendV2, score_layouts: bool = True
) -> tuple[list[int], int]:
    """The main function that generates the zigzag pattern with physical qubits that can be used
    as an `intial_layout` in a preset passmanager/transpiler.

    Args:
        num_orbitals (int): Number of orbitals.
        backend (BackendV2): A backend.
        score_layouts (bool): Optional. If `True`, it uses the `lightweight_layout_error_scoring`
            function to score the isomorphic layouts and returns the layout with less erroneous qubits.
            If `False`, returns the first isomorphic subgraph.

    Returns:
        A tuple of device compliant layout (list[int]) with zigzag pattern and an int representing
            number of alpha-beta-interactions.
    """
    backend_coupling_graph = _make_backend_cmap_pygraph(backend=backend)

    G, num_alpha_beta_qubits = create_lucj_zigzag_layout(
        num_orbitals=num_orbitals,
        backend_coupling_graph=backend_coupling_graph,
    )

    isomorphic_mappings = rustworkx.vf2_mapping(
        backend_coupling_graph, G, subgraph=True
    )
    isomorphic_mappings = list(isomorphic_mappings)

    edges = list(G.edge_list())

    layouts = []
    for mapping in isomorphic_mappings:
        initial_layout = [None] * (2 * num_orbitals + num_alpha_beta_qubits)
        for key, value in mapping.items():
            initial_layout[value] = key
        layouts.append(initial_layout)

    two_q_gate_name = IBM_TWO_Q_GATES.intersection(
        backend.configuration().basis_gates
    ).pop()

    if score_layouts:
        scores = lightweight_layout_error_scoring(
            backend=backend,
            virtual_edges=edges,
            physical_layouts=layouts,
            two_q_gate_name=two_q_gate_name,
        )

        return scores[0][0][:-num_alpha_beta_qubits], num_alpha_beta_qubits

    return layouts[0][:-num_alpha_beta_qubits], num_alpha_beta_qubits

initial_layout, _ = get_zigzag_physical_layout(num_orbitals, backend=backend)

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
 isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
 isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")

Gate counts (w/o pre-init passes): OrderedDict({'rz': 12438, 'sx': 12169, 'cz': 3986, 'x': 573, 'measure': 52, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'sx': 7059, 'rz': 6962, 'cz': 1876, 'measure': 52, 'x': 35, 'barrier': 1})

Schritt 3: Mit Qiskit-Primitiven ausführen

Nachdem wir den Schaltkreis für die Hardwareausführung optimiert haben, können wir'n auf der Zielhardware ausführen und Samples für die Grundzustandsenergieschätzung sammeln. Da wir nur an Schaltkreis haben, verwenden wir Qiskit Runtimes Job-Ausführungsmodus und führen unseren Schaltkreis aus.

sampler = Sampler(mode=backend)
job = sampler.run([isa_circuit], shots=100_000)

primitive_result = job.result()
pub_result = primitive_result[0]

Schritt 4: Nachbearbeiten und Ergebnis im gewünschten klassischen Format zurückgeben

Jetzt schätzen wir die Grundzustandsenergie des Hamiltonians mit der Funktion diagonalize_fermionic_hamiltonian. Diese Funktion führt die selbstkonsistente Konfigurations-Recovery-Prozedur aus, um die verrauschten Quantensamples iterativ zu verfeinern und die Energieschätzung zu verbessern. Wir übergeben a Callback-Funktion, damit wir die Zwischenergebnisse für die spätere Analyse speichern können. Schau in die API-Dokumentation für Erklärungen der Argumente von diagonalize_fermionic_hamiltonian.

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    pub_result.data.meas,
    samples_per_batch=samples_per_batch,
    norb=num_orbitals,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=12345,
)

Iteration 1
	Subsample 0
		Energy: -108.97781410104506
		Subspace dimension: 28561
	Subsample 1
		Energy: -108.97781410104506
		Subspace dimension: 28561
	Subsample 2
		Energy: -108.97781410104506
		Subspace dimension: 28561
Iteration 2
	Subsample 0
		Energy: -109.05150860754739
		Subspace dimension: 287296
	Subsample 1
		Energy: -109.08152283263908
		Subspace dimension: 264196
	Subsample 2
		Energy: -109.11636893067873
		Subspace dimension: 284089
Iteration 3
	Subsample 0
		Energy: -109.15878555367885
		Subspace dimension: 429025
	Subsample 1
		Energy: -109.16462831275786
		Subspace dimension: 473344
	Subsample 2
		Energy: -109.15895026995382
		Subspace dimension: 435600
Iteration 4
	Subsample 0
		Energy: -109.17784051223317
		Subspace dimension: 622521
	Subsample 1
		Energy: -109.1774651326829
		Subspace dimension: 657721
	Subsample 2
		Energy: -109.18085212360191
		Subspace dimension: 617796
Iteration 5
	Subsample 0
		Energy: -109.18636242518915
		Subspace dimension: 815409
	Subsample 1
		Energy: -109.18451014767518
		Subspace dimension: 837225
	Subsample 2
		Energy: -109.18333728638888
		Subspace dimension: 857476

Ergebnisse visualisieren

Der erste Plot zeigt, dass wir nach a paar Iterationen die Grundzustandsenergie auf etwa ~41 mH genau schätzen (chemische Genauigkeit wird üblicherweise mit 1 kcal/mol $\approx$ 1,6 mH angegeben). Die Energie kann durch mehr Konfigurationsrecovery-Iterationen oder durch Erhöhen der Samples pro Batch verbessert werden.

Der zweite Plot zeigt die durchschnittliche Besetzung jedes Raumorbitals nach der letzten Iteration. Wir sehen, dass sowohl die Spin-up- als auch die Spin-down-Elektronen in unseren Lösungen mit hoher Wahrscheinlichkeit die ersten fünf Orbitale besetzen.

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

Tutorial-Umfrage

Bitte nimm an dieser kurzen Umfrage teil, um Feedback zu dem Tutorial zu geben. Deine Einblicke helfen uns, unsere Inhalte und Nutzererfahrung zu verbessern.

Link zur Umfrage