Stichprobenbasierte Quantendiagonalisierung von an Chemie-Hamiltonian

Nutzungsschätzung: under a Minute auf an Heron r2 Prozessor (HINWEIS: Des is bloß a Schätzung. Ihri Laufzeit ko variieren.)

Hintergrund

In dem Tutorial zeign mia, wia verrauschte Quantenstichproben nachbearbeitet wern, um a Approximation vom Grundzustand vom Stickstoffmolekül $\text{N}_2$ bei Gleichgewichtsbindungslänge z'findn. Dabei verwend mia den stichprobenbasiertn Quantendiagonalisierungsalgorithmus (SQD) für Stichproben aus an 59-Qubit-Quantenschaltkreis (52 System-Qubits + 7 Ancilla-Qubits). A Qiskit-basierte Implementierung is verfügbar in de SQD Qiskit Addons, weitere Details findts in da entsprechenden Dokumentation mit an einfachn Beispiel zum Einstieg.

SQD is a Technik zum Findn von Eigenwerten und Eigenvektorn von Quantenoperatorn, wia dem Hamiltonian von an Quantensystem, unter Verwendung von Quanten- und verteiltem klassischem Computing mitanand. Des verteilte klassische Computing wird verwendet, um Stichproben von an Quantenprozessor z'verarbeiten und an Ziel-Hamiltonian in an durch sie aufgspanntn Unterraum z'projizieren und z'diagonalisieren. Des ermöglicht SQD, robust gegenüber durch Quantenrauschen verfälschte Stichproben z'sei und mit großen Hamiltonians umz'gehn, wia Chemie-Hamiltonians mit Millionen von Wechselwirkungstermen, de jenseits da Reichweite exakter Diagonalisierungsmethoden liegen. An SQD-basierter Workflow hat folgende Schritt:

1. Wähl an Schaltkreis-Ansatz und wend ihn auf an Quantencomputer auf an Referenzzustand o (in dem Fall den Hartree-Fock-Zustand).
2. Sampl Bitstrings aus dem resultierenden Quantenzustand.
3. Führ des selbstkonsistente Konfigurationswiederherstellungsverfahren auf de Bitstrings aus, um de Approximation vom Grundzustand z'erhaltn.

SQD funktioniert bekanntermaßen guad, wenn da Ziel-Eigenzustand dünn besetzt is: De Wellenfunktion wird von ana Menge von Basiszuständn $\mathcal{S} = \{|x\rangle \}$ getragn, deren Größe ned exponenziell mit da Problemgröße wachst.

Quantenchemie

De Eigenschaftn von Moleküln wern weitgehend durch de Struktur da Elektronen in ihnen bestimmt. Als fermionische Teilchen können Elektronen mit an mathematischen Formalismus namens Zweitquantisierung beschriebm wern. De Idee is, dass es a Anzahl von Orbitalen gibt, von denen jedes entweder leer is oder von an Fermion besetzt wird. A System von $N$ Orbitalen wird durch an Satz fermionischer Vernichtungsoperatorn $\{\hat{a}_p\}_{p=1}^N$ beschriebm, de de fermionischen Antikommutatorrelationen 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 Operator $\hat{a}_p^\dagger$ wird Erzeugungsoperator g'nannd.

Bis jetzt hat unsere Darstellung den Spin ned berücksichtigt, da a fundamentale Eigenschaft von Fermionen is. Bei da Berücksichtigung vom Spin komman de Orbitale in Paaren vor, de Raumorbitale g'nannd wern. Jedes Raumorbital besteht aus zwoa Spinorbitalen, von denen oas mit "Spin-$\alpha$" und oas mit "Spin-$\beta$" bezeichnet wird. Mia schreibm dann $\hat{a}_{p\sigma}$ für den Vernichtungsoperator, da mit dem Spinorbital mit Spin $\sigma$ ($\sigma \in \{\alpha, \beta\}$) im Raumorbital $p$ assoziiert is. Wenn mia $N$ als Anzahl da Raumorbitale nehman, gibt's insgesamt $2N$ Spinorbitale. Da Hilbert-Raum von dem System wird von $2^{2N}$ orthonormalen Basisvektorn aufg'spannt, de mit zweiteiligen Bitstrings bezeichnet wern: $\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$.

