화학 해밀토니안 에너지 추정을 위한 SQD

이번 레슨에서는 SQD를 적용하여 분자의 바닥 상태 에너지를 추정합니다.

특히, $4$단계 Qiskit 패턴 방식을 사용하여 다음 주제들을 다룹니다:

1. 1단계: 문제를 양자 회로와 연산자로 변환
   • $N_2$에 대한 분자 해밀토니안을 설정합니다.
   • 화학적 직관과 하드웨어 친화적인 로컬 유니터리 클러스터 재스트로(LUCJ) [1] 방법을 설명합니다.
2. 2단계: 타겟 하드웨어에 맞게 최적화
   • 하드웨어 실행을 위한 앤자츠의 게이트 수와 레이아웃을 최적화합니다.
3. 3단계: 타겟 하드웨어에서 실행
   • 최적화된 Circuit을 실제 QPU에서 실행하여 부분 공간의 샘플을 생성합니다.
4. 4단계: 결과 후처리
   • 자기 일관 구성 복구 루프 [2]를 소개합니다.
     • 입자 수에 대한 사전 지식과 가장 최근 반복에서 계산된 평균 궤도 점유율을 사용하여 비트스트링 샘플 전체를 후처리합니다.
     • 복구된 비트스트링에서 확률적으로 서브샘플 배치를 생성합니다.
     • 각 샘플링된 부분 공간에 대해 분자 해밀토니안을 투영하고 대각화합니다.
     • 모든 배치에서 찾은 최솟값 바닥 상태 에너지를 저장하고 평균 궤도 점유율을 업데이트합니다.

이 레슨 전반에 걸쳐 여러 소프트웨어 패키지를 사용합니다.

• PySCF: 분자를 정의하고 해밀토니안을 설정합니다.
• ffsim 패키지: LUCJ 앤자츠를 구성합니다.
• Qiskit: 하드웨어 실행을 위해 앤자츠를 트랜스파일합니다.
• Qiskit IBM Runtime: QPU에서 Circuit을 실행하고 샘플을 수집합니다.
• Qiskit addon SQD: 부분 공간 투영 및 행렬 대각화를 사용하여 구성 복구 및 바닥 상태 에너지를 추정합니다.

1. 문제를 양자 회로와 연산자로 변환

분자 해밀토니안

분자 해밀토니안은 다음과 같은 일반적인 형태를 갖습니다:

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \sum_{ \substack{prqs\\\sigma\tau} } \frac{(pr|qs)}{2} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma}$$

$\hat{a}^\dagger_{p\sigma}$/$\hat{a}_{p\sigma}$는 $p$번째 기저 집합 원소와 스핀 $\sigma$에 연관된 페르미온 생성/소멸 연산자입니다. $h_{pr}$과 $(pr|qs)$는 단체 및 이체 전자 적분입니다. pySCF를 사용하여 분자를 정의하고 6-31g 기저 집합에 대한 해밀토니안의 단체 및 이체 적분을 계산합니다.

# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy pyscf qiskit qiskit-addon-sqd qiskit-ibm-runtime

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# 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)]],  # Two N atoms 1 angstrom apart
    basis="6-31g",
    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)  # hcore: one-body integrals
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)  # eri: two-body integrals

# Compute exact energy for comparison
exact_energy = cas.run().e_tot

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

