Diagonalisasi kuantum berbasis sampel untuk Hamiltonian kimia

Estimasi penggunaan: kurang dari satu menit pada prosesor Heron r2 (CATATAN: Ini hanya estimasi. Waktu aktual bisa berbeda.)

Latar Belakang

Dalam tutorial ini, kita akan menunjukkan cara memproses sampel kuantum yang berisik untuk menemukan aproksimasi ke ground state molekul nitrogen $\text{N}_2$ pada panjang ikatan kesetimbangan, menggunakan algoritma diagonalisasi kuantum berbasis sampel (SQD), untuk sampel yang diambil dari sirkuit kuantum 59-Qubit (52 Qubit sistem + 7 Qubit ancilla). Implementasi berbasis Qiskit tersedia di SQD Qiskit addons, detail lebih lanjut dapat ditemukan di dokumentasi yang sesuai dengan contoh sederhana untuk memulai.

SQD adalah teknik untuk menemukan nilai eigen dan vektor eigen dari operator kuantum, seperti Hamiltonian sistem kuantum, menggunakan komputasi kuantum dan komputasi klasik terdistribusi bersama-sama. Komputasi klasik terdistribusi digunakan untuk memproses sampel yang diperoleh dari prosesor kuantum, serta untuk memproyeksikan dan mendiagonalisasi Hamiltonian target dalam subruang yang dibentangkan oleh sampel tersebut. Hal ini memungkinkan SQD untuk tahan terhadap sampel yang terkorupsi oleh noise kuantum dan menangani Hamiltonian besar, seperti Hamiltonian kimia dengan jutaan istilah interaksi, yang berada di luar jangkauan metode diagonalisasi eksak manapun. Alur kerja berbasis SQD memiliki langkah-langkah berikut:

1. Pilih ansatz sirkuit dan terapkan pada komputer kuantum ke sebuah state referensi (dalam kasus ini, state Hartree-Fock).
2. Sampel bitstring dari state kuantum yang dihasilkan.
3. Jalankan prosedur self-consistent configuration recovery pada bitstring untuk mendapatkan aproksimasi ke ground state.

SQD diketahui bekerja dengan baik ketika eigenstate target bersifat jarang: fungsi gelombang didukung dalam sekumpulan state basis $\mathcal{S} = \{|x\rangle \}$ yang ukurannya tidak bertambah secara eksponensial seiring dengan ukuran masalah.

Kimia kuantum

Sifat-sifat molekul sebagian besar ditentukan oleh struktur elektron di dalamnya. Sebagai partikel fermionic, elektron dapat dideskripsikan menggunakan formalisme matematika yang disebut kuantisasi kedua. Idenya adalah bahwa ada sejumlah orbital, yang masing-masing bisa kosong atau ditempati oleh sebuah fermion. Sistem $N$ orbital dideskripsikan oleh sekumpulan operator anihilasi fermionic $\{\hat{a}_p\}_{p=1}^N$ yang memenuhi relasi antikomutasi fermionic,

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

Adjoint $\hat{a}_p^\dagger$ disebut operator kreasi.

Sejauh ini, penjelasan kita belum memperhitungkan spin, yang merupakan sifat fundamental dari fermion. Ketika memperhitungkan spin, orbital hadir berpasangan yang disebut orbital spasial. Setiap orbital spasial terdiri dari dua orbital spin, satu berlabel "spin-$\alpha$" dan satu berlabel "spin-$\beta$". Kita kemudian menulis $\hat{a}_{p\sigma}$ untuk operator anihilasi yang terkait dengan orbital-spin dengan spin $\sigma$ ($\sigma \in \{\alpha, \beta\}$) dalam orbital spasial $p$. Jika kita ambil $N$ sebagai jumlah orbital spasial, maka total ada $2N$ orbital-spin. Ruang Hilbert sistem ini dibentangkan oleh $2^{2N}$ vektor basis ortonormal yang dilabeli dengan bitstring dua bagian $\lvert z \rangle = \lvert z_\beta z_\alpha \rangle = \lvert z_{\beta, N} \cdots z_{\beta, 1} z_{\alpha, N} \cdots z_{\alpha, 1} \rangle$.

