Diagonalisation quantique basée sur des échantillons d'un Hamiltonien de chimie

Estimation d'utilisation : moins d'une minute sur un processeur Heron r2 (REMARQUE : ceci est uniquement une estimation. Votre temps d'exécution peut varier.)

Contexte

Dans ce tutoriel, nous montrons comment post-traiter des échantillons quantiques bruités pour trouver une approximation de l'état fondamental de la molécule d'azote $\text{N}_2$ à la longueur de liaison d'équilibre, en utilisant l'algorithme de diagonalisation quantique basée sur des échantillons (SQD), pour des échantillons prélevés à partir d'un circuit quantique de 59 qubits (52 qubits système + 7 qubits ancillaires). Une implémentation basée sur Qiskit est disponible dans les extensions SQD de Qiskit ; plus de détails sont disponibles dans la documentation correspondante avec un exemple simple pour démarrer. SQD est une technique permettant de trouver les valeurs propres et les vecteurs propres d'opérateurs quantiques, tels qu'un Hamiltonien de système quantique, en combinant calcul quantique et calcul classique distribué. Le calcul classique distribué est utilisé pour traiter les échantillons obtenus à partir d'un processeur quantique, et pour projeter et diagonaliser un Hamiltonien cible dans un sous-espace engendré par ces échantillons. Cela permet à SQD d'être robuste face aux échantillons corrompus par le bruit quantique et de traiter de grands Hamiltoniens, tels que les Hamiltoniens de chimie avec des millions de termes d'interaction, hors de portée de toute méthode de diagonalisation exacte. Un workflow basé sur SQD comporte les étapes suivantes :

1. Choisir un ansatz de circuit et l'appliquer sur un ordinateur quantique à un état de référence (dans ce cas, l'état de Hartree-Fock).
2. Échantillonner des chaînes de bits à partir de l'état quantique résultant.
3. Exécuter la procédure de récupération auto-cohérente des configurations sur les chaînes de bits pour obtenir l'approximation de l'état fondamental.

SQD est reconnu pour bien fonctionner lorsque l'état propre cible est creux : la fonction d'onde est supportée dans un ensemble d'états de base $\mathcal{S} = \{|x\rangle \}$ dont la taille n'augmente pas exponentiellement avec la taille du problème.

Chimie quantique

Les propriétés des molécules sont largement déterminées par la structure des électrons qu'elles contiennent. En tant que particules fermioniques, les électrons peuvent être décrits à l'aide d'un formalisme mathématique appelé seconde quantification. L'idée est qu'il existe un certain nombre d'orbitales, chacune pouvant être soit vide, soit occupée par un fermion. Un système de $N$ orbitales est décrit par un ensemble d'opérateurs d'annihilation fermioniques $\{\hat{a}_p\}_{p=1}^N$ satisfaisant les relations d'anticommutation fermioniques,

$$\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'adjoint $\hat{a}_p^\dagger$ est appelé opérateur de création.

Jusqu'ici, notre exposé n'a pas pris en compte le spin, qui est une propriété fondamentale des fermions. En tenant compte du spin, les orbitales sont regroupées par paires appelées orbitales spatiales. Chaque orbitale spatiale est composée de deux orbitales de spin, l'une étiquetée « spin-$\alpha$ » et l'autre étiquetée « spin-$\beta$ ». On écrit alors $\hat{a}_{p\sigma}$ pour l'opérateur d'annihilation associé à l'orbitale de spin avec le spin $\sigma$ ($\sigma \in \{\alpha, \beta\}$) dans l'orbitale spatiale $p$. Si l'on prend $N$ comme le nombre d'orbitales spatiales, alors il y a au total $2N$ orbitales de spin. L'espace de Hilbert de ce système est engendré par $2^{2N}$ vecteurs de base orthonormés étiquetés par des chaînes de bits en deux parties $\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$.

Le Hamiltonien d'un système moléculaire peut s'écrire sous la forme

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

où les $h_{pr}$ et $h_{prqs}$ sont des nombres complexes appelés intégrales moléculaires qui peuvent être calculés à partir de la spécification de la molécule à l'aide d'un programme informatique. Dans ce tutoriel, nous calculons les intégrales à l'aide du logiciel PySCF.

Pour plus de détails sur la dérivation du Hamiltonien moléculaire, consultez un manuel de chimie quantique (par exemple, Modern Quantum Chemistry de Szabo et Ostlund). Pour une explication de haut niveau sur la façon dont les problèmes de chimie quantique sont transposés sur des ordinateurs quantiques, consultez le cours Mapping Problems to Qubits de la Qiskit Global Summer School 2024.

Ansatz local unitary cluster Jastrow (LUCJ)

SQD nécessite un ansatz de circuit quantique pour en tirer des échantillons. Dans ce tutoriel, nous utiliserons l'ansatz local unitary cluster Jastrow (LUCJ) en raison de sa combinaison de motivation physique et d'adéquation au matériel.

L'ansatz LUCJ est une forme spécialisée de l'ansatz général unitary cluster Jastrow (UCJ), qui a la forme

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

