Vai al contenuto principale

Diagonalizzazione quantistica basata su campioni di un Hamiltoniano chimico

Stima del tempo di utilizzo: meno di un minuto su un processore Heron r2 (NOTA: Questa è solo una stima. Il vostro tempo di esecuzione potrebbe variare.)

Contesto​

In questo tutorial, mostriamo come post-elaborare campioni quantistici rumorosi per trovare un'approssimazione allo stato fondamentale della molecola di azoto N2\text{N}_2 alla lunghezza di legame di equilibrio, utilizzando l'algoritmo di diagonalizzazione quantistica basata su campioni (SQD), per campioni prelevati da un circuito quantistico a 59 qubit (52 qubit di sistema + 7 qubit ancilla). Un'implementazione basata su Qiskit è disponibile negli SQD Qiskit addons, maggiori dettagli si possono trovare nella documentazione corrispondente con un semplice esempio per iniziare. SQD è una tecnica per trovare autovalori e autovettori di operatori quantistici, come l'Hamiltoniano di un sistema quantistico, utilizzando insieme calcolo quantistico e classico distribuito. Il calcolo classico distribuito viene utilizzato per elaborare i campioni ottenuti da un processore quantistico e per proiettare e diagonalizzare un Hamiltoniano target in un sottospazio da essi generato. Questo permette a SQD di essere robusto rispetto ai campioni corrotti dal rumore quantistico e di gestire Hamiltoniani di grandi dimensioni, come Hamiltoniani chimici con milioni di termini di interazione, oltre la portata di qualsiasi metodo di diagonalizzazione esatta. Un flusso di lavoro basato su SQD prevede i seguenti passaggi:

  1. Scegliete un ansatz di circuito e applicatelo su un computer quantistico a uno stato di riferimento (in questo caso, lo stato di Hartree-Fock).
  2. Campionate bitstring dallo stato quantistico risultante.
  3. Eseguite la procedura di recupero di configurazione auto-consistente sui bitstring per ottenere l'approssimazione allo stato fondamentale.

SQD è noto per funzionare bene quando l'autostato target è sparso: la funzione d'onda è supportata in un insieme di stati di base S={∣x⟩}\mathcal{S} = \{|x\rangle \} la cui dimensione non aumenta esponenzialmente con la dimensione del problema.

Chimica quantistica​

Le proprietà delle molecole sono in gran parte determinate dalla struttura degli elettroni al loro interno. Come particelle fermioniche, gli elettroni possono essere descritti utilizzando un formalismo matematico chiamato seconda quantizzazione. L'idea è che ci siano un certo numero di orbitali, ciascuno dei quali può essere vuoto o occupato da un fermione. Un sistema di NN orbitali è descritto da un insieme di operatori di annichilazione fermionici {a^p}p=1N\{\hat{a}_p\}_{p=1}^N che soddisfano le relazioni di anticommutazione fermioniche,

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

L'aggiunto a^p†\hat{a}_p^\dagger è chiamato operatore di creazione.

Finora, la nostra esposizione non ha tenuto conto dello spin, che è una proprietà fondamentale dei fermioni. Quando si tiene conto dello spin, gli orbitali vengono in coppie chiamate orbitali spaziali. Ciascun orbitale spaziale è composto da due orbitali di spin, uno etichettato come "spin-α\alpha" e uno etichettato come "spin-β\beta". Scriviamo quindi a^pσ\hat{a}_{p\sigma} per l'operatore di annichilazione associato all'orbitale di spin con spin σ\sigma (σ∈{α,β}\sigma \in \{\alpha, \beta\}) nell'orbitale spaziale pp. Se consideriamo NN come il numero di orbitali spaziali, allora ci sono in totale 2N2N orbitali di spin. Lo spazio di Hilbert di questo sistema è generato da 22N2^{2N} vettori di base ortonormali etichettati con bitstring a due parti ∣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.

L'Hamiltoniano di un sistema molecolare può essere scritto come

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

dove gli hprh_{pr} e gli hprqsh_{prqs} sono numeri complessi chiamati integrali molecolari che possono essere calcolati dalla specifica della molecola utilizzando un programma informatico. In questo tutorial, calcoliamo gli integrali utilizzando il pacchetto software PySCF.

Per dettagli su come viene derivato l'Hamiltoniano molecolare, consultate un libro di testo di chimica quantistica (ad esempio, Modern Quantum Chemistry di Szabo e Ostlund). Per una spiegazione ad alto livello di come i problemi di chimica quantistica vengono mappati su computer quantistici, consultate la lezione Mapping Problems to Qubits della Qiskit Global Summer School 2024.

Ansatz local unitary cluster Jastrow (LUCJ)​

SQD richiede un ansatz di circuito quantistico da cui estrarre campioni. In questo tutorial, useremo l'ansatz local unitary cluster Jastrow (LUCJ) per la sua combinazione di motivazione fisica e compatibilità con l'hardware.

L'ansatz LUCJ è una forma specializzata dell'ansatz generale unitary cluster Jastrow (UCJ), che ha la forma