Hamiltonian sistem molekuler dapat ditulis sebagai

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

di mana $h_{pr}$ dan $h_{prqs}$ adalah bilangan kompleks yang disebut integral molekuler yang dapat dihitung dari spesifikasi molekul menggunakan sebuah program komputer. Dalam tutorial ini, kita menghitung integral menggunakan paket perangkat lunak PySCF.

Untuk detail tentang cara Hamiltonian molekuler diturunkan, konsultasikan buku teks kimia kuantum (misalnya, Modern Quantum Chemistry oleh Szabo dan Ostlund). Untuk penjelasan tingkat tinggi tentang bagaimana masalah kimia kuantum dipetakan ke komputer kuantum, lihat kuliah Mapping Problems to Qubits dari Qiskit Global Summer School 2024.

Ansatz local unitary cluster Jastrow (LUCJ)

SQD memerlukan ansatz sirkuit kuantum untuk mengambil sampel. Dalam tutorial ini, kita akan menggunakan ansatz local unitary cluster Jastrow (LUCJ) karena kombinasinya antara motivasi fisik dan keramahan terhadap perangkat keras.

Ansatz LUCJ adalah bentuk khusus dari ansatz unitary cluster Jastrow (UCJ) umum, yang memiliki bentuk

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

di mana $\lvert \Phi_0 \rangle$ adalah state referensi, sering diambil sebagai state Hartree-Fock, dan $\hat{K}_\mu$ serta $\hat{J}_\mu$ memiliki bentuk

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

di mana kita telah mendefinisikan operator bilangan

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

Operator $e^{\hat{K}_\mu}$ adalah rotasi orbital, yang dapat diimplementasikan menggunakan algoritma yang dikenal dalam kedalaman linear dan menggunakan konektivitas linear. Mengimplementasikan suku $e^{i \mathcal{J}_k}$ dari ansatz UCJ memerlukan konektivitas all-to-all atau penggunaan jaringan fermionic swap, yang membuatnya menantang untuk prosesor kuantum pra-fault-tolerant yang berisik dengan konektivitas terbatas. Ide dari ansatz UCJ lokal adalah dengan menerapkan kendala kejarangan pada matriks $\mathbf{J}^{\alpha\alpha}$ dan $\mathbf{J}^{\alpha\beta}$ yang memungkinkan implementasinya dalam kedalaman konstan pada topologi Qubit dengan konektivitas terbatas. Kendala tersebut ditentukan oleh daftar indeks yang menunjukkan entri matriks mana di segitiga atas yang boleh bernilai bukan nol (karena matriksnya simetris, hanya segitiga atas yang perlu ditentukan). Indeks-indeks ini dapat ditafsirkan sebagai pasangan orbital yang diperbolehkan untuk berinteraksi.

Sebagai contoh, pertimbangkan topologi Qubit kisi persegi. Kita dapat menempatkan orbital $\alpha$ dan $\beta$ dalam garis sejajar pada kisi, dengan koneksi antar garis tersebut membentuk "rung" dari bentuk tangga, seperti ini:

Dengan pengaturan ini, orbital dengan spin yang sama terhubung dengan topologi garis, sedangkan orbital dengan spin berbeda terhubung jika mereka berbagi orbital spasial yang sama. Ini menghasilkan kendala indeks berikut pada matriks $\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*}$$

Dengan kata lain, jika matriks $\mathbf{J}$ bernilai bukan nol hanya pada indeks yang ditentukan di segitiga atas, maka suku $e^{i \mathcal{J}_k}$ dapat diimplementasikan pada topologi persegi tanpa menggunakan gate swap apapun, dalam kedalaman konstan. Tentu saja, menerapkan kendala seperti itu pada ansatz membuatnya kurang ekspresif, sehingga mungkin diperlukan lebih banyak pengulangan ansatz.

