화학 해밀토니언의 에너지 추정용 SQD

이 강의에서는 분자의 기저 상태 에너지를 추정하는 데 SQD를 적용해 보겠습니다.

특히 다음 주제는 $4$ -단계 키스킷 패턴 접근 방식을 사용하여 논의할 것입니다:

1. 1단계: 양자 회로 및 연산자에 문제 매핑하기
   • $N_2$ 에 대한 분자 해밀턴을 설정합니다.
   • 화학에서 영감을 받은 하드웨어 친화적인 로컬 유니티 클러스터 자스트로우(LUCJ) [1] 에 대해 설명하세요
2. 2단계: 대상 하드웨어에 최적화
   • 하드웨어 실행을 위한 앤사츠의 게이트 수 및 레이아웃 최적화
3. 3단계: 대상 하드웨어에서 실행
   • 실제 QPU에서 최적화된 회로를 실행하여 서브스페이스의 샘플을 생성합니다.
4. 4단계: 포스트 프로세스 결과
   • 자체 일관된 구성 복구 루프 도입 [2]
     • 가장 최근 반복에서 계산된 파티클 수와 평균 궤도 점유에 대한 사전 지식을 사용하여 전체 비트스트링 샘플 세트를 후처리합니다.
     • 복구된 비트스트링에서 하위 샘플 배치를 확률적으로 생성합니다.
     • 샘플링된 각 하위 공간에 분자 해밀턴을 투영하고 대각선화합니다.
     • 모든 배치에서 발견된 최소 지상 상태 에너지를 저장하고 평균 궤도 점유를 업데이트합니다.

강의 내내 여러 소프트웨어 패키지를 사용할 것입니다.

• PySCF 를 사용하여 분자를 정의하고 해밀턴을 설정합니다.
• ffsim 패키지를 사용하여 LUCJ ansatz를 구성합니다.
• Qiskit 를 사용하여 하드웨어 실행을 위한 트랜스파일링을 수행합니다.
• Qiskit IBM Runtime 를 사용하여 QPU에서 회로를 실행하고 샘플을 수집합니다.
• Qiskit addon SQD 하위 공간 투영 및 행렬 대각선화를 사용한 구성 복구 및 접지 상태 에너지 추정.

1. 문제를 양자 회로 및 연산자로 매핑한다

분자 해밀토니안

분자 해밀턴은 일반적인 형태를 취합니다:

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

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

Output:

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

이 강의에서는 조던-위그너(JW) 변환을 사용하여 페르미온 파동 함수를 큐비트 파동 함수에 매핑하여 양자 회로를 사용하여 준비할 수 있도록 합니다. JW 변환은 M개의 공간 오비탈을 갖는 페르미온의 포크 공간을 2M 개의 큐비트의 힐베르트 공간에 매핑합니다. 즉, 하나의 공간 오비탈이 두 개의 스핀 오비탈 로 나뉘는데, 하나는 스핀 업( $\alpha$ ) 전자와 연결되고 다른 하나는 스핀 다운( $\beta$ ) 전자와 연결됩니다. 스핀 오비탈은 점유되거나 비점유될 수 있습니다. 일반적으로 궤도 수를 언급할 때는 공간 궤도 수를 사용합니다. 스핀 궤도의 수는 두 배가 됩니다. 양자 회로에서는 각 스핀 궤도를 하나의 큐비트로 표현합니다. 따라서 큐비트 세트는 스핀업 또는 $\alpha$ -오비탈을 나타내고, 다른 세트는 스핀다운 또는 $\beta$ -오비탈을 나타냅니다. 예를 들어, 6-31g 기초 세트에 대한 $N_2$ 분자는 $16$ 공간 궤도를 갖습니다(즉, $16$ $\alpha$ + $16$ $\beta$ = $32$ 스핀 궤도). 따라서 $32$ -큐비트 양자 회로가 필요합니다(나중에 설명하는 것처럼 추가 앙실라 큐비트가 필요할 수도 있습니다). 큐비트는 전자 구성 또는 (슬레이터) 결정자를 나타내는 비트스트링을 생성하기 위해 계산적으로 측정됩니다. 이 단원에서는 비트스트링, 구성, 결정자라는 용어를 같은 의미로 사용할 것입니다. 비트 문자열은 스핀 궤도에서 전자의 점유를 알려줍니다. 비트 위치의 $1$ 은 해당 스핀 궤도가 점유되고 있음을 의미하며, $0$ 은 스핀 궤도가 비어 있음을 의미합니다. 전자 구조 문제는 입자 보존 문제이므로 고정된 수의 스핀 궤도만 점유해야 합니다. $N_2$ 분자는 $5$ 스핀업( $\alpha$ ) 및 $5$ 스핀다운( $\beta$ ) 전자를 가지고 있습니다. 따라서 $\alpha$ 및 $\beta$ 궤도를 나타내는 비트 문자열에는 $N_2$ 분자에 대해 각각 5개의 $1\text{s}$ 가 있어야 합니다.

