주 콘텐츠로 건너뛰기

바닥상태 추정을 위한 SqDRIFT 알고리즘

사용량 추정치: Heron r3 프로세서에서 180초 (참고: 이는 추정치일 뿐입니다. 실제 실행 시간은 다를 수 있습니다.)

C++ 버전을 찾고 계신가요?

이 튜토리얼은 Python을 사용합니다. 소스 코드와 빌드 지침을 포함한 C++ 구현은 C++ SqDRIFT 튜토리얼을 참조하세요.

학습 목표​

  • Trotterization과 비교하여 더 작은 깊이의 회로를 만드는 방법 배우기

  • qDRIFT와 SQD를 사용한 바닥상태 추정을 위한 종단 간 워크플로 따라하기

  • 이러한 워크플로를 구현하기 위해 qiskit-fermions를 다른 Qiskit 애드온과 함께 사용하는 방법 배우기

이 튜토리얼은 교육 목적으로 Python 노트북 형태로 제공됩니다.

사전 요구 사항​

배경​

SqDRIFT는 비트스트링을 샘플링할 안자츠를 선택해야 하는 필요성을, 목표 해밀토니안으로부터 직접 구성된 시간 발전 회로들의 앙상블로 대체하는 SKQD의 변형입니다. 이는 해밀토니안의 계수를 기반으로 더 작은 시간 발전 연산자를 해밀토니안으로부터 서브샘플링함으로써 달성되며, 이는 qDRIFT 트로터화 방법으로 알려져 있습니다.

이 튜토리얼은 qDRIFT 알고리즘을 위해 더 자연스러운 페르미온 회로를 만들기 위해 Qiskit Fermions를 사용하며, 그 뒤에 하드웨어 실행을 위한 전통적인 Qiskit 파이프라인에 회로를 연결하기 전에 페르미온 레이아웃과 합성 패스를 사용합니다.

해밀토니안이 다음과 같은 형태라고 하겠습니다.

H=∑i=1NcihiH = \sum_{i=1}^{N} c_i h_i

여기서, 일반성을 잃지 않고 ci>0c_i > 0이며 hih_i의 최대 고유값이 절댓값으로 11과 같다고 요구합니다. 부호나 복소수 형태의 앞자리 인수는 모두 hih_i에 흡수되므로, 계수 cic_i는 엄격하게 양의 가중치이고 hih_i는 각 항의 방향을 담당합니다. 여기서 NN은 해밀토니안 안의 항의 개수(또는 그룹화 후에는 그룹의 개수)이며, 이는 해밀토니안의 속성으로서 하나의 회로로 샘플링되는 연산자 개수와는 구별됩니다. 후자는 아래에서 nn으로 표기됩니다.

그러면 qDRIFT 알고리즘은 목표 시간 tt에 대해 어떤 연산자 VkV_k를 구현하며, 여기서 kk는 1⋯K1 \cdots K까지 진행하고 kk번째 SqDRIFT 회로를 나타내며, 다음과 같이 정의됩니다.

Vk=∏j=1ne−ihkjλt/nV_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }

여기서 nn은 회로당 샘플링된 연산자의 개수이고 KK는 앙상블 내 회로의 개수입니다. 이 곱은 전체 NN개의 해밀토니안 항이 아니라 nn번의 추출에 대해 실행되며, 항들이 복원 추출되기 때문에 동일한 hih_i가 하나의 VkV_k 안에 한 번 이상 나타날 수 있습니다.

다음 양은:

λ=∑i=1Nci\lambda = \sum_{i=1}^{N} c_i

계수들의 L1L_1 노름이며, 따라서 어떤 항이 추출되었는지와 무관하게 nn개의 각 단계는 동일한 시간 λt/n\lambda t / n 동안 발전합니다. 단계 각도의 균일성이 qDRIFT의 특징적인 성질입니다. 계수는 그 항이 얼마나 멀리 회전되는지가 아니라 그 항이 얼마나 자주 추출되는지를 통해 결과에 영향을 미칩니다. 인덱스는 다음 분포로부터 샘플링됩니다.

P[ki]=ciλP[k_i] = \frac{c_i}{\lambda}