Perangkat keras IBM® memiliki topologi Qubit kisi heavy-hex, di mana kita dapat mengadopsi pola "zig-zag", yang digambarkan di bawah. Dalam pola ini, orbital dengan spin yang sama dipetakan ke Qubit dengan topologi garis (lingkaran merah dan biru), dan koneksi antara orbital dengan spin berbeda hadir setiap 4 orbital spasial ke-4, dengan koneksi yang difasilitasi oleh Qubit ancilla (lingkaran ungu). Dalam kasus ini, kendala indeksnya adalah

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

Self-consistent configuration recovery

Prosedur self-consistent configuration recovery dirancang untuk mengekstrak sinyal sebanyak mungkin dari sampel kuantum yang berisik. Karena Hamiltonian molekuler mengkonservasi jumlah partikel dan spin Z, masuk akal untuk memilih ansatz sirkuit yang juga mengkonservasi simetri-simetri ini. Ketika diterapkan ke state Hartree-Fock, state yang dihasilkan memiliki jumlah partikel dan spin Z yang tetap dalam pengaturan tanpa noise. Oleh karena itu, bagian spin-$\alpha$ dan spin-$\beta$ dari setiap bitstring yang disampling dari state ini seharusnya memiliki Hamming weight yang sama seperti pada state Hartree-Fock. Karena adanya noise pada prosesor kuantum saat ini, beberapa bitstring yang diukur akan melanggar sifat ini. Bentuk sederhana dari postseleksi akan membuang bitstring-bitstring ini, tetapi hal itu boros karena bitstring tersebut mungkin masih mengandung sinyal. Prosedur self-consistent recovery berupaya memulihkan sebagian sinyal tersebut dalam post-processing. Prosedur ini bersifat iteratif dan memerlukan estimasi awal rata-rata penghunian setiap orbital dalam ground state sebagai input, yang pertama kali dihitung dari sampel mentah. Prosedur ini dijalankan dalam sebuah loop, dan setiap iterasi memiliki langkah-langkah berikut:

1. Untuk setiap bitstring yang melanggar simetri yang ditentukan, balik bit-nya dengan prosedur probabilistik yang dirancang untuk membawa bitstring lebih dekat ke estimasi saat ini dari rata-rata penghunian orbital, untuk mendapatkan bitstring baru.
2. Kumpulkan semua bitstring lama dan baru yang memenuhi simetri, dan subsample subset dari ukuran tetap yang dipilih sebelumnya.
3. Untuk setiap subset bitstring, proyeksikan Hamiltonian ke subruang yang dibentangkan oleh vektor basis yang sesuai (lihat bagian sebelumnya untuk deskripsi vektor basis ini), dan hitung estimasi ground state dari Hamiltonian yang diproyeksikan pada komputer klasik.
4. Perbarui estimasi rata-rata penghunian orbital dengan estimasi ground state yang memiliki energi terendah.

Diagram alur kerja SQD

Alur kerja SQD digambarkan dalam diagram berikut:

Persyaratan

Sebelum memulai tutorial ini, pastikan kamu telah menginstal yang berikut:

• Qiskit SDK v1.0 atau lebih baru, dengan dukungan visualisasi
• Qiskit Runtime v0.22 atau lebih baru (pip install qiskit-ibm-runtime)
• SQD Qiskit addon v0.11 atau lebih baru (pip install qiskit-addon-sqd)
• ffsim v0.0.58 atau lebih baru (pip install ffsim)

Pengaturan

# 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

Langkah 1: Petakan input klasik ke masalah kuantum

Dalam tutorial ini, kita akan mencari aproksimasi terhadap ground state molekul pada kesetimbangan dalam basis set cc-pVDZ. Pertama, kita tentukan molekul dan properti-propertinya.

# 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

Sebelum membangun Circuit ansatz LUCJ, kita terlebih dahulu melakukan perhitungan CCSD pada sel kode berikut. Amplitudo $t_1$ dan $t_2$ dari perhitungan ini akan digunakan untuk menginisialisasi parameter 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

