Přeskočit na hlavní obsah

Kvantová diagonalizace na základě vzorků pro chemický Hamiltonián

Odhadovaná spotřeba: méně než jedna minuta na procesoru Heron r2 (POZNÁMKA: Jde pouze o odhad. Skutečná doba běhu se může lišit.)

Pozadí​

V tomto tutoriálu ukážeme, jak post-procesovat hlučné kvantové vzorky a najít aproximaci základního stavu molekuly dusíku N2\text{N}_2 při rovnovážné délce vazby pomocí algoritmu kvantové diagonalizace na základě vzorků (SQD), pro vzorky získané z 59-qubitového kvantového obvodu (52 systémových Qubitů + 7 ancillárních Qubitů). Implementace v Qiskitu je dostupná v SQD Qiskit addons, více podrobností najdeš v příslušné dokumentaci s jednoduchým příkladem pro začátek. SQD je technika pro hledání vlastních hodnot a vlastních vektorů kvantových operátorů, jako je Hamiltonián kvantového systému, pomocí kombinace kvantového a distribuovaného klasického výpočtu. Klasický distribuovaný výpočet slouží ke zpracování vzorků získaných z kvantového procesoru a k projekci a diagonalizaci cílového Hamiltoniánu v podprostoru jimi generovaném. To umožňuje SQD být robustní vůči vzorkům poškozeným kvantovým šumem a pracovat s velkými Hamiltoniány, jako jsou chemické Hamiltoniány s miliony interakčních členů, které jsou mimo dosah jakýchkoliv metod přesné diagonalizace. Pracovní postup SQD se skládá z následujících kroků:

  1. Zvol Circuit ansatz a aplikuj ho na kvantovém počítači na referenční stav (v tomto případě stav Hartree-Fock).
  2. Vzorkuj bitové řetězce z výsledného kvantového stavu.
  3. Spusť proceduru self-consistent configuration recovery na bitových řetězcích, aby ses dostal k aproximaci základního stavu.

Je známo, že SQD funguje dobře, když je cílový vlastní stav řídký: vlnová funkce je podpořena na množině bázových stavů S={∣x⟩}\mathcal{S} = \{|x\rangle \}, jejíž velikost neroste exponenciálně s velikostí problému.

Kvantová chemie​

Vlastnosti molekul jsou z velké části určeny strukturou elektronů v nich. Jako fermionové částice lze elektrony popsat pomocí matematického formalismu zvaného druhá kvantizace. Myšlenka spočívá v tom, že existuje určitý počet orbitalů, z nichž každý může být buď prázdný, nebo obsazený fermionem. Systém NN orbitalů je popsán množinou fermionových annihilačních operátorů {a^p}p=1N\{\hat{a}_p\}_{p=1}^N, které splňují fermionové antikomutační relace,

a^pa^q+a^qa^p=0,a^pa^q†+a^q†a^p=δpq.\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*}

Adjungovaný operátor a^p†\hat{a}_p^\dagger se nazývá kreační operátor.

Dosud jsme v našem výkladu nebrali v úvahu spin, který je základní vlastností fermionů. Při zohlednění spinu se orbitaly sdružují do párů nazývaných prostorové orbitaly. Každý prostorový orbital se skládá ze dvou spinových orbitalů, z nichž jeden je označen „spin-α\alpha" a druhý „spin-β\beta". Píšeme pak a^pσ\hat{a}_{p\sigma} pro annihilační operátor přidružený ke spinovému orbitalu se spinem σ\sigma (σ∈{α,β}\sigma \in \{\alpha, \beta\}) v prostorovém orbitalu pp. Pokud NN označuje počet prostorových orbitalů, pak existuje celkem 2N2N spinových orbitalů. Hilbertův prostor tohoto systému je generován 22N2^{2N} ortonormálními bázovými vektory označenými dvoudílnými bitovými řetězci ∣z⟩=∣zβzα⟩=∣zβ,N⋯zβ,1zα,N⋯zα,1⟩\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.

Hamiltonián molekulárního systému lze zapsat jako

