주 콘텐츠로 건너뛰기

화학 해밀토니안의 샘플 기반 양자 대각화

사용량 추정: Heron r2 프로세서에서 1분 미만 (참고: 이는 추정치이며 실제 실행 시간은 다를 수 있습니다.)

학습 성과​

이 튜토리얼을 마친 후, 사용자는 다음을 이해할 수 있어야 합니다.

  • SQD Qiskit 애드온을 사용하여 양자 처리 장치(QPU)에서 샘플링된 비트스트링으로 분자 시스템의 기저 상태 에너지를 근사하는 방법.
  • ffsim을 사용하여 양자 화학 시뮬레이션을 위한 국소 유니터리 클러스터 야스트로프(LUCJ) Circuit을 구성하는 방법.

사전 요구 사항​

이 튜토리얼을 시작하기 전에 다음 주제에 익숙해지시길 권장합니다.

  • 양자 화학 및 2차 양자화
  • Sampler 프리미티브를 사용하여 양자 Circuit에서 샘플링하기

배경​

이 튜토리얼에서는 SQD Qiskit 애드온을 사용하여 샘플 기반 양자 대각화(SQD) 알고리즘을 구현함으로써, 잡음 있는 양자 샘플을 후처리하여 질소 분자 N2\text{N}_2의 평형 결합 길이에서의 기저 상태를 근사하는 방법을 소개합니다. 해당 소프트웨어에 대한 자세한 내용은 관련 문서와 시작을 위한 간단한 예제에서 확인할 수 있습니다.

이 튜토리얼은 양자 화학에 익숙한 사용자, 특히 분자의 기저 상태 에너지를 구하는 방법을 아는 분들에게 권장됩니다. 워크플로우에 대한 자세한 안내는 양자 대각화 알고리즘 강좌를 참조하세요.

SQD는 양자 및 분산 클래식 컴퓨팅을 함께 사용하여 양자 해밀토니안과 같은 양자 연산자의 고유값과 고유벡터를 찾는 기법입니다. 클래식 분산 컴퓨팅은 양자 프로세서에서 얻은 샘플을 처리하고, 샘플이 형성하는 부분 공간에서 목표 해밀토니안을 투영하여 대각화하는 데 사용됩니다. SQD 기반 워크플로우는 다음 단계로 구성됩니다.

  1. Circuit 앤사츠를 선택하고 양자 컴퓨터에서 기준 상태(이 경우 하트리-폭(Hartree-Fock) 상태)에 적용합니다.
  2. 결과 양자 상태에서 비트스트링을 샘플링합니다.
  3. 비트스트링에 대해 자기 일관적 구성 복구 절차를 실행하여 기저 상태 근사값을 구합니다.

SQD는 목표 고유상태가 희소할 때 잘 동작하는 것으로 알려져 있습니다. 즉, 파동 함수가 문제 크기에 따라 지수적으로 증가하지 않는 기저 상태 집합 S={∣x⟩}\mathcal{S} = \{|x\rangle \}에 의해 지지될 때 효과적입니다.

양자 화학​

분자 시스템의 해밀토니안은 다음과 같이 쓸 수 있습니다.

H^=∑prσhpr a^pσ†a^rσ+12∑prqsστhprqs a^pσ†a^qτ†a^sτa^rσ,\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},

여기서 hprh_{pr}과 hprqsh_{prqs}는 분자 적분이라고 불리는 복소수로, 컴퓨터 프로그램을 사용하여 분자의 사양으로부터 계산됩니다. 이 튜토리얼에서는 PySCF 소프트웨어 패키지를 사용하여 적분을 계산합니다.

분자 해밀토니안의 유도 과정에 대한 자세한 내용은 양자 화학 교재(예: Szabo와 Ostlund의 Modern Quantum Chemistry)를 참고하세요. 양자 화학 문제가 양자 컴퓨터에 어떻게 매핑되는지에 대한 고수준 설명은 Qiskit Global Summer School 2024의 강의 Mapping Problems to qubits를 확인하세요.

국소 유니터리 클러스터 야스트로프(LUCJ) 앤사츠​

SQD는 샘플을 추출하기 위한 양자 Circuit 앤사츠가 필요합니다. 이 튜토리얼에서는 물리적 동기와 하드웨어 친화성의 조합 덕분에 국소 유니터리 클러스터 야스트로프(LUCJ) 앤사츠를 사용합니다. 앤사츠 Circuit 구성에는 ffsim을 사용합니다.