Sekarang, kita gunakan ffsim untuk membuat Circuit ansatz. Karena molekul kita memiliki state Hartree-Fock closed-shell, kita gunakan varian spin-balanced dari ansatz UCJ, yaitu UCJOpSpinBalanced. Kita berikan pasangan interaksi yang sesuai untuk topologi Qubit heavy-hex lattice (lihat bagian latar belakang tentang ansatz LUCJ). Kita set optimize=True pada metode from_t_amplitudes untuk mengaktifkan double-factorization "terkompresi" dari amplitudo $t_2$ (lihat The local unitary cluster Jastrow (LUCJ) ansatz di dokumentasi ffsim untuk detailnya).

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

Langkah 2: Optimalkan masalah untuk eksekusi hardware kuantum

Selanjutnya, kita optimalkan Circuit untuk 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

Kami merekomendasikan langkah-langkah berikut untuk mengoptimalkan ansatz dan membuatnya kompatibel dengan hardware.

• Pilih Qubit fisik (initial_layout) dari hardware target yang mengikuti pola "zig-zag" yang telah dijelaskan sebelumnya. Menata Qubit dalam pola ini menghasilkan Circuit yang efisien dan kompatibel dengan hardware menggunakan lebih sedikit Gate. Di sini, kita sertakan beberapa kode untuk mengotomasi pemilihan pola "zig-zag" menggunakan heuristik penilaian untuk meminimalkan error yang terkait dengan layout yang dipilih.
• Buat staged pass manager menggunakan fungsi generate_preset_pass_manager dari Qiskit dengan pilihan backend dan initial_layout.
• Atur tahap pre_init dari staged pass manager ke ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT mencakup Transpiler pass Qiskit yang mengurai Gate menjadi rotasi orbital dan kemudian menggabungkan rotasi orbital tersebut, sehingga menghasilkan lebih sedikit Gate dalam Circuit akhir.
• Jalankan pass manager pada Circuit kamu.

Kode untuk pemilihan otomatis 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})

Langkah 3: Eksekusi menggunakan Qiskit primitives

Setelah mengoptimalkan Circuit untuk eksekusi hardware, kita siap menjalankannya di hardware target dan mengumpulkan sampel untuk estimasi energi ground state. Karena kita hanya memiliki satu Circuit, kita akan menggunakan mode eksekusi Job dari Qiskit Runtime dan mengeksekusi Circuit kita.

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

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

Langkah 4: Pasca-proses dan kembalikan hasil dalam format klasik yang diinginkan

Sekarang, kita estimasi energi ground state dari Hamiltonian menggunakan fungsi diagonalize_fermionic_hamiltonian. Fungsi ini melakukan prosedur configuration recovery yang self-consistent untuk secara iteratif menyempurnakan sampel kuantum yang berisik guna meningkatkan estimasi energi. Kita berikan fungsi callback agar bisa menyimpan hasil antara untuk analisis lebih lanjut. Lihat dokumentasi API untuk penjelasan argumen-argumen 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

Visualisasikan hasilnya

Plot pertama menunjukkan bahwa setelah beberapa iterasi kita mengestimasi energi ground state dalam rentang ~41 mH (akurasi kimia umumnya diterima sebagai 1 kcal/mol $\approx$ 1.6 mH). Energi dapat ditingkatkan dengan mengizinkan lebih banyak iterasi configuration recovery atau menambah jumlah sampel per batch.

Plot kedua menunjukkan rata-rata occupancy setiap orbital spasial setelah iterasi terakhir. Kita bisa melihat bahwa elektron spin-up maupun spin-down menempati lima orbital pertama dengan probabilitas tinggi dalam solusi kita.

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

Survei tutorial

Silakan isi survei singkat ini untuk memberikan masukan tentang tutorial ini. Pendapatmu akan membantu kami meningkatkan konten dan pengalaman pengguna.

Link ke survei

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.