Sample-basierte Quantendiagonalisierung von'm Chemie-Hamiltonian

Geschätzte Laufzeit: unta 'ne Minute uff'm Heron r2 Prozessor (HINWEIS: Dit is nur 'ne Schätzung. Deine Laufzeit kann vaschiedn sein.)

Hintajrund

In däm Tutorial zeiln wa, wie ma verrauschte Quantensamples nachbearbeitet, um 'ne Näherung an'n Grundzustand vom Stickstoffmolekül $\text{N}_2$ bei da Gleichjewichts-Bindungslänge zu findn, mithilfe vom Sample-basierten Quantendiagonalisierungs-Algorithmus (SQD), für Samples aus'm 59-Qubit-Quantenschaltkreis (52 System-Qubits + 7 Ancilla-Qubits). 'Ne Qiskit-basierte Implementierung findt sich in de SQD Qiskit Addons; mehr Details jibt's in de entsprechenden Docs mit 'nem einfachen Beispiel zum Loslärn.

SQD is 'ne Technik zum Findn von Eigenwerten un Eigenvektoren von Quantenoperatoren, zum Beispiel von'm Quantensystem-Hamiltonian, die Quantencomputing un verteiltes klassisches Computing zusammenbringt. Klassisches verteiltes Computing wird einjesetzt, um Samples zu vararb eiten, die von'm Quantenprozessor stammen, un um'n Ziel-Hamiltonian in'n von ihnen aufgespannten Unterraum zu projizieren un zu diagonalisieren. Dit macht SQD robust jejenüba Samples, die durch Quantenrausch verschmutzt sind, un ermöglicht et, mit jroßen Hamiltonians umzujehn — wie Chemie-Hamiltonians mit Millionen von Wechselwirkungstermen — die für jede exakte Diagonalisierungsmethode viel zu jroß wärn. 'N SQD-basierter Workflow hat folgende Schritte:

1. Wähl 'nen Circuit-Ansatz un wend ihn uff'm Quantencomputer uff'n Referenzzustand an (in däm Fall den Hartree-Fock-Zustand).
2. Sammel Bitstrings aus'm resultierenden Quantenzustand.
3. Führ dit Prozedur zur selbstkonsistenten Konfigurationswiederherstellung uff den Bitstrings durch, um de Näherung an'n Grundzustand zu kriegen.

SQD is bekannt dafür, jut zu funktionier'n, wenn da Ziel-Eigenzustand dünn besetzt is: Die Wellenfunktion is uff 'ne Menge von Basiszuständen $\mathcal{S} = \{|x\rangle \}$ gestützt, deren Größe nich exponentiell mit da Problemgröße wächst.

Quantenchemie

Die Eigenschaften von Moleküln werden jroßteils durch de Struktur da Elektronen in ihnen bestimmt. Als fermionische Teilchen können Elektronen mithilfe vom mathematischen Formalismus der zweiten Quantisierung beschrieben werden. Die Idee is, datt es 'ne Anzahl von Orbitalen jibt, die jeweiß entweda leer oda von'm Fermion besetzt sein können. 'N System von $N$ Orbitalen wird durch 'ne Menge fermionischer Vernichtungsoperatoren $\{\hat{a}_p\}_{p=1}^N$ beschrieben, die de fermionischen Antikommutatoren erfülln:

$$\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*}$$

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

Bis jetzt hat unsere Darstellung den Spin noch nich berücksichtigt, der 'ne fundamentale Eigenschaft von Fermionen is. Wenn ma den Spin mit einbezieht, kommen die Orbitale paarweise vor — man nennt se dann Raumorbitale. Jedes Raumorbital besteht aus zwei Spinorbitalen, eins mit da Bezeichnung "Spin-$\alpha$" un eins mit "Spin-$\beta$". Wa schreibn dann $\hat{a}_{p\sigma}$ für den Vernichtungsoperator, der zum Spinorbital mit Spin $\sigma$ ($\sigma \in \{\alpha, \beta\}$) im Raumorbital $p$ jehört. Wenn wa $N$ als Anzahl da Raumorbitale nehmen, jibt et insjesamt $2N$ Spinorbitale. Da Hilbertraum von däm System 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.

