{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "9e40af77-7f0f-4dd6-ab0a-420cf396050e",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"화학 Hamiltonian의 샘플 기반 양자 대각화\"\n",
        "description: \"샘플 기반 양자 대각화 알고리즘을 사용하여 잡음이 있는 양자 하드웨어로 질소 분자를 시뮬레이션합니다.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore fontdict fontsize milli LUCJ CCSD ccsd hcore pvdz */}\n",
        "\n",
        "<span id=\"sample-based-quantum-diagonalization-of-a-chemistry-hamiltonian\" />\n",
        "\n",
        "# 화학 Hamiltonian의 샘플 기반 양자 대각화\n",
        "\n",
        "*사용 예상 시간: Heron r2 프로세서에서 1분 미만(참고: 이는 예상치일 뿐입니다. 런타임은 다를 수 있습니다.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9701917b",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## 학습 성과\n",
        "\n",
        "* [SQD Qiskit 애드온](/docs/addons/qiskit-addon-sqd) 을 사용하여 양자 처리 장치(QPU)에서 샘플링한 비트스트링을 바탕으로 분자 계의 기저 상태 에너지를 근사하는 방법.\n",
        "* [ffsim을](https://github.com/qiskit-community/ffsim) 사용하여 양자 화학 시뮬레이션을 위한 국소 단일체 클러스터 Jastrow(LUCJ) 회로를 구축하는 방법.\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## 전제조건\n",
        "\n",
        "* 양자화학 및 제2 양자화\n",
        "* 샘플러 프리미티브를 사용하여 양자 회로에서 샘플링하기\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## 배경\n",
        "\n",
        "이 튜토리얼에서는 [SQD Qiskit 애드온](https://github.com/Qiskit/qiskit-addon-sqd) 을 사용하여 [샘플 기반 양자 대각화(SQD) 알고리즘을](https://arxiv.org/abs/2405.05068) 구현함으로써, 잡음이 포함된 양자 샘플을 후처리하여 평형 결합 길이에서 질소 분자 $\\text{N}_2$ 의 기저 상태를 근사하는 방법을 보여줍니다. 이 소프트웨어에 대한 자세한 내용은 [관련](/docs/addons/qiskit-addon-sqd) 설명서에서 확인할 수 있으며, 여기에는 시작에 도움이 되는 [간단한](/docs/addons/qiskit-addon-sqd/guides/quickstart) 예제도 포함되어 있습니다.\n",
        "\n",
        "이 튜토리얼은 양자화학, 특히 분자의 기저 상태 에너지를 구하는 방법에 익숙한 사용자에게 권장됩니다. 워크플로우에 대한 자세한 안내는 [양자 대각화 알고리즘 강좌](/learning/courses/quantum-diagonalization-algorithms) 를 참조하십시오.\n",
        "\n",
        "SQD는 양자 시스템의 해밀토니안과 같은 양자 연산자의 고유값과 고유벡터를 구하기 위해 양자 컴퓨팅과 분산형 고전 컴퓨팅을 함께 활용하는 기법입니다. 전통적인 분산 컴퓨팅은 양자 프로세서에서 얻은 샘플을 처리하고, 이들이 생성하는 부분공간 내에서 목표 해밀토니안을 투영 및 대각화하는 데 사용된다. SQD 기반 워크플로는 다음과 같은 단계로 구성됩니다:\n",
        "\n",
        "1. 회로 안사츠를 선택하고 양자 컴퓨터에서 이를 기준 상태(이 경우 [하트리-복](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) 상태)에 적용합니다.\n",
        "2. 결과 양자 상태에서 샘플 비트스트링을 추출합니다.\n",
        "3. 비트열에 대해 *자기일관성 구성 복원* 절차를 실행하여 기저 상태 근사치를 구한다.\n",
        "\n",
        "SQD는 목표 고유 상태가 희박할 때 잘 작동하는 것으로 알려져 있습니다. 파동 함수는 문제의 크기에 따라 기하급수적으로 증가하지 않는 기저 상태( $\\mathcal{S} = \\{|x\\rangle \\}$ )의 집합에서 지원됩니다.\n",
        "\n",
        "<span id=\"quantum-chemistry\" />\n",
        "\n",
        "### 양자화학\n",
        "\n",
        "분자 시스템의 해밀턴은 다음과 같이 쓸 수 있습니다\n",
        "\n",
        "$$\n",
        "\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n",
        "+ \\frac12\n",
        "\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n",
        "h_{prqs} \\,\n",
        "\\hat{a}^\\dagger_{p\\sigma}\n",
        "\\hat{a}^\\dagger_{q\\tau}\n",
        "\\hat{a}_{s\\tau}\n",
        "\\hat{a}_{r\\sigma},\n",
        "$$\n",
        "\n",
        "여기서 $h_{pr}$ 및 $h_{prqs}$ 은 컴퓨터 프로그램을 사용하여 분자의 사양으로부터 계산할 수 있는 분자 적분이라는 복소수입니다. 이 튜토리얼에서는 적분을 계산하기 위해 [PySCF](https://pyscf.org/) 소프트웨어 패키지를 사용하여 적분을 계산합니다.\n",
        "\n",
        "분자 해밀턴이 어떻게 유도되는지에 대한 자세한 내용은 양자 화학에 관한 교과서(예: Szabo와 Ostlund의 *현대 양자 화학* )를 참조하세요. 양자 화학 문제가 양자 컴퓨터에 어떻게 매핑되는지에 대한 자세한 설명은 2024년 키스킷 글로벌 서머 스쿨의 강연인 [*Mapping Problems to Qubits에서*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900) 확인해보세요.\n",
        "\n",
        "<span id=\"local-unitary-cluster-jastrow-lucj-ansatz\" />\n",
        "\n",
        "### 지역 단위 군집 자스트로(LUCJ) 가정\n",
        "\n",
        "SQD는 샘플을 추출하기 위해 양자 회로 안잭츠가 필요합니다. 이 튜토리얼에서는 물리적 근거와 하드웨어 친화성을 모두 갖춘 [Jastrow의 국소 단일 클러스터(LUCJ)](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) 가정을 사용할 것입니다. [ffsim을](https://qiskit-community.github.io/ffsim/) 사용하여 근사 회로를 구성할 것입니다.\n",
        "\n",
        "LUCJ 접근법은 큐비트 연결성이 제한된 QPU에 적용될 수 있다. 스핀 궤도는 큐비트에 매핑되어, 해당 안자츠를 적용할 때 SWAP 게이트를 이용한 라우팅이 필요하지 않게 됩니다. IBM® 하드웨어는 정사각 격자 또는 무거운 육각 격자 큐비트 토폴로지를 갖습니다. 헤비-헥스 격자의 경우, 아래에 묘사된 “지그재그” 패턴을 채택할 수 있습니다. 이 패턴에서, 동일한 스핀을 가진 궤도는 선형 토폴로지를 가진 큐비트(빨간색 및 파란색 원)에 매핑되며, 서로 다른 스핀을 가진 궤도 간에는 공간 궤도마다 4번째마다 연결이 존재하며, 이 연결은 보조 큐비트(보라색 원)에 의해 이루어집니다.\n",
        "\n",
        "![무거운 16진수 격자의 LUCJ 안사츠에 대한 큐비트 매핑 다이어그램](https://quantum.cloud.ibm.com/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/7e0ee7e1-2d24-417f-ac59-25c58db79aa9.avif)\n",
        "\n",
        "<span id=\"self-consistent-configuration-recovery\" />\n",
        "\n",
        "### 자기 일관성 있는 구성 복구\n",
        "\n",
        "자체 일관된 구성 복구 절차는 노이즈가 많은 양자 샘플에서 최대한 많은 신호를 추출하도록 설계되었습니다. 분자 해밀턴은 입자 수와 스핀 Z를 보존하기 때문에 이러한 대칭성을 보존하는 회로 안사츠를 선택하는 것이 합리적입니다. 하트트리-복 상태에 적용하면 결과 상태는 노이즈가 없는 설정에서 고정 파티클 수와 스핀 Z를 갖습니다. 따라서 이 상태에서 샘플링된 비트스트링의 스핀- $\\alpha$ 및 스핀- $\\beta$ 반쪽은 하트트리-복 상태와 동일한 [해머링 무게를](https://en.wikipedia.org/wiki/Hamming_weight) 가져야 합니다. 현재 양자 프로세서에는 노이즈가 존재하기 때문에 측정된 일부 비트 문자열은 이 속성을 위반하게 됩니다. 단순한 형태의 포스트셀렉션은 이러한 비트 문자열을 버리지만, 비트 문자열에 여전히 일부 신호가 포함되어 있을 수 있으므로 낭비입니다. 자체 일관된 복구 절차는 사후 처리에서 해당 신호의 일부를 복구하려고 시도합니다. 이 절차는 반복적이며, 먼저 원시 샘플에서 계산된 지상 상태의 각 궤도 평균 점유량 추정치를 입력으로 필요로 합니다. 절차는 반복으로 실행되며 각 반복에는 다음 단계가 있습니다:\n",
        "\n",
        "1. 지정된 대칭을 위반하는 각 비트 문자열에 대해 비트 문자열을 평균 궤도 점유에 대한 현재 추정치에 가깝게 만들도록 설계된 확률적 절차를 사용하여 비트를 뒤집어 새 비트 문자열을 얻습니다.\n",
        "2. 대칭을 만족하는 모든 이전 및 새 비트스트링과 미리 선택한 고정 크기의 하위 샘플 하위 집합을 수집합니다.\n",
        "3. 비트 문자열의 각 하위 집합에 대해 해밀턴을 해당 기저 벡터에 의해 펼쳐진 부분 공간에 투영하고(이러한 기저 벡터에 대한 설명은 [이전 섹션](#quantum-chemistry) 참조), 고전적인 컴퓨터에서 투영된 해밀턴의 기저 상태 추정치를 계산합니다.\n",
        "4. 평균 궤도 점유량 추정치를 가장 낮은 에너지의 지상 상태 추정치로 업데이트합니다.\n",
        "\n",
        "<span id=\"sqd-workflow-diagram\" />\n",
        "\n",
        "### SQD 워크플로 다이어그램\n",
        "\n",
        "SQD 워크플로우는 다음 다이어그램에 나와 있습니다:\n",
        "\n",
        "![SQD 알고리즘의 워크플로 다이어그램](https://quantum.cloud.ibm.com/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/fd7e816f-4e2e-4dd7-a7da-f71afb9ca68d.avif)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "88422c4b",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## 요구사항\n",
        "\n",
        "이 튜토리얼을 시작하기 전에 다음이 설치되어 있는지 확인하세요:\n",
        "\n",
        "* Qiskit SDK v1.0 또는 이후 버전, [시각화](/docs/api/qiskit/visualization) 지원 기능 포함\n",
        "* Qiskit Runtime v0.22 또는 이후 (`pip install qiskit-ibm-runtime`)\n",
        "* SQD Qiskit 애드온 v0.12 이상 (`pip install qiskit-addon-sqd`)\n",
        "* ffsim v0.0.75 이상 (`pip install ffsim`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c6e44a31",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## 설정\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "6e51c3d8",
      "metadata": {},
      "outputs": [],
      "source": [
        "import math\n",
        "\n",
        "import ffsim\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pyscf\n",
        "import pyscf.cc\n",
        "import pyscf.mcscf\n",
        "from qiskit import QuantumCircuit, QuantumRegister\n",
        "from qiskit.primitives import StatevectorSampler\n",
        "from qiskit.providers.fake_provider import GenericBackendV2\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "from qiskit_ibm_runtime import SamplerV2 as Sampler"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4bc6ee26-4371-4cd2-80a7-60752bf8775d",
      "metadata": {},
      "source": [
        "<span id=\"small-scale-simulator-example\" />\n",
        "\n",
        "## 소규모 시뮬레이터 예시\n",
        "\n",
        "이 튜토리얼에서는 질소 분자가 평형 결합 거리 근처에 있을 때의 기저 상태에 대한 근사값을 구해 보겠습니다. 먼저 실험을 시뮬레이션하고 제대로 작동하는지 확인하기 위해 작은 STO-6G 기반 집합을 사용합니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "afeb054c",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### 1단계: 고전적 입력을 양자 문제에 매핑하기\n",
        "\n",
        "먼저, 분자와 그 특성을 명시합니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "b821e660",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "converged SCF energy = -108.464957764796\n",
            "CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000\n",
            "norb = 8\n",
            "nelec = (5, 5)\n"
          ]
        }
      ],
      "source": [
        "# Specify molecule properties\n",
        "spin_sq = 0\n",
        "\n",
        "# Build N2 molecule\n",
        "mol = pyscf.gto.Mole()\n",
        "mol.build(\n",
        "    atom=[[\"N\", (0, 0, -0.5)], [\"N\", (0, 0, 0.5)]],\n",
        "    basis=\"sto-6g\",\n",
        "    symmetry=\"Dooh\",\n",
        ")\n",
        "\n",
        "# Define active space\n",
        "n_frozen = 2\n",
        "active_space = range(n_frozen, mol.nao_nr())\n",
        "\n",
        "# Get molecular integrals\n",
        "scf = pyscf.scf.RHF(mol).run()\n",
        "norb = len(active_space)\n",
        "n_electrons = int(sum(scf.mo_occ[active_space]))\n",
        "n_alpha = (n_electrons + mol.spin) // 2\n",
        "n_beta = (n_electrons - mol.spin) // 2\n",
        "nelec = (n_alpha, n_beta)\n",
        "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n",
        "mo = cas.sort_mo(active_space, base=0)\n",
        "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
        "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n",
        "\n",
        "# Compute exact energy using FCI\n",
        "reference_energy = cas.run().e_tot\n",
        "\n",
        "print(f\"norb = {norb}\")\n",
        "print(f\"nelec = {nelec}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "96bfe018",
      "metadata": {},
      "source": [
        "LUCJ 안사츠 회로를 구성하기 전에 먼저 다음 코드 셀에서 CCSD 계산을 수행합니다. 이 계산의 [$t_1$ 및 $t_2$ 진폭은](https://en.wikipedia.org/wiki/Coupled_cluster#Cluster_operator) ansatz의 매개변수를 초기화하는 데 사용됩니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "efe83d98",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "E(CCSD) = -108.5933309085008  E_corr = -0.1283731437052352\n"
          ]
        }
      ],
      "source": [
        "# Get CCSD t2 amplitudes for initializing the ansatz\n",
        "ccsd = pyscf.cc.CCSD(\n",
        "    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
        ").run()\n",
        "t1 = ccsd.t1\n",
        "t2 = ccsd.t2"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "f4d882fa",
      "metadata": {},
      "source": [
        "이제 [ffsim을](https://github.com/qiskit-community/ffsim) 사용하여 초기 가정 회로를 생성합니다. 우리 분자는 폐쇄 껍질 하트리-팍(Hartree-Fock) 상태를 가지므로, 스핀 균형형 UCJ 가설( [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced) )을 사용합니다. 우리는 $t_2$ 진폭의 \"압축된\" 이중 분해를 활성화하기 위해 메서드 `from_t_amplitudes` 내에서 를 설정했습니다 `optimize=True` (자세한 내용은 ffsim 문서의 ['국소 유니터리 클러스터 Jastrow(LUCJ) 가정' 항목을](https://qiskit-community.github.io/ffsim/explanations/lucj.html#Parameter-initialization-from-CCSD) 참조하십시오).\n",
        "\n",
        "LUCJ 안자츠는 QPU의 사용 가능한 연결성에 맞춰 조정되므로, 안자츠를 생성하기 전에 QPU 백엔드를 초기화해야 합니다. 우선, 정사각 격자 결합 맵과 LUCJ 안자츠가 자연스럽게 분해되는 게이트 세트를 갖춘 일반적인 백엔드를 만들어 보겠습니다. 그런 다음, 를 사용하여 “래더(ladder)” 레이아웃(ffsim 문서의 [설명](https://qiskit-community.github.io/ffsim/explanations/lucj.html#The-local-UCJ-\\(LUCJ\\)-ansatz) 참조)에 따라 LUCJ 안자츠를 주어진 백엔드로 트랜스파일링하는 데 특화된 패스 매니저를 생성할 `ffsim.qiskit.generate_lucj_pass_manager` 것입니다. 이 함수는 점수 기반 휴리스틱을 사용하여 선택된 레이아웃과 관련된 오류를 최소화하며, 이는 백엔드가 실제 QPU이거나 노이즈 모델이 적용된 시뮬레이터인 경우 중요합니다. 이 함수는 패스 매니저를 반환하는 것 외에도, 하드웨어에서 구현 가능한 알파-베타 결합 쌍을 반환합니다. 모든 쌍을 구현할 수 없는 경우, 경고가 표시됩니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "dd69a86c",
      "metadata": {},
      "outputs": [],
      "source": [
        "import warnings\n",
        "\n",
        "from qiskit.transpiler import CouplingMap\n",
        "\n",
        "warnings.formatwarning = lambda msg, *args, **kwargs: f\"Warning: {msg}\\n\"\n",
        "\n",
        "# Set ansatz properties\n",
        "n_reps = 4\n",
        "pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n",
        "\n",
        "# Let generate_lucj_pass_manager determine the alpha-beta interactions\n",
        "pairs_ab = None\n",
        "\n",
        "# Initialize backend\n",
        "coupling_map = CouplingMap.from_grid(norb, 2)\n",
        "backend = GenericBackendV2(\n",
        "    coupling_map.size(),\n",
        "    coupling_map=coupling_map,\n",
        "    basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\", \"swap\"],\n",
        ")\n",
        "\n",
        "# Create pass manager\n",
        "pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n",
        "    backend=backend,\n",
        "    norb=norb,\n",
        "    connectivity=\"square\",\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    optimization_level=3,\n",
        ")\n",
        "\n",
        "# Create the LUCJ ansatz operator\n",
        "ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
        "    t2=t2,\n",
        "    t1=t1,\n",
        "    n_reps=n_reps,\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    # Setting optimize=True enables the \"compressed\" factorization\n",
        "    optimize=True,\n",
        ")\n",
        "\n",
        "# create an empty quantum circuit\n",
        "qubits = QuantumRegister(2 * norb, name=\"q\")\n",
        "circuit = QuantumCircuit(qubits)\n",
        "\n",
        "# prepare Hartree-Fock state as the reference state and append it\n",
        "# to the quantum circuit\n",
        "circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n",
        "\n",
        "# apply the UCJ operator to the reference state\n",
        "circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
        "circuit.measure_all()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "db11bf6d",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-for-quantum-hardware-execution\" />\n",
        "\n",
        "### 2단계: 양자 하드웨어 실행에 최적화하기\n",
        "\n",
        "다음으로, 대상 하드웨어에 맞게 회로를 최적화합니다. 일반적으로 이 단계에서는 하드웨어 백엔드와 해당 백엔드의 패스 매니저를 초기화합니다. 그러나 LUCJ 접근법은 하드웨어 연결 환경에 맞춰져 있기 때문에, 우리는 이미 이전 단계에서 이러한 작업을 수행했습니다. 이제 남은 일은 회로에서 패스 매니저를 실행하여, QPU에서 직접 실행할 수 있는 ISA 회로로 변환하는 것뿐입니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "7d554aa5",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gate counts: OrderedDict({'xx_plus_yy': 254, 'cp': 88, 'p': 64, 'measure': 16, 'x': 10, 'barrier': 1})\n"
          ]
        }
      ],
      "source": [
        "isa_circuit = pass_manager.run(circuit)\n",
        "print(f\"Gate counts: {isa_circuit.count_ops()}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0cc1edef",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### 3단계: `Qiskit primitives` 명령어로 실행합니다\n",
        "\n",
        "하드웨어 실행을 위해 회로를 최적화한 후, 이제 대상 하드웨어에서 회로를 실행하고 기저 상태 에너지 추정을 위한 샘플을 수집할 준비가 되었습니다. 회로가 하나뿐이므로, ‘ IBM Quantum ’ 컴퓨트 서비스의 [작업 실행 모드를](/docs/guides/execution-modes) 사용하여 회로를 실행하겠습니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "93c1cef3-298e-4deb-8512-769fe94cd5a5",
      "metadata": {},
      "outputs": [],
      "source": [
        "rng = np.random.default_rng(12345)\n",
        "sampler = StatevectorSampler(seed=rng)\n",
        "job = sampler.run([isa_circuit], shots=100_000)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "332ecab3-77e6-473f-b0e7-af30f983393a",
      "metadata": {},
      "outputs": [],
      "source": [
        "primitive_result = job.result()\n",
        "pub_result = primitive_result[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6df05b6e",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### 4단계: 후처리 수행 및 원하는 클래식 형식으로 결과 반환\n",
        "\n",
        "QPU 출력의 품질을 판단하는 유용한 지표는 반환된 유효한 구성의 수입니다. 유효한 구성은 올바른 입자 수와 스핀 Z를 가지며, 이는 비트열의 오른쪽 절반이 스핀 업 전자의 수와 동일한 해밍 가중치를 가지며, 왼쪽 절반이 스핀 다운 전자의 수와 동일한 해밍 가중치를 가짐을 의미한다. 다음 셀은 유효한 샘플링된 구성의 비율을 계산합니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "718f8517",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Fraction of sampled configurations that are valid: 1.0\n"
          ]
        }
      ],
      "source": [
        "def is_valid_bitstring(\n",
        "    bitstring: str, norb: int, nelec: tuple[int, int]\n",
        ") -> bool:\n",
        "    n_alpha, n_beta = nelec\n",
        "    return (\n",
        "        len(bitstring) == 2 * norb\n",
        "        and bitstring[norb:].count(\"1\") == n_alpha\n",
        "        and bitstring[:norb].count(\"1\") == n_beta\n",
        "    )\n",
        "\n",
        "\n",
        "bit_array = pub_result.data.meas\n",
        "num_valid = sum(\n",
        "    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
        ")\n",
        "valid_fraction = num_valid / bit_array.num_shots\n",
        "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7b126f3a",
      "metadata": {},
      "source": [
        "우리는 노이즈가 없는 시뮬레이터에서 회로를 샘플링하고 있으므로 모든 비트열은 유효합니다. 노이즈가 많은 QPU에서 실행할 경우, 이 비율은 1보다 작겠지만, 다음 셀에서 계산된 것처럼 비트열을 균일한 확률로 무작위 표본 추출했을 때 예상되는 비율보다는 클 것으로 기대됩니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "6b3e4bca",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625\n"
          ]
        }
      ],
      "source": [
        "expected_fraction_random = (\n",
        "    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
        ")\n",
        "print(\n",
        "    f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n",
        "    f\"{expected_fraction_random}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "eb704101-0fe8-4d12-b572-b1d844e35a90",
      "metadata": {},
      "source": [
        "이제 `diagonalize_fermionic_hamiltonian` 함수를 사용하여 해밀턴의 기저 상태 에너지를 추정합니다. 이 기능은 자체 일관된 구성 복구 절차를 수행하여 노이즈가 있는 양자 샘플을 반복적으로 개선하여 에너지 추정치를 개선합니다. 나중에 분석할 수 있도록 중간 결과를 저장할 수 있도록 콜백 함수를 전달합니다. `diagonalize_fermionic_hamiltonian` 인자에 대한 설명은 [API 설명서를](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) 참조하세요.\n",
        "\n",
        "여기서 우리는 인자를 `initial_occupancies` 사용하여 `diagonalize_fermionic_hamiltonian` 기저 상태의 궤도 점유도에 대한 초기 추측값으로 하트리-포크(Hartree-Fock) 구성을 지정합니다. 이 접근법은 기저 상태가 하트리-포크 구성에 상당한 지지를 갖는 시스템에 합리적이지만, 다른 상황에서는 적합하지 않을 수 있다. 다만 더 발전된 계산 방법들이 그러한 경우에 더 나은 초기 추측값을 제공할 수 있다. 지정하면 유효한 구성이 샘플링되지 않은 `initial_occupancies` 경우에도 구성 복구가 실행될 수 있습니다. 이는 잡음이 많은 QPU에서 대규모 회로를 샘플링할 때 발생할 수 있는 상황입니다. 이 인수가 없으면 유효한 구성이 제공되지 않을 경우 구성 복구가 실패하고 오류가 발생합니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "2f32a352",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "Final energy: -108.59589540253083\n",
            "Final energy error: 9.19484552071026e-05\n"
          ]
        }
      ],
      "source": [
        "from functools import partial\n",
        "\n",
        "from qiskit_addon_sqd.fermion import (\n",
        "    SCIResult,\n",
        "    diagonalize_fermionic_hamiltonian,\n",
        "    solve_sci_batch,\n",
        ")\n",
        "\n",
        "# SQD options\n",
        "energy_tol = 1e-3\n",
        "occupancies_tol = 1e-3\n",
        "max_iterations = 5\n",
        "\n",
        "# Eigenstate solver options\n",
        "num_batches = 3\n",
        "samples_per_batch = 1000\n",
        "max_dim = 50\n",
        "symmetrize_spin = True\n",
        "carryover_threshold = 1e-4\n",
        "max_cycle = 200\n",
        "\n",
        "# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies\n",
        "initial_occupancies = (\n",
        "    np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
        "    np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
        ")\n",
        "\n",
        "# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
        "# you can omit this step, in which case you would not specify the sci_solver argument\n",
        "# in the call to diagonalize_fermionic_hamiltonian below.\n",
        "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n",
        "\n",
        "# List to capture intermediate results\n",
        "result_history = []\n",
        "\n",
        "\n",
        "def callback(results: list[SCIResult]):\n",
        "    result_history.append(results)\n",
        "    iteration = len(result_history)\n",
        "    print(f\"Iteration {iteration}\")\n",
        "    for i, result in enumerate(results):\n",
        "        print(f\"\\tSubsample {i}\")\n",
        "        print(f\"\\t\\tEnergy: {result.energy + nuclear_repulsion_energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    hcore,\n",
        "    eri,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=nelec,\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=energy_tol,\n",
        "    occupancies_tol=occupancies_tol,\n",
        "    max_iterations=max_iterations,\n",
        "    sci_solver=sci_solver,\n",
        "    symmetrize_spin=symmetrize_spin,\n",
        "    max_dim=max_dim,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    carryover_threshold=carryover_threshold,\n",
        "    callback=callback,\n",
        "    seed=rng,\n",
        ")\n",
        "\n",
        "final_energy = result.energy + nuclear_repulsion_energy\n",
        "energy_error = final_energy - reference_energy\n",
        "print(f\"Final energy: {final_energy}\")\n",
        "print(f\"Final energy error: {energy_error}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9d78906b-4759-4506-9c69-85d4e67766b3",
      "metadata": {},
      "source": [
        "<span id=\"visualize-the-results\" />\n",
        "\n",
        "#### 결과\n",
        "\n",
        "시각화\n",
        "\n",
        "첫 번째 그래프는 이번 시뮬레이션에서 첫 번째 반복 후 이미 정확한 해에 근접해 `1 mH` 있음을 보여줍니다(화학 정확도는 일반적으로 `1 kcal/mol`$\\approx$ 로 `1.6 mH`간주됩니다). 하지만 이는 작은 시스템이며, 샘플에 잡음이 없기 때문에 구성 복원이 필요하지 않습니다. 소음이 심한 QPU에서 실행되는 대규모 시스템의 경우, 여러 번의 구성 복구 반복이 필요할 수 있으며, 최종 정확도가 떨어질 수 있습니다. 일반적으로, 구성 복원 반복 횟수를 늘리거나 배치당 샘플 수를 늘림으로써 에너지 값을 개선할 수 있습니다.\n",
        "\n",
        "두 번째 플롯은 최종 반복 후 각 공간 궤도의 평균 점유도를 보여줍니다. 솔루션에서 스핀업과 스핀다운 전자가 모두 처음 5개의 궤도를 높은 확률로 점유하고 있음을 알 수 있습니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "caffd888-e89c-4aa9-8bae-4d1bb723b35e",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/caffd888-e89c-4aa9-8bae-4d1bb723b35e-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Data for energies plot\n",
        "x1 = range(len(result_history))\n",
        "min_e = [\n",
        "    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
        "    for result in result_history\n",
        "]\n",
        "e_diff = [abs(e - reference_energy) for e in min_e]\n",
        "\n",
        "# Chemical accuracy (+/- 1.6 milli-Hartree)\n",
        "chem_accuracy = 0.0016\n",
        "\n",
        "# Data for avg spatial orbital occupancy\n",
        "y2 = np.sum(result.orbital_occupancies, axis=0)\n",
        "x2 = range(len(y2))\n",
        "\n",
        "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
        "\n",
        "# Plot energies\n",
        "axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n",
        "axs[0].set_xticks(x1)\n",
        "axs[0].set_xticklabels(x1)\n",
        "axs[0].set_yscale(\"log\")\n",
        "axs[0].set_ylim(1e-5, 1e-2)\n",
        "axs[0].axhline(\n",
        "    y=chem_accuracy,\n",
        "    color=\"#BF5700\",\n",
        "    linestyle=\"--\",\n",
        "    label=\"chemical accuracy\",\n",
        ")\n",
        "axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
        "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
        "axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n",
        "axs[0].legend()\n",
        "\n",
        "# Plot orbital occupancy\n",
        "axs[1].bar(x2, y2, width=0.8)\n",
        "axs[1].set_xticks(x2)\n",
        "axs[1].set_xticklabels(x2)\n",
        "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
        "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
        "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ce0eecb3-8a23-4118-aa1e-a28afcec6334",
      "metadata": {},
      "source": [
        "<span id=\"large-scale-hardware-example\" />\n",
        "\n",
        "## 대규모 하드웨어 예시\n",
        "\n",
        "이제 실제 양자 하드웨어에서 더 큰 규모의 예제를 실행해 보겠습니다. 여기서는 cc-pVDZ 기저 세트를 바탕으로 질소 분자의 활성 공간을 도출해 보겠습니다. 또한 다음과 같은 변경 사항도 적용할 예정입니다:\n",
        "\n",
        "* 게이트 오류로 인한 노이즈를 줄이려면 `1` `n_reps` 값을 낮추십시오\n",
        "* 백엔드 프로세서 제품군에 따라 안자츠 연결 방식을 설정합니다(Heron의 경우 heavy-hex, Nighthawk의 경우 square)\n",
        "* 더 큰 부분공간 차원을 허용할 수 있도록 `1000` 를 `max_dim` 로 늘리십시오.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062",
      "metadata": {},
      "source": [
        "<span id=\"steps-1-4\" />\n",
        "\n",
        "### 1\\~4단계\n",
        "\n",
        "여기서는 모든 단계를 통합하여 대규모의 단일 워크플로우를 구성한 뒤, 이를 실제 양자 하드웨어에서 실행합니다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "3858949c-a55d-4ff8-a0fc-54fb53e131b5",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "converged SCF energy = -108.929838385609\n",
            "norb = 26\n",
            "nelec = (5, 5)\n",
            "E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450041\n",
            "Using backend ibm_boston\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].\n",
            "Removing interaction (24, 24) from the end.\n",
            "Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].\n",
            "Removing interaction (20, 20) from the end.\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gate counts: OrderedDict({'sx': 7036, 'rz': 6965, 'cz': 1858, 'measure': 52, 'x': 50, 'barrier': 1})\n",
            "Fraction of sampled configurations that are valid: 0.01677\n",
            "Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.17922180754668\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.17713541195594\n",
            "\t\tSubspace dimension: 962361\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.17407767345108\n",
            "\t\tSubspace dimension: 968256\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.1873108379026\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.1876970140745\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.18255107152576\n",
            "\t\tSubspace dimension: 1000000\n",
            "Iteration 3\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.19052148187927\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.19316375313544\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.19087348337145\n",
            "\t\tSubspace dimension: 1000000\n",
            "Iteration 4\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.19374847699348\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.19517039328333\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.19422416032815\n",
            "\t\tSubspace dimension: 1000000\n",
            "Iteration 5\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.19580768221108\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.19608553581007\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.19548131833729\n",
            "\t\tSubspace dimension: 1000000\n",
            "Final energy: -109.19608553581007\n",
            "Final energy error: 0.03194368084709254\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/3858949c-a55d-4ff8-a0fc-54fb53e131b5-3.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# ------------------------------ Step 1 ------------------------------\n",
        "# Build N2 molecule\n",
        "mol = pyscf.gto.Mole()\n",
        "mol.build(\n",
        "    atom=[[\"N\", (0, 0, -0.5)], [\"N\", (0, 0, 0.5)]],\n",
        "    basis=\"cc-pvdz\",\n",
        "    symmetry=\"Dooh\",\n",
        ")\n",
        "\n",
        "# Define active space\n",
        "n_frozen = 2\n",
        "active_space = range(n_frozen, mol.nao_nr())\n",
        "\n",
        "# Get molecular integrals\n",
        "scf = pyscf.scf.RHF(mol).run()\n",
        "norb = len(active_space)\n",
        "n_electrons = int(sum(scf.mo_occ[active_space]))\n",
        "n_alpha = (n_electrons + mol.spin) // 2\n",
        "n_beta = (n_electrons - mol.spin) // 2\n",
        "nelec = (n_alpha, n_beta)\n",
        "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n",
        "mo = cas.sort_mo(active_space, base=0)\n",
        "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
        "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n",
        "\n",
        "# Store reference energy from SCI calculation performed separately\n",
        "reference_energy = -109.22802921665716\n",
        "\n",
        "print(f\"norb = {norb}\")\n",
        "print(f\"nelec = {nelec}\")\n",
        "\n",
        "# Get CCSD t2 amplitudes for initializing the ansatz\n",
        "ccsd = pyscf.cc.CCSD(\n",
        "    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
        ").run()\n",
        "t1 = ccsd.t1\n",
        "t2 = ccsd.t2\n",
        "\n",
        "# Set ansatz properties\n",
        "n_reps = 1\n",
        "pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n",
        "\n",
        "# Let generate_lucj_pass_manager determine the alpha-beta interactions\n",
        "pairs_ab = None\n",
        "\n",
        "# Initialize backend\n",
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(operational=True, simulator=False)\n",
        "print(f\"Using backend {backend.name}\")\n",
        "\n",
        "# Create pass manager\n",
        "connectivity_map = {\"Heron\": \"heavy-hex\", \"Nighthawk\": \"square\"}\n",
        "pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n",
        "    backend=backend,\n",
        "    norb=norb,\n",
        "    connectivity=connectivity_map[\n",
        "        backend.configuration().processor_type[\"family\"]\n",
        "    ],\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    optimization_level=3,\n",
        ")\n",
        "\n",
        "# Create the LUCJ ansatz operator\n",
        "ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
        "    t2=t2,\n",
        "    t1=t1,\n",
        "    n_reps=n_reps,\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    # Setting optimize=True enables the \"compressed\" factorization\n",
        "    optimize=True,\n",
        ")\n",
        "\n",
        "# create an empty quantum circuit\n",
        "qubits = QuantumRegister(2 * norb, name=\"q\")\n",
        "circuit = QuantumCircuit(qubits)\n",
        "\n",
        "# prepare Hartree-Fock state as the reference state and append it\n",
        "# to the quantum circuit\n",
        "circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n",
        "\n",
        "# apply the UCJ operator to the reference state\n",
        "circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
        "circuit.measure_all()\n",
        "\n",
        "\n",
        "# ------------------------------ Step 2 ------------------------------\n",
        "\n",
        "isa_circuit = pass_manager.run(circuit)\n",
        "print(f\"Gate counts: {isa_circuit.count_ops()}\")\n",
        "\n",
        "\n",
        "# ------------------------------ Step 3 ------------------------------\n",
        "sampler = Sampler(mode=backend)\n",
        "sampler.options.environment.job_tags = [\"TUT_SQD\"]\n",
        "job = sampler.run([isa_circuit], shots=100_000)\n",
        "primitive_result = job.result()\n",
        "pub_result = primitive_result[0]\n",
        "\n",
        "\n",
        "# ------------------------------ Step 4 ------------------------------\n",
        "\n",
        "bit_array = pub_result.data.meas\n",
        "num_valid = sum(\n",
        "    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
        ")\n",
        "valid_fraction = num_valid / bit_array.num_shots\n",
        "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")\n",
        "expected_fraction_random = (\n",
        "    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
        ")\n",
        "print(\n",
        "    f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n",
        "    f\"{expected_fraction_random}\"\n",
        ")\n",
        "# SQD options\n",
        "energy_tol = 1e-3\n",
        "occupancies_tol = 1e-3\n",
        "max_iterations = 5\n",
        "\n",
        "# Eigenstate solver options\n",
        "num_batches = 3\n",
        "samples_per_batch = 1000\n",
        "max_dim = 1000\n",
        "symmetrize_spin = True\n",
        "carryover_threshold = 1e-4\n",
        "max_cycle = 200\n",
        "\n",
        "# Use the Hartree-Fock configuration as an initial guess for the\n",
        "# orbital occupancies\n",
        "initial_occupancies = (\n",
        "    np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
        "    np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
        ")\n",
        "\n",
        "# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
        "# you can omit this step, in which case you would not specify the\n",
        "# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.\n",
        "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n",
        "\n",
        "# List to capture intermediate results\n",
        "result_history = []\n",
        "\n",
        "\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    hcore,\n",
        "    eri,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=nelec,\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=energy_tol,\n",
        "    occupancies_tol=occupancies_tol,\n",
        "    max_iterations=max_iterations,\n",
        "    sci_solver=sci_solver,\n",
        "    symmetrize_spin=symmetrize_spin,\n",
        "    max_dim=max_dim,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    carryover_threshold=carryover_threshold,\n",
        "    callback=callback,\n",
        "    seed=rng,\n",
        ")\n",
        "\n",
        "final_energy = result.energy + nuclear_repulsion_energy\n",
        "energy_error = final_energy - reference_energy\n",
        "print(f\"Final energy: {final_energy}\")\n",
        "print(f\"Final energy error: {energy_error}\")\n",
        "\n",
        "# Data for energies plot\n",
        "x1 = range(len(result_history))\n",
        "min_e = [\n",
        "    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
        "    for result in result_history\n",
        "]\n",
        "e_diff = [abs(e - reference_energy) for e in min_e]\n",
        "\n",
        "# Chemical accuracy (+/- 1.6 milli-Hartree)\n",
        "chem_accuracy = 0.0016\n",
        "\n",
        "# Data for avg spatial orbital occupancy\n",
        "y2 = np.sum(result.orbital_occupancies, axis=0)\n",
        "x2 = range(len(y2))\n",
        "\n",
        "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
        "\n",
        "# Plot energies\n",
        "axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n",
        "axs[0].set_xticks(x1)\n",
        "axs[0].set_xticklabels(x1)\n",
        "axs[0].set_yscale(\"log\")\n",
        "axs[0].set_ylim(1e-3, 1e-1)\n",
        "axs[0].axhline(\n",
        "    y=chem_accuracy,\n",
        "    color=\"#BF5700\",\n",
        "    linestyle=\"--\",\n",
        "    label=\"chemical accuracy\",\n",
        ")\n",
        "axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
        "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
        "axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n",
        "axs[0].legend()\n",
        "\n",
        "# Plot orbital occupancy\n",
        "axs[1].bar(x2, y2, width=0.8)\n",
        "axs[1].set_xticks(x2)\n",
        "axs[1].set_xticklabels(x2)\n",
        "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
        "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
        "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "405e89ea-57da-4021-bb18-91e8d583d310",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## 다음 단계\n",
        "\n",
        "<Admonition type=\"tip\" title=\"권장사항\">\n",
        "  이 글이 흥미로웠다면, 다음 자료도 참고해 보시기 바랍니다:\n",
        "\n",
        "  * [페르미온 격자 모델에 대한 샘플 기반 크릴로프 양자 대각화](/docs/tutorials/sample-based-krylov-quantum-diagonalization) - 변분적 가설 대신 시간 진화 회로를 활용한 관련 튜토리얼\n",
        "  * [Dice 솔버를 활용한 SQD 화학 워크플로우 확장](/docs/addons/qiskit-addon-sqd/guides/integrate-dice-solver) — 대각화 처리에 더 효율적인 Dice 소프트웨어 사용법을 소개하는 페이지\n",
        "  * `diagonalize_fermionic_hamiltonian`[ SQD 애드온 API 문서](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) - 함수 참조\n",
        "  * [*양자 중심 슈퍼컴퓨터에서 정확한 대각화 기법의 한계를 넘어선 화학*](https://www.science.org/doi/10.1126/sciadv.adu9991) - 이 튜토리얼의 기초가 된 논문\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "hours": 1.5,
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "qpuSeconds": 60
  },
  "nbformat": 4,
  "nbformat_minor": 5
}