이 레슨에서는 조던-위그너(JW) 변환을 사용하여 페르미온 파동 함수를 큐비트 파동 함수로 변환하고, 이를 양자 회로로 준비합니다. JW 변환은 M개의 공간 궤도에서 페르미온의 폭 공간을 2M 큐비트의 힐베르트 공간으로 변환합니다. 즉, 공간 궤도 하나가 두 개의 _스핀 궤도_로 분리되는데, 하나는 스핀 업($\alpha$) 전자에, 다른 하나는 스핀 다운($\beta$) 전자에 연관됩니다. 스핀 궤도는 점유 또는 비점유 상태일 수 있습니다. 일반적으로 궤도 수를 언급할 때는 공간 궤도 수를 사용하며, 스핀 궤도 수는 그 두 배입니다. 양자 회로에서는 각 스핀 궤도를 Qubit 하나로 표현합니다. 따라서 일부 Qubit 집합은 스핀 업 또는 $\alpha$-궤도를 나타내고, 다른 집합은 스핀 다운 또는 $\beta$-궤도를 나타냅니다. 예를 들어, 6-31g 기저 집합에서 $N_2$ 분자는 $16$개의 공간 궤도(즉, $16$ $\alpha$ + $16$ $\beta$ = $32$개의 스핀 궤도)를 가집니다. 따라서 $32$-Qubit 양자 Circuit이 필요합니다(나중에 설명하는 것처럼 추가 앤실라 Qubit이 필요할 수 있습니다). Qubit은 계산 기저에서 측정되어 비트스트링을 생성하며, 이는 전자 구성 또는 (슬레이터) 행렬식을 나타냅니다. 이 레슨 전반에 걸쳐 비트스트링, 구성, 행렬식이라는 용어를 서로 바꿔 사용합니다. 비트스트링은 스핀 궤도의 전자 점유율을 나타냅니다: 비트 위치의 $1$은 해당 스핀 궤도가 점유됨을 의미하고, $0$은 스핀 궤도가 비어 있음을 의미합니다. 전자 구조 문제는 입자 수를 보존하므로, 정해진 수의 스핀 궤도만 점유되어야 합니다. $N_2$ 분자는 스핀 업($\alpha$) 전자 $5$개와 스핀 다운($\beta$) 전자 $5$개를 가집니다. 따라서 $N_2$ 분자의 $\alpha$ 및 $\beta$ 궤도를 나타내는 어떤 비트스트링이든 각각 다섯 개의 $1$을 포함해야 합니다.

1.1 샘플 생성을 위한 양자 Circuit: LUCJ 앤자츠

이 레슨에서는 양자 상태 준비 및 후속 샘플링을 위해 로컬 유니터리 결합 클러스터 재스트로(LUCJ) \[1\] 앤자츠를 사용합니다. 먼저, 전체 UCJ 앤자츠의 다양한 구성 요소와 로컬 버전에서 적용된 근사를 설명합니다. 다음으로 ffsim 패키지를 사용하여 LUCJ 앤자츠를 구성하고 Qiskit Transpiler로 하드웨어 실행을 위해 최적화합니다.

UCJ 앤자츠는 UCJ 연산자의 $L$개의 레이어 또는 반복의 곱에 대해 다음과 같은 형태를 갖습니다.

$$|\psi\rangle = \prod_{\mu=1}^{L}{(e^{K^{\mu}} \times {e^{iJ^{\mu}}} \times {e^{-K^{\mu}}})} |\Phi_{0}\rangle$$

여기서 $\vert \Phi_{0} \rangle$은 일반적으로 하트리-폭(HF) 상태로 취해지는 기준 상태입니다. 하트리-폭 상태는 가장 낮은 번호의 궤도가 점유된 상태로 정의되므로, HF 상태 준비는 점유된 궤도에 해당하는 Qubit을 1로 설정하기 위해 X Gate를 적용하는 것을 포함합니다. 예를 들어, 4개의 공간 궤도와 스핀 업 2개, 스핀 다운 2개의 경우 HF 상태 준비 블록은 다음과 같은 형태일 수 있습니다: UCJ 연산자 ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$의 단일 반복은 궤도 회전($e^{K^{(\mu)}}$과 $e^{-K^{(\mu)}}$)으로 둘러싸인 대각 쿨롱 진화($e^{iJ^{(\mu)}}$)로 구성됩니다. 궤도 회전 블록은 단일 스핀 종($\alpha$ (스핀 업)/$\beta$ (스핀 다운))에 작용합니다. 각 전자 종에 대해, 궤도 회전은 단일 Qubit $R_{z}$ Gate 레이어와 그 뒤를 잇는 2-Qubit 기번스 회전 Gate($XX + YY$ Gates)의 시퀀스로 구성됩니다.

2-Qubit Gate는 인접한 스핀 궤도(가장 가까운 이웃 Qubit)에 작용하므로, SWAP Gate 없이 IBM® QPU에서 구현할 수 있습니다. 대각 쿨롱 연산자라고도 알려진 $e^{iJ^{(\mu)}}$는 세 개의 블록으로 구성됩니다. 그 중 두 개는 같은 스핀 섹터에 작용($e^{iJ_{\alpha \alpha}^{(\mu)}}$와 $e^{iJ_{\beta \beta}^{(\mu)}}$)하고, 하나는 두 스핀 섹터 사이에 작용($e^{iJ_{\alpha \beta}^{(\mu)}}$)합니다.