LUCJ 앤사츠는 QPU의 제한된 qubit 연결성에 적응합니다. 스핀 오비탈은 앤사츠가 SWAP Gate를 사용한 라우팅 없이 구현될 수 있도록 Qubit에 매핑됩니다. IBM® 하드웨어는 헤비-헥스 격자 qubit 토폴로지를 사용하며, 이 경우 아래에 나타낸 "지그재그" 패턴을 채택할 수 있습니다. 이 패턴에서 같은 스핀을 가진 오비탈은 선형 토폴로지의 Qubit에 매핑되고(빨간색과 파란색 원), 다른 스핀의 오비탈 간 연결은 4번째 공간 오비탈마다 존재하며, 보조 qubit(보라색 원)에 의해 연결이 이루어집니다.

헤비-헥스 격자에서 LUCJ 앤사츠의 qubit 매핑 다이어그램

자기 일관적 구성 복구​

자기 일관적 구성 복구 절차는 잡음 있는 양자 샘플에서 가능한 한 많은 신호를 추출하도록 설계되었습니다. 분자 해밀토니안은 입자 수와 스핀 Z를 보존하므로, 이러한 대칭성을 보존하는 Circuit 앤사츠를 선택하는 것이 합리적입니다. 하트리-폭 상태에 적용하면, 잡음이 없는 경우 결과 상태는 고정된 입자 수와 스핀 Z를 가집니다. 따라서 이 상태에서 샘플링된 비트스트링의 스핀-α\alpha 및 스핀-β\beta 절반은 하트리-폭 상태와 동일한 해밍 무게를 가져야 합니다. 현재 양자 프로세서의 잡음으로 인해 일부 측정된 비트스트링은 이 특성을 위반합니다. 단순한 사후 선택 방법은 이러한 비트스트링을 버리겠지만, 이는 낭비적입니다. 왜냐하면 해당 비트스트링에도 일부 신호가 포함될 수 있기 때문입니다. 자기 일관적 복구 절차는 후처리에서 그 신호 일부를 복구하려고 시도합니다. 이 절차는 반복적이며, 기저 상태의 각 오비탈 평균 점유율 추정값을 입력으로 필요로 하고, 이 추정값은 먼저 원시 샘플에서 계산됩니다. 절차는 루프로 실행되며, 각 반복은 다음 단계로 구성됩니다.

  1. 지정된 대칭성을 위반하는 각 비트스트링에 대해, 비트스트링을 현재 평균 오비탈 점유율 추정값에 더 가깝게 만들도록 설계된 확률적 절차로 비트를 뒤집어 새로운 비트스트링을 얻습니다.
  2. 대칭성을 만족하는 이전 및 새 비트스트링을 모두 수집하고, 사전에 선택된 고정 크기의 부분 집합을 서브샘플링합니다.
  3. 각 비트스트링 부분 집합에 대해, 해당 기저 벡터가 생성하는 부분 공간으로 해밀토니안을 투영하고(이전 섹션에서 이 기저 벡터에 대한 설명 참조), 클래식 컴퓨터에서 투영된 해밀토니안의 기저 상태 추정값을 계산합니다.
  4. 가장 낮은 에너지를 가진 기저 상태 추정값으로 평균 오비탈 점유율 추정값을 업데이트합니다.

SQD 워크플로우 다이어그램​

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.75 이상 (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
import math

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.primitives import StatevectorSampler
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

소규모 시뮬레이터 예제​

이 튜토리얼에서는 평형 결합 거리 근처의 질소 분자 기저 상태 근사값을 구합니다. 먼저 소규모 STO-6G 기저 집합을 사용하여 실험을 시뮬레이션하고 올바르게 작동하는지 확인합니다.

1단계: 클래식 입력을 양자 문제로 매핑​

먼저 분자와 그 특성을 지정합니다.

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="sto-6g",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
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), norb)

# Compute exact energy using FCI
reference_energy = cas.run().e_tot

print(f"norb = {norb}")
print(f"nelec = {nelec}")
converged SCF energy = -108.464957764796
CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000
norb = 8
nelec = (5, 5)

LUCJ 앤사츠 Circuit을 구성하기 전에, 다음 코드 셀에서 먼저 CCSD 계산을 수행합니다. 이 계산에서 얻은 t1t_1 및 t2t_2 진폭은 앤사츠의 매개변수를 초기화하는 데 사용됩니다.

# 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) = -108.5933309085008  E_corr = -0.1283731437052354

이제 ffsim을 사용하여 앤사츠 Circuit을 만듭니다. 분자가 닫힌 껍질 하트리-폭 상태를 가지므로, UCJ 앤사츠의 스핀 균형 변형인 UCJOpSpinBalanced를 사용합니다. from_t_amplitudes 메서드에서 optimize=True를 설정하여 t2t_2 진폭의 "압축된" 이중 인수분해를 활성화합니다(자세한 내용은 ffsim 문서의 The local unitary cluster Jastrow (LUCJ) ansatz 참조).