Da Hamiltonian von an molekularen System ko g'schriebm wern 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 de $h_{pr}$ und de $h_{prqs}$ komplexe Zahlen san, de molekulare Integrale g'nannd wern und aus da Spezifikation vom Molekül mit an Computerprogramm berechnet wern können. In dem Tutorial berechnan mia de Integrale mit dem Softwarepaket PySCF.

Für Details darüber, wia da molekulare Hamiltonian hergeleitet wird, konsultiert a Lehrbuch zur Quantenchemie (zum Beispiel Modern Quantum Chemistry von Szabo und Ostlund). Für a übergeordnete Erklärung, wia Quantenchemieprobleme auf Quantencomputer abbildet wern, schaugts euch de Vorlesung Mapping Problems to Qubits von da Qiskit Global Summer School 2024 o.

Local Unitary Cluster Jastrow (LUCJ) Ansatz

SQD braucht an Quantenschaltkreis-Ansatz, aus dem Stichproben g'zogn wern können. In dem Tutorial verwend mia den Local Unitary Cluster Jastrow (LUCJ) Ansatz wegen seina Kombination aus physikalischer Motivation und Hardware-Freundlichkeit.

Da LUCJ-Ansatz is a spezialisierte Form vom allgemeinen Unitary Cluster Jastrow (UCJ) Ansatz, da de 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 da Hartree-Fock-Zustand, und de $\hat{K}_\mu$ und de $\hat{J}_\mu$ de Form ham

$$\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 mia den Teilchenzahloperator definiert ham

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

Da Operator $e^{\hat{K}_\mu}$ is a Orbitalrotation, de mit bekannten Algorithmen in linearer Tiefe und unter Verwendung linearer Konnektivität implementiert wern ko.

De Implementierung vom $e^{i \mathcal{J}_k}$ Term vom UCJ-Ansatz erfordert entweder vollständige Konnektivität oder de Verwendung von an fermionischen Swap-Netzwerk, was für verrauschte präfehlertolerante Quantenprozessorn mit begrenzter Konnektivität a Herausforderung is. De Idee vom lokalen UCJ-Ansatz is, Dünnbesetztheitsbedingungen auf de $\mathbf{J}^{\alpha\alpha}$- und $\mathbf{J}^{\alpha\beta}$-Matrizen aufz'erlegn, de es ermöglichen, sie in konstanter Tiefe auf Qubit-Topologien mit begrenzter Konnektivität z'implementieren. De Bedingungen wern durch a Liste von Indizes spezifiziert, de o'zeigen, welche Matrixeinträge im obern Dreieck von null verschieden sei dürfen (da de Matrizen symmetrisch san, muss nur des obere Dreieck spezifiziert wern). Dise Indizes können als Orbitalpaare interpretierd wern, de mitanand interagieren dürfen.

Betrachtets als Beispiel a quadratische Gitter-Qubit-Topologie. Mia können de $\alpha$- und $\beta$-Orbitale in parallelen Linien auf dem Gitter platzieren, wobei Verbindungen zwischen denen Linien "Sprossen" von ana Leiterform bilden, wie folgt:

Bei dem Setup san Orbitale mit demselbn Spin mit ana Linientopologie verbunden, während Orbitale mit unterschiedlichem Spin verbunden san, wenn sie dasselbe Raumorbital teilen. Des ergibt de folgenden Indexbedingungen für de $\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 Wortn: Wenn de $\mathbf{J}$-Matrizen nur an denen a'g'gebenen Indizes im obern Dreieck von null verschieden san, ko da $e^{i \mathcal{J}_k}$ Term auf ana quadratischen Topologie ohne Verwendung von Swap-Gates in konstanter Tiefe implementiert wern. Natürlich macht des Auferlegen solcher Bedingungen auf den Ansatz ihn weniger ausdrucksstark, sodass möglicherweise mehr Ansatz-Wiederholungen nötig san.