따라서 수열 (k1,…,kn)(k_1, \ldots, k_n)은 이 분포에서 추출된 항 인덱스들의 무작위 수열입니다. cic_i가 양수이고 합이 λ\lambda이므로, 이는 정규화된 확률 분포이며, 무작위 추출에 대한 결과 채널의 기댓값은 HH 아래에서의 발전을 근사하며, 그 오차는 nn이 커질수록 감소합니다. 근사 오차는 항의 개수 NN이 아니라 λ\lambda에 의존한다는 점에 유의하세요.

(SqDRIFT 논문에서는 항의 개수를 N\mathcal{N}으로, 수열 길이를 NN으로 표기합니다. 여기서는 두 가지를 명확히 구별하기 위해 NN과 nn을 사용합니다.)

이 튜토리얼은 이러한 무작위화된 회로들의 앙상블을 생성하는 방법을 보여줍니다. 이러한 회로를 만든 후, 서로 다른 연산자에 대해 Krylov 부분공간을 만드는 방식과 유사하게, 서로 다른 시간 매개변수를 가진 여러 연산자로부터 비트스트링을 샘플링합니다. 이는 바닥상태 벡터와 샘플링된 비트스트링 사이의 겹침을 더 높게 보장합니다.

요구 사항​

이 튜토리얼을 시작하기 전에 다음이 설치되어 있는지 확인하세요

  • Python(>=3.10) 가상 환경
  • pip>=25.1
  • qiskit ~= 2.5
  • qiskit-fermions==0.1.0 (이름이 복수형임에 유의하세요)
  • numpy
  • pyscf
  • qiskit-aer
  • qiskit-ibm-runtime
  • qiskit-addon-sqd

다음 명령으로 필요한 모든 패키지를 설치할 수 있습니다.

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

설정​

# Added by doQumentation — required packages for this notebook
!pip install -q numpy pyscf qiskit qiskit-addon-sqd qiskit-aer qiskit-fermions qiskit-ibm-runtime
# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

시뮬레이터 예제​

1단계: 고전 입력을 양자 문제로 매핑하기​

FCIDump 읽기 및 준비하기

이 튜토리얼에서는 질소(N2)의 전자 구조 해밀토니안을 불러오겠습니다. 페르미온 연산자를 만드는 다른 방법들도 있습니다. qiskit_fermions.operators.library 문서를 참조하세요.

이 FCIDump에 대하여. N2_sto_3g 파일은 실험적 평형 결합 길이인 1.09 A˚\AA의 원자 간 거리에서 최소 STO-3G 기저를 사용한 질소 분자(N2N_2)를 나타냅니다. 헤더에는 NORB=10, NELEC=14, MS2=0이 선언되어 있습니다. 즉 10개의 공간 오비탈(따라서 20개의 스핀 오비탈, Jordan-Wigner 변환 하에서 20개의 큐비트), 스핀 단일항 상태의 14개의 전자, 즉 일곱 개의 α\alpha 전자와 일곱 개의 β\beta 전자를 의미합니다. 모든 오비탈에는 대칭 라벨 1이 부여되어 있으며, 이는 점군 대칭이 활용되지 않는다는 것을 의미합니다. 전체 공간 STO-3G 덤프이므로 동결된 오비탈이 없으며, 상관 공간이 충분히 작아서 다음 셀에서 보여주듯이 비교를 위한 정확한 FCI 기준 에너지를 고전적으로 계산할 수 있습니다.

동일한 파일을 PySCF로 다시 생성할 수 있습니다.

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

적분값이 수렴된 SCF 오비탈에 의존하기 때문에, 다시 생성된 파일은 제공된 파일과 오비탈 위상이나 순서가 다를 수 있습니다. 그러나 총 에너지에는 영향을 미치지 않습니다.

파일 가져오기. FCIDump는 이 GitHub 저장소에서 찾을 수 있습니다. 아래 셀을 실행하여 튜토리얼의 나머지 부분이 예상하는 위치로 파일을 가져올 수 있습니다.

먼저 pyscf에서 제공하는 cisolver를 사용하여 기준 에너지를 구합니다. 이는 우리가 다루는 분자의 실제 바닥 상태 에너지입니다. 이를 위해 먼저 오비탈 수와 전자 수를 각각 나타내는 norb와 nelec을 선언합니다. 그런 다음 일전자 및 이전자 적분을 각각 나타내는 h1e와 h2e를 선언합니다. 이들은 모두 나중에 SQD에도 사용됩니다.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "assets/sqdrift/fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")
Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g
name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")
Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

