Diagonalizarea cuantică bazată pe eșantioane a unui Hamiltonian de chimie

Estimare de utilizare: sub un minut pe un procesor Heron r2 (NOTĂ: Aceasta este doar o estimare. Durata de rulare poate varia.)

Context

În acest tutorial, îți arătăm cum să post-procesezi eșantioane cuantice cu zgomot pentru a găsi o aproximare a stării fundamentale a moleculei de azot $\text{N}_2$ la lungimea de echilibru a legăturii, folosind algoritmul de diagonalizare cuantică bazată pe eșantioane (SQD), pentru eșantioane preluate dintr-un circuit cuantic de 59 de qubiți (52 qubiți de sistem + 7 qubiți ancilari). O implementare bazată pe Qiskit este disponibilă în addon-urile Qiskit pentru SQD; mai multe detalii se găsesc în documentația corespunzătoare, cu un exemplu simplu pentru a începe.

SQD este o tehnică de găsire a valorilor proprii și a vectorilor proprii ai operatorilor cuantici, cum ar fi Hamiltonianul unui sistem cuantic, folosind împreună calculul cuantic și calculul clasic distribuit. Calculul clasic distribuit este folosit pentru a procesa eșantioane obținute de la un procesor cuantic și pentru a proiecta și diagonaliza un Hamiltonian țintă într-un subspațiu generat de acestea. Aceasta permite SQD să fie robust față de eșantioanele afectate de zgomot cuantic și să gestioneze Hamiltonieni mari, cum ar fi Hamiltonenii de chimie cu milioane de termeni de interacțiune, dincolo de capacitățile oricărei metode exacte de diagonalizare. Un flux de lucru bazat pe SQD are următorii pași:

1. Alege un ansatz de circuit și aplică-l pe un calculator cuantic la o stare de referință (în acest caz, starea Hartree-Fock).
2. Eșantionează șiruri de biți din starea cuantică rezultată.
3. Rulează procedura de recuperare auto-consistentă a configurației pe șirurile de biți pentru a obține aproximarea stării fundamentale.

SQD funcționează bine atunci când starea proprie țintă este rară: funcția de undă este susținută de un set de stări de bază $\mathcal{S} = \{|x\rangle \}$ ale căror dimensiuni nu cresc exponențial cu dimensiunea problemei.

Chimie cuantică

Proprietățile moleculelor sunt determinate în mare parte de structura electronilor din interiorul lor. Ca particule fermionice, electronii pot fi descriși folosind un formalism matematic numit a doua cuantificare. Ideea este că există un număr de orbitale, fiecare putând fi fie goală, fie ocupată de un fermion. Un sistem de $N$ orbitale este descris de un set de operatori fermionici de anihilare $\{\hat{a}_p\}_{p=1}^N$ care satisfac relațiile fermionice de anticomutare,

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

Adjunctul $\hat{a}_p^\dagger$ se numește operator de creare.

Până acum, expunerea noastră nu a ținut cont de spin, care este o proprietate fundamentală a fermionilor. Când se ia în considerare spinul, orbitalele apar în perechi numite orbitale spațiale. Fiecare orbital spațial este compus din două orbitale de spin, unul etichetat „spin-$\alpha$" și unul etichetat „spin-$\beta$". Scriem $\hat{a}_{p\sigma}$ pentru operatorul de anihilare asociat cu orbitalul de spin cu spinul $\sigma$ ($\sigma \in \{\alpha, \beta\}$) din orbitalul spațial $p$. Dacă luăm $N$ ca numărul de orbitale spațiale, atunci există în total $2N$ orbitale de spin. Spațiul Hilbert al acestui sistem este generat de $2^{2N}$ vectori de bază ortonormali etichetați cu șiruri de biți în două părți $\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$.

Hamiltonianul unui sistem molecular poate fi scris ca

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

unde $h_{pr}$ și $h_{prqs}$ sunt numere complexe numite integrale moleculare, care pot fi calculate din specificația moleculei folosind un program de calculator. În acest tutorial, calculăm integralele folosind pachetul software PySCF.

Pentru detalii despre derivarea Hamiltonianului molecular, consultă un manual de chimie cuantică (de exemplu, Modern Quantum Chemistry de Szabo și Ostlund). Pentru o explicație de nivel înalt a modului în care problemele de chimie cuantică sunt mapate pe calculatoare cuantice, urmărește prelegerea Mapping Problems to Qubits de la Qiskit Global Summer School 2024.

Ansatzul local unitar cluster Jastrow (LUCJ)

SQD necesită un ansatz de circuit cuantic pentru a genera eșantioane. În acest tutorial, vom folosi ansatzul local unitar cluster Jastrow (LUCJ) datorită combinației dintre motivația fizică și compatibilitatea cu hardware-ul.

Ansatzul LUCJ este o formă specializată a ansatzului general unitar cluster Jastrow (UCJ), care are forma

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