LUCJ 앤사츠는 QPU의 사용 가능한 연결성에 적응하므로, 앤사츠를 생성하기 전에 QPU Backend를 초기화해야 합니다. 여기서는 헤비-헥스 커플링 맵과 LUCJ 앤사츠가 자연스럽게 분해되는 Gate 집합을 갖춘 범용 Backend를 생성합니다. 그런 다음 ffsim.qiskit.generate_lucj_pass_manager를 사용하여 LUCJ 앤사츠 배경 섹션에서 설명한 "지그재그" 레이아웃에 따라 주어진 Backend에 LUCJ 앤사츠를 트랜스파일하는 데 특화된 패스 매니저를 생성합니다. 이 함수는 점수 휴리스틱을 사용하여 선택한 레이아웃과 관련된 오류를 최소화하는데, 이는 실제 QPU나 잡음 모델이 있는 시뮬레이터를 사용할 때 중요합니다. 패스 매니저를 반환하는 것 외에도, 하드웨어에서 구현 가능한 알파-베타 커플링 쌍도 반환합니다. 일부 쌍을 구현할 수 없는 경우 경고를 출력합니다.

import warnings

from qiskit.transpiler import CouplingMap

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
coupling_map = CouplingMap.from_heavy_hex(3)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # 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),
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, 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(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

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

다음으로, 대상 하드웨어에 맞게 Circuit를 최적화합니다. 일반적으로 이 단계에서는 하드웨어 Backend와 해당 Backend용 패스 매니저를 초기화합니다. 그러나 LUCJ 앤사츠가 하드웨어 연결성에 적응하므로, 이미 이전 단계에서 이 작업을 수행했습니다. 이제 패스 매니저를 Circuit에 실행하여 QPU에서 직접 실행할 수 있는 ISA Circuit으로 트랜스파일하기만 하면 됩니다.

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")
Gate counts: OrderedDict({'xx_plus_yy': 86, 'p': 16, 'measure': 16, 'cp': 15, 'x': 10, 'swap': 2, 'barrier': 1})

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

하드웨어 실행을 위한 Circuit 최적화가 완료되면, 대상 하드웨어에서 실행하여 바닥 상태 에너지 추정을 위한 샘플을 수집할 준비가 됩니다. Circuit가 하나뿐이므로 Qiskit Runtime의 작업 실행 모드를 사용하여 Circuit를 실행합니다.

rng = np.random.default_rng()
sampler = StatevectorSampler(seed=rng)
job = sampler.run([isa_circuit], shots=100_000)
Warning: Trying to add QuantumRegister to a QuantumCircuit having a layout
primitive_result = job.result()
pub_result = primitive_result[0]

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

QPU 출력의 품질을 판단하는 유용한 지표는 반환된 유효한 구성의 수입니다. 유효한 구성은 올바른 입자 수와 스핀 Z를 가지며, 이는 비트스트링의 오른쪽 절반이 스핀-업 전자 수와 같은 해밍 무게를 가지고, 왼쪽 절반이 스핀-다운 전자 수와 같은 해밍 무게를 가져야 함을 의미합니다. 다음 셀은 유효한 샘플링된 구성의 비율을 계산합니다.

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
    )

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")
Fraction of sampled configurations that are valid: 1.0

잡음 없는 시뮬레이터에서 Circuit를 샘플링하기 때문에 모든 비트스트링이 유효합니다. 잡음 있는 QPU에서 실행하면 비율이 1보다 작아지지만, 다음 셀에서 계산한 균일 무작위 비트스트링에서 예상되는 비율보다는 크길 바랍니다.

expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)
Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

이제 diagonalize_fermionic_hamiltonian 함수를 사용하여 해밀토니안의 바닥 상태 에너지를 추정합니다. 이 함수는 자기 일관적 구성 복구 절차를 수행하여 노이즈가 있는 양자 샘플을 반복적으로 개선함으로써 에너지 추정값을 향상시킵니다. 나중에 분석을 위해 중간 결과를 저장할 수 있도록 콜백 함수를 전달합니다. diagonalize_fermionic_hamiltonian의 인수 설명은 API 문서를 참조하세요.

여기서는 diagonalize_fermionic_hamiltonian의 initial_occupancies 인수를 사용하여 오비탈 점유율의 초기 추정값으로 하트리-폭 구성을 지정합니다. 이 방법은 기저 상태가 하트리-폭 구성에 상당한 지지를 가지는 시스템에 적합하지만, 다른 상황에서는 적합하지 않을 수 있습니다. 이 경우 더 고급 계산 방법을 통해 더 나은 초기 추정값을 얻을 수 있습니다. initial_occupancies를 지정하면 유효한 구성이 샘플링되지 않은 경우에도 구성 복구를 실행할 수 있습니다. 이는 잡음 있는 QPU에서 대형 Circuit를 샘플링할 때 발생할 수 있습니다. 이 인수 없이는 유효한 구성이 제공되지 않으면 구성 복구가 실패하고 오류가 발생합니다.

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

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# 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,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    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,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")
Iteration 1
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Iteration 2
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Final energy: -108.59275573641656
Final energy error: 0.0032316145694579745

