Pular para o conteúdo principal

Diagonalização quântica baseada em amostragem de um Hamiltoniano de química

Estimativa de uso: menos de um minuto em um processador Heron r2 (NOTA: Isto é apenas uma estimativa. Seu tempo de execução pode variar.)

Contexto​

Neste tutorial, mostramos como pós-processar amostras quânticas ruidosas para encontrar uma aproximação do estado fundamental da molécula de nitrogênio N2\text{N}_2 no comprimento de ligação de equilíbrio, usando o algoritmo de diagonalização quântica baseada em amostragem (SQD), para amostras obtidas de um circuito quântico de 59 qubits (52 qubits de sistema + 7 qubits ancilla). Uma implementação baseada em Qiskit está disponível nos addons Qiskit SQD, mais detalhes podem ser encontrados na documentação correspondente com um exemplo simples para começar. SQD é uma técnica para encontrar autovalores e autovetores de operadores quânticos, como um Hamiltoniano de sistema quântico, usando computação quântica e computação clássica distribuída juntas. A computação clássica distribuída é usada para processar amostras obtidas de um processador quântico e para projetar e diagonalizar um Hamiltoniano alvo em um subespaço gerado por elas. Isso permite que o SQD seja robusto a amostras corrompidas por ruído quântico e lide com Hamiltonianos grandes, como Hamiltonianos de química com milhões de termos de interação, além do alcance de quaisquer métodos de diagonalização exata. Um fluxo de trabalho baseado em SQD tem as seguintes etapas:

  1. Escolher um ansatz de circuito e aplicá-lo em um computador quântico a um estado de referência (neste caso, o estado de Hartree-Fock).
  2. Amostrar bitstrings do estado quântico resultante.
  3. Executar o procedimento de recuperação de configuração autoconsistente nas bitstrings para obter a aproximação do estado fundamental.

O SQD é conhecido por funcionar bem quando o autoestado alvo é esparso: a função de onda é suportada em um conjunto de estados de base S={∣x⟩}\mathcal{S} = \{|x\rangle \} cujo tamanho não aumenta exponencialmente com o tamanho do problema.

Química quântica​

As propriedades das moléculas são amplamente determinadas pela estrutura dos elétrons dentro delas. Como partículas fermiônicas, os elétrons podem ser descritos usando um formalismo matemático chamado segunda quantização. A ideia é que existem vários orbitais, cada um dos quais pode estar vazio ou ocupado por um férmion. Um sistema de NN orbitais é descrito por um conjunto de operadores de aniquilação fermiônicos {a^p}p=1N\{\hat{a}_p\}_{p=1}^N que satisfazem as relações de anticomutação fermiônicas,

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

O adjunto a^p†\hat{a}_p^\dagger é chamado de operador de criação.

Até agora, nossa exposição não considerou o spin, que é uma propriedade fundamental dos férmions. Ao considerar o spin, os orbitais vêm em pares chamados orbitais espaciais. Cada orbital espacial é composto por dois orbitais de spin, um rotulado como "spin-α\alpha" e um rotulado como "spin-β\beta". Então escrevemos a^pσ\hat{a}_{p\sigma} para o operador de aniquilação associado ao orbital de spin com spin σ\sigma (σ∈{α,β}\sigma \in \{\alpha, \beta\}) no orbital espacial pp. Se tomarmos NN como o número de orbitais espaciais, então há um total de 2N2N orbitais de spin. O espaço de Hilbert deste sistema é gerado por 22N2^{2N} vetores de base ortonormais rotulados com bitstrings de duas partes ∣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.

O Hamiltoniano de um sistema molecular pode ser escrito como

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

onde os hprh_{pr} e hprqsh_{prqs} são números complexos chamados integrais moleculares que podem ser calculados a partir da especificação da molécula usando um programa de computador. Neste tutorial, calculamos as integrais usando o pacote de software PySCF.

Para detalhes sobre como o Hamiltoniano molecular é derivado, consulte um livro-texto sobre química quântica (por exemplo, Modern Quantum Chemistry de Szabo e Ostlund). Para uma explicação de alto nível sobre como problemas de química quântica são mapeados em computadores quânticos, confira a palestra Mapping Problems to Qubits da Qiskit Global Summer School 2024.

Ansatz local unitary cluster Jastrow (LUCJ)​

O SQD requer um ansatz de circuito quântico para extrair amostras. Neste tutorial, usaremos o ansatz local unitary cluster Jastrow (LUCJ) devido à sua combinação de motivação física e compatibilidade com hardware.

O ansatz LUCJ é uma forma especializada do ansatz geral unitary cluster Jastrow (UCJ), que tem a 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

