Stichprobenbasierte Quantendiagonalisierung eines Chemie-Hamiltonians

Nutzungsschätzung: unter einer Minute auf einem Heron r2 Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann variieren.)

Hintergrund

In diesem Tutorial zeigen wir, wie verrauschte Quantenstichproben nachbearbeitet werden, um eine Approximation des Grundzustands des Stickstoffmoleküls $\text{N}_2$ bei Gleichgewichtsbindungslänge zu finden. Dabei verwenden wir den stichprobenbasierten Quantendiagonalisierungsalgorithmus (SQD) für Stichproben aus einem 59-Qubit-Quantenschaltkreis (52 System-Qubits + 7 Ancilla-Qubits). Eine Qiskit-basierte Implementierung ist verfügbar in den SQD Qiskit Addons, weitere Details findest du in der entsprechenden Dokumentation mit einem einfachen Beispiel für den Einstieg. SQD ist eine Technik zum Finden von Eigenwerten und Eigenvektoren von Quantenoperatoren, wie dem Hamiltonian eines Quantensystems, unter Verwendung von Quanten- und verteiltem klassischem Computing gemeinsam. Das verteilte klassische Computing wird verwendet, um Stichproben von einem Quantenprozessor zu verarbeiten und einen Ziel-Hamiltonian in einen durch sie aufgespannten Unterraum zu projizieren und zu diagonalisieren. Dies ermöglicht SQD, robust gegenüber durch Quantenrauschen verfälschten Stichproben zu sein und mit großen Hamiltonians umzugehen, wie Chemie-Hamiltonians mit Millionen von Wechselwirkungstermen, die jenseits der Reichweite exakter Diagonalisierungsmethoden liegen. Ein SQD-basierter Workflow hat folgende Schritte:

1. Wähle einen Schaltkreis-Ansatz und wende ihn auf einem Quantencomputer auf einen Referenzzustand an (in diesem Fall den Hartree-Fock-Zustand).
2. Sample Bitstrings aus dem resultierenden Quantenzustand.
3. Führe das selbstkonsistente Konfigurationswiederherstellungsverfahren auf den Bitstrings aus, um die Approximation des Grundzustands zu erhalten.

SQD funktioniert bekanntermaßen gut, wenn der Ziel-Eigenzustand dünn besetzt ist: Die Wellenfunktion wird von einer Menge von Basiszuständen $\mathcal{S} = \{|x\rangle \}$ getragen, deren Größe nicht exponentiell mit der Problemgröße wächst.

Quantenchemie

Die Eigenschaften von Molekülen werden weitgehend durch die Struktur der Elektronen in ihnen bestimmt. Als fermionische Teilchen können Elektronen mit einem mathematischen Formalismus namens zweite Quantisierung beschrieben werden. Die Idee ist, dass es eine Anzahl von Orbitalen gibt, von denen jedes entweder leer ist oder von einem Fermion besetzt wird. Ein System von $N$ Orbitalen wird durch einen Satz fermionischer Vernichtungsoperatoren $\{\hat{a}_p\}_{p=1}^N$ beschrieben, die die fermionischen Antikommutatorrelationen 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*}$$

Der adjungierte Operator $\hat{a}_p^\dagger$ wird Erzeugungsoperator genannt.

Bisher hat unsere Darstellung den Spin nicht berücksichtigt, der eine fundamentale Eigenschaft von Fermionen ist. Bei Berücksichtigung des Spins kommen die Orbitale in Paaren vor, die Raumorbitale genannt werden. Jedes Raumorbital besteht aus zwei Spinorbitalen, von denen eines mit "Spin-$\alpha$" und eines mit "Spin-$\beta$" bezeichnet wird. Wir schreiben dann $\hat{a}_{p\sigma}$ für den Vernichtungsoperator, der mit dem Spinorbital mit Spin $\sigma$ ($\sigma \in \{\alpha, \beta\}$) im Raumorbital $p$ assoziiert ist. Wenn wir $N$ als Anzahl der Raumorbitale nehmen, gibt es insgesamt $2N$ Spinorbitale. Der Hilbert-Raum dieses Systems wird von $2^{2N}$ orthonormalen Basisvektoren aufgespannt, die mit zweiteiligen Bitstrings bezeichnet werden: $\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$.

Der Hamiltonian eines molekularen Systems kann geschrieben werden 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 sind, die molekulare Integrale genannt werden und aus der Spezifikation des Moleküls mit einem Computerprogramm berechnet werden können. In diesem Tutorial berechnen wir die Integrale mit dem Softwarepaket PySCF.

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

Local Unitary Cluster Jastrow (LUCJ) Ansatz

SQD benötigt einen Quantenschaltkreis-Ansatz, aus dem Stichproben gezogen werden können. In diesem Tutorial verwenden wir den Local Unitary Cluster Jastrow (LUCJ) Ansatz aufgrund seiner Kombination aus physikalischer Motivation und Hardware-Freundlichkeit.