H^=∑prσhpr a^pσ†a^rσ+12∑prqsστhprqs a^pσ†a^qτ†a^sτa^rσ,\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},

kde hprh_{pr} a hprqsh_{prqs} jsou komplexní čísla nazývaná molekulární integrály, která lze vypočítat ze specifikace molekuly pomocí počítačového programu. V tomto tutoriálu počítáme integrály pomocí softwarového balíku PySCF.

Podrobnosti o tom, jak se odvozuje molekulární Hamiltonián, najdeš v učebnici kvantové chemie (například Modern Quantum Chemistry od Szaba a Ostlunda). Přehledné vysvětlení, jak se problémy kvantové chemie mapují na kvantové počítače, nabízí přednáška Mapping Problems to Qubits z Qiskit Global Summer School 2024.

Ansatz lokálního unitárního shlukového Jastrow (LUCJ)​

SQD vyžaduje ansatz kvantového obvodu, z něhož se vzorkuje. V tomto tutoriálu použijeme ansatz lokálního unitárního shlukového Jastrow (LUCJ) díky jeho kombinaci fyzikální motivace a přátelskosti k hardwaru.

Ansatz LUCJ je specializovanou formou obecného ansatzu unitárního shlukového Jastrow (UCJ), který má tvar

∣Ψ⟩=∏μ=1LeK^μeiJ^μe−K^μ∣Φ0⟩ \lvert \Psi \rangle = \prod_{\mu=1}^L e^{\hat{K}_\mu} e^{i \hat{J}_\mu} e^{-\hat{K}_\mu} | \Phi_0 \rangle

kde ∣Φ0⟩\lvert \Phi_0 \rangle je referenční stav, často volený jako stav Hartree-Fock, a K^μ\hat{K}_\mu a J^μ\hat{J}_\mu mají tvar

K^μ=∑pq,σKpqμ a^pσ†a^qσ†  ,  J^μ=∑pq,στJpq,στμ n^pσn^qτ  ,\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} \;,

kde jsme definovali operátor počtu částic

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

Operátor eK^μe^{\hat{K}_\mu} je orbitální rotace, která lze implementovat pomocí známých algoritmů v lineární hloubce a s využitím lineárního propojení. Implementace členu eiJke^{i \mathcal{J}_k} ansatzu UCJ vyžaduje buď propojení všech se všemi, nebo použití fermionové swap sítě, což je náročné pro hlučné pre-fault-tolerantní kvantové procesory s omezeným propojením. Myšlenka lokálního ansatzu UCJ spočívá v zavedení podmínek řídkosti na matice Jαα\mathbf{J}^{\alpha\alpha} a Jαβ\mathbf{J}^{\alpha\beta}, které umožňují jejich implementaci v konstantní hloubce na qubitových topologiích s omezeným propojením. Podmínky jsou specifikovány seznamem indexů označujících, které položky matice v horním trojúhelníku smějí být nenulové (protože matice jsou symetrické, stačí specifikovat pouze horní trojúhelník). Tyto indexy lze interpretovat jako dvojice orbitalů, kterým je povoleno interagovat.

Jako příklad uvažujme topologii Qubitů čtvercové mřížky. Orbitaly α\alpha a β\beta můžeme umístit do rovnoběžných přímek v mřížce, přičemž spojení mezi těmito přímkami tvoří „příčky" žebříku, takto:

Diagram mapování Qubitů pro ansatz LUCJ na čtvercové mřížce

Při tomto uspořádání jsou orbitaly stejného spinu propojeny topologií přímky, zatímco orbitaly různého spinu jsou propojeny, pokud sdílejí stejný prostorový orbital. To vede k následujícím indexovým omezením na matice J\mathbf{J}:

Jαα:{ (p,p+1)  ,  p=0,…,N−2 }Jαβ:{ (p,p)  ,  p=0,…,N−1 }\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*}

Jinými slovy, pokud jsou matice J\mathbf{J} nenulové pouze na uvedených indexech v horním trojúhelníku, pak lze člen eiJke^{i \mathcal{J}_k} implementovat na čtvercové topologii bez použití swap Gate, v konstantní hloubce. Samozřejmě, zavedení takových omezení na ansatz snižuje jeho expresivitu, takže může být zapotřebí více opakování ansatzu.