unde $\lvert \Phi_0 \rangle$ este o stare de referință, adesea aleasă a fi starea Hartree-Fock, iar $\hat{K}_\mu$ și $\hat{J}_\mu$ au forma

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

unde am definit operatorul de număr

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

Operatorul $e^{\hat{K}_\mu}$ este o rotație orbitală, care poate fi implementată folosind algoritmi cunoscuți în adâncime liniară și cu conectivitate liniară. Implementarea termenului $e^{i \mathcal{J}_k}$ al ansatzului UCJ necesită fie conectivitate totală, fie utilizarea unei rețele de swap fermionic, ceea ce o face dificilă pentru procesoarele cuantice pre-tolerante la erori cu conectivitate limitată. Ideea ansatzului UCJ local este de a impune constrângeri de raritate asupra matricelor $\mathbf{J}^{\alpha\alpha}$ și $\mathbf{J}^{\alpha\beta}$, care permit implementarea lor în adâncime constantă pe topologii de qubiți cu conectivitate limitată. Constrângerile sunt specificate printr-o listă de indici care indică ce intrări din triunghiul superior al matricei pot fi nenule (deoarece matricele sunt simetrice, este necesar să se specifice doar triunghiul superior). Acești indici pot fi interpretați ca perechi de orbitale care pot interacționa.

Ca exemplu, consideră o topologie de qubiți cu rețea pătrată. Putem plasa orbitalele $\alpha$ și $\beta$ pe linii paralele ale rețelei, cu conexiuni între aceste linii formând „treptele" unei scări, astfel:

Cu această configurare, orbitalele cu același spin sunt conectate printr-o topologie liniară, în timp ce orbitalele cu spinuri diferite sunt conectate când împărtășesc același orbital spațial. Aceasta produce următoarele constrângeri de indici asupra matricelor $\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*}$$

Cu alte cuvinte, dacă matricele $\mathbf{J}$ sunt nenule doar la indicii specificați în triunghiul superior, atunci termenul $e^{i \mathcal{J}_k}$ poate fi implementat pe o topologie pătrată fără a folosi porți swap, în adâncime constantă. Desigur, impunerea unor astfel de constrângeri asupra ansatzului îl face mai puțin expresiv, astfel că pot fi necesare mai multe repetări ale ansatzului.

Hardware-ul IBM® are o topologie de qubiți cu rețea hex grea, caz în care putem adopta un model „zig-zag", reprezentat mai jos. În acest model, orbitalele cu același spin sunt mapate pe qubiți cu o topologie liniară (cercuri roșii și albastre), iar o conexiune între orbitale cu spin diferit există la fiecare al 4-lea orbital spațial, conexiunea fiind facilitată de un qubit ancilar (cercuri violet). În acest caz, constrângerile de indici sunt

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

Recuperarea auto-consistentă a configurației

Procedura de recuperare auto-consistentă a configurației este concepută pentru a extrage cât mai mult semnal posibil din eșantioanele cuantice cu zgomot. Deoarece Hamiltonianul molecular conservă numărul de particule și spinul Z, are sens să alegi un ansatz de circuit care conservă și aceste simetrii. Când este aplicat la starea Hartree-Fock, starea rezultată are un număr fix de particule și spin Z în cazul fără zgomot. Prin urmare, jumătățile spin-$\alpha$ și spin-$\beta$ ale oricărui șir de biți eșantionat din această stare ar trebui să aibă aceeași greutate Hamming ca în starea Hartree-Fock. Din cauza zgomotului din procesoarele cuantice actuale, unele șiruri de biți măsurate vor încălca această proprietate. O formă simplă de postselecție ar elimina aceste șiruri de biți, dar aceasta este risipitoare deoarece acele șiruri de biți ar putea conține totuși un anumit semnal. Procedura de recuperare auto-consistentă încearcă să recupereze o parte din acel semnal în post-procesare. Procedura este iterativă și necesită ca intrare o estimare a ocupanțelor medii ale fiecărui orbital în starea fundamentală, care este calculată mai întâi din eșantioanele brute. Procedura rulează în buclă, iar fiecare iterație are următorii pași:

1. Pentru fiecare șir de biți care încalcă simetriile specificate, inversează biții printr-o procedură probabilistică concepută pentru a aduce șirul de biți mai aproape de estimarea curentă a ocupanțelor medii ale orbitalelor, pentru a obține un nou șir de biți.
2. Colectează toate șirurile de biți vechi și noi care satisfac simetriile și eșantionează subseturi de dimensiune fixă, aleasă în prealabil.
3. Pentru fiecare subset de șiruri de biți, proiectează Hamiltonianul în subspațiul generat de vectorii de bază corespunzători (consulta secțiunea anterioară pentru o descriere a acestor vectori de bază) și calculează o estimare a stării fundamentale a Hamiltonianului proiectat pe un calculator clasic.
4. Actualizează estimarea ocupanțelor medii ale orbitalelor cu estimarea stării fundamentale cu cea mai mică energie.