결과 시각화​

첫 번째 플롯은 이 시뮬레이션에서 첫 번째 반복 후 이미 정확한 답의 1 mH 이내에 있음을 보여줍니다 (화학적 정확도는 일반적으로 1 kcal/mol ≈\approx 1.6 mH로 허용됩니다). 그러나 이것은 소규모 시스템이며 샘플이 잡음 없이 생성되었기 때문에 구성 복구가 필요하지 않습니다. 잡음 있는 QPU에서 실행되는 더 큰 시스템의 경우, 여러 번의 구성 복구 반복이 필요할 수 있으며 최종 정확도가 더 낮을 수 있습니다. 일반적으로 구성 복구 반복 횟수를 늘리거나 배치당 샘플 수를 늘리면 에너지를 개선할 수 있습니다.

두 번째 플롯은 마지막 반복 후 각 공간 궤도의 평균 점유율을 보여줍니다. 스핀-업 전자와 스핀-다운 전자 모두 우리의 솔루션에서 높은 확률로 처음 다섯 개의 궤도를 점유하고 있음을 알 수 있습니다.

# 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 - reference_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

Output of the previous code cell

대규모 하드웨어 예제​

이제 실제 양자 하드웨어에서 더 큰 예제를 실행합니다. 여기서는 cc-pVDZ 기저 집합에서 질소 분자의 활성 공간을 도출합니다.

1-4단계​

여기서는 모든 단계를 더 큰 규모의 단일 워크플로우로 통합하여 실제 양자 하드웨어에서 실행합니다.

# ------------------------------ Step 1 ------------------------------
# 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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
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), norb)

# Store reference energy from SCI calculation performed separately
reference_energy = -109.22802921665716

print(f"norb = {norb}")
print(f"nelec = {nelec}")

# 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

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)
print(f"Using backend {backend.name}")

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # 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),
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, 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(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

# ------------------------------ Step 2 ------------------------------

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

# ------------------------------ Step 3 ------------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_SQD"]
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]

# ------------------------------ Step 4 ------------------------------

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")
expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)
# 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

# Use the Hartree-Fock configuration as an initial guess for the
# orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# 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 = []

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    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,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

# 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 - reference_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()
converged SCF energy = -108.929838385609
norb = 26
nelec = (5, 5)
E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450045
Using backend ibm_boston
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].
Removing interaction (24, 24) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].
Removing interaction (20, 20) from the end.
Gate counts: OrderedDict({'sx': 7039, 'rz': 6990, 'cz': 1858, 'x': 61, 'measure': 52, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.02124
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.13889134249762
		Subspace dimension: 120409
	Subsample 1
		Energy: -109.11785470455858
		Subspace dimension: 110889
	Subsample 2
		Energy: -109.13234360554011
		Subspace dimension: 130321
Iteration 2
	Subsample 0
		Energy: -109.16392179579177
		Subspace dimension: 223729
	Subsample 1
		Energy: -109.16281938332986
		Subspace dimension: 223729
	Subsample 2
		Energy: -109.16955816711932
		Subspace dimension: 233289
Iteration 3
	Subsample 0
		Energy: -109.17905772999075
		Subspace dimension: 324900
	Subsample 1
		Energy: -109.17532445048462
		Subspace dimension: 357604
	Subsample 2
		Energy: -109.1733168689756
		Subspace dimension: 348100
Iteration 4
	Subsample 0
		Energy: -109.18437778820451
		Subspace dimension: 474721
	Subsample 1
		Energy: -109.18450164209159
		Subspace dimension: 476100
	Subsample 2
		Energy: -109.18493571190754
		Subspace dimension: 487204
Iteration 5
	Subsample 0
		Energy: -109.18616522497996
		Subspace dimension: 622521
	Subsample 1
		Energy: -109.18652868888333
		Subspace dimension: 644809
	Subsample 2
		Energy: -109.18753326484406
		Subspace dimension: 585225
Final energy: -109.18753326484406
Final energy error: 0.040495951813099396

Output of the previous code cell

다음 단계​

권장 사항

이 작업에 관심이 있으시다면 다음 자료를 참고해 보세요.

출처: IBM Quantum docs — 2026년 6월 10일 업데이트
doQumentation의 영어 버전 — 2026년 6월 10일 업데이트
이 번역은 2026년 4월 9일의 영어 버전을 기반으로 합니다
사이트 또는 번역 문제?콘텐츠 문제? IBM Quantum docs에서 편집