해밀토니안 불러오기

필요한 데이터가 준비되었으므로, qiskit-fermions와 호환되는 형식으로 FCI 파일에서 해밀토니안을 읽어옵니다

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

qiskit-fermions를 이용한 페르미온 워크플로

먼저 페르미온 회로에 특화된 Transpiler 패스와 게이트를 제공하는 qiskit-fermions를 사용하여 해밀토니안을 페르미온 회로 모델로 매핑합니다. 이들은 나중에 이 워크플로에서 Qiskit의 기존 Transpiler 패스보다 먼저 사용됩니다.

항 그룹화

결과의 재현성을 보장하기 위해 먼저 canonical_order를 사용하여 항을 오직 그 구조에 따라서만 정렬합니다. 따라서 canon 목록 내 연산자의 순서는 고정됩니다. 앞으로 사용할 QDriftTrotterization 패스가 qDRIFT 연산자를 생성하기 위해 무작위 인덱스를 샘플링하기 때문에, 이 과정은 생성된 연산자의 재현성을 보장합니다.

이 단계에서는 동일한 계수를 가진 관련 항들을 그룹화하여 전자 구조 해밀토니안에 존재하는 많은 대칭성을 활용합니다. 이렇게 하면 qDRIFT 프로토콜이 샘플링하는 연산자 계수 분포가 변경되지만, 이는 수렴 보장에는 영향을 미치지 않습니다. 결정적으로, 대칭성에 의해 관련된 항을 그룹화하면 파울리 항의 유리한 상쇄가 일어나고, 이러한 항의 작용 하에 상태를 시간 전개할 때 전체적으로 더 짧은 회로 깊이를 얻을 수 있습니다.

qiskit-fermions는 이러한 그룹화를 수행해 주는 group_terms_by_electronic_structure 함수를 제공합니다.

group_terms_by_electronic_structure는 정규 순서화된 항을 가정한다는 점에 유의하세요.

대각 항 필터링

회로를 생성하는 데 사용되는 해밀토니안에서 대각 항을 제거하여, nn개의 qDRIFT 샘플링 슬롯이 구성 간 인구(population)를 이동시키는 항에 사용되도록 합니다. 이러한 항은 다음 단계에서 Evolution 게이트가 구성되기 전, 바로 이 시점에서 해밀토니안으로부터 걸러내는 것이 가장 좋습니다.

여기서 말하는 항은 점유수 기저(occupation-number basis)에서 대각인 항, 즉 수 연산자 ai†aia^\dagger_i a_i의 곱입니다. 이 범주에는 세 가지 종류의 항이 있습니다.

  • 수 연산자가 하나도 없는 곱인 상수 에너지 오프셋으로, 그 시간 전개는 전역 위상만을 기여합니다.

  • 개별 수 연산자 nin_i로, 그 시간 전개는 단일 큐비트 ZZ 회전으로 귀결됩니다.

  • ninjn_i n_j와 같은 고차 곱입니다.

이들 자체만으로는 점유수 구성 간에 인구를 이동시키지 않습니다. 이미 존재하는 구성의 위상에만 작용할 뿐입니다. 하지만 이들이 아무 영향도 없는 것은 아닙니다. 이러한 상대 위상은 회로 뒷부분의 여기(excitation) 항이 만들어내는 간섭에 영향을 주기 때문에, 이를 필터링하면 실제로 생성되는 전개가 달라지고 샘플링 분포도 바뀔 수 있습니다. 이는 샘플링된 분포를 그대로 유지하는 단계가 아니라, 샘플링을 여기 항에 집중시키기 위해 회로 생성 단계에서 의도적으로 도입한 근사입니다. qDRIFT 수렴 보장을 그대로 유지하는 위의 대칭성 그룹화와는 달리, 이 필터는 전개되는 연산자 자체를 변경합니다. 따라서 회로는 더 이상 전체 해밀토니안 하의 전개를 근사하지 않으며, qDRIFT 오차 범위는 원래 연산자가 아닌 필터링된 연산자에 적용됩니다. 이는 회로가 구성을 제안하는 데 사용되는 샘플링 휴리스틱에 불과하기 때문에 여기서는 허용됩니다. 필터는 회로를 구성하는 데 사용되는 해밀토니안에만 적용되고, 이후의 고전적 대각화는 대각 항을 포함한 전체 해밀토니안을 사용하므로 에너지 추정치 자체에서 어떤 항도 손실되지 않습니다. SQD의 정확도는 이 고전적 단계에 달려 있으며, 이 단계는 구성이 어떻게 제안되었는지와 무관하게 샘플링된 부분공간 내에서 변분적으로 남아 있습니다.

