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