∣Ψ⟩=∏μ=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

dove ∣Φ0⟩\lvert \Phi_0 \rangle è uno stato di riferimento, spesso preso come stato di Hartree-Fock, e gli K^μ\hat{K}_\mu e J^μ\hat{J}_\mu hanno la forma

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

dove abbiamo definito l'operatore numero

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

L'operatore eK^μe^{\hat{K}_\mu} è una rotazione orbitale, che può essere implementata utilizzando algoritmi noti in profondità lineare e utilizzando connettività lineare. L'implementazione del termine eiJke^{i \mathcal{J}_k} dell'ansatz UCJ richiede connettività all-to-all o l'uso di una rete di swap fermionico, rendendola impegnativa per i processori quantistici rumorosi pre-fault-tolerant che hanno connettività limitata. L'idea dell'ansatz UCJ locale è imporre vincoli di sparsità sulle matrici Jαα\mathbf{J}^{\alpha\alpha} e Jαβ\mathbf{J}^{\alpha\beta} che permettono loro di essere implementate in profondità costante su topologie di qubit con connettività limitata. I vincoli sono specificati da un elenco di indici che indica quali elementi della matrice nel triangolo superiore possono essere diversi da zero (poiché le matrici sono simmetriche, è necessario specificare solo il triangolo superiore). Questi indici possono essere interpretati come coppie di orbitali che possono interagire.

Come esempio, considerate una topologia di qubit a reticolo quadrato. Possiamo posizionare gli orbitali α\alpha e β\beta in linee parallele sul reticolo, con connessioni tra queste linee che formano "pioli" di una forma a scala, in questo modo:

Qubit mapping diagram for the LUCJ ansatz on a square lattice

Con questa configurazione, gli orbitali con lo stesso spin sono connessi con una topologia lineare, mentre gli orbitali con spin diversi sono connessi quando condividono lo stesso orbitale spaziale. Questo produce i seguenti vincoli di indice sulle matrici 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*}

In altre parole, se le matrici J\mathbf{J} sono diverse da zero solo agli indici specificati nel triangolo superiore, allora il termine eiJke^{i \mathcal{J}_k} può essere implementato su una topologia quadrata senza utilizzare alcun gate di swap, in profondità costante. Naturalmente, imporre tali vincoli sull'ansatz lo rende meno espressivo, quindi potrebbero essere necessarie più ripetizioni dell'ansatz.

L'hardware IBM® ha una topologia di qubit a reticolo heavy-hex, nel qual caso possiamo adottare un pattern "zig-zag", raffigurato qui sotto. In questo pattern, gli orbitali con lo stesso spin sono mappati su qubit con topologia lineare (cerchi rossi e blu), e una connessione tra orbitali di spin diverso è presente ogni 4° orbitale spaziale, con la connessione facilitata da un qubit ancilla (cerchi viola). In questo caso, i vincoli di indice sono

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

Qubit mapping diagram for the LUCJ ansatz on a heavy-hex lattice

Recupero di configurazione auto-consistente​

La procedura di recupero di configurazione auto-consistente è progettata per estrarre il maggior segnale possibile da campioni quantistici rumorosi. Poiché l'Hamiltoniano molecolare conserva il numero di particelle e lo spin Z, ha senso scegliere un ansatz di circuito che conservi anch'esso queste simmetrie. Quando applicato allo stato di Hartree-Fock, lo stato risultante ha un numero di particelle e uno spin Z fissi nell'impostazione senza rumore. Pertanto, le metà con spin-α\alpha e spin-β\beta di qualsiasi bitstring campionato da questo stato dovrebbero avere lo stesso peso di Hamming dello stato di Hartree-Fock. A causa della presenza di rumore negli attuali processori quantistici, alcuni bitstring misurati violeranno questa proprietà. Una semplice forma di postselezione scarterebbe questi bitstring, ma questo è uno spreco perché quei bitstring potrebbero ancora contenere del segnale. La procedura di recupero auto-consistente tenta di recuperare parte di quel segnale nella post-elaborazione. La procedura è iterativa e richiede come input una stima delle occupazioni medie di ciascun orbitale nello stato fondamentale, che viene prima calcolata dai campioni grezzi. La procedura viene eseguita in un ciclo, e ogni iterazione ha i seguenti passaggi:

  1. Per ciascun bitstring che viola le simmetrie specificate, capovolgete i suoi bit con una procedura probabilistica progettata per avvicinare il bitstring alla stima attuale delle occupazioni orbitali medie, per ottenere un nuovo bitstring.
  2. Raccogliete tutti i vecchi e nuovi bitstring che soddisfano le simmetrie, e sottocampionate sottoinsiemi di dimensione fissa, scelta in anticipo.
  3. Per ciascun sottoinsieme di bitstring, proiettate l'Hamiltoniano nel sottospazio generato dai vettori di base corrispondenti (vedere la sezione precedente per una descrizione di questi vettori di base), e calcolate una stima dello stato fondamentale dell'Hamiltoniano proiettato su un computer classico.
  4. Aggiornate la stima delle occupazioni orbitali medie con la stima dello stato fondamentale con l'energia più bassa.