Da Hamiltonian von'm Molekülsystem lässt sich schreibn 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}$ un $h_{prqs}$ komplexe Zahlen sind, die man Molekülintegrale nennt un die aus da Spezifikation vom Molekül mithilfe von'm Computerprogramm berechnet werden können. In däm Tutorial berechnen wa die Integrale mithilfe vom PySCF-Softwarepaket.

Für Details darüba, wie da Molekül-Hamiltonian hergeleitet wird, kiek dir'n Lehrbuch zur Quantenchemie an (zum Beispiel Modern Quantum Chemistry von Szabo un Ostlund). Für 'ne Hochebenen-Erklärung davon, wie Quantenchemieprobleme uff Quantencomputa abjebildet werden, kiek dir de Vorlesung Mapping Problems to Qubits von da Qiskit Global Summer School 2024 an.

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

SQD braucht 'nen Quantenschaltkreis-Ansatz, aus dem Samples jezogn werden können. In däm Tutorial nutzen wa den lokalen unitären Cluster-Jastrow-(LUCJ-)Ansatz, weil er physikalische Motivation un Hardware-Freundlichkeit jut kombiniert.

Da LUCJ-Ansatz is 'ne spezialisierte Form vom allgemeinen unitären Cluster-Jastrow-(UCJ-)Ansatz, 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$ 'n Referenzzustand is — oft der Hartree-Fock-Zustand — un die $\hat{K}_\mu$ un $\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 wa den Teilchenzahloperator definiert haben als

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

Da Operator $e^{\hat{K}_\mu}$ is 'ne Orbitalrotation, die mithilfe bekannter Algorithmen in linearer Tiefe un mit linearer Konnektivität implementiert werden kann. Da $e^{i \mathcal{J}_k}$-Term vom UCJ-Ansatz erfordert entweder All-to-All-Konnektivität oda die Verwendung von'm fermionischen Swap-Netzwerk, wat für verrauschte, vor-fehlertolerante Quantenprozessoren mit bejrenzter Konnektivität schwierig is. Die Idee vom lokalen UCJ-Ansatz is, Dünnheits-Einschränkungen uff die $\mathbf{J}^{\alpha\alpha}$- un $\mathbf{J}^{\alpha\beta}$-Matrizen zu verhängen, die deren Implementierung in konstanter Tiefe uff Qubit-Topologien mit bejrenzter Konnektivität ermöglichen. Die Einschränkungen werden durch 'ne Liste von Indizes angegeben, die anzeigen, welche Matrixeinträge im oberen Dreieck ungleich null sein dürfen (da die Matrizen symmetrisch sind, muss nur das obere Dreieck angegeben werden). Diese Indizes können als Paare von Orbitalen interpretiert werden, die miteinander wechselwirken dürfen.

Als Beispiel betrachten wa 'ne quadratische Gitter-Qubit-Topologie. Wa können die $\alpha$- un $\beta$-Orbitale in parallelen Linien uff dem Gitter anordnen, mit Verbindungen zwischen diesen Linien, die "Sprossen" einer Leiterform bilden — so wie hier:

Mit dieser Anordnung sind Orbitale mit gleichem Spin durch 'ne Linientopologie verbunden, während Orbitale mit verschiedenem Spin verbunden sind, wenn se dasselbe Raumorbital teilen. Dit jibt folgende Indexeinschrä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*}$$

Mit anderen Worten: Wenn die $\mathbf{J}$-Matrizen nur an den angegebenen Indizes im oberen Dreieck ungleich null sind, dann kann da $e^{i \mathcal{J}_k}$-Term uff 'ner quadratischen Topologie ohne Swap-Gates in konstanter Tiefe implementiert werden. Natürlich macht dit Verhängen von solchen Einschränkungen den Ansatz weniger ausdrucksstark, daher können mehr Ansatz-Wiederholungen erforderlich sein.