filter_diagonal_terms() 함수는 연산자에서 이러한 항을 제자리에서 제거합니다. 이 함수는 정규 순서화된 구조, 즉 생성 모드의 다중집합이 소멸 모드의 다중집합과 일치하는지를 통해 이러한 항을 식별하므로, 이미 정규 순서화된 연산자에만 유효합니다. 이 가정은 실행 시점에 검사되지 않습니다.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))
5060

해밀토니안의 항을 그룹화했으므로, 이제 회로 앙상블을 생성하기 위해 다음 매개변수를 결정합니다.

  • 생성할 회로 수: num_circuits
  • 여기 그룹 기준으로 본 각 회로의 길이: num_exc
  • 서로 다른 전개 시간에 대한 인자: times

페르미온 회로 생성

이제 각 시간 단계에 대해 페르미온 회로를 생성합니다. 각 회로는 앞서 선언한 전개 시간을 갖는 단일 전개 게이트로 구성됩니다. 전개 연산자는 해밀토니안입니다. 이후 이 회로들에 Transpiler 패스를 실행하여 qDRIFT 회로를 생성합니다.

Ansatz 준비

InitializeModes 클래스를 사용하여 Hartree-Fock 상태를 준비합니다. 질소의 경우, 이 과정은 단순히 처음 num_elec_a개의 큐비트에 X 게이트를 적용한 다음 num_elec_b개의 큐비트에 X 게이트를 적용하는 것으로, 질소에서는 두 값 모두 7입니다. 이 상태는 질소의 일곱 개 α\alpha 전자와 일곱 개 β\beta 전자를 나타냅니다.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

2단계: 양자 하드웨어 실행을 위한 문제 최적화​

회로가 준비되었으므로, 먼저 qiskit-fermions에서 제공되는 패스를 사용하여 페르미온 수준의 최적화를 수행한 다음, 선택한 Backend에 맞게 회로를 Transpiler로 변환합니다. 이는 시뮬레이터 실험이므로 먼저 AerSimulator에 대해 이를 수행합니다. 각 그룹의 가중치 계산

이 단계에서는 해밀토니안 내 계수에 비례하는 확률로 항을 확률적으로 qDRIFT 샘플링합니다. qDRIFT Transpiler 패스가 이를 대신 수행합니다. 이제 해밀토니안이 장거리 결합 및 2차보다 높은 차수의 항을 포함하더라도, 제한된 큐비트 연결성에도 불구하고 하드웨어에서 더 효율적으로 실행할 수 있는 더 얕은 회로를 생성할 수 있습니다. 항 그룹화 후, 가중치를 기반으로 연산자를 샘플링합니다. 각 연산자 hih_i에 대해 가중치 WhiW_{h_i}는 다음과 같이 정의됩니다.

Whi=∣ci∣/λW_{h_i} = |c_i| / \lambda

페르미온 및 하드웨어 네이티브 최적화

generate_preset_jw_pass_manager() 함수는 FermionicCircuit을 입력받아 하드웨어에서 실행하기 위해 Transpiler로 변환할 수 있는 최적화된 최종 회로를 생성하는 MultiStagePassManager를 반환합니다. 기본 최적화 단계를 QDriftTrotterization 패스를 포함하는 FermionicPassManager로 대체합니다.

  • QDriftTrotterization 패스는 샘플링에 사용할 회로를 생성하기 위해 내부적으로 가중치 계산 및 샘플링을 사용합니다

  • RelabelModes 패스는 큐비트 간 연결성을 최적화하고 게이트 깊이를 줄이기 위해 페르미온 모드를 순열할 수 있는 또 다른 최적화 패스입니다. 자세한 내용은 API 참조를 참고하세요