$e^{iJ^{(\mu)}}$의 모든 블록은 수-수 Gate $U_{nn}(\phi)$ [1]로 구성됩니다. $U_{nn}(\phi)$ Gate는 $R_{ZZ}(\frac{\phi}{2})$ Gate와 두 개의 별도 Qubit에 작용하는 두 개의 단일 Qubit $Rz(-\frac{\phi}{2})$ Gate로 더 분해될 수 있습니다.

같은 스핀 성분($J_{\alpha \alpha}$와 $J_{\beta \beta}$)은 모든 가능한 Qubit 쌍 사이에 $U_{nn}$ Gate를 가집니다. 그러나 초전도 QPU는 연결성이 제한적이므로, 인접하지 않은 Qubit 사이의 Gate를 구현하려면 Qubit을 스왑해야 합니다.

예를 들어, $N = 4$개의 공간 궤도에 대한 다음 $e^{iJ_{\alpha \alpha}^{(\mu)}}$(또는 $e^{iJ_{\beta \beta}^{(\mu)}}$) 블록을 고려해 보세요. 선형 Qubit 연결성에서, 마지막 세 개의 Gate는 인접하지 않은 Qubit 사이에 작용하므로(예: Q0과 Q2는 직접 연결되어 있지 않음) 직접 구현할 수 없습니다. 따라서 인접하게 만들기 위해 SWAP Gate가 필요합니다(다음 그림은 $3$개의 SWAP Gate를 사용한 예를 보여줍니다). 다음으로, $J_{\alpha \beta}$는 다른 스핀 섹터에서 같은 인덱스의 궤도 사이에 Gate를 구현합니다(예: $0\alpha$와 $0\beta$ 사이). 마찬가지로, Qubit이 QPU에서 물리적으로 인접하지 않은 경우, 이러한 Gate도 SWAP을 필요로 합니다. 위의 논의에서, UCJ 앤자츠는 인접하지 않은 Qubit 상호작용으로 인해 SWAP Gate가 필요하기 때문에 하드웨어 실행에 몇 가지 어려움이 있습니다. UCJ 앤자츠의 로컬 변형인 LUCJ는 대각 쿨롱 연산자에서 일부 $U_{nn}$을 제거하여 이 문제를 해결합니다.

같은 전자 종 블록($J_{\alpha \alpha}$와 $J_{\beta \beta}$)에서, LUCJ 버전은 최근접 이웃 연결성과 호환되는 $U_{nn}$ Gate만 유지하고 인접하지 않은 Qubit 사이의 Gate를 제거합니다. 다음 그림은 인접하지 않은 Gate를 제거한 후의 LUCJ 블록을 보여줍니다. 다음으로, 다른 전자 종 사이에 작용하는 $J_{\alpha \beta}$ 블록의 LUCJ 버전은 장치 토폴로지에 따라 다른 형태를 취할 수 있습니다.

여기서도 LUCJ 버전은 호환되지 않는 Gate를 제거합니다. 아래 그림은 격자, 육각형, 헤비-헥스, 선형 등 다양한 Qubit 토폴로지에 대한 $J_{\alpha \beta}$ 블록의 변형을 보여줍니다.

• 정사각형: 모든 $\alpha$와 $\beta$ 궤도 사이에 SWAP 없이 $U_{nn}$ Gate를 가질 수 있으므로, $U_{nn}$ Gate를 제거할 필요가 없습니다.
• 헤비-헥스: $\alpha$-$\beta$ 상호작용은 매 $4$번째 인덱스(예: 0번째, 4번째, 8번째) 스핀 궤도 사이에 유지되며 앤실라 매개 방식, 즉 $\alpha$와 $\beta$ 궤도를 나타내는 선형 체인 사이에 앤실라 Qubit이 필요합니다. 이 배치는 제한된 수의 SWAP만 필요로 합니다.
• 육각형: $\alpha$와 $\beta$가 두 인접한 선형 체인에 배치될 때, 0번째, 2번째, 4번째와 같이 하나 건너 하나씩의 궤도가 최근접 이웃이 됩니다.
• 선형: 하나의 $\alpha$와 하나의 $\beta$ 궤도만 연결되므로, $J_{\alpha \beta}$ 블록에는 Gate가 하나만 있습니다. LUCJ 버전을 구성하기 위해 UCJ 앤자츠에서 Gate를 제거하면 하드웨어 호환성이 향상되지만, 앤자츠의 표현력이 줄어듭니다. 따라서 LUCJ 앤자츠를 사용할 때는 수정된 UCJ 연산자의 반복 횟수($L$)를 더 늘려야 할 수 있습니다.