où $\lvert \Phi_0 \rangle$ est un état de référence, souvent pris comme l'état de Hartree-Fock, et les $\hat{K}_\mu$ et $\hat{J}_\mu$ ont la forme

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

où l'on a défini l'opérateur nombre

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

L'opérateur $e^{\hat{K}_\mu}$ est une rotation orbitale, qui peut être implémentée à l'aide d'algorithmes connus en profondeur linéaire et avec une connectivité linéaire. L'implémentation du terme $e^{i \mathcal{J}_k}$ de l'ansatz UCJ nécessite soit une connectivité tous-à-tous, soit l'utilisation d'un réseau d'échange fermionique, ce qui le rend difficile pour les processeurs quantiques bruités pré-tolérants aux fautes ayant une connectivité limitée. L'idée de l'ansatz UCJ local est d'imposer des contraintes de parcimonie sur les matrices $\mathbf{J}^{\alpha\alpha}$ et $\mathbf{J}^{\alpha\beta}$ qui permettent de les implémenter en profondeur constante sur des topologies de qubits à connectivité limitée. Les contraintes sont spécifiées par une liste d'indices indiquant quelles entrées de matrice dans le triangle supérieur sont autorisées à être non nulles (puisque les matrices sont symétriques, seul le triangle supérieur doit être spécifié). Ces indices peuvent être interprétés comme des paires d'orbitales autorisées à interagir.

À titre d'exemple, considérons une topologie de qubits en réseau carré. Nous pouvons placer les orbitales $\alpha$ et $\beta$ sur des lignes parallèles du réseau, avec des connexions entre ces lignes formant les « barreaux » d'une forme en échelle, comme ceci :

Avec cette configuration, les orbitales de même spin sont connectées selon une topologie linéaire, tandis que les orbitales de spin différent sont connectées lorsqu'elles partagent la même orbitale spatiale. Cela donne les contraintes d'indices suivantes sur les matrices $\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*}$$

En d'autres termes, si les matrices $\mathbf{J}$ sont non nulles uniquement aux indices spécifiés dans le triangle supérieur, alors le terme $e^{i \mathcal{J}_k}$ peut être implémenté sur une topologie carrée sans utiliser de portes d'échange, en profondeur constante. Bien entendu, imposer de telles contraintes sur l'ansatz le rend moins expressif, de sorte que davantage de répétitions de l'ansatz peuvent être nécessaires.

Le matériel IBM® possède une topologie de qubits en réseau heavy-hex, auquel cas nous pouvons adopter un motif en « zigzag », illustré ci-dessous. Dans ce motif, les orbitales de même spin sont associées à des qubits ayant une topologie linéaire (cercles rouges et bleus), et une connexion entre les orbitales de spin différent est présente tous les 4 orbitales spatiales, la connexion étant facilitée par un qubit ancillaire (cercles violets). Dans ce cas, les contraintes d'indices sont

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

Récupération auto-cohérente des configurations

La procédure de récupération auto-cohérente des configurations est conçue pour extraire autant de signal que possible des échantillons quantiques bruités. Étant donné que le Hamiltonien moléculaire conserve le nombre de particules et le spin Z, il est logique de choisir un ansatz de circuit qui conserve également ces symétries. Lorsqu'il est appliqué à l'état de Hartree-Fock, l'état résultant a un nombre de particules et un spin Z fixes dans le cas sans bruit. Par conséquent, les moitiés spin-$\alpha$ et spin-$\beta$ de toute chaîne de bits échantillonnée à partir de cet état devraient avoir le même poids de Hamming que dans l'état de Hartree-Fock. En raison de la présence de bruit dans les processeurs quantiques actuels, certaines chaînes de bits mesurées violeront cette propriété. Une forme simple de post-sélection consisterait à écarter ces chaînes de bits, mais cela serait un gaspillage car elles pourraient encore contenir du signal. La procédure de récupération auto-cohérente tente de récupérer une partie de ce signal en post-traitement. La procédure est itérative et nécessite en entrée une estimation des occupations moyennes de chaque orbitale dans l'état fondamental, qui est d'abord calculée à partir des échantillons bruts. La procédure est exécutée en boucle, et chaque itération comporte les étapes suivantes :

1. Pour chaque chaîne de bits qui viole les symétries spécifiées, inverser ses bits selon une procédure probabiliste conçue pour rapprocher la chaîne de bits de l'estimation actuelle des occupations orbitales moyennes, afin d'obtenir une nouvelle chaîne de bits.
2. Rassembler toutes les anciennes et nouvelles chaînes de bits satisfaisant les symétries, et sous-échantillonner des sous-ensembles d'une taille fixe, choisie à l'avance.
3. Pour chaque sous-ensemble de chaînes de bits, projeter le Hamiltonien dans le sous-espace engendré par les vecteurs de base correspondants (voir la section précédente pour une description de ces vecteurs de base), et calculer une estimation de l'état fondamental du Hamiltonien projeté sur un ordinateur classique.
4. Mettre à jour l'estimation des occupations orbitales moyennes avec l'estimation de l'état fondamental ayant l'énergie la plus basse.

