Diagonalizacja kwantowa oparta na próbkowaniu hamiltonianów chemicznych

Szacowane zużycie zasobów: poniżej jednej minuty na procesorze Heron r2 (UWAGA: To jest jedynie szacunek. Twój czas wykonania może się różnić.)

Tło

W tym samouczku pokazujemy, jak poddać przetwarzaniu końcowemu zaszumione próbki kwantowe w celu znalezienia przybliżenia stanu podstawowego cząsteczki azotu $\text{N}_2$ przy równowagowej długości wiązania — z użyciem algorytmu diagonalizacji kwantowej opartej na próbkowaniu (SQD) — dla próbek pobranych z 59-qubitowego Circuit kwantowego (52 qubity systemowe + 7 qubitów ancilla). Implementacja oparta na Qiskit jest dostępna w dodatkach SQD do Qiskit; więcej szczegółów można znaleźć w odpowiedniej dokumentacji oraz w prostym przykładzie na początek.

SQD to technika wyznaczania wartości własnych i wektorów własnych operatorów kwantowych, takich jak Hamiltonian układu kwantowego, przy użyciu obliczeń kwantowych i rozproszonych obliczeń klasycznych łącznie. Klasyczne obliczenia rozproszone służą do przetwarzania próbek uzyskanych z procesora kwantowego oraz do rzutowania i diagonalizacji docelowego Hamiltonianu w podprzestrzeni rozpinanej przez te próbki. Dzięki temu SQD jest odporny na próbki zaburzone przez szum kwantowy i radzi sobie z dużymi Hamiltonianami — takimi jak Hamiltoniany chemiczne zawierające miliony wyrazów oddziaływania — które wykraczają poza zasięg jakichkolwiek metod dokładnej diagonalizacji. Przepływ pracy oparty na SQD składa się z następujących kroków:

1. Wybierz ansatz Circuit i zastosuj go na komputerze kwantowym do stanu referencyjnego (w tym przypadku stanu Hartree-Focka).
2. Pobierz ciągi bitów z wynikowego stanu kwantowego.
3. Uruchom procedurę samospójnego odtwarzania konfiguracji na ciągach bitów, aby uzyskać przybliżenie stanu podstawowego.

Wiadomo, że SQD działa dobrze, gdy docelowy stan własny jest rzadki: funkcja falowa jest oparta na zbiorze stanów bazowych $\mathcal{S} = \{|x\rangle \}$, którego rozmiar nie rośnie wykładniczo wraz z rozmiarem problemu.

Chemia kwantowa

Właściwości cząsteczek są w dużej mierze zdeterminowane przez strukturę elektronów w ich wnętrzu. Jako cząstki fermionowe elektrony można opisać za pomocą formalizmu matematycznego zwanego drugą kwantyzacją. Idea polega na tym, że istnieje pewna liczba orbitali, z których każdy może być albo pusty, albo zajęty przez fermion. Układ $N$ orbitali opisuje zbiór fermionowych operatorów anihilacji $\{\hat{a}_p\}_{p=1}^N$ spełniających fermionowe relacje antykomutacji:

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

Sprzężony $\hat{a}_p^\dagger$ nosi nazwę operatora kreacji.

Jak dotąd nasze rozważania nie uwzględniają spinu, który jest fundamentalną właściwością fermionów. Uwzględniając spin, orbitale tworzą pary zwane orbitalami przestrzennymi. Każdy orbital przestrzenny składa się z dwóch orbitali spinowych: jednego oznaczonego jako „spin-$\alpha$" i jednego oznaczonego jako „spin-$\beta$". Piszemy więc $\hat{a}_{p\sigma}$ dla operatora anihilacji związanego z orbitalem spinowym o spinie $\sigma$ ($\sigma \in \{\alpha, \beta\}$) w orbitaly przestrzennym $p$. Jeżeli przez $N$ oznaczymy liczbę orbitali przestrzennych, to łączna liczba orbitali spinowych wynosi $2N$. Przestrzeń Hilberta tego układu jest rozpięta przez $2^{2N}$ ortonormalnych wektorów bazowych oznaczonych dwuczęściowymi ciągami bitów $\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$.