MultiStagePassManager의 나머지 단계는 자동으로 실행되며 페르미온-큐비트 매핑 전체를 처리합니다.

  • F2QLayout: 프리셋 PassManager는 nn개의 페르미온 비트를 nn개의 큐비트에 단순하게 매핑하는 TrivialF2QLayout 패스를 적용합니다.

  • F2QSynth: 페르미온 기반 회로 명령어를 큐비트 기반 명령어로 매핑하는 Transpiler 패스입니다.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))
400

페르미온 수준 최적화가 완료되었으므로, 이제 시뮬레이터에서 실행하기 위해 회로를 Transpiler로 변환할 수 있습니다.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

3단계: Qiskit 프리미티브를 사용하여 실행​

회로가 준비되었으므로, 이제 AerSimulator에서 Qiskit 프리미티브를 사용하여 회로를 실행할 수 있습니다. 서로 다른 회로들의 모든 카운트를 결합합니다. 이를 불리언 벡터로 변환한 다음 마지막으로 SQD로 후처리합니다.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")
Executing 400 circuits with 100 shots each...
400 length before post processing

4단계: 후처리 및 원하는 고전적 형식으로 결과 반환​

SQD에 비트스트링 사용하기

이제 선택된 비트스트링에 대해 대각화 방식을 실행하여 분자의 바닥 상태 에너지에 해당하는 최소 고유값을 찾을 수 있습니다. 콜백 함수를 생성하고, 초기 점유 상태를 선언하고, 매개변수를 설정한 다음 마지막으로 대각화 방식을 실행합니다. 콜백 함수는 각 반복에서 현재 반복 횟수와 현재 고유값 추정치를 출력하는 데 사용됩니다.

마지막으로, 바닥 상태 추정치를 얻기 위해 결과 에너지에 nuclear_repulsion_energy를 더합니다.

참고: 부분공간의 차원은 노이즈가 없는 시뮬레이터에서도 반복마다 고정되지 않습니다. 각 서브샘플은 서로 다른 구성 집합을 추출하고, 복구 단계가 반복 사이에 풀(pool)을 재구성하기 때문에, 보고되는 차원은 서브샘플마다 달라집니다. 노이즈 없는 샘플링 자체만으로는 선택된 부분공간의 차원을 고정하지 않습니다. 반면 하드웨어 실행은 체계적으로 더 큰 부분공간을 얻는 경향이 있는데, 이는 노이즈가 있는 샷이 입자 수 대칭성을 깨뜨리고 구성 복구가 이를 추가적인 기저 벡터로 바꾸기 때문입니다. 이러한 이유로, 하드웨어 섹션에서는 비트스트링을 가지치기(pruning)하는 또 다른 단계도 소개할 것입니다.

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")
40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

하드웨어 예제​

이 예제에서는 20개의 큐비트(10개의 공간 오비탈)를 사용합니다. 이 선택은 빠르게 실행되어야 하는 튜토리얼을 위한 편의상의 선택일 뿐, 이 방법에 대한 엄격한 상한선은 아닙니다.

고전적 단계의 비용은 큐비트 수에 의해 직접 결정되지 않습니다. SQD는 샘플링된 구성들이 생성하는 부분공간에 투영된 해밀토니안을 대각화하므로, 고전적 비용을 좌우하는 것은 이 선택된 부분공간의 차원입니다 — 이는 여기서 samples_per_batch, num_batches, 그리고 회로가 실제로 생성하는 서로 다른 구성의 수에 의해 결정되며, 투영된 해밀토니안을 적용하는 데 필요한 희소 선형대수 연산도 함께 작용합니다. 전체 CI 공간은 오비탈과 전자 수에 따라 조합적으로 증가하지만, 선택된 부분공간은 그중 작고 조정 가능한 일부이며, 우리는 그 크기를 직접 제어합니다. 결과적으로 큐비트 수와 고전적 난이도는 어느 정도 독립적으로 조절할 수 있습니다. 즉, 더 넓은 오비탈 공간을 적당한 부분공간으로 샘플링하는 것이 매우 큰 부분공간에 대해 대각화하는 더 작은 시스템보다 더 저렴할 수 있습니다.