Hardware IBM® má topologii Qubitů heavy-hex mřížky, v tomto případě můžeme použít vzor „cik-cak", znázorněný níže. V tomto vzoru jsou orbitaly stejného spinu mapovány na Qubity s topologií přímky (červené a modré kroužky) a spojení mezi orbitaly různého spinu se vyskytuje každé 4. prostorové orbitaly, přičemž spojení je zprostředkováno ancilárním Qubitem (fialové kroužky). V tomto případě jsou indexová omezení

Jαα:{ (p,p+1)  ,  p=0,…,N−2 }Jαβ:{ (p,p)  ,  p=0,4,8,…(p≤N−1) }\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*}

Diagram mapování Qubitů pro ansatz LUCJ na heavy-hex mřížce

Self-consistent configuration recovery​

Procedura self-consistent configuration recovery je navržena k extrakci co největšího signálu z hlučných kvantových vzorků. Protože molekulární Hamiltonián zachovává počet částic a spin Z, dává smysl volit Circuit ansatz, který tyto symetrie také zachovává. Při aplikaci na stav Hartree-Fock má výsledný stav v bezhlučném nastavení pevný počet částic a spin Z. Proto by spin-α\alpha a spin-β\beta části každého bitového řetězce vzorkovaného z tohoto stavu měly mít stejnou Hammingovu váhu jako ve stavu Hartree-Fock. Kvůli přítomnosti šumu v současných kvantových procesorech budou některé naměřené bitové řetězce tuto vlastnost porušovat. Jednoduchá forma postselektion by tyto bitové řetězce zahodila, ale to je plýtvání, protože mohou stále obsahovat určitý signál. Procedura self-consistent recovery se pokouší část tohoto signálu získat zpět v post-processingu. Procedura je iterativní a jako vstup vyžaduje odhad průměrného obsazení každého orbitalu v základním stavu, který se nejprve vypočítá ze surových vzorků. Procedura běží ve smyčce a každá iterace má tyto kroky:

  1. Pro každý bitový řetězec, který porušuje zadané symetrie, překlop jeho bity pravděpodobnostní procedurou navrženou tak, aby bitový řetězec přiblížila k aktuálnímu odhadu průměrného obsazení orbitalů, a získej nový bitový řetězec.
  2. Shromáždi všechny staré a nové bitové řetězce, které splňují symetrie, a vzorkuj podmnožiny pevné velikosti, zvolené předem.
  3. Pro každou podmnožinu bitových řetězců promítni Hamiltonián do podprostoru generovaného odpovídajícími bázovými vektory (popis těchto bázových vektorů viz předchozí sekci) a vypočítej odhad základního stavu promítnutého Hamiltoniánu na klasickém počítači.
  4. Aktualizuj odhad průměrného obsazení orbitalů pomocí odhadu základního stavu s nejnižší energií.

Diagram pracovního postupu SQD​

Pracovní postup SQD je znázorněn v následujícím diagramu:

Diagram pracovního postupu algoritmu SQD

Požadavky​

Před zahájením tohoto tutoriálu se ujisti, že máš nainstalováno následující:

  • Qiskit SDK v1.0 nebo novější, s podporou vizualizace
  • Qiskit Runtime v0.22 nebo novější (pip install qiskit-ibm-runtime)
  • SQD Qiskit addon v0.11 nebo novější (pip install qiskit-addon-sqd)
  • ffsim v0.0.58 nebo novější (pip install ffsim)

Nastavení​

# 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

Krok 1: Mapování klasických vstupů na kvantový problém​

V tomto tutoriálu najdeme aproximaci základního stavu molekuly při rovnovážné délce vazby v bázi cc-pVDZ. Nejprve specifikujeme molekulu a její vlastnosti.

# 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

Před sestavením Circuit ansatzu LUCJ nejprve provedeme výpočet CCSD v následující buňce kódu. Amplitudy t1t_1 a t2t_2 z tohoto výpočtu budou použity k inicializaci parametrů ansatzu.

