주 콘텐츠로 건너뛰기

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

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

학습 목표​

  • Trotter 분해에 비해 깊이가 더 작은 회로를 만드는 방법을 배웁니다

  • qDRIFT와 SQD를 사용한 바닥 상태 추정의 전체 워크플로를 따라가 봅니다

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

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

사전 지식​

배경​

SqDRIFT는 SKQD의 변형으로, 비트열을 샘플링할 앤사츠를 선택해야 하는 필요를 목표 해밀토니안에서 직접 구성한 시간 발전 회로의 앙상블로 대체합니다. 이는 해밀토니안의 계수에 따라 더 작은 시간 발전 연산자를 서브샘플링하여 이루어지며, 이를 qDRIFT Trotter 분해 방법이라고 합니다.

이 튜토리얼은 Qiskit Fermions를 사용해 qDRIFT 알고리즘을 위한 더 자연스러운 페르미온 회로를 만든 다음, 페르미온 레이아웃 및 합성 패스를 사용한 후 그 회로를 하드웨어 실행을 위한 기존 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까지 가고 kthk_{th} 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이므로 이는 정규화된 확률 분포이며, 무작위 추출에 대한 결과 채널의 기댓값은 nn이 커질수록 오차가 줄어들면서 HH에 따른 발전을 근사합니다. 근사 오차는 항의 수 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 — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")
# 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개 전자이며, 7개의 α\alpha 전자와 7개의 β\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 오비탈에 의존하므로, 다시 생성한 파일은 제공된 파일과 오비탈 위상이나 순서가 다를 수 있습니다. 총 에너지는 영향을 받지 않습니다.

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

먼저 pyscf에서 제공하는 cisolver로 기준 에너지를 구합니다. 이것은 우리가 다루는 분자의 실제 바닥 상태 에너지입니다. 이를 위해 먼저 각각 오비탈 수와 전자 수인 norb와 nelec를 선언합니다. 그런 다음 각각 1전자 및 2전자 적분인 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 = "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 = "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 샘플링 슬롯이 구성 사이에서 점유를 이동시키는 항에 쓰이도록 합니다. 이러한 항은 다음 단계에서 Evolution 게이트를 구성하기 전, 이 시점에 해밀토니안에서 걸러내는 것이 가장 좋습니다.

해당 항은 점유수 기저에서 대각인 항, 즉 수 연산자의 곱 ai†aia^\dagger_i a_i입니다. 이 설명에 해당하는 항은 세 종류입니다.

  • 상수 에너지 오프셋: 수 연산자가 0개인 곱이며, 그 시간 발전은 전역 위상만 기여합니다.

  • 개별 수 연산자 nin_i: 그 시간 발전은 단일 큐비트 ZZ 회전으로 환원됩니다.

  • ninjn_i n_j 같은 고차 곱.

이들 중 어느 것도 단독으로는 점유수 구성 사이에서 점유를 이동시키지 않으며, 이미 존재하는 구성의 위상에만 작용합니다. 그러나 이들이 아무 영향이 없는 것은 아닙니다. 이 상대 위상이 회로 뒤쪽의 들뜸 항이 만드는 간섭에 반영되므로, 이를 걸러내면 실제로 생성되는 발전이 바뀌고 샘플링 분포가 바뀔 수 있습니다. 이것은 샘플링 분포를 건드리지 않는 단계가 아니라, 들뜸 항에 샘플링을 집중하기 위해 회로 생성 단계에서 의도적으로 도입한 근사입니다. 위의 대칭 그룹화는 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 회로를 만듭니다.

앤사츠 준비

InitializeModes 클래스를 사용해 Hartree-Fock 상태를 준비합니다. 질소의 경우, 처음 num_elec_a개의 큐비트와 이어서 num_elec_b개의 큐비트에 X 게이트를 적용하기만 하면 되며, 질소에서는 둘 다 7입니다. 이 상태는 질소의 7개의 α\alpha 전자와 7개의 β\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에서 제공하는 패스로 페르미온 수준의 최적화를 수행한 다음, 선택한 백엔드에 맞게 회로를 트랜스파일합니다. 시뮬레이터 실험이므로 먼저 AerSimulator에 대해 수행합니다. 각 그룹의 가중치 계산

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

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