Die IBM®-Hardware hat 'ne Heavy-Hex-Gitter-Qubit-Topologie — in däm Fall können wa 'n "Zickzack"-Muster verwenden, das unten dargestellt is. In däm Muster werden Orbitale mit gleichem Spin uff Qubits mit Linientopologie abjebildet (rote un blaue Kreise), un 'ne Verbindung zwischen Orbitalen mit verschiedenem Spin is bei jedem 4. Raumorbital vorhanden, wobei die Verbindung durch'n Ancilla-Qubit (lila Kreise) hergestellt wird. In däm Fall lauten die Indexeinschränkungen:

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

Die selbstkonsistente Konfigurationswiederherstellungs-Prozedur is dafür ausjelegt, so viel Signal wie möglich aus verrauschten Quantensamples rauszuholen. Weil da Molekül-Hamiltonian die Teilchenzahl un den Spin Z erhält, macht et Sinn, 'nen Schaltkreis-Ansatz zu wähln, der diese Symmetrien ooch erhält. Wenn er uff den Hartree-Fock-Zustand anjewendet wird, hat da resultierende Zustand im rauschfreien Fall 'ne feste Teilchenzahl un festes Spin Z. Daher sollten die Spin-$\alpha$- un Spin-$\beta$-Hälften jedes aus däm Zustand gesampelten Bitstrings dasselbe Hamming-Gewicht wie im Hartree-Fock-Zustand haben. Wegen da Anwesenheit von Rausch in heutichn Quantenprozessoren werden manche gemessenen Bitstrings diese Eigenschaft verletzen. 'Ne einfache Form von Nachauswahl würde diese Bitstrings wegjeschmeißen, aber dit is Verschwendung, weil die Bitstrings noch 'n bisschen Signal enthalten könnten. Die selbstkonsistente Wiederherstellungsprozedur versucht, etwas von dem Signal in da Nachbearbeitung wiederzuholen. Die Prozedur is iterativ un braucht als Eingabe 'ne Schätzung der durchschnittlichen Besetzungen von jedem Orbital im Grundzustand, die zuerst aus den Roh-Samples berechnet wird. Die Prozedur läuft in 'ner Schleife, un jede Iteration hat folgende Schritte:

1. Für jeden Bitstring, der die angegebenen Symmetrien verletzt, kippe seine Bits mit 'nem probabilistischen Verfahren, das darauf ausgelegt is, den Bitstring näher an die aktuelle Schätzung der durchschnittlichen Orbitalbesetzungen zu bringen, um 'nen neuen Bitstring zu erhalten.
2. Sammel alle alten un neuen Bitstrings, die die Symmetrien erfüllen, un ziehe Teilmengen von 'ner festen Größe, die vorher festjelegt wurde.
3. Für jede Teilmenge von Bitstrings projiziere den Hamiltonian in den von den entsprechenden Basisvektoren aufgespannten Unterraum (kiek dir den vorherigen Abschnitt für 'ne Beschreibung dieser Basisvektoren an) un berechne 'ne Grundzustandsschätzung vom projizierten Hamiltonian uff'm klassischen Computer.
4. Aktualisiere die Schätzung der durchschnittlichen Orbitalbesetzungen mit der Grundzustandsschätzung, die die niedrigste Energie hat.

SQD-Workflow-Diagramm

Da SQD-Workflow is im folgenden Diagramm dargestellt:

Voraussetzungen

Bevor du mit däm Tutorial anfängst, stell sicher, datt du folgendes installiert hast:

• Qiskit SDK v1.0 oda neuer, mit Visualisierungs-Unterstützung
• Qiskit Runtime v0.22 oda neuer (pip install qiskit-ibm-runtime)
• SQD Qiskit Addon v0.11 oda neuer (pip install qiskit-addon-sqd)
• ffsim v0.0.58 oda 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 uff'n Quantenproblem abbilden

In däm Tutorial werden wa 'ne Näherung an'n Grundzustand vom Molekül bei da Gleichjewichtsstruktur im cc-pVDZ-Basissatz findn. Erstmal spezifizieren wa det Molekül un 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

Bevor wa den LUCJ-Ansatz-Schaltkreis aufbauen, führen wa im folgenden Code-Block erstmal 'ne CCSD-Rechnung durch. Die $t_1$- un $t_2$-Amplituden aus dieser Rechnung werden verwendet, um die Parameter vom Ansatz 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 wa ffsim, um den Ansatz-Schaltkreis zu erstellen. Da unser Molekül 'nen closed-shell-Hartree-Fock-Zustand hat, verwenden wa die spinbalancierte Variante vom UCJ-Ansatz, UCJOpSpinBalanced. Wa übergeben Wechselwirkungspaare, die für 'ne Heavy-Hex-Gitter-Qubit-Topologie passend sind (kiek dir den Hintajrund-Abschnitt zum LUCJ-Ansatz an). Wa setzen optimize=True in da from_t_amplitudes-Methode, um die "komprimierte" Doppel-Faktorisierung da $t_2$-Amplituden zu aktivieren (kiek dir The local unitary cluster Jastrow (LUCJ) ansatz in ffsims Dokumentation für Details an).

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 optimiern

Als nächstes optimieren wa den Schaltkreis für 'ne 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

Wa empfehlen die folgenden Schritte, um den Ansatz zu optimieren un hardware-kompatibel zu machen.

• Wähl physikalische Qubits (initial_layout) von da Ziel-Hardware aus, die dem vorher beschriebenen "Zickzack"-Muster entsprechen. Das Anlegen von Qubits in däm Muster führt zu 'nem effizienten hardware-kompatiblen Schaltkreis mit weniger Gates. Hier haben wa 'n bisschen Code, der die Auswahl vom "Zickzack"-Muster automatisiert un 'ne Heuristik zur Fehlerbewertung verwendet, um die Fehler der gewählten Anordnung zu minimieren.
• Erzeug 'nen mehrstufigen Pass-Manager mit da Funktion generate_preset_pass_manager von Qiskit mit deiner Wahl von backend un initial_layout.
• Setz die pre_init-Stufe von deinem mehrstufigen Pass-Manager uff ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT enthält Qiskit-Transpiler-Passes, die Gates in Orbitalrotationen zerlegen un dann die Orbitalrotationen zusammenführen, was zu weniger Gates im finalen Schaltkreis führt.
• Führ den Pass-Manager uff deinem Schaltkreis aus.

Code für die automatische Auswahl vom "Zickzack"-Layout

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ührn

Nachdem wa den Schaltkreis für die Hardware-Ausführung optimiert haben, können wa ihn uff da Ziel-Hardware ausführn un Samples für die Grundzustandsenergieschätzung sammeln. Da wa nur 'nen Schaltkreis haben, verwenden wa Qiskit Runtimes Job-Ausführungsmodus un 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 un Ergebnis im gewünschten klassischen Format zurückliefan

Jetzt schätzen wa die Grundzustandsenergie vom Hamiltonian mit da diagonalize_fermionic_hamiltonian-Funktion. Diese Funktion führt die selbstkonsistente Konfigurationswiederherstellungs-Prozedur durch, um die verrauschten Quantensamples iterativ zu verfeinern un die Energieschätzung zu vabessern. Wa übergeben 'ne Callback-Funktion, damit wa die Zwischenergebnisse für die spätere Analyse speichern können. Kiek dir die API-Dokumentation für Erklärungen zu den Argumenten von diagonalize_fermionic_hamiltonian an.

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 visualisiern

Da erste Plot zeigt, datt wa nach 'n paar Iterationen die Grundzustandsenergie inna ~41 mH schätzen (chemische Jenauigkeit wird üblicherweise als 1 kcal/mol $\approx$ 1.6 mH akzeptiert). Die Energie kann verbessert werden, indem mehr Iterationen von da Konfigurationswiederherstellung erlaubt oda die Anzahl der Samples pro Batch erhöht werden.

Da zweite Plot zeigt die durchschnittliche Besetzung von jedem Raumorbital nach da finalen Iteration. Wa können sehn, datt sowohl die Spin-up- als ooch 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

Mach bitte diese kurze Umfrage, um Feedback zu däm Tutorial zu geben. Deine Einblicke helfen uns, unsere Inhalte un die Benutzererfahrung zu vabessern.

Link zur Umfrage