1.1 양자 회로에 의한 샘플 생성: LUCJ 가정

이 강의에서는 양자 상태 준비와 후속 샘플링을 위해 국소 단일 결합 클러스터 자스트로우(LUCJ) [1] 안사츠를 사용하겠습니다. 먼저 전체 UCJ 안사츠의 다양한 구성 요소와 로컬 버전에서 만들어진 근사치를 설명하겠습니다. 다음으로 ffsim 패키지를 사용하여 LUCJ 안사츠를 구성하고 하드웨어 실행을 위해 키스킷 트랜스파일러를 사용하여 최적화합니다.

UCJ 안사츠의 형식은 다음과 같습니다( $L$ 레이어 또는 UCJ 연산자의 반복 제품의 경우)

여기서 $\vert \Phi_{0} \rangle$ 는 참조 상태이며, 일반적으로 하트트리-복(HF) 상태로 사용됩니다. 하트리-복 상태는 가장 낮은 번호의 궤도가 점유된 것으로 정의되므로, HF 상태 준비에는 점유된 궤도에 해당하는 큐비트를 1로 설정하기 위해 X 게이트를 적용합니다. 예를 들어, 4개의 공간 궤도와 2개의 업 스핀 및 2개의 다운 스핀에 대한 HF 상태 준비 블록은 다음과 같이 보일 수 있습니다:

8개의 큐비트를 나타내는

UCJ 연산자 ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ 의 단일 반복은 대각선 쿨롱 진화( $e^{iJ^{(\mu)}}$ )와 궤도 회전( $e^{K^{(\mu)}}$ 및 $e^{-K^{(\mu)}}$ )으로 구성됩니다.

궤도 회전 블록은 단일 스핀 종에서 작동합니다( $\alpha$ (업스핀)/ $\beta$ (다운스핀)). 각 전자 종에 대해 궤도 회전은 단일 큐비트 $R_{z}$ 게이트 층과 2큐비트 주어진 회전 게이트( $XX + YY$ 게이트) 시퀀스로 구성됩니다.

2큐비트 게이트는 인접한 스핀 오비탈(가장 가까운 이웃 큐비트)에 작용하므로 SWAP 게이트 없이 IBM® QPU에서 구현할 수 있습니다.

4개의 알파 궤도 큐비트와 4개의 베타 궤도 큐비트를 나타낸

대각선 쿨롱 연산자라고도 하는 $e^{iJ^{(\mu)}}$ 는 세 개의 블록으로 구성됩니다. 그 중 두 개는 동일한 스핀 섹터에서 작동하고( $e^{iJ_{\alpha \alpha}^{(\mu)}}$ 및 $e^{iJ_{\beta \beta}^{(\mu)}}$ ), 하나는 두 스핀 섹터 사이에서 작동합니다( $e^{iJ_{\alpha \beta}^{(\mu)}}$ ).

$e^{iJ^{(\mu)}}$ 의 모든 블록은 숫자-번호 게이트 $U_{nn}(\phi)$ [1] 로 구성됩니다. $U_{nn}(\phi)$ 게이트는 $R_{ZZ}(\frac{\phi}{2})$ 게이트와 두 개의 개별 큐비트에서 작동하는 두 개의 단일 큐비트 $Rz(-\frac{\phi}{2})$ 게이트로 더 세분화할 수 있습니다.