1.2 LUCJ 앤자츠 초기화

LUCJ는 매개변수화된 앤자츠이며, 하드웨어 실행 전에 매개변수를 초기화해야 합니다. 앤자츠를 초기화하는 한 가지 방법은 고전적인 결합 클러스터 단일 및 이중(CCSD) 방법의 t1 및 t2 진폭을 사용하는 것입니다. 여기서 t1 진폭은 단일 여기 연산자의 계수이고, t2 진폭은 이중 여기 연산자에 대한 것입니다.

LUCJ 앤자츠를 t1 및 t2 진폭으로 초기화하면 적절한 결과를 생성하지만, 앤자츠 매개변수에 추가 최적화가 필요할 수 있습니다.

# 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]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 ffsim을 사용한 LUCJ 앤자츠 구성

위에서 계산한 t1 및 t2 진폭으로 앤자츠를 생성하고 초기화하기 위해 ffsim 패키지를 사용합니다. 우리 분자는 닫힌 껍질 하트리-폭 상태를 가지므로, UCJ 앤자츠의 스핀 균형 변형인 UCJOpSpinBalanced를 사용합니다.

IBM 하드웨어는 헤비-헥스 토폴로지를 가지므로, Qubit 상호작용을 위해 [1]에서 사용되고 위에서 설명한 지그재그 패턴을 채택합니다. 이 패턴에서 같은 스핀을 가진 궤도(Qubit)는 선형 토폴로지로 연결됩니다(빨간색과 파란색 원). 헤비-헥스 토폴로지로 인해, 다른 스핀의 궤도는 매 4번째 궤도, 즉 0번째, 4번째, 8번째 등 사이에 연결됩니다(보라색 원).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
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),
)

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()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

반복 레이어가 있는 LUCJ 앤자츠는 일부 인접 블록을 병합하여 최적화할 수 있습니다. n_reps=2의 경우를 생각해 보세요. 중간에 있는 두 개의 궤도 회전 블록은 단일 궤도 회전 블록으로 병합될 수 있습니다. ffsim 패키지에는 이러한 인접 블록을 병합하여 Circuit을 최적화하는 ffsim.qiskit.PRE_INIT이라는 패스 매니저가 있습니다.

2. 타겟 하드웨어에 맞게 최적화

먼저, 원하는 Backend를 가져옵니다. 이 Backend에 맞게 Circuit을 최적화하고, 최적화된 Circuit을 동일한 Backend에서 실행하여 부분 공간에 대한 샘플을 생성합니다.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

다음으로, 앤자츠를 최적화하여 하드웨어 호환 가능하게 만드는 다음 단계를 권장합니다.

• 위에서 설명한 지그재그 패턴(앤실라 Qubit을 중간에 두고 두 개의 선형 체인)을 따르는 타겟 하드웨어에서 물리적 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 qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]

initial_layout = spin_a_layout + spin_b_layout

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': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. 타겟 하드웨어에서 실행

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

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. 결과 후처리

SQD 워크플로의 후처리 단계는 다음 다이어그램으로 요약할 수 있습니다.

계산 기저에서 LUCJ 앤자츠를 샘플링하면 잡음이 있는 구성(configuration)의 풀인 $\tilde{\mathcal{\chi}}$가 생성되며, 이는 후처리 루틴에 사용됩니다. 이 루틴에는 _구성 복구(configuration recovery)_라는 방법(세부 사항은 나중에 설명)이 포함되어 있으며, 전자 수가 올바르지 않은 구성을 확률적으로 보정합니다. 올바른 전자 수를 가진 구성인 $\tilde{\mathcal{\chi}}_{R}$만 부분 샘플링되어, 각 고유 구성의 출현 빈도를 기반으로 여러 배치에 분산됩니다. 각 샘플 배치는 부분 공간($\mathcal{S^{(k)}}$)을 정의합니다. 다음으로 분자 해밀토니안 $H$가 부분 공간에 투영됩니다:

$$H_{\mathcal{S}^{(k)}} = P_{\mathcal{S}^{(k)}} H _{\mathcal{S}^{(k)}} \text{ with } P_{\mathcal{S}^{(k)}} = \sum_{x \in \mathcal{S}^{(k)}} \vert x \rangle \langle x \vert$$

각 투영된 해밀토니안 $H_{\mathcal{S}^{(k)}}$는 고유값 풀이기(Eigensolver)에 입력되어 대각화를 통해 고유값과 고유벡터를 계산하고 고유 상태를 재구성합니다. 이 레슨에서는 PySCF의 Davidson 방법을 사용하여 대각화를 수행하는 qiskit-addon-sqd 패키지로 해밀토니안을 투영하고 대각화합니다.

$$H_{\mathcal{S}^{(k)}} \vert \psi^{(k)} \rangle = E^{(k)} \vert \psi^{(k)} \rangle$$

그런 다음 배치에서 최저 고유값(에너지)을 수집하고, 평균 궤도 점유율 $\text{n}$도 계산합니다. 평균 점유율 정보는 구성 복구 단계에서 잡음 구성을 확률적으로 보정하는 데 사용됩니다.

다음으로 자기 일관 구성 복구 루프를 자세히 설명하고, $N_2$ 해밀토니안의 기저 상태 에너지를 추정하기 위한 위의 단계를 구현하는 구체적인 코드 예제를 제시합니다.

4.1 구성 복구: 개요

비트스트링(Slater 행렬식)의 각 비트는 스핀 궤도를 나타냅니다. 비트스트링의 오른쪽 절반은 스핀업 궤도를, 왼쪽 절반은 스핀다운 궤도를 나타냅니다. 1은 해당 궤도가 전자로 채워져 있음을, 0은 비어 있음을 의미합니다. 올바른 입자 수(스핀업 전자와 스핀다운 전자 모두)는 사전에 알 수 있습니다. 비트스트링 $x$에 $N_x$개의 전자(즉, 비트스트링에 $N_x$개의 1이 있음)가 있다고 가정합니다. 올바른 입자 수는 $N$입니다. $N_x \neq N$이면, 비트스트링이 잡음에 의해 손상된 것임을 알 수 있습니다. 자기 일관 구성 루틴은 평균 궤도 점유율 정보를 활용하여 $|N_x - N|$개의 비트를 확률적으로 뒤집어 비트스트링을 보정하려 합니다. 평균 궤도 점유율($n$)은 궤도가 전자로 채워질 확률을 알려줍니다. $N_x < N$이면 전자가 부족하므로 일부 0을 1로 바꿔야 하며, 그 반대도 마찬가지입니다.

뒤집기 확률은 i번째 스핀 궤도에 대해 $|x[i] - avg\_occupancy[i]|$가 될 수 있습니다. [2]에서 저자들은 수정된 ReLU 함수를 사용한 가중 뒤집기 확률을 사용했습니다.

$$\begin{align} w(y) = \begin{cases} \delta \frac{y}{h} & \text{if } y \leq h\\ \nonumber \delta + (1 - \delta) \frac{y - h}{1 - h} & \text{if } y > h \end{cases} \end{align}$$

여기서 $h$는 ReLU 함수의 "꺾임" 위치를 정의하고, 매개변수 $\delta$는 꺾임 점에서 ReLU 함수의 값을 정의합니다. $\delta = 0$이면 $w$는 진짜 ReLU 함수가 되고, $\delta > 0$이면 수정된 ReLU가 됩니다. 논문에서 저자들은 $\delta = 0.01$과 $h =$ 알파(또는 베타) 입자 수/알파(또는 베타) 스핀 궤도 수 $= N/M$(채움 인수)를 사용했습니다.

평균 궤도 점유율($n$)은 사전에 알 수 없습니다. 기저 상태 추정의 첫 번째 반복은 두 스핀 종 모두에서 올바른 입자 수를 가진 구성만으로 시작합니다. 첫 번째 반복 후 기저 상태의 추정값을 얻고, 이 추정값을 이용해 $n$의 초기 추측값을 구성할 수 있습니다. 이 $n$의 추측값은 구성을 복구하고, 다음 기저 상태 추정 반복을 실행하여 자기 일관적으로 $n$의 추측값을 정제하는 데 사용됩니다. 이 과정은 중단 기준이 충족될 때까지 반복됩니다.