# 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

Nyní použijeme ffsim k vytvoření Circuit ansatzu. Protože naše molekula má uzavřený stav Hartree-Fock, použijeme spin-vyvážený variant ansatzu UCJ, UCJOpSpinBalanced. Předáme interakční páry vhodné pro topologii Qubitů heavy-hex mřížky (viz sekci pozadí o ansatzu LUCJ). V metodě from_t_amplitudes nastavíme optimize=True, aby bylo povoleno „komprimované" dvojité faktorizování amplitud t2t_2 (viz The local unitary cluster Jastrow (LUCJ) ansatz v dokumentaci ffsim pro podrobnosti).

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()

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru​

Dále optimalizujeme Circuit pro cílový 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

Doporučujeme následující kroky k optimalizaci ansatzu a jeho kompatibilitě s hardwarem.

  • Vyber fyzické Qubity (initial_layout) z cílového hardwaru, které odpovídají vzoru „cik-cak" popsanému výše. Rozložení Qubitů v tomto vzoru vede k efektivnímu Circuit kompatibilnímu s hardwarem s méně Gate. Zde uvádíme kód pro automatizaci výběru vzoru „cik-cak", který používá heuristiku skóre k minimalizaci chyb spojených s vybraným rozložením.
  • Vygeneruj staged pass manager pomocí funkce generate_preset_pass_manager z Qiskitu s tvým výběrem backend a initial_layout.
  • Nastav fázi pre_init svého staged pass manageru na ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT zahrnuje průchody Transpileru Qiskitu, které rozkládají Gate na orbitální rotace a následně slučují orbitální rotace, čímž výsledný Circuit obsahuje méně Gate.
  • Spusť pass manager na svém Circuit.
Kód pro automatizovaný výběr rozložení „cik-cak"
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

Krok 3: Spuštění pomocí Qiskit primitiv​

Po optimalizaci Circuit pro spuštění na hardwaru jsme připraveni spustit ho na cílovém hardwaru a shromáždit vzorky pro odhad energie základního stavu. Protože máme pouze jeden Circuit, použijeme režim spuštění úlohy Qiskit Runtime a spustíme náš Circuit.

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})

Krok 4: Post-processing a vrácení výsledku v požadovaném klasickém formátu​

Nyní odhadneme energii základního stavu Hamiltoniánu pomocí funkce diagonalize_fermionic_hamiltonian. Tato funkce provádí postup self-konzistentního obnovení konfigurací, který iterativně zpřesňuje zašuměné kvantové vzorky za účelem zlepšení odhadu energie. Předáváme callback funkci, abychom mohli uložit mezivýsledky pro pozdější analýzu. Vysvětlení argumentů funkce diagonalize_fermionic_hamiltonian najdeš v dokumentaci API.

sampler = Sampler(mode=backend)
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]
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

Vizualizace výsledků​

První graf ukazuje, že po několika iteracích odhadujeme energii základního stavu s přesností ~41 mH (chemická přesnost je typicky považována za 1 kcal/mol ≈\approx 1,6 mH). Energii lze zlepšit povolením více iterací obnovy konfigurace nebo zvýšením počtu vzorků na dávku.

Druhý graf zobrazuje průměrnou obsazenost každého prostorového orbitalu po poslední iteraci. Vidíme, že elektrony se spinem nahoru i se spinem dolů obsazují s vysokou pravděpodobností prvních pět orbitalů v našich řešeních.

# 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()

Output of the previous code cell

Průzkum k tutoriálu​

Prosím, vyplň tento krátký průzkum a poskytni nám zpětnou vazbu k tomuto tutoriálu. Tvoje postřehy nám pomohou zlepšit náš obsah a uživatelský zážitek.

Odkaz na průzkum

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.

Zdroj: IBM Quantum docs — aktualizováno 15. 5. 2026
Anglická verze na doQumentation — aktualizováno 15. 5. 2026
Tento překlad vychází z anglické verze ze dne 9. 4. 2026
Problém s webem nebo překladem?Problém v obsahu? Upravit v IBM Quantum docs