Diagramme du workflow SQD

Le workflow SQD est représenté dans le diagramme suivant :

Prérequis

Avant de commencer ce tutoriel, assure-toi d'avoir installé les éléments suivants :

• Qiskit SDK v1.0 ou ultérieur, avec le support de visualisation
• Qiskit Runtime v0.22 ou ultérieur (pip install qiskit-ibm-runtime)
• Extension SQD de Qiskit v0.11 ou ultérieur (pip install qiskit-addon-sqd)
• ffsim v0.0.58 ou ultérieur (pip install ffsim)

Configuration

# 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

Étape 1 : Traduire les entrées classiques en un problème quantique

Dans ce tutoriel, nous trouverons une approximation de l'état fondamental de la molécule à l'équilibre dans l'ensemble de base cc-pVDZ. Tout d'abord, nous spécifions la molécule et ses propriétés.

# 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

Avant de construire le circuit d'ansatz LUCJ, nous effectuons d'abord un calcul CCSD dans la cellule de code suivante. Les amplitudes $t_1$ et $t_2$ issues de ce calcul seront utilisées pour initialiser les paramètres de l'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

Maintenant, nous utilisons ffsim pour créer le circuit d'ansatz. Puisque notre molécule possède un état de Hartree-Fock à couche fermée, nous utilisons la variante à spin équilibré de l'ansatz UCJ, UCJOpSpinBalanced. Nous passons des paires d'interaction appropriées pour une topologie de qubits en réseau heavy-hex (voir la section de contexte sur l'ansatz LUCJ). Nous définissons optimize=True dans la méthode from_t_amplitudes pour activer la double factorisation « compressée » des amplitudes $t_2$ (voir The local unitary cluster Jastrow (LUCJ) ansatz dans la documentation de ffsim pour plus de détails).

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

Étape 2 : Optimiser le problème pour l'exécution sur matériel quantique

Ensuite, nous optimisons le circuit pour un matériel cible.

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

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

Using backend ibm_fez

Nous recommandons les étapes suivantes pour optimiser l'ansatz et le rendre compatible avec le matériel.

• Sélectionner les qubits physiques (initial_layout) du matériel cible qui respectent le motif en « zigzag » décrit précédemment. Disposer les qubits selon ce motif conduit à un circuit efficace compatible avec le matériel, comportant moins de portes. Ici, nous incluons du code pour automatiser la sélection du motif en « zigzag » qui utilise une heuristique de notation pour minimiser les erreurs associées à la disposition sélectionnée.
• Générer un gestionnaire de passes par étapes à l'aide de la fonction generate_preset_pass_manager de Qiskit avec ton choix de backend et d'initial_layout.
• Définir l'étape pre_init de ton gestionnaire de passes par étapes sur ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT inclut des passes de transpilation Qiskit qui décomposent les portes en rotations orbitales puis fusionnent les rotations orbitales, résultant en moins de portes dans le circuit final.
• Exécuter le gestionnaire de passes sur ton circuit.

Code pour la sélection automatisée de la disposition en « zigzag »

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

Étape 3 : Exécuter à l'aide des primitives Qiskit

Après avoir optimisé le circuit pour l'exécution sur le matériel, nous sommes prêts à l'exécuter sur le matériel cible et à collecter des échantillons pour l'estimation de l'énergie de l'état fondamental. Comme nous n'avons qu'un seul circuit, nous utiliserons le mode d'exécution de tâches de Qiskit Runtime pour exécuter notre circuit.

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

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

Étape 4 : Post-traiter et retourner le résultat dans le format classique souhaité

Maintenant, nous estimons l'énergie de l'état fondamental du Hamiltonien à l'aide de la fonction diagonalize_fermionic_hamiltonian. Cette fonction effectue la procédure de récupération auto-cohérente des configurations pour affiner itérativement les échantillons quantiques bruités et améliorer l'estimation de l'énergie. Nous passons une fonction de rappel afin de pouvoir sauvegarder les résultats intermédiaires pour une analyse ultérieure. Consulte la documentation de l'API pour les explications des arguments de 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

Visualiser les résultats

Le premier graphique montre qu'après quelques itérations, nous estimons l'énergie de l'état fondamental à environ ~41 mH près (la précision chimique est généralement acceptée comme étant de 1 kcal/mol $\approx$ 1.6 mH). L'énergie peut être améliorée en autorisant davantage d'itérations de récupération de configurations ou en augmentant le nombre d'échantillons par lot.

Le second graphique montre l'occupation moyenne de chaque orbitale spatiale après l'itération finale. Nous pouvons observer que les électrons de spin-up et de spin-down occupent les cinq premières orbitales avec une forte probabilité dans nos solutions.

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

Enquête sur le tutoriel

Réponds à cette courte enquête pour donner ton avis sur ce tutoriel. Tes retours nous aideront à améliorer nos contenus et l'expérience utilisateur.

Lien vers l'enquête