onde ∣Φ0⟩\lvert \Phi_0 \rangle é um estado de referência, frequentemente tomado como o estado de Hartree-Fock, e os K^μ\hat{K}_\mu e J^μ\hat{J}_\mu têm a 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} \;,

onde definimos o operador de número

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

O operador eK^μe^{\hat{K}_\mu} é uma rotação orbital, que pode ser implementada usando algoritmos conhecidos em profundidade linear e usando conectividade linear. Implementar o termo eiJke^{i \mathcal{J}_k} do ansatz UCJ requer conectividade todos-para-todos ou o uso de uma rede de swap fermiônica, tornando-o desafiador para processadores quânticos pré-tolerância a falhas ruidosos que têm conectividade limitada. A ideia do ansatz UCJ local é impor restrições de esparsidade nas matrizes Jαα\mathbf{J}^{\alpha\alpha} e Jαβ\mathbf{J}^{\alpha\beta} que permitem que sejam implementadas em profundidade constante em topologias de qubit com conectividade limitada. As restrições são especificadas por uma lista de índices indicando quais entradas de matriz no triângulo superior podem ser diferentes de zero (já que as matrizes são simétricas, apenas o triângulo superior precisa ser especificado). Esses índices podem ser interpretados como pares de orbitais que podem interagir.

Como exemplo, considere uma topologia de qubit em rede quadrada. Podemos colocar os orbitais α\alpha e β\beta em linhas paralelas na rede, com conexões entre essas linhas formando "degraus" de uma forma de escada, assim:

Qubit mapping diagram for the LUCJ ansatz on a square lattice

Com essa configuração, orbitais com o mesmo spin são conectados com uma topologia de linha, enquanto orbitais com spins diferentes são conectados quando compartilham o mesmo orbital espacial. Isso produz as seguintes restrições de índice nas matrizes 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*}

Em outras palavras, se as matrizes J\mathbf{J} são diferentes de zero apenas nos índices especificados no triângulo superior, então o termo eiJke^{i \mathcal{J}_k} pode ser implementado em uma topologia quadrada sem usar nenhum gate swap, em profundidade constante. É claro que impor tais restrições ao ansatz o torna menos expressivo, então mais repetições de ansatz podem ser necessárias.

O hardware IBM® tem uma topologia de qubit em rede heavy-hex, caso em que podemos adotar um padrão "zig-zag", representado abaixo. Neste padrão, orbitais com o mesmo spin são mapeados para qubits com uma topologia de linha (círculos vermelhos e azuis), e uma conexão entre orbitais de spin diferente está presente a cada 4º orbital espacial, com a conexão sendo facilitada por um qubit ancilla (círculos roxos). Neste caso, as restrições de índice são

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

Recuperação de configuração autoconsistente​

O procedimento de recuperação de configuração autoconsistente é projetado para extrair o máximo de sinal possível de amostras quânticas ruidosas. Como o Hamiltoniano molecular conserva o número de partículas e o spin Z, faz sentido escolher um ansatz de circuito que também conserve essas simetrias. Quando aplicado ao estado de Hartree-Fock, o estado resultante tem um número fixo de partículas e spin Z no cenário sem ruído. Portanto, as metades de spin-α\alpha e spin-β\beta de qualquer bitstring amostrada deste estado devem ter o mesmo peso de Hamming que no estado de Hartree-Fock. Devido à presença de ruído nos processadores quânticos atuais, algumas bitstrings medidas violarão essa propriedade. Uma forma simples de pós-seleção descartaria essas bitstrings, mas isso é desperdiçador porque essas bitstrings ainda podem conter algum sinal. O procedimento de recuperação autoconsistente tenta recuperar parte desse sinal no pós-processamento. O procedimento é iterativo e requer como entrada uma estimativa das ocupações médias de cada orbital no estado fundamental, que é primeiro calculada a partir das amostras brutas. O procedimento é executado em loop, e cada iteração tem as seguintes etapas:

  1. Para cada bitstring que viola as simetrias especificadas, inverter seus bits com um procedimento probabilístico projetado para aproximar a bitstring da estimativa atual das ocupações orbitais médias, para obter uma nova bitstring.
  2. Coletar todas as bitstrings antigas e novas que satisfazem as simetrias e subamostrar subconjuntos de um tamanho fixo, escolhido antecipadamente.
  3. Para cada subconjunto de bitstrings, projetar o Hamiltoniano no subespaço gerado pelos vetores de base correspondentes (veja a seção anterior para uma descrição desses vetores de base), e calcular uma estimativa de estado fundamental do Hamiltoniano projetado em um computador clássico.
  4. Atualizar a estimativa das ocupações orbitais médias com a estimativa de estado fundamental com a menor energia.