Hamiltonian układu molekularnego można zapisać jako

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

gdzie $h_{pr}$ i $h_{prqs}$ to liczby zespolone zwane całkami molekularnymi, które można obliczyć na podstawie specyfikacji cząsteczki przy użyciu programu komputerowego. W tym samouczku wyznaczamy całki za pomocą pakietu oprogramowania PySCF.

Szczegółowe informacje na temat wyprowadzenia Hamiltonianu molekularnego można znaleźć w podręczniku chemii kwantowej (na przykład Modern Quantum Chemistry Szabo i Ostlunda). Ogólne omówienie sposobu, w jaki problemy chemii kwantowej są odwzorowywane na komputery kwantowe, zawiera wykład Mapping Problems to Qubits z Qiskit Global Summer School 2024.

Ansatz lokalnego unitarnego klastra Jastrow (LUCJ)

SQD wymaga Circuit ansatz do pobierania próbek. W tym samouczku użyjemy ansatzu lokalnego unitarnego klastra Jastrow (LUCJ) ze względu na jego połączenie fizycznej motywacji i przyjazności dla sprzętu.

Ansatz LUCJ jest wyspecjalizowaną postacią ogólnego ansatzu unitarnego klastra Jastrow (UCJ), który ma postać

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

gdzie $\lvert \Phi_0 \rangle$ jest stanem referencyjnym, często przyjmowanym jako stan Hartree-Focka, a $\hat{K}_\mu$ i $\hat{J}_\mu$ mają postać

$$\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} \;,$$

gdzie zdefiniowaliśmy operator liczby cząstek

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

Operator $e^{\hat{K}_\mu}$ jest rotacją orbitalną, którą można zaimplementować przy użyciu znanych algorytmów w liniowej głębokości i przy liniowej łączności. Implementacja członu $e^{i \mathcal{J}_k}$ ansatzu UCJ wymaga łączności każdy-z-każdym lub użycia sieci fermionic swap, co jest trudne do zrealizowania na zaszumionych, przed-odpornych na błędy procesorach kwantowych o ograniczonej łączności. Idea lokalnego ansatzu UCJ polega na nałożeniu ograniczeń rzadkości na macierze $\mathbf{J}^{\alpha\alpha}$ i $\mathbf{J}^{\alpha\beta}$, które pozwalają na ich implementację w stałej głębokości na topologiach qubitów o ograniczonej łączności. Ograniczenia są określone przez listę indeksów wskazujących, które wpisy macierzy w górnym trójkącie mogą być niezerowe (ponieważ macierze są symetryczne, wystarczy podać tylko górny trójkąt). Indeksy te można interpretować jako pary orbitali, którym wolno oddziaływać.

Rozważmy jako przykład topologię qubitów w postaci sieci kwadratowej. Możemy umieścić orbitale $\alpha$ i $\beta$ w równoległych liniach na siatce, a połączenia między tymi liniami tworzą „szczeble" drabiny, jak pokazano poniżej:

Przy takim układzie orbitale o tym samym spinie są połączone topologią liniową, a orbitale o różnych spinach są połączone, gdy należą do tego samego orbitalu przestrzennego. Daje to następujące ograniczenia indeksów na macierze $\mathbf{J}$:

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

Innymi słowy, jeżeli macierze $\mathbf{J}$ są niezerowe tylko w określonych indeksach górnego trójkąta, to człon $e^{i \mathcal{J}_k}$ można zaimplementować na topologii kwadratowej bez użycia jakichkolwiek bramek swap, w stałej głębokości. Oczywiście nałożenie takich ograniczeń na ansatz zmniejsza jego wyrazistość, przez co może być wymagana większa liczba powtórzeń ansatzu.