Diagramma del flusso di lavoro SQD​

Il flusso di lavoro SQD è raffigurato nel seguente diagramma:

Workflow diagram of the SQD algorithm

Requisiti​

Prima di iniziare questo tutorial, assicuratevi di avere quanto segue installato:

  • Qiskit SDK v1.0 o successivo, con supporto per la visualizzazione
  • Qiskit Runtime v0.22 o successivo (pip install qiskit-ibm-runtime)
  • SQD Qiskit addon v0.11 o successivo (pip install qiskit-addon-sqd)
  • ffsim v0.0.58 o successivo (pip install ffsim)

Configurazione​

# 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

Passo 1: Mappare gli input classici su un problema quantistico​

In questo tutorial, troveremo un'approssimazione allo stato fondamentale della molecola all'equilibrio nel set di base cc-pVDZ. Prima, specifichiamo la molecola e le sue proprietà.

# 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

Prima di costruire il circuito ansatz LUCJ, eseguiamo prima un calcolo CCSD nella seguente cella di codice. Le ampiezze t1t_1 e t2t_2 da questo calcolo verranno utilizzate per inizializzare i parametri dell'ansatz.

# 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

Ora, utilizziamo ffsim per creare il circuito ansatz. Poiché la nostra molecola ha uno stato di Hartree-Fock a shell chiuso, utilizziamo la variante bilanciata in spin dell'ansatz UCJ, UCJOpSpinBalanced. Passiamo coppie di interazione appropriate per una topologia di qubit a reticolo heavy-hex (vedere la sezione di background sull'ansatz LUCJ). Impostiamo optimize=True nel metodo from_t_amplitudes per abilitare la "compressa" doppia fattorizzazione delle ampiezze t2t_2 (vedere The local unitary cluster Jastrow (LUCJ) ansatz nella documentazione di ffsim per i dettagli).

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

Passo 2: Ottimizzare il problema per l'esecuzione sull'hardware quantistico​

Successivamente, ottimizziamo il circuito per un hardware target.

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

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

Raccomandiamo i seguenti passaggi per ottimizzare l'ansatz e renderlo compatibile con l'hardware.

  • Selezionate qubit fisici (initial_layout) dall'hardware target che aderisce al pattern "zig-zag" descritto in precedenza. Disporre i qubit in questo pattern porta a un circuito efficiente e compatibile con l'hardware con meno gate. Qui, includiamo del codice per automatizzare la selezione del pattern "zig-zag" che utilizza un'euristica di punteggio per minimizzare gli errori associati al layout selezionato.
  • Generate un pass manager a stadi utilizzando la funzione generate_preset_pass_manager da qiskit con la vostra scelta di backend e initial_layout.
  • Impostate lo stadio pre_init del vostro pass manager a stadi su ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT include pass del transpiler di qiskit che decompongono i gate in rotazioni orbitali e poi fondono le rotazioni orbitali, risultando in meno gate nel circuito finale.
  • Eseguite il pass manager sul vostro circuito.
Codice per la selezione automatizzata del layout "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})

Passo 3: Eseguire utilizzando le primitive Qiskit​

Dopo aver ottimizzato il circuito per l'esecuzione hardware, siamo pronti ad eseguirlo sull'hardware di destinazione e raccogliere campioni per la stima dell'energia dello stato fondamentale. Poiché abbiamo un solo circuito, utilizzeremo la modalità di esecuzione Job di Qiskit Runtime ed eseguiremo il nostro circuito.

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

Passo 4: Post-elaborare e restituire il risultato nel formato classico desiderato​

Ora stimiamo l'energia dello stato fondamentale dell'Hamiltoniana utilizzando la funzione diagonalize_fermionic_hamiltonian. Questa funzione esegue la procedura di recupero della configurazione autoconsistente per raffinare iterativamente i campioni quantistici rumorosi e migliorare la stima dell'energia. Passiamo una funzione di callback in modo da poter salvare i risultati intermedi per un'analisi successiva. Consultate la documentazione API per spiegazioni degli argomenti di 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

Visualizzare i risultati​

Il primo grafico mostra che dopo un paio di iterazioni stimiamo l'energia dello stato fondamentale entro ~41 mH (l'accuratezza chimica è tipicamente accettata come 1 kcal/mol ≈\approx 1.6 mH). L'energia può essere migliorata consentendo più iterazioni di recupero della configurazione o aumentando il numero di campioni per batch.

Il secondo grafico mostra l'occupazione media di ciascun orbitale spaziale dopo l'iterazione finale. Possiamo vedere che sia gli elettroni con spin-up che quelli con spin-down occupano i primi cinque orbitali con alta probabilità nelle nostre soluzioni.

# 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

Sondaggio sul tutorial​

Vi preghiamo di rispondere a questo breve sondaggio per fornire feedback su questo tutorial. Le vostre opinioni ci aiuteranno a migliorare la nostra offerta di contenuti e l'esperienza utente.

Link to survey