Diagrama de fluxo de trabalho SQD​

O fluxo de trabalho SQD é representado no seguinte diagrama:

Workflow diagram of the SQD algorithm

Requisitos​

Antes de iniciar este tutorial, certifique-se de ter o seguinte instalado:

  • Qiskit SDK v1.0 ou posterior, com suporte a visualização
  • Qiskit Runtime v0.22 ou posterior (pip install qiskit-ibm-runtime)
  • SQD Qiskit addon v0.11 ou posterior (pip install qiskit-addon-sqd)
  • ffsim v0.0.58 ou posterior (pip install ffsim)

Configuração​

# 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: Mapear entradas clássicas para um problema quântico​

Neste tutorial, encontraremos uma aproximação do estado fundamental da molécula no equilíbrio no conjunto de base cc-pVDZ. Primeiro, especificamos a molécula e suas propriedades.

# 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

Antes de construir o circuito ansatz LUCJ, primeiro realizamos um cálculo CCSD na célula de código a seguir. As amplitudes t1t_1 e t2t_2 deste cálculo serão usadas para inicializar os parâmetros do 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

Agora, usamos ffsim para criar o circuito ansatz. Como nossa molécula tem um estado de Hartree-Fock de camada fechada, usamos a variante com spin balanceado do ansatz UCJ, UCJOpSpinBalanced. Passamos pares de interação apropriados para uma topologia de qubit em rede heavy-hex (veja a seção de contexto sobre o ansatz LUCJ). Definimos optimize=True no método from_t_amplitudes para habilitar a "compressão" de dupla fatoração das amplitudes t2t_2 (veja The local unitary cluster Jastrow (LUCJ) ansatz na documentação do ffsim para detalhes).

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: Otimizar o problema para execução em hardware quântico​

Em seguida, otimizamos o circuito para um hardware alvo.

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

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

Recomendamos os seguintes passos para otimizar o ansatz e torná-lo compatível com hardware.

  • Selecionar qubits físicos (initial_layout) do hardware alvo que adere ao padrão "zig-zag" descrito anteriormente. Dispor qubits neste padrão leva a um circuito eficiente e compatível com hardware com menos gates. Aqui, incluímos algum código para automatizar a seleção do padrão "zig-zag" que usa uma heurística de pontuação para minimizar os erros associados ao layout selecionado.
  • Gerar um gerenciador de passes em estágios usando a função generate_preset_pass_manager do qiskit com sua escolha de backend e initial_layout.
  • Definir o estágio pre_init do seu gerenciador de passes em estágios como ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT inclui passes do transpilador qiskit que decompõem gates em rotações orbitais e depois mesclam as rotações orbitais, resultando em menos gates no circuito final.
  • Executar o gerenciador de passes no seu circuito.
Código para seleção automatizada do 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: Executar usando primitivos do Qiskit​

Após otimizar o circuito para execução em hardware, estamos prontos para executá-lo no hardware de destino e coletar amostras para estimativa da energia do estado fundamental. Como temos apenas um circuito, usaremos o modo de execução de Job do Qiskit Runtime e executaremos nosso circuito.

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

Passo 4: Pós-processar e retornar resultado no formato clássico desejado​

Agora, estimamos a energia do estado fundamental do Hamiltoniano usando a função diagonalize_fermionic_hamiltonian. Esta função executa o procedimento de recuperação de configuração auto-consistente para refinar iterativamente as amostras quânticas ruidosas e melhorar a estimativa de energia. Passamos uma função de callback para que possamos salvar os resultados intermediários para análise posterior. Consulte a documentação da API para explicações dos argumentos 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

Visualizar os resultados​

O primeiro gráfico mostra que após algumas iterações estimamos a energia do estado fundamental dentro de ~41 mH (precisão química é tipicamente aceita como sendo 1 kcal/mol ≈\approx 1.6 mH). A energia pode ser melhorada permitindo mais iterações de recuperação de configuração ou aumentando o número de amostras por lote.

O segundo gráfico mostra a ocupação média de cada orbital espacial após a iteração final. Podemos ver que tanto os elétrons spin-up quanto spin-down ocupam os primeiros cinco orbitais com alta probabilidade em nossas soluções.

# 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

Pesquisa do tutorial​

Por favor, responda a esta breve pesquisa para fornecer feedback sobre este tutorial. Suas percepções nos ajudarão a melhorar nossas ofertas de conteúdo e experiência do usuário.

Link to survey