Sprzęt IBM® posiada topologię qubitów w postaci sieci heavy-hex; w tym przypadku możemy przyjąć wzorzec „zygzak", przedstawiony poniżej. W tym wzorcu orbitale o tym samym spinie są odwzorowywane na qubitach z topologią liniową (czerwone i niebieskie kółka), a połączenie między orbitalami różnych spinów pojawia się co 4. orbital przestrzenny, przy czym połączenie jest realizowane przez qubit ancilla (fioletowe kółka). W tym przypadku ograniczenia indeksów wynoszą:

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

Samospójne odtwarzanie konfiguracji

Procedura samospójnego odtwarzania konfiguracji jest zaprojektowana tak, aby wydobyć jak najwięcej sygnału z zaszumionych próbek kwantowych. Ponieważ Hamiltonian molekularny zachowuje liczbę cząstek i składową Z spinu, sensowne jest wybranie Circuit ansatz, który również zachowuje te symetrie. Kiedy zastosujemy go do stanu Hartree-Focka, wynikowy stan będzie miał ustaloną liczbę cząstek i składową Z spinu w warunkach bezszumowych. Dlatego obie połówki (spin-$\alpha$ i spin-$\beta$) każdego ciągu bitów próbkowanego z tego stanu powinny mieć taką samą wagę Hamminga jak stan Hartree-Focka. Ze względu na obecność szumu w obecnych procesorach kwantowych niektóre zmierzone ciągi bitów naruszają tę właściwość. Prosta postać postelekcji odrzucałaby te ciągi bitów, ale jest to marnotrawstwo, ponieważ mogą one nadal zawierać pewien sygnał. Procedura samospójnego odtwarzania próbuje odzyskać część tego sygnału w przetwarzaniu końcowym. Procedura jest iteracyjna i wymaga jako danych wejściowych oszacowania średnich obsadzeń każdego orbitalu w stanie podstawowym, które jest najpierw obliczane na podstawie surowych próbek. Procedura jest uruchamiana w pętli, a każda iteracja składa się z następujących kroków:

1. Dla każdego ciągu bitów naruszającego określone symetrie odwróć jego bity za pomocą probabilistycznej procedury zaprojektowanej tak, aby przybliżyć ciąg bitów do bieżącego oszacowania średnich obsadzeń orbitalnych, uzyskując nowy ciąg bitów.
2. Zbierz wszystkie stare i nowe ciągi bitów spełniające symetrie i pobierz z nich podzbiory o ustalonym z góry rozmiarze.
3. Dla każdego podzbioru ciągów bitów zrzutuj Hamiltonian na podprzestrzeń rozpinaną przez odpowiadające wektory bazowe (opis tych wektorów bazowych znajdziesz w poprzedniej sekcji) i oblicz klasycznym komputerem przybliżenie stanu podstawowego zrzutowanego Hamiltonianu.
4. Zaktualizuj oszacowanie średnich obsadzeń orbitalnych na podstawie przybliżenia stanu podstawowego o najniższej energii.

Diagram przepływu pracy SQD

Przepływ pracy SQD przedstawiono na poniższym diagramie:

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane:

• Qiskit SDK w wersji 1.0 lub nowszej z obsługą wizualizacji
• Qiskit Runtime w wersji 0.22 lub nowszej (pip install qiskit-ibm-runtime)
• Dodatek SQD do Qiskit w wersji 0.11 lub nowszej (pip install qiskit-addon-sqd)
• ffsim w wersji 0.0.58 lub nowszej (pip install ffsim)

Konfiguracja

# 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: Odwzorowanie wejść klasycznych na problem kwantowy

W tym samouczku znajdziemy przybliżenie stanu podstawowego cząsteczki przy równowadze w bazie cc-pVDZ. Najpierw określamy cząsteczkę i jej właściwości.

# 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

Przed skonstruowaniem Circuit ansatzu LUCJ najpierw przeprowadzamy obliczenie CCSD w następującej komórce kodu. Amplitudy $t_1$ i $t_2$ z tego obliczenia zostaną użyte do inicjalizacji parametrów 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

