{ "cells": [ { "cell_type": "markdown", "id": "1682996d", "metadata": {}, "source": [ "---\n", "title: \"Qiskit Serverless 를 이용한 암시적 용매 계산\"\n", "description: \"SQD 및 IEF-PCM을 사용하여 양자 하드웨어에서 암시적 용매 효과를 계산하는 방법을 배워보세요.\"\n", "---\n", "\n" ] }, { "cell_type": "markdown", "id": "d2c31ae8", "metadata": {}, "source": [ "\n", "\n", "# Qiskit Serverless 를 이용한 암시적 용매 계산\n", "\n", "*예상 소요 시간: Heron r2 프로세서 기준 2분 (참고: 이는 예상치에 불과합니다.) (실제 실행 시간은 다를 수 있습니다.)*\n", "\n", "{/* cspell:ignore avas AVAS hcore textit TRIC dmas mocore ncore ncas mocas fermilevel ecore orbts iiter IITER edup textcoords xytext fontsize fontweight frameon */}\n", "\n" ] }, { "cell_type": "markdown", "id": "8bf80006", "metadata": {}, "source": [ "\n", "\n", "## 학습 성과\n", "\n", "* [Qiskit Serverless](/docs/guides/serverless) 를 사용하여 원격 워크플로를 구성하고 실행하는 방법\n", "* 양자 컴퓨터를 사용하여 암시적 용매 효과를 계산하는 방법\n", "\n", "\n", "\n", "## 전제조건\n", "\n", "* [샘플 기반 양자 대각화](/learning/courses/quantum-diagonalization-algorithms/sqd-overview)\n", "\n", "\n", "\n", "## 배경\n", "\n", "암시적 용매 계산은 계산 생물물리학 분야에서 자주 사용된다. 이 모델들은 용매 계를 직접 모델링하지 않은 채, 용질 화합물이 용매와 어떻게 상호작용하는지를 설명합니다. 대신, 용질 시스템 모델을 경험적으로 특성화된 유전체 매체의 수학적 표현으로 감싸는 근사법을 적용한다. 이러한 유전체 근사치는 용질과 상호작용하며, 용질 자체는 직접 모델링된다. 유전체는 전자장과 상호작용함으로써 용질계의 기저 상태 에너지와 같은 특성에 영향을 미친다. 이는 화합물이 다양한 유전율 환경에서 서로 다른 거동을 보이기 때문에, 예를 들어 신약 개발 등에 활용되는 생물물리학적 모델에 있어 중요한 요소입니다. 공기 중(진공 상태)에서 화합물을 모델링하면 물 속에서 모델링했을 때와는 다른 거동을 나타냅니다. 의약 화합물은 대부분 물로 이루어진 인체 내부로 흡수되어야 하므로, 진공 상태가 아닌 물과 같은 용액 속에서 화합물을 모델링하는 것이 유용합니다. 암시적 용매 모델을 사용하면 이러한 동작을 적은 계산 비용으로 구현할 수 있지만, 최종 결과는 일반적으로 용질과 용매 분자를 모두 직접 표현하는, 계산 비용이 더 많이 드는 명시적 용매 모델에 비해 근사적인 수준에 그치게 됩니다.\n", "\n", "이 튜토리얼에서는 양자 알고리즘인 샘플 기반 양자 대각화(SQD)를 계산 비용이 상대적으로 적은 암시적 용매 모델에 어떻게 적용할 수 있는지 보여드립니다. 이 예시에서는 메틸아민이 물에 녹을 때 어떤 반응을 보이는지 설명합니다. 우리는 양자 알고리즘을 CASCI라고 불리는 최신의 고전적 비교 방법과 비교하고, 두 계산 결과가 매우 유사함을 보여준다. 우리는 양자 중심의 초고성능 컴퓨팅 아키텍처를 소형화하여 구현했으며, 루틴 내 양자 샘플링 단계에서 발생하는 계산 집약적인 고전적 후처리 작업을 Qiskit Serverless 내의 클라우드 기반 환경으로 오프로드합니다. 또한 이 코드는 계산 시간을 단축하기 위해 원격으로 접근 가능한 CPU 코어 간에 병렬 처리를 수행합니다.\n", "\n", "Qiskit Serverless 인프라를 관리할 필요 없이 분산형 양자 및 고전 워크로드를 실행할 수 있는 프레임워크입니다. 서버 프로비저닝( EC2s, 클러스터, Docker 컨테이너의 생성)이 필요 없으며, 오케스트레이션 도구( Kubernetes, Docker Swarm)도 필요 없고, 모니터링이나 유지보수도 필요하지 않습니다. 각 서버리스 작업은 깨끗한 컨테이너에서 실행되어 코드를 처리한 후 종료됩니다. 작업 간에는 데이터가 유지되지 않습니다. 코드를 작성한 다음 작업을 제출하기만 하면 됩니다. 서버리스 작업 내에서 프로그램은 IBM Quantum® 백엔드에 원활하게 액세스하여, 해당 백엔드에서 결과를 처리할 수 있습니다. Qiskit Serverless 를 통해 사용자는 상시 가동되는 원격 CPU 코어와 메모리에 액세스할 수 있으며, 이를 통해 특정 기존 워크로드를 원격 리소스에 분산할 수 있습니다. 또한 사용자는 프로그램의 병렬 처리에서 일정한 이점을 얻을 수 있을 뿐만 아니라, 실행 도중 장치가 종료되는 데서 비롯되는 흔한 문제점도 피할 수 있습니다. Qiskit Serverless 에 대한 자세한 내용은 [해당](/docs/guides/serverless) 문서와 [GitHub](https://qiskit.github.io/qiskit-serverless/index.html) 에 게시된 추가 자료를 참조하십시오.\n", "\n", "이 튜토리얼에서는 다음 내용의 실제 적용 사례를 보여줍니다:\n", "\n", "* 샘플 기반 양자 대각화\n", "* 양자 컴퓨팅을 위한 클라이언트-서버 계산 모델\n", "\n", "이 튜토리얼은 다음 문헌에 기술된 클리블랜드 클리닉의 연구에서 영감을 받아 이를 바탕으로 작성되었습니다 [Kaliakin, Danil 외. \"암시적 용매 샘플 기반 양자 대각화.\" 『Journal of Physical Chemistry B』 129.23 (2025): 5788-5796](https://pubs.acs.org/doi/10.1021/acs.jpcb.5c01030), 이 논문은 암시적 용매 계산을 위한 전체 워크플로를 제시하고, 이를 반복적 용매 자기일관성(\"The Heartwood Algorithm\", M. Motta, T. Pellegrini, 2025), 기하학적 최적화, 그리고 큐비트 배열의 자동 선택. 암시적 용매 계산을 실행하기 위한 간소화된 블랙박스 인터페이스에 대해서는, 클리블랜드 클리닉(Cleveland Clinic)의 연구를 바탕으로 클리블랜드 클리닉과 퀀텀 [컴퓨팅 연구](/docs/guides/function-template-chemistry-workflow) 소( IBM® )가 공동으로 개발한 ‘SQD IEF-PCM Qiskit 함수 템플릿’을 참고하시기 바랍니다.\n", "\n" ] }, { "cell_type": "markdown", "id": "55b94021", "metadata": {}, "source": [ "\n", "\n", "## 요구사항\n", "\n", "이 튜토리얼을 시작하기 전에 다음 항목이 설치되어 있는지 확인하십시오:\n", "\n" ] }, { "cell_type": "markdown", "id": "9456ea67", "metadata": {}, "source": [ "* Qiskit SDK v2.0 또는 그 이후 버전( [시각화](/docs/api/qiskit/visualization) 기능 지원)\n", "* Qiskit Runtime v0.40 또는 그 이후 (`pip install qiskit-ibm-runtime`)\n", "* Qiskit IBM 카탈로그 `pip install qiskit_ibm_catalog`\n", "* Qiskit IBM 서버리스 `pip install qiskit_serverless`\n", "* Qiskit 애드온: 샘플 기반 양자 대각화(SQD) v0.12.0 `pip install qiskit_addon_sqd`\n", "* PySCF `pip install pyscf`\n", "* FFSIM `pip install ffsim`\n", "* Matplotlib `pip install matplotlib`\n", "* 지오메트릭 `pip install geometric`\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "7db2e559", "metadata": {}, "source": [ "\n", "\n", "## 설정\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2c88b910", "metadata": {}, "outputs": [], "source": [ "# Establish Quantum Resource connection\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "\n", "service = QiskitRuntimeService()\n", "\n", "backend = service.least_busy()\n", "print(f\"Using backend {backend.name}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "af9286cf", "metadata": {}, "outputs": [], "source": [ "# Establish Classical HPC Resource connection\n", "from qiskit_ibm_catalog import QiskitFunction, QiskitServerless\n", "\n", "client = QiskitServerless()" ] }, { "cell_type": "markdown", "id": "aee647af", "metadata": {}, "source": [ "메인 노트북 프로그램 바로 옆에. `source_files`라는 이름의 디렉터리를 생성합니다. 원격 컴퓨팅 환경과 공유하려는 `Python` 파일을 이 디렉터리에 넣어 두십시오. 다음 두 개의 파일을 만들어야 합니다:\n", "\n", "* `source_files\\diagonalization_engine.py`\n", "* `source_files\\classical_simulation.py`\n", "\n", "아래 각 스크립트의 텍스트를 클릭하여 펼친 다음, 해당 내용을 복사하여 다음 경로 이름의 로컬 파일에 붙여넣으세요.\n", "\n", "\n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " 자세한 내용은 앞서 언급된 [SQD IEF-PCM Qiskit 함수 템플릿 가이드](/docs/guides/function-template-chemistry-workflow) (클리블랜드 클리닉과 IBM 이 공동으로 개발)를 참조하시기 바랍니다. 라이브러리도 [`qiskit_addon_sqd`](/docs/addons/qiskit-addon-sqd) 함께 참조하십시오.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7347c7e5", "metadata": { "tags": [ "id-diagonalization" ] }, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "import numpy as np\n", "from json.encoder import JSONEncoder\n", "from json.decoder import JSONDecoder\n", "from functools import partial\n", "import os\n", "\n", "from qiskit_serverless import (\n", " distribute_task,\n", " get_arguments,\n", " get,\n", " save_result,\n", " get_runtime_service,\n", ")\n", "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", " solve_sci,\n", ")\n", "\n", "\n", "### Argument retrieval\n", "args = get_arguments()\n", "\n", "data = args[\"data\"] # Chemistry Data\n", "energy_tol = args[\"energy_tol\"] # SQD option\n", "occupancies_tol = args[\"occupancies_tol\"] # SQD option\n", "max_iterations = args[\"max_iterations\"] # SQD option\n", "symmetrize_spin = args[\"symmetrize_spin\"] # Eigenstate solver option\n", "carryover_threshold = args[\"carryover_threshold\"] # Eigenstate solver option\n", "num_batches = args[\"num_batches\"] # Eigenstate solver option\n", "samples_per_batch = args[\"samples_per_batch\"] # Eigenstate solver option\n", "max_cycle = args[\"max_cycle\"] # Eigenstate solver option\n", "mem = args[\"mem\"] # Memory per Worker\n", "\n", "\n", "# --- fan‑out target: 1 CPU + mem GB RAM per call -------------\n", "@distribute_task(target={\"cpu\": 1, \"mem\": mem * 1024**3})\n", "def _solve_sci_worker(\n", " ix, ci_strs, one_body_tensor, two_body_tensor, norb, nelec, spin_sq\n", "):\n", " print(f\">>>>> WORKER {ix} INITIATED\")\n", " res = solve_sci(\n", " ci_strs,\n", " one_body_tensor,\n", " two_body_tensor,\n", " norb=norb,\n", " nelec=nelec,\n", " spin_sq=spin_sq,\n", " )\n", "\n", " print(f\">>>>> WORKER {ix} COMPLETE\")\n", " return res\n", "\n", "\n", "def distribute_solve_sci_batch(\n", " ci_strings: list[tuple[np.ndarray, np.ndarray]],\n", " one_body_tensor: np.ndarray,\n", " two_body_tensor: np.ndarray,\n", " norb: int,\n", " nelec: tuple[int, int],\n", " *,\n", " spin_sq: float | None = None,\n", " **kwargs,\n", ") -> list[SCIResult]:\n", " \"\"\"Diagonalize Hamiltonian in subspaces, parallelizing across\n", " vCPUs in the Serverless environment.\n", "\n", " Args:\n", " ci_strings: List of pairs (strings_a, strings_b) of arrays of\n", " spin-alpha CI strings and spin-beta CI strings whose Cartesian\n", " product gives the basis of the subspace in which to perform a\n", " diagonalization.\n", " one_body_tensor: The one-body tensor of the Hamiltonian.\n", " two_body_tensor: The two-body tensor of the Hamiltonian.\n", " norb: The number of spatial orbitals.\n", " nelec: The numbers of alpha and beta electrons.\n", " spin_sq: Target value for the total spin squared for the ground state.\n", " If ``None``, no spin will be imposed.\n", " **kwargs: Keyword arguments to pass to\n", " `pyscf.fci.selected_ci.kernel_fixed_space`\n", " (https://pyscf.org/pyscf_api_docs/pyscf.fci.html#pyscf.fci.selected_ci.kernel_fixed_space\n", "\n", " Returns:\n", " The results of the diagonalizations in the subspaces given by ci_strings.\n", " \"\"\"\n", " inputs = [\n", " (ix, ci_strs, one_body_tensor, two_body_tensor, norb, nelec, spin_sq)\n", " for ix, ci_strs in enumerate(ci_strings)\n", " ]\n", "\n", " # fan‑out: spawn one worker per input tuple\n", " print(\">>>>> ENTERING WORKER FAN-OUT\")\n", " refs = [_solve_sci_worker(*input_) for input_ in inputs]\n", " print(\">>>>> WAITING ON WORKERS TO FINISH TASKS\")\n", "\n", " # fan‑in: block until every worker finishes\n", " results = get(refs)\n", " print(\">>>>> DISTRIBUTED JOBS COMPLETED\")\n", "\n", " return results\n", "\n", "\n", "# A caveat of executing a Python program remotely is\n", "# that the inputs to the remote program must be passed\n", "# over an internet network. Similarly, the outputs\n", "# must be passed back to the local program via the same\n", "# structure. Python objects are not always able to be\n", "# passed over a network, and must be encoded in a\n", "# JSON serializable format.\n", "i_data = JSONDecoder().decode(data)\n", "\n", "# i_data has all of the information needed from the\n", "# local program to pick up where the computation left off\n", "# after its submission to the remote environment.\n", "[\n", " job_id,\n", " hcore,\n", " eri,\n", " num_orbitals,\n", " nuclear_repulsion_energy,\n", " num_elec_a,\n", " num_elec_b,\n", "] = i_data\n", "\n", "# Re-convert data back into numpy format, after serialization\n", "hcore = np.array(hcore)\n", "eri = np.array(eri)\n", "nuclear_repulsion_energy = np.float64(nuclear_repulsion_energy)\n", "\n", "# Instantiate Runtime Service to retrieve the\n", "# bitstrings from the QPU job. We provided these\n", "# credentials upon Serverless setup.\n", "service = get_runtime_service()\n", "\n", "# retrieving the QPU job data from the Serverless side\n", "job = service.job(job_id)\n", "primitive_result = job.result()\n", "pub_result = primitive_result[0]\n", "bit_array = pub_result.data.meas # Getting the bitstrings\n", "\n", "# Pass options to the built-in eigensolver\n", "sci_solver = partial(\n", " distribute_solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle\n", ")\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\">>>>> SQD ITERATION {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\">>>>> SUBSAMPLE {i}\")\n", " print(f\">>>>> \\tENERGY: {result.energy + nuclear_repulsion_energy}\")\n", " print(\n", " f\">>>>> \\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=num_orbitals,\n", " nelec=(num_elec_a, num_elec_b),\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", " carryover_threshold=carryover_threshold,\n", " callback=callback,\n", " seed=12345,\n", ")\n", "\n", "print(\">>>>> EXACT DIAGONALIZATION COMPLETE. CLEANING UP, SERIALIZING DATA.\")\n", "# Numpy arrays are not JSON serializable.\n", "# Convert them to List objects before using the JSONEncoder\n", "o_data = JSONEncoder().encode(\n", " [\n", " result.energy + nuclear_repulsion_energy,\n", " result.energy,\n", " result.rdm1.tolist(),\n", " result.rdm2.tolist(),\n", " [x.tolist() for x in result.orbital_occupancies],\n", " [\n", " result.sci_state.nelec,\n", " result.sci_state.norb,\n", " [x.tolist() for x in result.sci_state.orbital_occupancies()],\n", " [x.tolist() for x in result.sci_state.rdm()],\n", " ],\n", " ]\n", ")\n", "\n", "# JSON-safe package\n", "save_result({\"outputs\": o_data}) # single JSON blob returned to client" ] }, { "cell_type": "code", "execution_count": null, "id": "41807b3a", "metadata": { "tags": [ "id-classical" ] }, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "from json.encoder import JSONEncoder\n", "from json.decoder import JSONDecoder\n", "\n", "from qiskit_serverless import get_arguments, save_result\n", "\n", "import pyscf\n", "from pyscf import gto, scf\n", "from pyscf.solvent import pcm\n", "from pyscf.mcscf import avas\n", "\n", "import psutil\n", "\n", "mem_info = (\n", " psutil.virtual_memory()\n", ") # Get information about virtual memory (RAM)\n", "total_ram_gb = mem_info.total / (1024**3) # Convert bytes to GB\n", "print(f\">>>>> SERVERLESS TOTAL RAM: {total_ram_gb:.2f} GB\")\n", "\n", "### Argument retrieval\n", "args = get_arguments()\n", "data = args[\"data\"] # Chemistry Data\n", "\n", "i_data = JSONDecoder().decode(data)\n", "[mol_geo, eps, ao_labels] = i_data\n", "\n", "print(\">>>>> DEFINING MOLECULE\")\n", "mol = gto.M()\n", "mol.atom = mol_geo\n", "mol.basis = \"cc-pVDZ\"\n", "mol.unit = \"Ang\"\n", "mol.charge = 0\n", "mol.spin = 0\n", "mol.verbose = 0\n", "\n", "print(\">>>>> BUILDING MOLECULE\")\n", "mol.build()\n", "\n", "print(\">>>>> DEFINING PCM\")\n", "cm = pcm.PCM(mol)\n", "cm.eps = eps # for water\n", "cm.method = \"IEF-PCM\"\n", "\n", "print(\">>>>> BUILDING RESTRICTED HARTREE FOCK\")\n", "mf = scf.RHF(mol).PCM(cm) # This is the Final SCF object\n", "mf.kernel(verbose=0)\n", "\n", "print(\">>>>> RUNNING AVAS\")\n", "avas_ = avas.AVAS(mf, ao_labels, with_iao=True, canonicalize=True, verbose=0)\n", "avas_.kernel()\n", "norb, ne_act, mo_avas = avas_.ncas, avas_.nelecas, avas_.mo_coeff\n", "\n", "print(\">>>>> STARTING CASCI\")\n", "mc_pcm = pyscf.mcscf.CASCI(mf, norb, ne_act).PCM(\n", " cm\n", ") # Make sure to decorate the CASCI object with PCM\n", "mc_pcm.mo_coeff = mo_avas\n", "# mc_pcm.max_memory = 140000\n", "\n", "(CASCI_E, _, _, _, _) = mc_pcm.kernel(verbose=0)\n", "\n", "print(f\">>>>> CASCI_E: {CASCI_E}\")\n", "o_data = JSONEncoder().encode([float(CASCI_E)])\n", "\n", "# JSON-safe package\n", "save_result({\"outputs\": o_data}) # single JSON blob returned to client" ] }, { "cell_type": "markdown", "id": "c94aebd8", "metadata": {}, "source": [ "클라우드 환경에서 실행될 프로그램을 공유해야 하며, 소스 코드를 수정할 때마다 이를 다시 업로드해야 합니다:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cc00669f", "metadata": {}, "outputs": [], "source": [ "client.upload(\n", " QiskitFunction(\n", " title=\"diagonalization_engine\",\n", " entrypoint=\"diagonalization_engine.py\", # lives in ./source_files\n", " working_dir=\"source_files\",\n", " )\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "ef726eee", "metadata": {}, "outputs": [], "source": [ "client.upload(\n", " QiskitFunction(\n", " title=\"classical_simulation\",\n", " entrypoint=\"classical_simulation.py\", # lives in ./source_files\n", " working_dir=\"source_files\",\n", " )\n", ")" ] }, { "cell_type": "markdown", "id": "0cd9185e", "metadata": {}, "source": [ "\n", "\n", "## 소규모 시뮬레이터 예시\n", "\n" ] }, { "cell_type": "markdown", "id": "b0db693d", "metadata": {}, "source": [ "이 튜토리얼은 시뮬레이터 탐구 단계를 넘어 확장 가능한 양자 애플리케이션을 보여주기 위한 것이므로, 소규모 시뮬레이터를 사용하지 않습니다. 대신, 이 글의 뒷부분에서 CASCI라고 불리는 기존의 최첨단 비교 기법을 사용하여 이 방법을 어떻게 구현할 수 있는지 보여드리겠습니다.\n", "\n" ] }, { "cell_type": "markdown", "id": "431a5bd2-e6ed-471b-ad9e-c4edd27784a8", "metadata": {}, "source": [ "\n", "\n", "## 대규모 하드웨어 예시\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5df67bda", "metadata": {}, "outputs": [], "source": [ "# This is a useful helper function that displays\n", "# remote job execution details to the user's local machine\n", "def feedback_serverless(serverless_job):\n", " import time\n", "\n", " # Wait for the job to execute\n", " print(f\">>>>> Serverless status: {serverless_job.job_id}\")\n", " timer = 0\n", " while timer < 10000:\n", " if (\n", " serverless_job.status() == \"QUEUED\"\n", " or serverless_job.status() == \"INITIALIZING\"\n", " or serverless_job.status() == \"RUNNING\"\n", " ):\n", " print(f\">>>>> [{timer}s] Serverless job {serverless_job.job_id}: \\\n", " {serverless_job.status()}\")\n", " time.sleep(10)\n", " timer += 10\n", "\n", " elif serverless_job.status() == \"ERROR\":\n", " print(\n", " f\">>>>> Serverless job {serverless_job.job_id}: {serverless_job.status()}\"\n", " )\n", " print(\">>>>> Logs:\")\n", " print(serverless_job.logs())\n", " break\n", "\n", " elif serverless_job.status() == \"DONE\":\n", " print(\n", " f\">>>>> Serverless job {serverless_job.job_id}: {serverless_job.status()}\"\n", " )\n", " break\n", "\n", " else:\n", " break\n", "\n", " return" ] }, { "attachments": {}, "cell_type": "markdown", "id": "988ee237", "metadata": {}, "source": [ "\n", "\n", "## 1단계: 고전적 입력을 양자 문제에 매핑하기\n", "\n" ] }, { "cell_type": "markdown", "id": "cf4bbd59", "metadata": {}, "source": [ "\n", "\n", "### 1.1: 알려진 $\\textit{a priori}$ 분자 기하 구조를 사용하여 분자 객체를 초기화합니다\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "566e06b4", "metadata": {}, "outputs": [], "source": [ "# Reference guide for building molecule structures:\n", "# https://pyscf.org/user/gto.html\n", "# Video tutorial on building molecular objects in PySCF:\n", "# https://www.youtube.com/watch?v=cNC2cY9E9j0\n", "\n", "molecule_name = \"Methylamine\"\n", "\n", "methylamine_geo = \"\"\"\n", " N -0.7154 0.0000 0.0000;\n", " C 0.7154 0.0000 0.0000;\n", " H 1.1069 0.0916 1.0174;\n", " H 1.0996 0.8349 -0.5930;\n", " H 1.0996 -0.9274 -0.4345;\n", " H -1.0625 0.8564 0.4294;\n", " H -1.0625 -0.7661 0.5753;\n", "\"\"\"" ] }, { "cell_type": "code", "execution_count": null, "id": "4a7aec02", "metadata": {}, "outputs": [], "source": [ "# Imports\n", "import pyscf\n", "from pyscf import gto # Deals with molecular initialization\n", "from pyscf import scf # Solvation methods\n", "\n", "# Explicitly defining the Methylamine molecule\n", "mol = gto.M()\n", "mol.atom = methylamine_geo\n", "mol.basis = \"cc-pVDZ\"\n", "mol.unit = \"Ang\"\n", "mol.charge = 0\n", "mol.spin = 0\n", "mol.verbose = 0\n", "\n", "mol.build()" ] }, { "cell_type": "markdown", "id": "99359d80", "metadata": {}, "source": [ "\n", "\n", "### 1.2: 극화 연속체 모델(PCM)을 사용하여 용매화 효과를 정의한다\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9294b528", "metadata": {}, "outputs": [], "source": [ "# You can explore other solvents (such as methanol) by\n", "# retrieving other dielectric parameters from:\n", "# https://gaussian.com/scrf/\n", "from pyscf.solvent import pcm\n", "\n", "eps_water = 78.3553 # If solvating in a different medium,\n", "# set this constant appropriately using a known value" ] }, { "cell_type": "code", "execution_count": null, "id": "ffccfb92", "metadata": {}, "outputs": [], "source": [ "cm = pcm.PCM(mol)\n", "cm.eps = eps_water # PySCF defaults to water solvation,\n", "# but here we show this solvation parameter explicitly\n", "\n", "cm.method = (\n", " \"IEF-PCM\" # Alternative solvation models include C-PCM, SS(V)PE, COSMO\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "bfef3c9b", "metadata": {}, "outputs": [], "source": [ "# Create a \"Restricted Hartree-Fock\" object for the solute,\n", "# then wrap the SCF object with a Polarizable Continuum Model\n", "mf_pcm0 = scf.RHF(mol).PCM(\n", " cm\n", ") # Restricted Hartree-Fock misses instantaneous correlations,\n", "# post-HF methods like CCSD, CI, MP2 might be worth exploring" ] }, { "cell_type": "markdown", "id": "b5e6c20a", "metadata": {}, "source": [ "\n", "\n", "### 1.3: TRIC을 이용한 기하학적 최적화\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "493402b3", "metadata": {}, "outputs": [], "source": [ "# Geometry optimization with geomeTRIC\n", "from pyscf.geomopt.geometric_solver import (\n", " optimize,\n", ") # GeomeTRIC under the hood, for geometry optimization\n", "\n", "mol_opt = optimize(\n", " mf_pcm0, tol_grad=3e-4, verbose=0\n", ") # Use geomeTRIC/TRIC under the hood" ] }, { "cell_type": "markdown", "id": "58c168bf", "metadata": {}, "source": [ "\n", "\n", "### 1.4: 관련 변수를 사용하여 연속체 모델과 평균장 객체를 준비합니다\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1834cb22", "metadata": {}, "outputs": [], "source": [ "from pyscf.mcscf import avas\n", "\n", "# Re-define PCM\n", "cm = pcm.PCM(mol_opt)\n", "cm.eps = eps_water # for water\n", "cm.method = \"IEF-PCM\"\n", "\n", "# Re-build Restricted Hartree Fock object\n", "mf_opt = scf.RHF(mol_opt).PCM(cm)\n", "mf_opt.kernel(verbose=0)\n", "\n", "# Run AVAS\n", "ao_labels = [\"C 2s\", \"C 2p\", \"N 2s\", \"N 2p\", \"H 1s\"]\n", "avas_ = avas.AVAS(\n", " mf_opt, ao_labels, with_iao=True, canonicalize=True, verbose=0\n", ")\n", "avas_.kernel()\n", "norb, ne_act, mo_avas = avas_.ncas, avas_.nelecas, avas_.mo_coeff\n", "\n", "num_elec_a = (ne_act + mol_opt.spin) // 2\n", "num_elec_b = (ne_act - mol_opt.spin) // 2" ] }, { "cell_type": "markdown", "id": "ac6f36e3", "metadata": {}, "source": [ "\n", "\n", "## 2단계: 양자 하드웨어 실행을 위해 문제 최적화하기\n", "\n" ] }, { "cell_type": "markdown", "id": "52506835", "metadata": {}, "source": [ "여기서 소개한 헬퍼 함수에 대한 자세한 내용은 [‘화학 해밀토니안의 샘플 기반 양자 대각화’](/docs/tutorials/sample-based-quantum-diagonalization) 튜토리얼을 참조하십시오.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dc27c42e", "metadata": {}, "outputs": [], "source": [ "# Standard SQD helper functions (From SQD Tutorial)\n", "\n", "from typing import Sequence\n", "import rustworkx\n", "from qiskit.providers import BackendV2\n", "from qiskit import QuantumCircuit, QuantumRegister\n", "\n", "from rustworkx import NoEdgeBetweenNodes, PyGraph\n", "\n", "IBM_TWO_Q_GATES = {\"cx\", \"ecr\", \"cz\"}\n", "\n", "\n", "def create_linear_chains(num_orbitals: int) -> PyGraph:\n", " \"\"\"In zig-zag layout, there are two linear chains (with connecting\n", " qubits between the chains). This function creates those two linear\n", " chains: a rustworkx PyGraph with two disconnected linear chains.\n", " Each chain contains `num_orbitals` number of nodes, that is, in the\n", " final graph there are `2 * num_orbitals` number of nodes.\n", "\n", " Args:\n", " num_orbitals (int): Number orbitals or nodes in each linear chain.\n", " They are also known as alpha-alpha interaction qubits.\n", "\n", " Returns:\n", " A rustworkx.PyGraph with two disconnected linear chains each with\n", " `num_orbitals` number of nodes.\n", " \"\"\"\n", " G = rustworkx.PyGraph()\n", "\n", " for n in range(num_orbitals):\n", " G.add_node(n)\n", "\n", " for n in range(num_orbitals - 1):\n", " G.add_edge(n, n + 1, None)\n", "\n", " for n in range(num_orbitals, 2 * num_orbitals):\n", " G.add_node(n)\n", "\n", " for n in range(num_orbitals, 2 * num_orbitals - 1):\n", " G.add_edge(n, n + 1, None)\n", "\n", " return G\n", "\n", "\n", "def create_lucj_zigzag_layout(\n", " num_orbitals: int, backend_coupling_graph: PyGraph\n", ") -> tuple[PyGraph, int]:\n", " \"\"\"This function creates the complete zigzag graph that 'can be mapped'\n", " to an IBM QPU with heavy-hex connectivity (the zigzag must be an\n", " isomorphic sub-graph to the QPU/backend coupling graph for it to be\n", " mapped). The zigzag pattern includes both linear chains (alpha-alpha\n", " interactions) and connecting qubits between the linear chains\n", " (alpha-beta interactions).\n", "\n", " Args:\n", " num_orbitals (int): Number of orbitals, that is, number of nodes in\n", " each alpha-alpha linear chain.\n", " backend_coupling_graph (PyGraph): The coupling graph of the backend\n", " on which the LUCJ ansatz will be mapped and run. This function takes\n", " the coupling graph as a undirected `rustworkx.PyGraph` where there\n", " is only one 'undirected' edge between two nodes, that is, qubits.\n", " Usually, the coupling graph of an IBM backend is directed (for\n", " example, Eagle devices such as ibm_brisbane) or may have two edges\n", " between two nodes (for example, Heron `ibm_torino`). A user\n", " needs to make such graphs undirected or remove duplicate edges\n", " (or do both) to make them compatible with this function.\n", "\n", " Returns:\n", " G_new (PyGraph): The graph with IBM backend compliant zigzag pattern.\n", " num_alpha_beta_qubits (int): Number of connecting qubits between the\n", " linear chains in the zigzag pattern. While we want as many\n", " connecting (alpha-beta) qubits between the linear (alpha-alpha)\n", " chains, we cannot accommodate all due to qubit and connectivity\n", " constraints of backends. This is the maximum number of connecting\n", " qubits the zigzag pattern can have while being backend compliant\n", " (that is, isomorphic to backend coupling graph).\n", " \"\"\"\n", " isomorphic = False\n", " G = create_linear_chains(num_orbitals=num_orbitals)\n", "\n", " num_iters = num_orbitals\n", " while not isomorphic:\n", " G_new = G.copy()\n", " num_alpha_beta_qubits = 0\n", " for n in range(num_iters):\n", " if n % 4 == 0:\n", " new_node = 2 * num_orbitals + num_alpha_beta_qubits\n", " G_new.add_node(new_node)\n", " G_new.add_edge(n, new_node, None)\n", " G_new.add_edge(new_node, n + num_orbitals, None)\n", " num_alpha_beta_qubits = num_alpha_beta_qubits + 1\n", " isomorphic = rustworkx.is_subgraph_isomorphic(\n", " backend_coupling_graph, G_new\n", " )\n", " num_iters -= 1\n", "\n", " return G_new, num_alpha_beta_qubits\n", "\n", "\n", "def lightweight_layout_error_scoring(\n", " backend: BackendV2,\n", " virtual_edges: Sequence[Sequence[int]],\n", " physical_layouts: Sequence[int],\n", " two_q_gate_name: str,\n", ") -> list[list[list[int], float]]:\n", " \"\"\"Lightweight and heuristic function to score isomorphic layouts. There\n", " can be many zigzag patterns, each with different set of physical qubits,\n", " that can be mapped to a backend. Some of them might include fewer noise\n", " qubits and couplings than others. This function computes a simple error\n", " score for each such layout. It sums up 2Q gate error for all couplings\n", " in the zigzag pattern (layout) and measurement of errors of physical\n", " qubits in the layout to compute the error score.\n", "\n", " Note:\n", " This lightweight scoring can be refined using concepts such as\n", " mapomatic.\n", "\n", " Args:\n", " backend (BackendV2): A backend.\n", " virtual_edges (Sequence[Sequence[int]]): Edges in the device-\n", " compliant zigzag pattern where nodes are numbered from 0 to (2 *\n", " num_orbitals + num_alpha_beta_qubits).\n", " physical_layouts (Sequence[int]): All physical layouts of the zigzag\n", " pattern that are isomorphic to each other and to the larger backend\n", " coupling map.\n", " two_q_gate_name (str): The name of the two-qubit gate of the\n", " backend. The name is used for fetching two-qubit gate error from\n", " backend properties.\n", "\n", " Returns:\n", " scores (list): A list of lists where each sublist contains two\n", " items. First item is the layout, and second item is a float\n", " representing error score of the layout. The layouts in the `scores`\n", " are sorted in the ascending order of error score.\n", " \"\"\"\n", " props = backend.properties()\n", " scores = []\n", " for layout in physical_layouts:\n", " total_2q_error = 0\n", " for edge in virtual_edges:\n", " physical_edge = (layout[edge[0]], layout[edge[1]])\n", " try:\n", " ge = props.gate_error(two_q_gate_name, physical_edge)\n", " except Exception:\n", " ge = props.gate_error(two_q_gate_name, physical_edge[::-1])\n", " total_2q_error += ge\n", " total_measurement_error = 0\n", " for qubit in layout:\n", " meas_error = props.readout_error(qubit)\n", " total_measurement_error += meas_error\n", " scores.append([layout, total_2q_error + total_measurement_error])\n", " return sorted(scores, key=lambda x: x[1])\n", "\n", "\n", "def _make_backend_cmap_pygraph(backend: BackendV2) -> PyGraph:\n", " graph = backend.coupling_map.graph\n", " if not graph.is_symmetric():\n", " graph.make_symmetric()\n", " backend_coupling_graph = graph.to_undirected()\n", "\n", " edge_list = backend_coupling_graph.edge_list()\n", " removed_edge = []\n", " for edge in edge_list:\n", " if set(edge) in removed_edge:\n", " continue\n", " try:\n", " backend_coupling_graph.remove_edge(edge[0], edge[1])\n", " removed_edge.append(set(edge))\n", " except NoEdgeBetweenNodes:\n", " pass\n", "\n", " return backend_coupling_graph\n", "\n", "\n", "def get_zigzag_physical_layout(\n", " num_orbitals: int, backend: BackendV2, score_layouts: bool = True\n", ") -> tuple[list[int], int]:\n", " \"\"\"The main function that generates the zigzag pattern\n", " with physical qubits that can be used as an `intial_layout` in a\n", " preset passmanager/transpiler.\n", "\n", " Args:\n", " num_orbitals (int): Number of orbitals.\n", " backend (BackendV2): A backend.\n", " score_layouts (bool): Optional. If `True`, it uses the\n", " `lightweight_layout_error_scoring` function to score the\n", " isomorphic layouts and returns the layout with\n", " fewer erroneous qubits.\n", " If `False`, returns the first isomorphic subgraph.\n", "\n", " Returns:\n", " A tuple of device compliant layout (list[int]) with zigzag pattern\n", " and an int representing number of alpha-beta-interactions.\n", " \"\"\"\n", " backend_coupling_graph = _make_backend_cmap_pygraph(backend=backend)\n", "\n", " G, num_alpha_beta_qubits = create_lucj_zigzag_layout(\n", " num_orbitals=num_orbitals,\n", " backend_coupling_graph=backend_coupling_graph,\n", " )\n", "\n", " isomorphic_mappings = rustworkx.vf2_mapping(\n", " backend_coupling_graph, G, subgraph=True\n", " )\n", " isomorphic_mappings = list(isomorphic_mappings)\n", "\n", " edges = list(G.edge_list())\n", "\n", " layouts = []\n", " for mapping in isomorphic_mappings:\n", " initial_layout = [None] * (2 * num_orbitals + num_alpha_beta_qubits)\n", " for key, value in mapping.items():\n", " initial_layout[value] = key\n", " layouts.append(initial_layout)\n", "\n", " two_q_gate_name = IBM_TWO_Q_GATES.intersection(\n", " backend.configuration().basis_gates\n", " ).pop()\n", "\n", " if score_layouts:\n", " scores = lightweight_layout_error_scoring(\n", " backend=backend,\n", " virtual_edges=edges,\n", " physical_layouts=layouts,\n", " two_q_gate_name=two_q_gate_name,\n", " )\n", "\n", " return scores[0][0][:-num_alpha_beta_qubits], num_alpha_beta_qubits\n", "\n", " return layouts[0][:-num_alpha_beta_qubits], num_alpha_beta_qubits" ] }, { "cell_type": "code", "execution_count": null, "id": "87f8e3ec", "metadata": {}, "outputs": [], "source": [ "from qiskit.transpiler import generate_preset_pass_manager\n", "import ffsim\n", "\n", "# Initial LUCJ ansatz layout\n", "initial_layout, _ = get_zigzag_physical_layout(norb, backend=backend)\n", "\n", "# Initialize a pass manager\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend, initial_layout=initial_layout\n", ")\n", "\n", "pass_manager.pre_init = ffsim.qiskit.PRE_INIT" ] }, { "attachments": {}, "cell_type": "markdown", "id": "b4d480b3", "metadata": {}, "source": [ "\n", "\n", "## 3단계 및 4단계: Qiskit을 사용하여 실행하고 Qiskit Serverless 를 사용하여 후처리하기\n", "\n", "여기서는 3단계(실행)와 4단계(후처리)를 결합합니다. 이는 암시적 용매 모델의 적용 맥락상, 최종 계산 결과를 개선하기 위해 실행과 후처리 과정을 여러 차례 반복하는 반복적 피드백 루프가 필요하기 때문입니다.\n", "\n" ] }, { "cell_type": "markdown", "id": "8715bb6b", "metadata": {}, "source": [ "\n", "\n", "### 3.1 제한 하트리-팍 에너지를 계산하라\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3aa879ad", "metadata": {}, "outputs": [], "source": [ "# Run the kernel to get the RHF energy\n", "mf_opt = scf.RHF(mol_opt).PCM(cm)\n", "hf_e = float(mf_opt.kernel())" ] }, { "cell_type": "code", "execution_count": null, "id": "480dd476", "metadata": {}, "outputs": [], "source": [ "print(f\"Restricted Hartree-Fock Energy: {hf_e}\")" ] }, { "cell_type": "markdown", "id": "fd0145c6", "metadata": {}, "source": [ "\n", "\n", "### 3.2: CASCI를 통해 표준 기준 에너지를 설정합니다\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f629119d", "metadata": {}, "outputs": [], "source": [ "# Setup a Serverless Client\n", "worker = client.load(\"classical_simulation\")" ] }, { "cell_type": "code", "execution_count": null, "id": "bb9d2407", "metadata": {}, "outputs": [], "source": [ "from json.encoder import JSONEncoder\n", "\n", "ao_labels = [\"C 2s\", \"C 2p\", \"N 2s\", \"N 2p\", \"H 1s\"]\n", "data_e = JSONEncoder().encode([mol_opt.tostring(), eps_water, ao_labels])" ] }, { "cell_type": "code", "execution_count": null, "id": "93874c58", "metadata": {}, "outputs": [], "source": [ "serverless_job = worker.run(data=data_e)" ] }, { "cell_type": "code", "execution_count": null, "id": "756f0fdf", "metadata": {}, "outputs": [], "source": [ "# Optionally, check the Serverless status feedback\n", "# Don't sit here and stare at the feedback unless debugging.\n", "# You can go develop something else while the Serverless job runs.\n", "_ = feedback_serverless(serverless_job)" ] }, { "cell_type": "code", "execution_count": null, "id": "d3ede43c", "metadata": {}, "outputs": [], "source": [ "# If you make a mistake and need to cancel something\n", "\n", "# for job in client.jobs():\n", "# job.cancel()" ] }, { "cell_type": "code", "execution_count": null, "id": "7379c49f", "metadata": {}, "outputs": [], "source": [ "from json.decoder import JSONDecoder\n", "\n", "CASCI_E = JSONDecoder().decode(serverless_job.result()[\"outputs\"])[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "0a5f4b88", "metadata": {}, "outputs": [], "source": [ "# We have approximated the red, classical baseline from\n", "# Figure 1 for Methanol (North-West panel)\n", "print(f\"CASCI/IEF-PCM(cc-pVDZ): E={CASCI_E}\")" ] }, { "cell_type": "markdown", "id": "b22a1b00", "metadata": {}, "source": [ "\n", "\n", "## 애플리케이션 매개변수 설정\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "618375af", "metadata": {}, "outputs": [], "source": [ "# Systematically vary these parameters to improve hardware results\n", "\n", "# Set to \"True\" to run on real hardware\n", "use_hardware = True\n", "\n", "# Error suppression/mitigation options\n", "# >> Configure within Sampler primitive\n", "\n", "# Transpiler Options\n", "optimization_level = 3\n", "\n", "# Heartwood algorithm options\n", "n_iter = 15 # How many update loops to run\n", "resample = 1 # (resample=1 -> resample the QPU after every\n", "# update loop; resample=n_iter -> sample QPU only once)\n", "shots = 10000\n", "\n", "# SQD options\n", "energy_tol = 1e-4\n", "occupancies_tol = 1e-3\n", "max_iterations = 12\n", "\n", "# Eigenstate solver options\n", "num_batches = 5\n", "samples_per_batch = 300\n", "symmetrize_spin = True\n", "carryover_threshold = 1e-5\n", "max_cycle = 200\n", "\n", "# Classical post-processing options\n", "local = (\n", " False # Remote, Serverless (False) versus Local Post-Processing (True)\n", ")\n", "mem = 16 # Memory allocated to each diagonalization worker (Gb)" ] }, { "cell_type": "code", "execution_count": null, "id": "ce75a126", "metadata": {}, "outputs": [], "source": [ "# Heartwood algorithm subroutines\n", "import numpy as np\n", "import pyscf\n", "from pyscf import ao2mo, cc\n", "from functools import reduce\n", "import ffsim\n", "\n", "from json.encoder import JSONEncoder\n", "from json.decoder import JSONDecoder\n", "import time\n", "\n", "from functools import partial\n", "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", " solve_sci_batch,\n", ")\n", "\n", "\n", "def update_rdm(casci_object, dmas):\n", " \"\"\"\n", " Inputs:\n", " mc -> CASCI object\n", " dmas -> Spin-summed 1-particle reduced density matrix\n", "\n", " This function returns the CASCI/SQD one-body density matrix in\n", " the full basis of atomic orbitals, written as the sum (last line)\n", " of two terms:\n", " - a contribution from the core orbitals,\n", " np.dot(mocore, mocore.conj().T) * 2, (core = inactive and doubly-occupied)\n", " - a contribution from the active-space orbitals and electrons (dmas)\n", " rotated from the active-space to the AO basis (the reduce operation)\n", "\n", " Outputs:\n", " rho_approximation: The CASCI/SQD one-body density matrix\n", " in the full basis of atomic orbitals\n", "\n", " \"\"\"\n", " mo_coeff = casci_object.mo_coeff\n", " ncore = casci_object.ncore\n", " ncas = casci_object.ncas\n", " mocore = mo_coeff[:, :ncore]\n", " mocas = mo_coeff[:, ncore : ncore + ncas]\n", " dm1 = np.dot(mocore, mocore.conj().T) * 2\n", "\n", " rho_approximation = dm1 + reduce(np.dot, (mocas, dmas, mocas.conj().T))\n", " return rho_approximation\n", "\n", "\n", "def run_active_space_calculation(\n", " h1e_cas, h2e_cas, norb, ne_act, orbs, fermilevel, ecore\n", "):\n", " # ----- perform an HF and a CCSD calculation in the active space\n", " from pyscf import tools\n", " from datetime import datetime\n", "\n", " now = datetime.now().strftime(\"%H:%M:%S\")\n", " print(\">>>>> ACTIVE SPACE CALCULATIONS \")\n", " tools.fcidump.from_integrals(\n", " f\"as_fcidump_{now}.txt\",\n", " h1e_cas,\n", " h2e_cas,\n", " norb,\n", " ne_act,\n", " ms=0,\n", " nuc=ecore,\n", " ) # Forcefully represents the active space in the correct structure\n", " mf_as = tools.fcidump.to_scf(f\"as_fcidump_{now}.txt\")\n", " os.remove(f\"as_fcidump_{now}.txt\")\n", " mf_as.kernel()\n", " print(\">>>>> RUNNING CCSD\")\n", "\n", " mf_cc = cc.CCSD(mf_as)\n", " mf_cc.kernel()\n", " orbts = mf_as.mo_coeff\n", " t1, t2 = mf_cc.t1, mf_cc.t2\n", " print(\">>>>> UPDATED t1, t2 PARAMETERS\")\n", "\n", " # ----- update the HF orbitals\n", " active = list(\n", " range(fermilevel - ne_act // 2, fermilevel - ne_act // 2 + norb)\n", " )\n", " orbs[:, active] = np.dot(orbs[:, active], orbts)\n", " return orbs, t1, t2\n", "\n", "\n", "def get_lucj(norb, num_elec_a, num_elec_b, t1, t2, n_reps=1):\n", " print(\">>>>> CONSTRUCTING LUCJ CIRCUIT\")\n", "\n", " alpha_alpha_indices = [(p, p + 1) for p in range(norb - 1)]\n", " alpha_beta_indices = [(p, p) for p in range(0, norb, 4)]\n", "\n", " ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n", " t1=t1, # <---- Update t1 each loop\n", " t2=t2, # <---- Update t2 each loop\n", " n_reps=n_reps,\n", " interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),\n", " )\n", " nelec = (num_elec_a, num_elec_b)\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\n", " # and append it 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", " return circuit\n", "\n", "\n", "# Classical diagonalization engine sent to HPC\n", "def classically_diagonalize(\n", " bit_array=None, # Bit string array (only needed if locally processing data)\n", " nuclear_repulsion_energy=None, # Electronic energy from the core orbitals\n", " hcore=None, # 1-electron hamiltonian integrals\n", " eri=None, # 2-electron hamiltonian integrals\n", " num_orbitals=None, # Number of spatial orbitals\n", " nelec=None, # Number of electrons\n", " num_elec_a=None, # Alpha orbitals\n", " num_elec_b=None, # Beta orbitals\n", " job_id=None, # QPU bitstring Job ID\n", " client=None, # Diagonalization engine worker\n", " energy_tol=1e-4, # SQD option\n", " occupancies_tol=1e-3, # SQD option\n", " max_iterations=12, # SQD option\n", " num_batches=8, # Eigenstate solver option\n", " samples_per_batch=300, # Eigenstate solver option\n", " symmetrize_spin=False, # Eigenstate solver option\n", " carryover_threshold=1e-5, # Eigenstate solver option\n", " max_cycle=200, # Eigenstate solver option\n", " local=True, # Remote vs Local Diagonalization\n", " mem=16.0, # Memory per Serverless Worker (Gb)\n", "):\n", " print(\">>>>> STARTING DIAGONALIZATION ENGINE \")\n", " # Pass options to the built-in eigensolver. If you just want to use\n", " # the defaults, you can omit this step, in which case you would not\n", " # specify the sci_solver argument in the call to\n", " # diagonalize_fermionic_hamiltonian below.\n", " if local:\n", " sci_solver = partial(\n", " solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle\n", " )\n", "\n", " # List to capture intermediate results\n", " result_history = []\n", "\n", " def callback(results: list[SCIResult]):\n", " result_history.append(results)\n", " iteration = len(result_history)\n", " print(f\">>>>> SQD ITERATION {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\">>>>> SUBSAMPLE {i}\")\n", " print(\n", " f\">>>>> \\tENERGY: {result.energy + nuclear_repulsion_energy}\"\n", " )\n", " print(\n", " f\">>>>> \\tSUBSPACE DIMENSION: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", " result = diagonalize_fermionic_hamiltonian(\n", " hcore,\n", " eri,\n", " bit_array,\n", " samples_per_batch=samples_per_batch,\n", " norb=num_orbitals,\n", " nelec=(nelec // 2, nelec // 2),\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", " carryover_threshold=carryover_threshold,\n", " callback=callback,\n", " seed=12345,\n", " )\n", "\n", " result = (result.energy, result.rdm1, result.rdm2)\n", "\n", " else:\n", " # Serverless Logic\n", " print(\n", " f\">>>>> SENDING QISKIT RUNTIME JOB {job_id} TO QISKIT SERVERLESS\"\n", " )\n", "\n", " data = [\n", " job_id,\n", " hcore.tolist(),\n", " eri.tolist(),\n", " int(num_orbitals),\n", " float(nuclear_repulsion_energy),\n", " int(num_elec_a),\n", " int(num_elec_b),\n", " ]\n", "\n", " # Encode the execution dependencies with the JSONEncoder\n", " data_e = JSONEncoder().encode(data)\n", "\n", " # Send to Serverless\n", " worker = client.load(\"diagonalization_engine\")\n", " serverless_job = worker.run(\n", " data=data_e,\n", " energy_tol=energy_tol, # SQD option\n", " occupancies_tol=occupancies_tol, # SQD option\n", " max_iterations=max_iterations, # SQD option\n", " symmetrize_spin=symmetrize_spin, # Eigenstate solver option\n", " carryover_threshold=carryover_threshold, # Eigenstate solver option\n", " num_batches=num_batches, # Eigenstate solver option\n", " samples_per_batch=samples_per_batch, # Eigenstate solver option\n", " max_cycle=max_cycle, # Eigenstate solver option\n", " mem=mem, # Memory per Worker (Gb)\n", " )\n", "\n", " # Wait for the job to execute\n", " _ = feedback_serverless(serverless_job)\n", "\n", " o_data = JSONDecoder().decode(serverless_job.result()[\"outputs\"])\n", " result = (o_data[1], np.array(o_data[2]), np.array(o_data[3]))\n", "\n", " print(f\">>>>>>>>>> Active Space Energy: {o_data[1]}\")\n", " print(f\">>>>>>>>>> rdm1: {o_data[2]}\")\n", " print(f\">>>>>>>>>> rdm2: {o_data[3]}\")\n", "\n", " return result" ] }, { "cell_type": "code", "execution_count": null, "id": "d1e23956", "metadata": {}, "outputs": [], "source": [ "# The Heartwood algorithm\n", "import numpy as np\n", "\n", "from qiskit_ibm_runtime import SamplerV2 as Sampler\n", "from qiskit_addon_sqd.counts import generate_bit_array_uniform\n", "\n", "mc = pyscf.mcscf.CASCI(mf_opt, ncas=norb, nelecas=ne_act).PCM(cm)\n", "mc.with_solvent.method = mf_opt.with_solvent.method # Here we make sure\n", "# that mc is also using the same solvent method defined earlier (IEF-PCM)\n", "mc.with_solvent.eps = mf_opt.with_solvent.eps # Set the dielectric parameters\n", "mc.mo_coeff = mo_avas.copy() # Update the molecular orbitals to include\n", "# those computed in the presence of the solvent\n", "\n", "h1e_cas, ecore = (\n", " mc.get_h1eff()\n", ") # <-- h1eff is the 1-electron hamiltonian integrals. h1e_cas is\n", "# a common alias. ecore is the electronic energy from the core orbitals.\n", "h2e_cas = ao2mo.restore(\n", " 1, mc.get_h2eff(), norb\n", ") # <-- get the 2-electron hamiltonian integrals\n", "\n", "mc.mo_coeff, t1, t2 = run_active_space_calculation(\n", " h1e_cas,\n", " h2e_cas,\n", " norb,\n", " ne_act,\n", " mo_avas.copy(),\n", " mf_opt.mol.nelectron // 2,\n", " ecore,\n", ")\n", "\n", "# Sampler primitive options\n", "sampler = Sampler(mode=backend)\n", "\n", "# Explore error suppression techniques and see if they can improve result quality\n", "sampler.options.dynamical_decoupling.enable = True\n", "sampler.options.dynamical_decoupling.sequence_type = \"XY4\"\n", "sampler.options.twirling.enable_measure = True\n", "sampler.options.environment.job_tags = [\"TUT_ISC\"]\n", "# sampler.options.twirling.enable_gates = False\n", "# sampler.options.twirling.num_randomizations = 10\n", "# sampler.options.twirling.shots_per_randomization = 1024\n", "\n", "# initial approximation for rdm1\n", "with_solvent_e, with_solvent_v = None, None # Don't touch\n", "data = []\n", "for iiter in range(n_iter):\n", " print(f\">>>>> IMPLICIT SOLVENT ITERATION {iiter+1}/{n_iter}\")\n", " if with_solvent_v is not None:\n", " # Subsequent update loops enter here\n", " mc.get_hcore = lambda *args: mc._scf.get_hcore() + with_solvent_v\n", " else:\n", " # First update loop starts here\n", " # hcore is the CAS space (classically computed) 1-electron\n", " # hamiltonian, which we default to at the start of the routine.\n", " mc.get_hcore = (\n", " lambda *args: mc._scf.get_hcore()\n", " ) # REF: https://pyscf.org/pyscf_api_docs/pyscf.mcscf.html#pyscf.scf.hf.CASBase.get_h1cas\n", "\n", " # Alias mapping\n", " # hcore : h1e_cas : h1e_eff\n", " # nuclear_repulsion_energy : ecore\n", " # eri : h2e_cas : h2e_eff\n", "\n", " h1e_cas, ecore = (\n", " mc.get_h1eff()\n", " ) # <-- h1eff is the 1-electron hamiltonian integrals. h1e_cas is a\n", " # common alias. ecore is the electronic energy from the core orbitals.\n", " h2e_cas = ao2mo.restore(\n", " 1, mc.get_h2eff(), norb\n", " ) # <-- get the 2-electron hamiltonian integrals\n", "\n", " mc.mo_coeff, t1, t2 = run_active_space_calculation(\n", " h1e_cas,\n", " h2e_cas,\n", " norb,\n", " ne_act,\n", " mo_avas.copy(),\n", " mf_opt.mol.nelectron // 2,\n", " ecore,\n", " )\n", "\n", " if use_hardware:\n", " if (\n", " iiter % resample == 0\n", " ): # <-- Toggle how often you refresh your bitstrings here. The\n", " # developer suggests that you do it every time, but benevolently\n", " # provides the freedom to disagree with him via the resample\n", " # control variable.\n", " # The \"Quantum-Centric\" part\n", " print(\">>>>> GENERATING BITSTRINGS USING QUANTUM HARDWARE\")\n", " # LUCJ Ansatz construction\n", " circuit = get_lucj(norb, num_elec_a, num_elec_b, t1, t2, n_reps=1)\n", "\n", " print(f\">>>>> TRANSPILING LUCJ TO {backend.name}\")\n", " isa_circuit = pass_manager.run(circuit)\n", " print(f\">>>>> SUBMITTING ISA_CIRCUIT TO {backend.name}\")\n", " job = sampler.run(\n", " [isa_circuit], shots=shots\n", " ) # <----- Error Suppression/Mitigation configured performed upstream\n", " job_id = str(job.job_id())\n", "\n", " timer = 0\n", " while job.status() != \"DONE\":\n", " timer += 10\n", " print(\n", " f\">>>>> [{timer}s] RUNTIME JOB {job_id}: {job.status()}\"\n", " )\n", " time.sleep(10)\n", "\n", " primitive_result = job.result()\n", " print(f\">>>>> RETRIEVED {job_id} FROM {backend.name}\")\n", "\n", " pub_result = primitive_result[0]\n", " bit_array = pub_result.data.meas\n", "\n", " else:\n", " print(\">>>>> GENERATING BITSTRINGS CLASSICALLY\")\n", " rng = np.random.default_rng(24)\n", " bit_array = generate_bit_array_uniform(\n", " 100_000, 2 * norb, rand_seed=rng\n", " ) # <-- Sample bitstrings from a uniform distribution. This is\n", " # useful for debug, but runs out of steam on large systems\n", " job_id = float(\n", " \"nan\"\n", " ) # <-- we will check that valid job_id's were passed during grading\n", " local = True\n", "\n", " # The \"Classical Post-processing\" part\n", " result = classically_diagonalize(\n", " bit_array=bit_array,\n", " nuclear_repulsion_energy=ecore, # Electronic energy from the core orbitals\n", " hcore=h1e_cas, # 1-electron hamiltonian integrals\n", " eri=h2e_cas, # 2-electron hamiltonian integrals\n", " num_orbitals=norb, # Number of spatial orbitals\n", " nelec=ne_act, # Number of electrons\n", " num_elec_a=ne_act // 2, # Alpha orbitals\n", " num_elec_b=ne_act // 2, # Beta orbitals\n", " job_id=job_id, # QPU bitstring Job ID\n", " client=client, # Diagonalization engine worker\n", " energy_tol=energy_tol, # SQD option\n", " occupancies_tol=occupancies_tol, # SQD option\n", " max_iterations=max_iterations, # SQD option\n", " num_batches=num_batches, # Eigenstate solver option\n", " samples_per_batch=samples_per_batch, # Eigenstate solver option\n", " symmetrize_spin=symmetrize_spin, # Eigenstate solver option\n", " carryover_threshold=carryover_threshold, # Eigenstate solver option\n", " max_cycle=max_cycle, # Eigenstate solver option\n", " local=local, # Remote vs Local Diagonalization\n", " mem=mem, # Memory per Worker (Gb)\n", " )\n", "\n", " # e : SQD-based estimate of the energy\n", " # rdm1: Spin-summed 1-particle reduced density matrix\n", " e, rdm1 = result[0], result[1]\n", " rho_approximation = update_rdm(\n", " mc, rdm1\n", " ).copy() # <--- Reconstruct the one-body density matrix in the\n", " # atomic orbital basis to update the external potential due to\n", " # the solvent\n", "\n", " if with_solvent_e is not None:\n", " # Subsequent update loops enter here\n", " edup = np.einsum(\n", " \"ij,ji->\", with_solvent_v, rho_approximation\n", " ) # <-- edup: Incrementing the energy calculation with\n", " # subsequent iterations\n", " e += ecore + with_solvent_e - edup\n", "\n", " else:\n", " # First update loop enters here\n", " e += (\n", " ecore # Pulled from the CAS space object (molecule's core energy)\n", " )\n", "\n", " # Outputs:\n", " # with_solvent_e : scalar energy correction due to solvent polarization\n", " # with_solvent_v : Fock-like matrix to be added to the core Hamiltonian in SCF\n", " with_solvent_e, with_solvent_v = mc.with_solvent._get_vind(\n", " rho_approximation\n", " )\n", " data.append((iiter, float(e), job_id))\n", " print(f\">>>>> END IITER {iiter}\")\n", " print(f\">>>>> TOTAL ENERGY: {e}\\n\")" ] }, { "cell_type": "code", "execution_count": null, "id": "3905922f", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from matplotlib.ticker import ScalarFormatter\n", "\n", "\n", "def plot_data(data, baseline=0, name=None, save=False):\n", " x_vals, y_vals, job_ids = zip(*data)\n", " fig, ax = plt.subplots(figsize=(10, 6))\n", "\n", " # Plot line + markers\n", " ax.plot(\n", " x_vals,\n", " y_vals,\n", " color=\"navy\",\n", " linewidth=2,\n", " marker=\"o\",\n", " markersize=5,\n", " label=\"Energy trajectory\",\n", " )\n", " ax.axhline(\n", " baseline,\n", " color=\"red\",\n", " linestyle=\"--\",\n", " linewidth=1.5,\n", " label=\"Reference energy\",\n", " )\n", "\n", " # Force plain formatting\n", " ax.yaxis.set_major_formatter(ScalarFormatter(useMathText=True))\n", " ax.ticklabel_format(style=\"plain\", axis=\"y\")\n", "\n", " # Annotate each point with its exact value\n", " for x, y, job_id in zip(x_vals, y_vals, job_ids):\n", " ax.annotate(\n", " f\"{y:.8f}, ID: {job_id}\",\n", " (x, y),\n", " textcoords=\"offset points\",\n", " xytext=(0, 8), # vertical offset\n", " ha=\"center\",\n", " fontsize=8,\n", " rotation=25,\n", " color=\"navy\",\n", " )\n", "\n", " # Annotate the Classical Reference line\n", " for x, y in zip([0.5], [baseline]):\n", " ax.annotate(\n", " f\"{y:.5f}\",\n", " (x, y),\n", " textcoords=\"offset points\",\n", " xytext=(0, 8), # vertical offset\n", " ha=\"center\",\n", " fontsize=8,\n", " rotation=25,\n", " color=\"red\",\n", " )\n", "\n", " # Titles, labels, etc\n", " ax.set_title(\n", " f\"SQD/IEF-PCM(cc-pVDZ) - {name}\\nEnergy Convergence\",\n", " fontsize=14,\n", " fontweight=\"bold\",\n", " pad=15,\n", " )\n", " ax.set_xlabel(\"Update Iterations\", fontsize=12)\n", " ax.set_ylabel(\"Total Energy (Hartrees)\", fontsize=12)\n", " ax.grid(True, linestyle=\"--\", linewidth=0.6, alpha=0.7)\n", " ax.legend(frameon=True, loc=\"best\")\n", " plt.tight_layout()\n", "\n", " if save:\n", " plt.savefig(f\"./results/{name}_energy_convergence.png\")\n", " return fig, ax" ] }, { "cell_type": "code", "execution_count": null, "id": "6fdb639b", "metadata": {}, "outputs": [], "source": [ "# Plot your data\n", "\n", "fig, ax = plot_data(data, baseline=CASCI_E, name=molecule_name, save=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "75f48e6a-c7e4-46f3-9d39-a7a877427a04", "metadata": {}, "source": [ "\n", "\n", "## 다음 단계\n", "\n", "\n", " 이 글이 흥미로웠다면, 다음 자료들도 참고해 보시기 바랍니다:\n", "\n", " * [양자 리소스 관리 인터페이스(QRMI)](https://github.com/qiskit-community/qrmi) — 양자 및 고전 리소스를 Slurm과 같은 HPC 워크로드 관리자에 통합하여, 이 튜토리얼에서 시연한 분산 컴퓨팅 패턴을 확장합니다\n", " * [Qiskit Serverless 가이드](/docs/guides/serverless)\n", " * [Qiskit Serverless GitHub](https://qiskit.github.io/qiskit-serverless/index.html)\n", " * [화학 해밀토니안의 샘플 기반 양자 대각화](/docs/tutorials/sample-based-quantum-diagonalization) 튜토리얼\n", " * [샘플 기반 양자 대각화(SQD) 개요](/docs/addons/qiskit-addon-sqd)\n", " * [암시적 용매 모델을 사용한 전자 구조 시뮬레이션 템플릿 배포 및 실행](/docs/guides/function-template-chemistry-workflow) (클리블랜드 클리닉과 IBM 이 공동 개발한 SQD IEF-PCM Qiskit 함수 템플릿)\n", "\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": 120 }, "nbformat": 4, "nbformat_minor": 5 }