1단계에서 항을 그룹화했기 때문에, 여기서 각 hih_i는 하나의 그룹 전체입니다. cic_i는 그룹 ii에 속한 항들의 평균 절댓값 계수이며, 그룹의 각 항은 계수를 부호만 남긴 값으로 줄여서 시간 전개합니다.

페르미온 및 하드웨어 고유 최적화

generate_preset_jw_pass_manager() 함수는 FermionicCircuit을 받아 하드웨어에서 실행하도록 Transpile할 수 있는 최적화된 최종 회로를 만들어 주는 MultiStagePassManager를 반환합니다. 기본 최적화 단계를 QDriftTrotterization 패스가 포함된 FermionicPassManager로 교체합니다.

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

  • RelabelModes 패스는 페르미온 모드를 치환하여 큐비트 간 연결성을 최적화하고 게이트 깊이를 줄이는 데 사용할 수 있는 또 다른 최적화 패스입니다. 자세한 내용은 API 레퍼런스를 참고하세요

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

  • F2QLayout: 프리셋 패스 매니저는 TrivialF2QLayout 패스를 적용하며, 이 패스는 nn개의 페르미온 비트를 nn개의 큐비트에 그대로 대응시킵니다.

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

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

페르미온 수준의 최적화를 마쳤으므로, 이제 시뮬레이터에서 실행할 회로를 Transpile할 수 있습니다.

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를 더합니다.

참고: 잡음이 없는 시뮬레이터에서도 부분공간 차원은 반복마다 고정되지 않습니다. 각 부분 샘플이 서로 다른 구성 집합을 뽑고, 복구 단계가 반복 사이에 풀을 다시 구성하기 때문에 보고되는 차원이 부분 샘플마다 달라집니다. 잡음 없는 샘플링만으로는 선택된 부분공간의 차원이 정해지지 않습니다. 반면 하드웨어 실행에서는 대체로 부분공간이 체계적으로 더 커지는데, 잡음이 있는 샷이 입자 수 대칭성을 깨뜨리고 구성 복구가 이를 추가 기저 벡터로 바꾸기 때문입니다. 이 때문에 하드웨어 섹션에서는 비트열을 가지치기하는 단계도 추가로 소개합니다.

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 전자 7개와 β\beta 전자 7개만 가질 수 있으므로, 출력의 앞쪽 절반과 뒤쪽 절반에서 1의 개수가 7개보다 많거나 적은 비트열은 버릴 수 있습니다. 비트열이 유효한지 확인하고 유효하지 않으면 버리는 함수를 정의합니다. 불필요한 비트열을 걸러내면 나머지는 대각화 방식으로 전달됩니다. 아래의 PRUNE 플래그로 두 동작을 전환하세요.

가지치기는 회로 수, 시간 전개 값의 집합, 대각항 필터링과 함께 최종 부분공간을 결정하는 여러 선택 중 하나일 뿐임을 염두에 두세요. 가지치기를 한 실행과 하지 않은 실행을 비교하는 것은 나머지 모든 조건이 고정된 경우에만 의미가 있습니다. 이 튜토리얼의 C++ 버전은 복구하는 대신 사후 선택을 하고 그 밖의 매개변수도 다르기 때문에 이 내용을 더 자세히 다룹니다.

name = "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

다음 단계​

추천 자료

이 내용이 흥미로웠다면 다음 자료도 살펴보세요.

출처: IBM Quantum docs — 2026년 10월 2일 업데이트
doQumentation의 영어 버전 — 2026년 10월 2일 업데이트
사이트 또는 번역 문제?콘텐츠 문제? IBM Quantum docs에서 편집