동일한 스핀 컴포넌트( $J_{\alpha \alpha}$ 및 $J_{\beta \beta}$ )는 가능한 모든 큐비트 쌍 사이에 $U_{nn}$ 게이트가 있습니다. 그러나 초전도 QPU는 연결성이 제한적이기 때문에 인접하지 않은 큐비트 간의 게이트를 구현하려면 큐비트를 교체해야 합니다.

예를 들어 $N = 4$ 공간 궤도에 대해 다음 $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (또는 $e^{iJ_{\beta \beta}^{(\mu)}}$ ) 블록을 고려합니다. 선형 큐비트 연결의 경우, 마지막 세 개의 게이트는 인접하지 않은 큐비트 사이에서 작동하므로 직접 구현할 수 없습니다(예: Q0 와 Q2 은 직접 연결되지 않음). 따라서 이들을 인접하게 만들기 위해 SWAP 게이트가 필요합니다(다음 그림은 $3$ SWAP 게이트를 사용한 예시).

다음으로 $J_{\alpha \beta}$ 은 서로 다른 스핀 섹터의 동일한 인덱싱된 오비탈 간에 게이트를 구현합니다(예: $0\alpha$ 과 $0\beta$ 사이). 마찬가지로 큐비트가 QPU에서 물리적으로 인접하지 않은 경우 이러한 게이트에도 SWAP이 필요합니다.

위의 논의에서 알 수 있듯이 UCJ 안사츠는 인접하지 않은 큐비트 상호작용으로 인해 SWAP 게이트가 필요하기 때문에 HW 실행에 몇 가지 장애물에 직면해 있습니다. UCJ ansatz의 로컬 변형인 LUCJ는 대각선 쿨롱 연산자에서 $U_{nn}$ 일부를 제거하여 이 문제를 해결합니다.

동일한 전자 종 블록( $J_{\alpha \alpha}$ 및 $J_{\beta \beta}$ )에서 가장 가까운 이웃 연결과 호환되는 $U_{nn}$ 게이트만 유지하고 인접하지 않은 큐비트 사이의 게이트를 제거한 LUCJ 버전입니다. 다음 그림은 인접하지 않은 게이트를 제거한 후의 LUCJ 블록을 보여줍니다.

다음으로, 서로 다른 전자 종 사이에서 작동하는 $J_{\alpha \beta}$ 블록의 LUCJ 버전은 디바이스 토폴로지에 따라 다른 모양을 취할 수 있습니다.

여기서도 LUCJ 버전은 호환되지 않는 게이트를 제거합니다. 아래 그림은 그리드, 육각형, 무거운 육각형, 선형 등 다양한 큐비트 토폴로지에 대한 $J_{\alpha \beta}$ 블록의 변형을 보여줍니다.

• 정사각형 : 모든 $\alpha$ 및 $\beta$ 궤도 사이에 $U_{nn}$ 게이트를 스왑하지 않고도 가질 수 있으므로 $U_{nn}$ 게이트를 제거할 필요가 없습니다.
• 헤비헥스 : $\alpha$ - $\beta$ 상호작용은 모든 $4$ -번째 인덱스(예: 0번째, 4번째, 8번째) 스핀 궤도 사이에 유지되며, $\alpha$ 및 $\beta$ 궤도를 나타내는 선형 체인 사이에 안실라 큐비트가 필요한 안실라 매개체입니다. 이 방식은 제한된 수의 스왑이 필요합니다.
• 육각형 : 0번째, 2번째, 4번째 색인 궤도와 같은 다른 모든 궤도는 $\alpha$ 및 $\beta$ 을 두 개의 인접한 선형 체인으로 배치할 때 가장 가까운 이웃이 됩니다.
• 선형 : $\alpha$ 궤도와 $\beta$ 궤도가 하나만 연결되어 있으므로 $J_{\alpha \beta}$ 블록에는 게이트가 하나만 있습니다.