Der LUCJ-Ansatz ist eine spezialisierte Form des allgemeinen Unitary 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$ ein Referenzzustand ist, 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 Teilchenzahloperator definiert haben

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

Der Operator $e^{\hat{K}_\mu}$ ist eine Orbitalrotation, die mit bekannten Algorithmen in linearer Tiefe und unter Verwendung linearer Konnektivität implementiert werden kann. Die Implementierung des $e^{i \mathcal{J}_k}$ Terms des UCJ-Ansatzes erfordert entweder vollständige Konnektivität oder die Verwendung eines fermionischen Swap-Netzwerks, was für verrauschte präfehlertolerante Quantenprozessoren mit begrenzter Konnektivität herausfordernd ist. Die Idee des lokalen UCJ-Ansatzes besteht darin, Dünnbesetztheitsbedingungen 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 begrenzter Konnektivität zu implementieren. Die Bedingungen werden durch eine Liste von Indizes spezifiziert, die angeben, welche Matrixeinträge im oberen Dreieck von null verschieden sein dürfen (da die Matrizen symmetrisch sind, muss nur das obere Dreieck spezifiziert werden). Diese Indizes können als Orbitalpaare interpretiert werden, die interagieren dürfen.

Betrachte als Beispiel eine quadratische Gitter-Qubit-Topologie. Wir können die $\alpha$- und $\beta$-Orbitale in parallelen Linien auf dem Gitter platzieren, wobei Verbindungen zwischen diesen Linien "Sprossen" einer Leiterform bilden, wie folgt:

Bei diesem Setup sind Orbitale mit demselben Spin mit einer Linientopologie verbunden, während Orbitale mit unterschiedlichen Spins verbunden sind, wenn sie dasselbe Raumorbital teilen. Dies ergibt die folgenden Indexbedingungen 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*}$$

Mit anderen Worten: Wenn die $\mathbf{J}$-Matrizen nur an den angegebenen Indizes im oberen Dreieck von null verschieden sind, kann der $e^{i \mathcal{J}_k}$ Term auf einer quadratischen Topologie ohne Verwendung von Swap-Gates in konstanter Tiefe implementiert werden. Natürlich macht die Auferlegung solcher Bedingungen auf den Ansatz ihn weniger ausdrucksstark, sodass möglicherweise mehr Ansatz-Wiederholungen erforderlich sind.

Die IBM® Hardware hat eine Heavy-Hex-Gitter-Qubit-Topologie, in welchem Fall wir ein "Zickzack"-Muster verwenden können, das unten dargestellt ist. In diesem Muster werden Orbitale mit demselben Spin auf Qubits mit einer Linientopologie abgebildet (rote und blaue Kreise), und eine Verbindung zwischen Orbitalen unterschiedlichen Spins ist an jedem 4. Raumorbital vorhanden, wobei die Verbindung durch ein Ancilla-Qubit (violette Kreise) ermöglicht wird. In diesem Fall sind die Indexbedingungen

$$\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 Konfigurationswiederherstellung

Das selbstkonsistente Konfigurationswiederherstellungsverfahren ist darauf ausgelegt, so viel Signal wie möglich aus verrauschten Quantenstichproben zu extrahieren. Da der molekulare Hamiltonian die Teilchenzahl und Spin-Z erhält, ist es sinnvoll, einen Schaltkreis-Ansatz zu wählen, der diese Symmetrien ebenfalls erhält. Wenn er auf den Hartree-Fock-Zustand angewendet wird, hat der resultierende Zustand im rauschfreien Fall eine feste Teilchenzahl und Spin-Z. Daher sollten die Spin-$\alpha$- und Spin-$\beta$-Hälften jedes aus diesem Zustand gesampelten Bitstrings dasselbe Hamming-Gewicht wie im Hartree-Fock-Zustand haben. Aufgrund des Vorhandenseins von Rauschen in aktuellen Quantenprozessoren werden einige gemessene Bitstrings diese Eigenschaft verletzen. Eine einfache Form der Postselektion würde diese Bitstrings verwerfen, aber das ist verschwenderisch, da diese Bitstrings möglicherweise noch etwas Signal enthalten. Das selbstkonsistente Wiederherstellungsverfahren versucht, einen Teil dieses Signals in der Nachbearbeitung wiederherzustellen. Das Verfahren ist iterativ und benötigt als Eingabe eine Schätzung der durchschnittlichen Besetzungen jedes Orbitals im Grundzustand, die zuerst aus den Rohstichproben berechnet wird. Das Verfahren wird in einer Schleife ausgeführt, und jede Iteration hat die folgenden Schritte:

1. Für jeden Bitstring, der die spezifizierten Symmetrien verletzt, werden seine Bits mit einem probabilistischen Verfahren geflippt, das darauf ausgelegt ist, den Bitstring näher an die aktuelle Schätzung der durchschnittlichen Orbitalbesetzungen zu bringen, um einen neuen Bitstring zu erhalten.
2. Sammle alle alten und neuen Bitstrings, die die Symmetrien erfüllen, und entnehme Teilmengen einer im Voraus gewählten festen Größe.
3. Für jede Teilmenge von Bitstrings wird der Hamiltonian in den durch die entsprechenden Basisvektoren aufgespannten Unterraum projiziert (siehe den vorherigen Abschnitt für eine Beschreibung dieser Basisvektoren) und eine Grundzustandsschätzung des projizierten Hamiltonians auf einem klassischen Computer berechnet.
4. Aktualisiere die Schätzung der durchschnittlichen Orbitalbesetzungen mit der Grundzustandsschätzung mit der niedrigsten Energie.

SQD-Workflow-Diagramm

Der SQD-Workflow ist im folgenden Diagramm dargestellt:

Anforderungen

Stelle vor Beginn dieses Tutorials sicher, dass du Folgendes installiert hast:

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

Setup

# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy pyscf qiskit qiskit-addon-sqd qiskit-ibm-runtime rustworkx

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 ein Quantenproblem abbilden

In diesem Tutorial finden wir eine Approximation des Grundzustands des Moleküls im Gleichgewicht im cc-pVDZ-Basissatz. Zuerst spezifizieren wir das Molekül und seine Eigenschaften.

# 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

Vor dem Konstruieren des LUCJ-Ansatz-Schaltkreises führen wir zunächst eine CCSD-Berechnung in der folgenden Code-Zelle 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 einen geschlossenschaligen Hartree-Fock-Zustand hat, verwenden wir die spin-balancierte Variante des UCJ-Ansatzes, UCJOpSpinBalanced. Wir übergeben Wechselwirkungspaare, die für eine Heavy-Hex-Gitter-Qubit-Topologie geeignet sind (siehe den Hintergrundabschnitt zum LUCJ-Ansatz). Wir setzen optimize=True in der from_t_amplitudes-Methode, um die "komprimierte" Doppelfaktorisierung der $t_2$-Amplituden zu ermöglichen (siehe The local unitary cluster Jastrow (LUCJ) ansatz in der ffsim-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 Quantenhardware-Ausführung optimieren

Als Nächstes optimieren wir den Schaltkreis für eine Ziel-Hardware.

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 hardware-kompatibel zu machen.

• Wähle physikalische Qubits (initial_layout) aus der Ziel-Hardware aus, die dem zuvor beschriebenen "Zickzack"-Muster entspricht. Das Anlegen von Qubits in diesem Muster führt zu einem effizienten hardware-kompatiblen Schaltkreis mit weniger Gates. Hier fügen wir Code ein, um die Auswahl des "Zickzack"-Musters zu automatisieren, der eine Bewertungsheuristik verwendet, um die mit dem ausgewählten Layout verbundenen Fehler zu minimieren.
• Generiere einen gestuften Pass-Manager mit der Funktion generate_preset_pass_manager von Qiskit mit deiner Wahl von backend und initial_layout.
• Setze die pre_init-Stufe deines gestuften Pass-Managers auf ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT enthält Qiskit-Transpiler-Pässe, die Gates in Orbitalrotationen zerlegen und dann die Orbitalrotationen zusammenführen, was zu weniger Gates im endgültigen Schaltkreis führt.
• Führe den Pass-Manager auf deinem Schaltkreis aus.

Code für automatisierte 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: Ausführen mit Qiskit-Primitiven

Nach der Optimierung des Schaltkreises für die Hardware-Ausführung sind wir bereit, ihn auf der Ziel-Hardware auszuführen und Stichproben für die Grundzustandsenergieabschätzung zu sammeln. Da wir nur einen Schaltkreis haben, verwenden wir den Job-Ausführungsmodus von Qiskit Runtime 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: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Jetzt schätzen wir die Grundzustandsenergie des Hamiltonians mit der Funktion diagonalize_fermionic_hamiltonian. Diese Funktion führt das selbstkonsistente Konfigurationswiederherstellungsverfahren aus, um die verrauschten Quantenstichproben iterativ zu verfeinern und die Energieabschätzung zu verbessern. Wir übergeben eine Callback-Funktion, damit wir die Zwischenergebnisse für eine spätere Analyse speichern können. Siehe 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

Visualisierung der Ergebnisse

Der erste Plot zeigt, dass wir nach ein paar Iterationen die Grundzustandsenergie innerhalb von ~41 mH schätzen (chemische Genauigkeit wird typischerweise als 1 kcal/mol $\approx$ 1.6 mH akzeptiert). Die Energie kann verbessert werden, indem mehr Iterationen der Konfigurationswiederherstellung erlaubt werden oder die Anzahl der Stichproben pro Batch erhöht wird.

Der zweite Plot zeigt die durchschnittliche Besetzung jedes Raumorbitals nach der letzten Iteration. Wir können sehen, dass sowohl die Spin-Up- als auch die Spin-Down-Elektronen die ersten fünf Orbitale mit hoher Wahrscheinlichkeit in unseren Lösungen 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 diesem Tutorial zu geben. Deine Einblicke helfen uns, unsere Inhaltsangebote und Benutzererfahrung zu verbessern.

Link zur Umfrage