De IBM® Hardware hat a Heavy-Hex-Gitter-Qubit-Topologie, in welchem Fall mia a "Zickzack"-Muster verwenden können, des untn dargestellt is. In dem Muster wern Orbitale mit demselbn Spin auf Qubits mit ana Linientopologie abbildet (rote und blaue Kreise), und a Verbindung zwischen Orbitalen unterschiedlicher Spins is an jedem 4. Raumorbital vorhanden, wobei de Verbindung durch an Ancilla-Qubit (violette Kreise) ermöglicht wird. In dem Fall san de 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

Des selbstkonsistente Konfigurationswiederherstellungsverfahren is dazu ausgelegt, so viel Signal wia möglich aus verrauschten Quantenstichproben rauszholn. Da da molekulare Hamiltonian de Teilchenzahl und Spin-Z erhält, is es sinnvoll, an Schaltkreis-Ansatz z'wähln, da dise Symmetrien a erhält. Wenn er auf den Hartree-Fock-Zustand o'g'wendet wird, hat da resultierend Zustand im rauschfreien Fall a feste Teilchenzahl und Spin-Z. Daher sollatn de Spin-$\alpha$- und Spin-$\beta$-Hälften von jedem aus dem Zustand g'sampltn Bitstring dasselbe Hamming-Gewicht wia im Hartree-Fock-Zustand ham. Wegen dem Vorhandensein von Rauschen in aktuellen Quantenprozessorn wern manche g'messene Bitstrings dise Eigenschaft verletzen. A einfache Form da Postselektion würde dise Bitstrings verwerfen, aber des is verschwenderisch, weil dise Bitstrings vielleicht no a weng Signal enthaltn. Des selbstkonsistente Wiederherstellungsverfahren versucht, an Teil von dem Signal in da Nachbearbeitung wiederherzustelln. Des Verfahren is iterativ und braucht als Eingabe a Schätzung da durchschnittlichn Besetzungen von jedem Orbital im Grundzustand, de zuerst aus de Rohstichproben berechnet wird. Des Verfahren wird in ana Schleife ausführt, und jede Iteration hat de folgenden Schritt:

1. Für jeden Bitstring, da de spezifiziertn Symmetrien verletzt, flipp seine Bits mit an probabilistischen Verfahren, des dazu ausgelegt is, den Bitstring näher an de aktuelle Schätzung da durchschnittlichn Orbitalbesetzungen z'bringn, um an neuen Bitstring z'erhaltn.
2. Samml alli altn und neuen Bitstrings, de de Symmetrien erfülln, und entnimm Teilmengen von ana im Voraus g'wähltn festen Größe.
3. Für jede Teilmenge von Bitstrings projizier den Hamiltonian in den durch de entsprechenden Basisvektorn aufg'spanntn Unterraum (schaugts in den vorigen Abschnitt für a Beschreibung von denen Basisvektorn) und berechne a Grundzustandsschätzung vom projiziertn Hamiltonian auf an klassischn Computer.
4. Aktualisier de Schätzung da durchschnittlichn Orbitalbesetzungen mit da Grundzustandsschätzung mit da niedrigstn Energie.

SQD-Workflow-Diagramm

Da SQD-Workflow is im folgenden Diagramm dargestellt:

Anforderungen

Stellt's vor dem Beginn von dem Tutorial sicher, dass's Folgendes installierd ham:

• 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 a Quantenproblem abbilden

In dem Tutorial findn mia a Approximation vom Grundzustand vom Molekül im Gleichgewicht im cc-pVDZ-Basissatz. Zuerst spezifizieren mia des Molekül und seine Eigenschaftn.