다양한 큐비트 배열에 대한

UCJ 안사츠에서 게이트를 제거하여 LUCJ 버전을 구성하면 HW 호환성이 높아지지만, 안사츠의 표현력이 다소 떨어집니다. 따라서 LUCJ 안사츠를 사용할 때는 수정된 UCJ 연산자의 반복( $L$ )이 더 필요할 수 있습니다.

1.2 LUCJ 접근법 초기화

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

t1 및 t2 진폭으로 LUCJ ansatz를 초기화하면 괜찮은 결과를 얻을 수 있지만, ansatz 매개변수는 추가 최적화가 필요할 수 있습니다.

# 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

Output:

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 LUCJ 접근법을 구축하기 위해 ffsim

Ffsim 패키지를 사용하여 위에서 계산한 t1 및 t2 진폭으로 ansatz를 생성하고 초기화합니다. 우리 분자는 닫힌 쉘 하트리-복크 상태이므로, 스핀 균형이 잡힌 UCJ 안사츠의 변형을 사용합니다, UCJOpSpinBalanced.

IBM 하드웨어는 16진수 토폴로지를 가지고 있으므로, [1] 에서 사용하고 위에서 설명한 지그재그 패턴을 큐비트 상호작용에 적용합니다. 이 패턴에서 동일한 스핀을 가진 궤도(큐비트)는 선 토폴로지(빨간색과 파란색 원)로 연결됩니다. 무거운 16진수 토폴로지로 인해 서로 다른 스핀의 궤도는 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 패키지에는 이러한 인접 블록을 병합하여 회로를 최적화하는 ffsim.qiskit.PRE_INIT 이라는 패스 관리자가 있습니다.

2. 대상 하드웨어에 최적화

먼저 원하는 백엔드를 가져옵니다. 백엔드에 맞게 회로를 최적화한 다음, 동일한 백엔드에서 최적화된 회로를 실행하여 하위 공간에 대한 샘플을 생성합니다.

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

다음으로, 다음 단계를 통해 ansatz를 최적화하고 하드웨어와 호환되도록 하는 것이 좋습니다.

• 위에서 설명한 지그재그 패턴(두 개의 선형 체인 사이에 앙실라 큐비트가 있는 두 개의 선형 체인)을 준수하는 대상 하드웨어에서 물리적 큐비트(initial_layout)를 선택합니다. 이 패턴으로 큐비트를 배치하면 더 적은 게이트로 효율적인 하드웨어 호환 회로를 만들 수 있습니다.
• 키스킷의 generate_preset_pass_manager 함수를 사용하여 단계별 패스 관리자 생성 backend 와 initial_layout.
• 스테이지 패스 관리자의 pre_init 단계를 ffsim.qiskit.PRE_INIT 으로 설정합니다. ffsim.qiskit.PRE_INIT 에는 게이트를 궤도 회전으로 분해한 다음 궤도 회전을 병합하여 최종 회로에서 더 적은 수의 게이트를 생성하는 키스킷 트랜스파일러 패스가 포함되어 있습니다.
• 서킷에서 패스 관리자를 실행하세요.

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

Output:

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. 대상 하드웨어에서 실행

하드웨어 실행을 위해 회로를 최적화한 후, 이제 대상 하드웨어에서 회로를 실행하고 기저 상태 에너지 추정을 위한 샘플을 수집할 준비가 되었습니다. 회로가 하나뿐이므로, 를 qiskit-ibm-runtime작업 실행 모드 사용하여 회로를 실행하겠습니다.

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 안사츠를 샘플링하면 노이즈 구성 풀( $\tilde{\mathcal{\chi}}$ )이 생성되며, 이 풀은 후처리 루틴에 사용됩니다. 여기에는 전자 번호가 잘못된 구성을 확률적으로 수정하는 구성 복구라는 방법(자세한 내용은 나중에 설명)이 포함됩니다. 올바른 전자 번호( $\tilde{\mathcal{\chi}}_{R}$ )를 가진 구성만 하위 샘플링하여 각 고유 구성의 출현 빈도에 따라 여러 배치로 배포합니다. 각 샘플 배치는 하위 공간( $\mathcal{S^{(k)}}$ )을 정의합니다. 다음으로 분자 해밀턴( $H$ )을 하위 공간에 투영합니다:

투영된 각 해밀토니안 $H_{\mathcal{S}^{(k)}}$ 은 고유값과 고유 벡터를 계산하여 고유 상태를 재구성하기 위해 대각선화되는 고유 솔버에 공급됩니다. 이 단원에서는 대각선화를 위해 PySCF 의 데이비슨 방법을 사용하는 qiskit-addon-sqd 패키지를 사용하여 해밀턴을 투영하고 대각선화합니다.

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

다음으로, 자체 일관된 구성 복구 루프에 대해 자세히 설명하고 위에서 언급한 단계를 구현하여 $N_2$ 해밀턴의 접지 상태 에너지를 추정하는 구체적인 코드 예제를 보여드립니다.

4.1 구성 복구: 개요

비트스트링의 각 비트(슬레이터 결정자)는 스핀 궤도를 나타냅니다. 비트스트링의 오른쪽 절반은 스핀업 궤도를, 왼쪽 절반은 스핀다운 궤도를 나타냅니다. 1 은 전자가 궤도를 점유하고 있음을 의미하고 0 은 궤도가 비어 있음을 의미합니다. 우리는 정확한 입자 수(업스핀 전자와 다운스핀 전자 모두)를 선험적으로 알고 있습니다. $N_x$ 전자가 있는 $x$ 행렬이 있다고 가정해 보겠습니다(즉, 비트 문자열에 $N_x$ 개의 $1$ s가 있습니다). 올바른 파티클 수는 $N$ 입니다. $N_x \neq N$ 이면 비트 문자열이 노이즈에 의해 손상된 것입니다. 자체 일관된 구성 루틴은 평균 궤도 점유 정보를 활용하여 $|N_x - N|$ 비트를 확률적으로 뒤집어 비트스트링을 수정하려고 시도합니다. 평균 궤도 점유도( $n$ )는 전자가 궤도를 얼마나 많이 차지하고 있는지를 알려줍니다. $N_x < N$ 의 경우 전자가 더 적으므로 일부 $0$ 를 $1$ 로, 또는 그 반대로 뒤집어야 합니다.

뒤집기 확률은 i-번째 스핀 궤도의 경우 $|x[i] - avg\_occupancy[i]|$ 입니다. [2] 에서 저자는 수정된 ReLU 함수를 사용하여 가중치를 부여한 뒤집기 확률을 사용했습니다.

여기서 $h$ 은 ReLU 함수의 "모서리" 위치를 정의하고, $\delta$ 매개변수는 모서리에 있는 ReLU 함수의 값을 정의합니다. $\delta = 0$ 의 경우 $w$ 은 true ReLU 함수가 되고, $\delta >0$ 의 경우 수정된 ReLU 함수가 됩니다. 이 논문에서 저자는 $\delta = 0.01$ 및 $h =$ 알파(또는 베타) 입자 수/알파(또는 베타) 스핀 궤도 수 $= N/M$ (충전 계수)를 사용했습니다.

평균 궤도 점유량( $n$ )은 선험적으로 알려져 있지 않습니다. 접지 상태 추정의 첫 번째 반복은 두 스핀 종 모두에서 정확한 입자 수만 있는 구성으로 시작됩니다. 첫 번째 반복이 끝나면 접지 상태의 추정치를 갖게 되고, 이 추정치를 사용하여 $n$ 의 첫 번째 추측을 만들 수 있습니다. 이 추측( $n$ )은 구성을 복구하고, 다음 접지 상태 추정 반복을 실행하고, $n$ 의 추측을 자체적으로 개선하는 데 사용됩니다. 이 과정은 중지 기준이 충족될 때까지 반복됩니다.