$N = 2$이고 $x = |1000\rangle$ ($N_x = 1$)인 다음 예제를 살펴보세요. 입자 수를 보정하기 위해 0 중 하나를 1로 뒤집어야 하며, 선택지는 1100, 1010, 1001입니다. 뒤집기 확률에 따라 선택지 중 하나가 복구된 구성(즉, 올바른 수의 입자를 가진 비트스트링)으로 선택됩니다.

첫 번째 반복에서 두 배치를 실행하고 추정된 기저 상태가 다음과 같다고 가정합니다:

$$\begin{align}\nonumber \text{Batch0: } \vert \psi \rangle &= 0.8 \times \vert 1001 \rangle + 0.6 \times \vert 0110 \rangle \\ \nonumber \text{Batch1: } \vert \psi \rangle &= \frac{1}{\sqrt{3}} \left( \vert 1001 \rangle + \vert 0101 \rangle + \vert 0110 \rangle \right) \nonumber \end{align}$$

계산 기저 상태와 그 진폭을 이용하여 각 스핀 궤도(Qubit)별 전자 점유 확률(줄여서 점유율)을 계산할 수 있습니다 (확률 = |진폭|$^2$). 아래는 추정 기저 상태에 나타나는 각 비트스트링의 Qubit 별 점유율과 배치 별 총 궤도 점유율을 나타낸 표입니다. Qiskit 순서 규칙에 따라 가장 오른쪽 비트가 Qubit-0(Q0)을, 가장 왼쪽 비트가 Q3을 나타냅니다.

점유율 (Batch0):

	Q3	Q2	Q1	Q0
1001	0.64	0.0	0.0	0.64
0110	0.0	0.36	0.36	0.0
n (Batch0)	0.64	0.36	0.36	0.64

점유율 (Batch1)

	Q3	Q2	Q1	Q0
1001	0.33	0.00	0.00	0.33
0101	0.0	0.33	0.00	0.33
0110	0.0	0.33	0.33	0.00
n (Batch1)	0.33	0.66	0.33	0.66

점유율 (배치 평균)

	Q3	Q2	Q1	Q0
n (Batch0)	0.64	0.36	0.36	0.64
n (Batch1)	0.33	0.66	0.33	0.66
n (average)	0.49	0.51	0.35	0.65

위에서 계산한 평균 궤도 점유율을 이용하여 구성 $x = \vert 1000 \rangle$에서 서로 다른 궤도에 대한 뒤집기 확률을 구할 수 있습니다. Q3이 나타내는 궤도는 이미 채워져 있어 뒤집을 필요가 없으므로 p(뒤집기)를 $0$으로 설정합니다. 나머지 비어 있는 궤도에 대한 뒤집기 확률은 각각 $\vert x[i] - \text{n}[i] \vert$입니다. p(뒤집기)와 함께, 위에서 설명한 수정된 ReLU 함수를 사용한 뒤집기 확률 가중치도 계산합니다.

뒤집기 확률 ($x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$)

	Q3	Q2	Q1	Q0
p(뒤집기) ($\vert x[i] - \text{n}[i] \vert$)	0	0.51	0.35	0.65
w(p(뒤집기))	0	0.03	0.007	0.31

마지막으로, 위의 가중 확률을 이용하여 채워지지 않은 Q2, Q1, Q0 궤도 중 하나를 뒤집을 수 있습니다. 위의 값에 따르면 Q0이 가장 뒤집힐 가능성이 높으며, 가능한 복구 구성은 $\vert \text{1001} \rangle$이 됩니다. 완전한 자기 일관 구성 복구 과정은 다음과 같이 요약할 수 있습니다:

첫 번째 반복: 양자 컴퓨터가 생성한 비트스트링(구성 또는 Slater 행렬식)이 집합 $\widetilde{\chi}$를 형성한다고 가정합니다. 이 집합에는 각 스핀 섹터에서 올바른($\widetilde{\chi}_{correct}$) 입자 수와 올바르지 않은($\widetilde{\chi}_{incorrect}$) 입자 수를 가진 구성이 모두 포함됩니다.