Diagrama fluxului de lucru SQD

Fluxul de lucru SQD este reprezentat în diagrama următoare:

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalat următoarele:

• Qiskit SDK v1.0 sau o versiune ulterioară, cu suport pentru vizualizare
• Qiskit Runtime v0.22 sau o versiune ulterioară (pip install qiskit-ibm-runtime)
• Addon-ul SQD pentru Qiskit v0.11 sau o versiune ulterioară (pip install qiskit-addon-sqd)
• ffsim v0.0.58 sau o versiune ulterioară (pip install ffsim)

Configurare

# 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

Pasul 1: Maparea intrărilor clasice la o problemă cuantică

În acest tutorial, vom găsi o aproximare a stării fundamentale a moleculei la echilibru în setul de baze cc-pVDZ. Mai întâi, specificăm molecula și proprietățile sale.

# 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

Înainte de a construi Circuit-ul ansatz LUCJ, efectuăm mai întâi un calcul CCSD în celula de cod următoare. Amplitudinile $t_1$ și $t_2$ din acest calcul vor fi folosite pentru a inițializa parametrii ansatzului.

# 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

Acum, folosim ffsim pentru a crea Circuit-ul ansatz. Deoarece molecula noastră are o stare Hartree-Fock cu coajă închisă, folosim varianta echilibrată de spin a ansatzului UCJ, UCJOpSpinBalanced. Transmitem perechi de interacțiune potrivite pentru o topologie de qubiți cu rețea hex grea (consulta secțiunea de fundal despre ansatzul LUCJ). Setăm optimize=True în metoda from_t_amplitudes pentru a activa factorizarea dublă „comprimată" a amplitudinilor $t_2$ (consulta The local unitary cluster Jastrow (LUCJ) ansatz din documentația ffsim pentru detalii).

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

Pasul 2: Optimizează problema pentru execuția pe hardware quantum

În continuare, optimizăm circuitul pentru un hardware țintă.

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

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

Using backend ibm_fez

Recomandăm următorii pași pentru a optimiza ansatz-ul și a-l face compatibil cu hardware-ul.

• Selectează qubiți fizici (initial_layout) de pe hardware-ul țintă care respectă modelul „zig-zag" descris anterior. Dispunerea qubiților în acest model duce la un circuit eficient, compatibil cu hardware-ul, cu mai puține gate-uri. Aici includem cod pentru automatizarea selecției modelului „zig-zag", care folosește o euristică de punctaj pentru a minimiza erorile asociate layout-ului selectat.
• Generează un pass manager etapizat folosind funcția generate_preset_pass_manager din Qiskit, cu backend-ul și initial_layout-ul alese de tine.
• Setează etapa pre_init a pass manager-ului etapizat la ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT include pase de Transpiler Qiskit care descompun gate-urile în rotații orbitale și apoi le combină, rezultând mai puține gate-uri în circuitul final.
• Rulează pass manager-ul pe circuitul tău.

Cod pentru selecția automată a layout-ului „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})

Pasul 3: Execuție folosind primitivele Qiskit

După optimizarea circuitului pentru execuția pe hardware, suntem pregătiți să îl rulăm pe hardware-ul țintă și să colectăm mostre pentru estimarea energiei stării fundamentale. Deoarece avem un singur circuit, vom folosi modul de execuție Job din Qiskit Runtime și vom executa circuitul nostru.

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

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

Pasul 4: Post-procesare și returnarea rezultatului în formatul clasic dorit

Acum estimăm energia stării fundamentale a Hamiltonianului folosind funcția diagonalize_fermionic_hamiltonian. Această funcție realizează procedura de recuperare auto-consistentă a configurației pentru a rafina iterativ mostrele cuantice zgomotoase și a îmbunătăți estimarea energiei. Transmitem o funcție callback pentru a putea salva rezultatele intermediare în vederea analizei ulterioare. Consultă documentația API pentru explicații privind argumentele funcției 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

Vizualizarea rezultatelor

Primul grafic arată că după câteva iterații estimăm energia stării fundamentale cu o eroare de aproximativ ~41 mH (precizia chimică este acceptată în general ca fiind 1 kcal/mol $\approx$ 1.6 mH). Energia poate fi îmbunătățită permițând mai multe iterații de recuperare a configurației sau mărind numărul de mostre per lot.

Al doilea grafic arată ocupanța medie a fiecărei orbite spațiale după iterația finală. Putem observa că atât electronii cu spin-up, cât și cei cu spin-down ocupă primele cinci orbite cu probabilitate ridicată în soluțiile noastre.

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

Sondaj tutorial

Te rog să completezi acest scurt sondaj pentru a oferi feedback despre acest tutorial. Opiniile tale ne vor ajuta să îmbunătățim conținutul și experiența utilizatorilor.

Link către sondaj

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.