# 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 vom LUCJ-Ansatz-Schaltkreis führn mia zunächst a CCSD-Berechnung in da folgenden Code-Zelle durch. De $t_1$- und $t_2$-Amplituden aus dera Berechnung wern verwendet, um de Parameter vom Ansatz z'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 verwend mia ffsim, um den Ansatz-Schaltkreis z'erstelln. Da unser Molekül an geschlossenschaligen Hartree-Fock-Zustand hat, verwend mia de spin-balancierte Variante vom UCJ-Ansatz, UCJOpSpinBalanced. Mia übergeben Wechselwirkungspaare, de für a Heavy-Hex-Gitter-Qubit-Topologie geeignet san (schaugts in den Hintergrundabschnitt zum LUCJ-Ansatz). Mia setzn optimize=True in da from_t_amplitudes-Methode, um de "komprimierte" Doppelfaktorisierung da $t_2$-Amplituden z'ermöglichn (schaugts The local unitary cluster Jastrow (LUCJ) ansatz in da 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 mia den Schaltkreis für a 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

Mia empfehlen de folgenden Schritt, um den Ansatz z'optimieren und hardware-kompatibel z'machn.

• Wähl physikalische Qubits (initial_layout) aus da Ziel-Hardware aus, de dem zuvor beschriebenen "Zickzack"-Muster entsprechn. Des Anlegen von Qubits in dem Muster führt zu an effizientn hardware-kompatibln Schaltkreis mit weniger Gates. Do fügen mia Code ei, um de Auswahl vom "Zickzack"-Muster z'automatisieren, da a Bewertungsheuristik verwendet, um de mit dem ausgewähltn Layout verbundenen Fehler z'minimieren.
• Generier an gestuftn Pass-Manager mit da Funktion generate_preset_pass_manager von Qiskit mit Eana Wahl von backend und initial_layout.
• Setz de pre_init-Stufe von Eurem gestuftn Pass-Manager auf ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT enthält Qiskit-Transpiler-Pässe, de Gates in Orbitalrotationen zerlegn und dann de Orbitalrotationen zusammenführn, was zu weniger Gates im endgültigem Schaltkreis führt.
• Führ den Pass-Manager auf Eurem Schaltkreis aus.

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

Nach da Optimierung vom Schaltkreis für de Hardware-Ausführung san mia bereit, ihn auf da Ziel-Hardware auszuführn und Stichproben für de Grundzustandsenergieabschätzung z'sammln. Da mia nur an Schaltkreis ham, verwend mia den Job-Ausführungsmodus von Qiskit Runtime und führn unsern 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 vom Ergebnis im g'wünschtn klassischn Format

Jetzt schätzn mia de Grundzustandsenergie vom Hamiltonian mit da Funktion diagonalize_fermionic_hamiltonian. Dise Funktion führt des selbstkonsistente Konfigurationswiederherstellungsverfahren aus, um de verrauschten Quantenstichproben iterativ z'verfeinern und de Energieabschätzung z'verbessern. Mia übergeben a Callback-Funktion, damit mia de Zwischenergebnisse für a spätere Analyse speichern können. Schaugts in de API-Dokumentation für Erklärungen da 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 da Ergebnisse

Da erste Plot zeigt, dass mia nach a paar Iterationen de Grundzustandsenergie innerhalb von ~41 mH schätzn (chemische Genauigkeit wird typischerweise als 1 kcal/mol $\approx$ 1.6 mH akzeptierd). De Energie ko verbessert wern, indem mehr Iterationen da Konfigurationswiederherstellung erlaubt wern oder de Anzahl da Stichproben pro Batch erhöht wird.

Da zweite Plot zeigt de durchschnittliche Besetzung von jedem Raumorbital nach da letzten Iteration. Mia können segn, dass sowohl de Spin-Up- als a de Spin-Down-Elektronen de 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 nehmt's an dera kurzen Umfrage teil, um Feedback zu dem Tutorial z'gebm. Eure Einblicke helfn uns, unsere Inhaltsangebote und Benutzererfahrung z'verbessern.

Link zur Umfrage