$N = 2$ 와 $x = |1000\rangle$ ( $N_x = 1$ )에 대한 다음 예시를 살펴봅시다. 입자 번호에 맞게 수정하려면 0s 중 하나를 1로 뒤집어야 하며, 선택 가능한 값은 1100, 1010, 1001 입니다. 뒤집기 확률에 따라 선택 항목 중 하나가 복구된 구성 (또는 정확한 파티클 수를 가진 비트스트링)으로 선택됩니다.

첫 번째 반복에서 두 개의 배치를 실행하고 그 결과 예상되는 지상 상태는 다음과 같다고 가정합니다:

계산 기준 상태와 그 진폭을 사용하여 스핀 궤도(큐비트) 당 전자의 점유 확률(간단히 점유 )을 계산할 수 있습니다(확률 = |진폭| $^2$ ). 아래에서는 추정 기저 상태에 나타나는 각 비트 문자열에 대한 큐비트별 점유율을 표로 작성하고 배치에 대한 총 궤도 점유율을 계산합니다. 키스킷 순서 규칙에 따라, 가장 오른쪽 비트는 qubit-0 ( Q0 ), 가장 왼쪽 비트는 Q3 를 나타냅니다.

점유 ( Batch0 ):

점유 ( Batch1 )

점유율(배치 전체 평균)

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

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

마지막으로 위의 가중치를 사용하여 비어 있는 Q2, Q1, Q0 궤도 중 하나를 뒤집을 수 있습니다. 위의 값에 따르면 Q0 이 뒤집힐 가능성이 가장 높으며, 복구 가능한 구성은 $\vert \text{1001} \rangle$ 입니다.

전체 자체 일관된 구성 복구 프로세스는 다음과 같이 요약할 수 있습니다:

첫 번째 반복: 양자 컴퓨터에서 생성된 비트스트링(구성 또는 슬레이터 결정자)이 각 스핀 섹터의 입자 수가 올바른( $\widetilde{\chi}_{correct}$ ) 구성과 잘못된( $\widetilde{\chi}_{incorrect}$ ) 구성을 모두 포함하는 $\widetilde{\chi}$ 집합을 형성한다고 가정합니다.

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 기저 상태 추정

먼저 후처리를 위해 개수를 비트 문자열 행렬과 확률 배열로 변환합니다.

행렬의 각 행은 하나의 고유한 비트스트링을 나타냅니다. 큐비트는 키스킷의 비트 문자열 오른쪽부터 인덱싱되므로 0 열은 큐비트 N-1 를 나타내고, N-1 열은 큐비트 0 를 나타내며, 여기서 N 는 큐비트 수입니다.

알파 궤도는 열 인덱스 범위 (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: 각 아이겐솔버가 실행하는 최대 데이비슨 사이클 수

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 postselect 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_cycle=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)

Output:

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 (화학적 정확도는 일반적으로 1kcal/mol로 인정됨 $\approx$ 1.6 mH ) 이내의 기저 상태 에너지를 추정하는 것을 보여줍니다. 두 번째 플롯은 최종 반복 후 각 공간 궤도의 평균 점유도를 보여줍니다. 솔루션에서 스핀업과 스핀다운 전자가 모두 처음 5개의 궤도를 높은 확률로 점유하고 있음을 알 수 있습니다.

추정된 지상 상태 에너지는 괜찮지만 화학적 정확도 한계( $\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()

Output:

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

독자를 위한 연습

samples_per_batch 파라미터를 점진적으로 증가시키고(예: $1000$ 에서 $10000$ 단계로 $1000$; 내 컴퓨터의 메모리가 허용됨) 예상 접지 상태 에너지를 비교합니다.

참조

[1] M. 모타 외.., "상관된 전자 상태를 위한 물리적 직관과 하드웨어 효율성의 연결: 전자 구조를 위한 로컬 단일 클러스터 자스트로우 안사츠"(2023). Chem. Sci.., 2023, 14, 11213.

[2] J. 로블레도-모레노 외.., "양자 중심 슈퍼컴퓨터의 정확한 솔루션을 넘어서는 화학"(2024). arXiv:quant-ph/2405.05068.