화학 해밀토니안의 샘플 기반 양자 대각화
사용량 추정: Heron r2 프로세서에서 1분 미만 (참고: 이는 추정치이며 실제 실행 시간은 다를 수 있습니다.)
배경
이 튜토리얼에서는 잡음 있는 양자 샘플을 후처리하여 질소 분자 의 평형 결합 길이에서의 기저 상태 근사값을 구하는 방법을 소개합니다. 59-Qubit 양자 Circuit(52개 시스템 Qubit + 7개 보조 Qubit)에서 얻은 샘플을 사용하여 샘플 기반 양자 대각화(SQD) 알고리즘을 적용합니다. Qiskit 기반 구현은 SQD Qiskit 애드온에서 제공되며, 관련 문서와 시작을 위한 간단한 예제에서 자세한 내용을 확인할 수 있습니다. SQD는 양자 및 분산 클래식 컴퓨팅을 함께 사용하여 양자 해밀토니안과 같은 양자 연산자의 고유값과 고유벡터를 찾는 기법입니다. 클래식 분산 컴퓨팅은 양자 프로세서에서 얻은 샘플을 처리하고, 샘플이 형성하는 부분 공간에서 목표 해밀토니안을 투영하여 대각화하 는 데 사용됩니다. 이를 통해 SQD는 양자 잡음으로 손상된 샘플에도 견고하게 동작하며, 수백만 개의 상호작용 항을 가진 화학 해밀토니안처럼 정확한 대각화 방법으로는 다룰 수 없는 대형 해밀토니안도 처리할 수 있습니다. SQD 기반 워크플로우는 다음 단계로 구성됩니다.
- Circuit 앤사츠를 선택하고 양자 컴퓨터에서 기준 상태(이 경우 하트리-폭(Hartree-Fock) 상태)에 적용합니다.
- 결과 양자 상태에서 비트스트링을 샘플링합니다.
- 비트스트링에 대해 자기 일관적 구성 복구 절차를 실행하여 기저 상태의 근사값을 구합니다.
SQD는 목표 고유상태가 희소할 때 잘 동작하는 것으로 알려져 있습니다. 즉, 파동 함수가 문제 크기에 따라 지수적으로 증가하지 않는 기저 상태 집합 에 의해 지지될 때 효과적입니다.
양자 화학
분자의 특성은 주로 전자의 구조에 의해 결정됩니다. 페르미온 입자인 전자는 2차 양자화라는 수학적 형식으로 기술할 수 있습니다. 이 형식에서는 여러 개의 오비탈이 존재하며, 각 오비탈은 비어 있거나 페르미온으로 점유될 수 있습니다. 개의 오비탈로 구성된 시스템은 페르미온 반교환 관계를 만족하는 소멸 연산자 집합 으로 기술됩니다.
수반 연산자 는 생성 연산자라고 합니다.
지금까지의 설명에서는 페르미온의 기본 특성인 스핀을 고려하지 않았습니다. 스핀을 고려하면 오비탈은 공간 오비탈이라고 불리는 쌍으로 구성됩니다. 각 공간 오비탈은 "스핀-"와 "스핀-"로 레이블된 두 개의 스핀 오비탈로 이루어집니다. 공간 오비탈 에서 스핀 ()를 가진 스핀 오비탈과 연관된 소멸 연산자를 로 씁니다. 을 공간 오비탈의 수라 하면 총 개의 스핀 오비탈이 존재합니다. 이 시스템의 힐베르트 공간은 두 부분으로 구성된 비트스트링으로 레이블된 개의 정규 직교 기저 벡터 로 생성됩니다.
분자 시스템의 해밀토니안은 다음과 같이 쓸 수 있습니다.
여기서 과 는 분자 적분이라고 불리는 복소수로, 컴퓨터 프로그램을 사용하여 분자의 사양으로부터 계산됩니다. 이 튜토리얼에서는 PySCF 소프트웨어 패키지를 사용하여 적분을 계산합니다.
분자 해밀토니안의 유도 과정에 대한 자세한 내용은 양자 화학 교재(예: Szabo와 Ostlund의 Modern Quantum Chemistry)를 참고하세요. 양자 화학 문제가 양자 컴퓨터에 어떻게 매핑되는지에 대한 고수준 설명은 Qiskit Global Summer School 2024의 강의 Mapping Problems to Qubits를 확인하세요.
국소 유니터리 클러스터 야스트로프(LUCJ) 앤사츠
SQD는 샘플을 추출하기 위한 양자 Circuit 앤사츠가 필요합니다. 이 튜토리얼에서는 물리적 동기와 하드웨어 친화성의 조합 덕분에 국소 유니터리 클러스터 야스트로프(LUCJ) 앤사츠를 사용합니다.
LUCJ 앤사츠는 일반 유니터리 클러스터 야스트로프(UCJ) 앤사츠의 특수화된 형태로, 다음과 같은 형식을 가집니다.
여기서 은 기준 상태(보통 하트리-폭 상태로 선택)이며, 와 는 다음과 같은 형식입니다.
여기서 수 연산자를 다음과 같이 정의합니다.
연산자 는 오비탈 회전으로, 선형 깊이와 선형 연결성을 사용하는 알려진 알고리즘으로 구현할 수 있습니다. UCJ 앤사츠의 항을 구현하려면 전역 연결성이나 페르미온 스왑 네트워크가 필요하여, 제한된 연결성을 가진 잡음 있는 내결함성 이전 양자 프로세서에서는 구현이 어렵습니다. 국소 UCJ 앤사츠의 아이디어는 및 행렬에 희소성 제약 조건을 부과하여 제한된 연결성을 가진 Qubit 토폴로지에서 상수 깊이로 구현할 수 있도록 하는 것입니다. 제약 조건은 상삼각 행렬에서 0이 아닌 값을 가질 수 있는 행렬 항목의 인덱스 목록으로 지정됩니다(행렬이 대칭이므로 상삼각만 지정하면 됩니다). 이 인덱스는 상호작용이 허용되는 오비탈 쌍으로 해석할 수 있습니다.
예를 들어, 정사각형 격자 Qubit 토폴로지를 생각해봅시다. 및 오비탈을 격자의 평행한 두 줄에 배치하고, 두 줄 사이의 연결이 사다리 모양의 "가로대"를 형성하도록 다음과 같이 배치할 수 있습니다.
이 설정에서 같은 스핀을 가진 오비탈은 선형 토폴로지로 연결되고, 다른 스핀을 가진 오비탈은 같은 공간 오비탈을 공유할 때 연결됩니다. 이를 통해 행렬에 대한 다음과 같은 인덱스 제약 조건이 생깁니다.
즉, 행렬이 상삼각의 지정된 인덱스에서만 0이 아닌 값을 가지면, 항을 스왑 Gate 없이 상수 깊이로 정사각형 토폴로지에서 구현할 수 있습니다. 물론, 앤사츠에 이러한 제약을 부과하면 표현력이 줄어들어 더 많은 앤사츠 반복이 필요할 수 있습니다.
IBM® 하드웨어는 헤비-헥스 격자 Qubit 토폴로지를 사용하며, 이 경우 아래에 나타낸 "지그재그" 패턴을 채택할 수 있습니다. 이 패턴에서 같은 스핀을 가진 오비탈은 선형 토폴로지의 Qubit에 매핑되고(빨간색과 파란색 원), 다른 스핀의 오비탈 간 연결은 4번째 공간 오비탈마다 존재하며, 보조 Qubit(보라색 원)에 의해 연결이 이루어집니다. 이 경우 인덱스 제약 조건은 다음과 같습니다.
자기 일관적 구성 복구
자기 일관적 구성 복구 절차는 잡음 있는 양자 샘플에서 가능한 한 많은 신호를 추출하도록 설계되었습니다. 분자 해밀토니안은 입자 수와 스핀 Z를 보존하므로, 이러한 대칭성을 보존하는 Circuit 앤사츠를 선택하는 것이 합리적입니다. 하트리-폭 상태에 적용하면, 잡음이 없는 경우 결과 상태는 고정된 입자 수와 스핀 Z를 가집니다. 따라서 이 상태에서 샘플링된 비트스트링의 스핀- 및 스핀- 절반은 하트리-폭 상태와 동일한 해밍 무게를 가져야 합니다. 현재 양자 프로세서의 잡음으로 인해 일부 측정된 비트스트링은 이 특성을 위반합니다. 단순한 사후 선택 방법은 이러한 비트스트링을 버리겠지만, 이는 낭비적입니다. 왜냐하면 해당 비트스트링에도 일부 신호가 포함될 수 있기 때문입니다. 자기 일관적 복구 절차는 후처리에서 그 신호 일부를 복구하려고 시도합니다. 이 절차는 반복적이며, 기저 상태의 각 오비탈 평균 점유율 추정값을 입력으로 필요로 하고, 이 추정값은 먼저 원시 샘플에서 계산됩니다. 절차는 루프로 실행되며, 각 반복은 다음 단계로 구성됩니다.
- 지정된 대칭성을 위반하는 각 비트스트링에 대해, 비트스트링을 현재 평균 오비탈 점유율 추정값에 더 가깝게 만들도록 설계된 확률적 절차로 비트를 뒤집어 새로운 비트스트링을 얻습니다.
- 대칭성을 만족하는 이전 및 새 비트스트 링을 모두 수집하고, 사전에 선택된 고정 크기의 부분 집합을 서브샘플링합니다.
- 각 비트스트링 부분 집합에 대해, 해당 기저 벡터가 생성하는 부분 공간으로 해밀토니안을 투영하고(이전 섹션에서 이 기저 벡터에 대한 설명 참조), 클래식 컴퓨터에서 투영된 해밀토니안의 기저 상태 추정값을 계산합니다.
- 가장 낮은 에너지를 가진 기저 상태 추정값으로 평균 오비탈 점유율 추정값을 업데이트합니다.
SQD 워크플로우 다이어그램
SQD 워크플로우는 다음 다이어그램에 나타나 있습니다.
요구 사항
이 튜토리얼을 시작하기 전에 다음이 설치되어 있는지 확인하세요.
- Qiskit SDK v1.0 이상, 시각화 지원 포함
- Qiskit Runtime v0.22 이상 (
pip install qiskit-ibm-runtime) - SQD Qiskit 애드온 v0.11 이상 (
pip install qiskit-addon-sqd) - ffsim v0.0.58 이상 (
pip install ffsim)
설정
# 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
1단계: 클래식 입력을 양자 문제로 매핑
이 튜토리얼에서는 cc-pVDZ 기저 집합에서 평형 상태의 분자 기저 상태 근사값을 구합니다. 먼저 분자와 그 특성을 지정합니다.
# 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
LUCJ 앤사츠 Circuit을 구성하기 전에, 다음 코드 셀에서 먼저 CCSD 계산을 수행합니다. 이 계산에서 얻은 및 진폭은 앤사츠의 매개변수를 초기화하는 데 사용됩니다.
# 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
이제 ffsim을 사용하여 앤사츠 Circuit을 만듭니다. 분자가 닫힌 껍질 하트리-폭 상태를 가지므로, UCJ 앤사츠의 스핀 균형 변형인 UCJOpSpinBalanced를 사용합니다. 헤비-헥스 격자 Qubit 토폴로지에 적합한 상호작용 쌍을 전달합니다(LUCJ 앤사츠 배경 섹션 참조). from_t_amplitudes 메서드에서 optimize=True를 설정하여 진폭의 "압축된" 이중 인수분해를 활성화합니다(자세한 내용은 ffsim 문서의 The local unitary cluster Jastrow (LUCJ) ansatz 참조).
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()
Step 2: 양자 하드웨어 실행을 위한 문제 최적화
다음으로, 대상 하드웨어에 맞게 Circuit를 최적화합니다.
service = QiskitRuntimeService()
backend = service.least_busy(
operational=True, simulator=False, min_num_qubits=133
)
print(f"Using backend {backend.name}")
Using backend ibm_fez
ansatz를 최적화하고 하드웨어 호환 가능하게 만들기 위해 다음 단계를 권장합니다.
- 앞서 설명한 "지그재그" 패턴을 따르는 물리적 Qubit(
initial_layout)을 대상 하드웨어에서 선택합니다. 이 패턴으로 Qubit를 배치하면 Gate 수가 적은 효율적인 하드웨어 호환 Circuit를 생성할 수 있습니다. 여기서는 선택 된 레이아웃과 관련된 오류를 최소화하는 점수 휴리스틱을 사용하여 "지그재그" 패턴 선택을 자동화하는 코드를 포함합니다. - qiskit의 generate_preset_pass_manager 함수를 사용하여 선택한
backend와initial_layout으로 단계별 패스 매니저를 생성합니다. - 단계별 패스 매니저의
pre_init단계를ffsim.qiskit.PRE_INIT으로 설정합니다.ffsim.qiskit.PRE_INIT에는 Gate를 궤도 회전으로 분해한 다음 궤도 회전을 병합하는 Qiskit Transpiler 패스가 포함되어 있어 최종 Circuit의 Gate 수가 줄어듭니다. - Circuit에 패스 매니저를 실행합니다.
"지그재그" 레이아웃 자동 선택 코드
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})
Step 3: Qiskit 프리미티브를 사용하여 실행
하드웨어 실행을 위한 Circuit 최적화가 완료되면, 대상 하드웨어에서 실행하여 바닥 상태 에너지 추정을 위한 샘플을 수집할 준비가 됩니다. Circuit가 하나뿐이므로 Qiskit Runtime의 작업 실행 모드를 사용하여 Circuit를 실행합니다.
sampler = Sampler(mode=backend)
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]
Step 4: 후처리 및 원하는 고전 형식으로 결과 반환
이제 diagonalize_fermionic_hamiltonian 함수를 사용하여 해밀토니안의 바닥 상태 에너지를 추정합니다. 이 함수는 자기 일관적 구성 복구 절차를 수행하여 노이즈가 있는 양자 샘플을 반복적으로 개선함으로써 에너지 추정값을 향상시킵니다. 나중에 분석을 위해 중간 결과를 저장할 수 있도록 콜백 함수를 전달합니다. diagonalize_fermionic_hamiltonian의 인수 설명은 API 문서를 참조하세요.
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
결과 시각화
첫 번째 플롯은 몇 번의 반복 후에 바닥 상태 에너지를 ~41 mH 이내로 추정함을 보여줍니다 (화학적 정확도는 일반적으로 1 kcal/mol 1.6 mH로 허용됩니다). 구성 복구 반복 횟수를 늘리거나 배치당 샘플 수를 늘리면 에너지를 개선할 수 있습니다.
두 번째 플롯은 마지막 반복 후 각 공간 궤도의 평균 점유율을 보여줍니다. 스핀-업 전자와 스핀-다운 전자 모두 우리의 솔루션에서 높은 확률로 처음 다섯 개의 궤도를 점유하고 있음을 알 수 있습니다.
# 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()
튜토리얼 설문조사
이 튜토리얼에 대한 피드백을 제공하기 위해 간단한 설문조사에 참여해 주세요. 여러분의 의견은 콘텐츠 및 사용자 경험을 개선하는 데 도움이 됩니다.
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.