따라서 실제로는 실행 가능한 시스템 크기가 원하는 정확도에 필요한 부분공간 차원과, 고유값 해석기가 사용할 수 있는 메모리 및 코어 수에 달려 있습니다. 일반적으로 더 넓은 오비탈 공간은 화학적 정확도에 도달하기 위해 더 큰 부분공간을 필요로 하며, 이것이 결국 분산 자원을 필요로 하게 만드는 이유입니다. 이 단계를 확장하는 방법에 대해서는 qiskit-addon-sqd-hpc를 참고하세요. 고정된 상한을 가정하기보다는, 반복에 걸쳐 보고되는 부분공간 차원과 에너지 수렴을 관찰하면서 에너지가 더 이상 개선되지 않거나 사용 가능한 메모리를 모두 소진할 때까지 부분공간 크기를 늘리는 것이 실용적인 접근 방식입니다.

참고: 하드웨어의 노이즈로 인한 샘플링 오차 때문에, 하드웨어 실행에서 대각화를 위해 생성되는 부분공간은 시뮬레이터를 사용할 때보다 더 커집니다. 이는 대각화하려는 부분공간의 차원을 늘리지만, SQD가 노이즈에 대해 강건하기 때문에 이 워크플로는 여전히 정확한 답을 제공합니다.

허위 문자열 가지치기

여기서는 추가 단계를 수행할지 선택할 수 있습니다. 회로 실행으로부터 모든 비트스트링을 얻은 후, SQD를 실행하기 전에 유효하지 않은 비트스트링을 걸러내거나, 가지치기 없이 진행할 수 있습니다. 하드웨어 실행에서는 일반적으로 가지치기를 건너뛰는 것이 더 낫습니다. 대칭성이 깨진 샷을 그대로 두면 구성 복구가 이를 유효한 구성으로 복원할 수 있어, 해당 샷을 완전히 버리는 대신 부분공간을 넓힐 수 있기 때문입니다.

질소는 일곱 개의 α\alpha 전자와 일곱 개의 β\beta 전자만 가질 수 있으므로, 출력의 앞쪽 절반과 뒤쪽 절반에서 1의 개수가 7개보다 많거나 적은 비트스트링은 버릴 수 있습니다. 비트스트링이 유효한지 확인하고, 유효하지 않으면 버리는 함수를 정의합니다. 허위 비트스트링을 걸러낸 후, 나머지는 대각화 방식으로 전달됩니다. 아래의 PRUNE 플래그를 사용하여 두 동작을 전환할 수 있습니다.

가지치기는 회로 수, 전개 시간 집합, 대각 항 필터링과 더불어 최종 부분공간을 형성하는 여러 선택 사항 중 하나일 뿐이라는 점을 기억하세요. 가지치기를 수행한 실행과 수행하지 않은 실행을 비교하는 것은 다른 모든 요소가 고정되어 있을 때만 의미가 있습니다. C++ 동반 자료는 이 점을 더 자세히 다루는데, 이는 복구가 아니라 후선택(postselect)을 수행하며 다른 매개변수들도 다르기 때문입니다.

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")
Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400
Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

다음 단계​

권장 사항

이 작업이 흥미로웠다면 다음 자료에도 관심이 있을 수 있습니다.

  • 페르미온 격자 모델의 샘플 기반 Krylov 양자 대각화 - 변분 Ansatz 대신 시간 전개 회로를 사용하는 관련 튜토리얼입니다.
  • 화학 해밀토니안의 샘플 기반 양자 대각화 - 양자 화학 시뮬레이션을 위한 국소 유니터리 클러스터 야스트로프(LUCJ) 회로를 구성하는 방법에 대한 튜토리얼입니다.
  • SqDRIFT 논문 - 이 튜토리얼의 기반이 되는 문헌입니다. (이 논문에서 논의된 일부 최적화는 현재 진행 중인 작업이며, 이 튜토리얼은 사용된 라이브러리의 발전에 따라 향후 변경될 수 있다는 점에 유의하세요.)
출처: IBM Quantum docs — 2026년 9월 22일 업데이트
doQumentation의 영어 버전 — 2026년 9월 22일 업데이트
사이트 또는 번역 문제?콘텐츠 문제? IBM Quantum docs에서 편집