Teraz używamy ffsim do stworzenia Circuit ansatzu. Ponieważ nasza cząsteczka ma stan Hartree-Focka z zamkniętą powłoką, używamy zbilansowanej spinowo odmiany ansatzu UCJ: UCJOpSpinBalanced. Przekazujemy pary oddziaływań odpowiednie dla topologii qubitów w postaci sieci heavy-hex (patrz sekcja tła dotycząca ansatzu LUCJ). Ustawiamy optimize=True w metodzie from_t_amplitudes, aby włączyć „skompresowaną" podwójną faktoryzację amplitud $t_2$ (szczegóły znajdziesz w The local unitary cluster Jastrow (LUCJ) ansatz w dokumentacji ffsim).

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: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

Następnie optymalizujemy Circuit pod kątem docelowego sprzętu.

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

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

Using backend ibm_fez

Zalecamy następujące kroki, aby zoptymalizować ansatz i uczynić go kompatybilnym ze sprzętem.

• Wybierz fizyczne Qubit (initial_layout) z docelowego sprzętu, które są zgodne ze wzorcem „zig-zag" opisanym wcześniej. Rozmieszczenie Qubit w tym wzorcu prowadzi do wydajnego, kompatybilnego ze sprzętem Circuit z mniejszą liczbą Gate. Poniżej zamieszczono kod automatyzujący dobór wzorca „zig-zag", który używa heurystyki punktacji w celu minimalizacji błędów związanych z wybranym układem.
• Wygeneruj etapowy menedżer przebiegów (pass manager) przy użyciu funkcji generate_preset_pass_manager z Qiskit, podając wybrany backend i initial_layout.
• Ustaw etap pre_init swojego etapowego menedżera przebiegów na ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT zawiera przejścia Transpiler Qiskit, które rozkładają Gate na rotacje orbitalne, a następnie łączą je ze sobą, co skutkuje mniejszą liczbą Gate w końcowym Circuit.
• Uruchom menedżer przebiegów na swoim Circuit.

Kod do automatycznego doboru układu „zig-zag"

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

Krok 3: Wykonanie za pomocą prymitywów Qiskit

Po zoptymalizowaniu obwodu do wykonania na sprzęcie jesteśmy gotowi uruchomić go na docelowym urządzeniu i zebrać próbki do estymacji energii stanu podstawowego. Ponieważ mamy tylko jeden obwód, użyjemy trybu wykonania zadań Qiskit Runtime i wykonamy nasz Circuit.

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

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

Krok 4: Przetwarzanie końcowe i zwracanie wyniku w żądanym formacie klasycznym

Teraz szacujemy energię stanu podstawowego hamiltonianu za pomocą funkcji diagonalize_fermionic_hamiltonian. Funkcja ta wykonuje procedurę samospójnego odtwarzania konfiguracji w celu iteracyjnego udoskonalania zaszumionych próbek kwantowych i poprawy estymacji energii. Przekazujemy funkcję zwrotną (callback), aby móc zapisywać wyniki pośrednie do późniejszej analizy. Wyjaśnienia argumentów funkcji diagonalize_fermionic_hamiltonian znajdziesz w dokumentacji API.

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

Wizualizacja wyników

Pierwszy wykres pokazuje, że po kilku iteracjach szacujemy energię stanu podstawowego z dokładnością do ~41 mH (dokładność chemiczna jest zazwyczaj przyjmowana jako 1 kcal/mol $\approx$ 1.6 mH). Energię można poprawić, zezwalając na więcej iteracji odtwarzania konfiguracji lub zwiększając liczbę próbek w każdej partii.

Drugi wykres pokazuje średnie obsadzenie każdego orbitalnego przestrzennego po ostatniej iteracji. Widać, że zarówno elektrony ze spinem w górę, jak i w dół obsadzają pierwsze pięć orbitali z dużym prawdopodobieństwem w naszych rozwiązaniach.

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

Ankieta o samouczku

Wypełnij tę krótką ankietę, aby podzielić się opinią na temat tego samouczka. Twoje uwagi pomogą nam ulepszyć treści i doświadczenie użytkownika.

Link do ankiety

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.