1. ($\widetilde{\chi}_{correct}$)의 구성이 무작위로 샘플링되어 부분 공간 투영을 위한 벡터 배치 $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$를 생성합니다. 배치 수와 각 배치의 샘플 수는 사용자 정의 매개변수입니다. 배치당 샘플 수가 많을수록 부분 공간 차원이 커지고 대각화 계산 비용도 높아집니다. 반면, 샘플 수가 너무 적으면 기저 상태 지지 벡터를 놓쳐 잘못된 추정이 나올 수 있습니다.
2. 배치에 대해 고유 상태 풀이기(부분 공간 투영 및 대각화)를 실행하고 근사 고유 상태를 얻습니다. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. 근사 고유 상태에서 $n$의 초기 추측값을 구성합니다.

이후 반복:

1. $n$을 이용하여 $\widetilde{\chi}_{incorrect}$에서 잘못된 입자 수를 가진 구성을 보정합니다. 이를 $\widetilde{\chi}_{correct\_new}$라 하면, $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$가 올바른 입자 수를 가진 새로운 구성 집합을 형성합니다.
2. $\widetilde{\chi}_{R}$을 샘플링하여 배치 $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$를 생성합니다.
3. 고유 상태 풀이기가 새 배치로 실행되어 기저 상태의 새 추정값 $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$을 생성합니다.
4. 근사 고유 상태에서 $n$의 정제된 추측값을 구성합니다.
5. 중단 기준이 충족되지 않으면 2.1 단계로 돌아갑니다.

4.2 기저 상태 추정

먼저, 후처리를 위해 카운트를 비트스트링 행렬과 확률 배열로 변환합니다.

행렬의 각 행은 하나의 고유 비트스트링을 나타냅니다. Qiskit에서는 비트스트링의 오른쪽부터 Qubit 인덱스가 시작되므로, 열 0은 Qubit N-1을, 열 N-1은 Qubit 0을 나타내며, N은 Qubit 수입니다.

알파 궤도는 열 인덱스 범위 (N, N/2](오른쪽 절반)에, 베타 궤도는 열 범위 (N/2, 0](왼쪽 절반)에 표현됩니다.

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

이 기법에서 중요한 사용자 제어 옵션이 몇 가지 있습니다:

• iterations: 자기 일관 구성 복구 반복 횟수
• n_batches: 고유 상태 풀이기의 다양한 호출에 사용되는 구성 배치 수
• samples_per_batch: 각 배치에 포함할 고유 구성 수
• max_davidson_cycles: 각 고유값 풀이기가 실행하는 최대 Davidson 사이클 수

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine the full set of noisy configurations
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We post-select here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_davidson=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 결과 논의

첫 번째 플롯은 몇 번의 반복 후 기저 상태 에너지를 약 24 mH 이내로 추정함을 보여줍니다 (화학적 정확도는 일반적으로 1 kcal/mol $\approx$ 1.6 mH로 허용됩니다). 두 번째 플롯은 최종 반복 후 각 공간 궤도의 평균 점유율을 보여줍니다. 우리의 해에서 스핀업과 스핀다운 전자 모두 처음 다섯 궤도를 높은 확률로 점유하고 있음을 알 수 있습니다.

추정된 기저 상태 에너지는 괜찮지만, 화학적 정확도 한계($\pm \approx 1.6$ mH) 이내는 아닙니다. 이 차이는 투영 및 대각화에 사용한 작은 부분 공간 차원에 기인합니다. samples_per_batch=500을 사용했기 때문에 부분 공간이 최대 $500$개의 벡터로 구성되어 기저 상태 지지에 필요한 벡터가 누락되어 있습니다. samples_per_batch 매개변수를 늘리면 더 많은 고전 계산 자원과 실행 시간을 대가로 정확도가 향상될 수 있습니다.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

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

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

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

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Exact energy: -109.04667 Ha
SQD energy: -109.02234 Ha
Absolute error: 0.02434 Ha

독자를 위한 연습 문제

samples_per_batch 매개변수를 점진적으로 늘리며(예: $1000$에서 $10000$까지 $1000$ 단위로; 컴퓨터 메모리가 허용하는 범위 내에서) 추정된 기저 상태 에너지를 비교해 보세요.

References

[1] M. Motta et al., "Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure" (2023). Chem. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer" (2024). arXiv:quant-ph/2405.05068.