{
"cells": [
{
"cell_type": "markdown",
"id": "1682996d",
"metadata": {},
"source": [
"---\n",
"title: \"Calculs de solvant implicite à l'aide d' Qiskit Serverless\"\n",
"description: \"Apprenez à calculer les effets du solvant implicite sur le matériel quantique à l'aide des méthodes SQD et IEF-PCM.\"\n",
"---\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d2c31ae8",
"metadata": {},
"source": [
"\n",
"\n",
"# Calculs de solvant implicite à l'aide d' Qiskit Serverless\n",
"\n",
"*Estimation de la durée d'exécution : 2 minutes sur un processeur Heron r2 (REMARQUE : il s'agit uniquement d'une estimation. (La durée d'exécution peut varier.)*\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",
"## Résultats d'apprentissage\n",
"\n",
"À l'issue de ce tutoriel, les utilisateurs devraient avoir compris :\n",
"\n",
"* Comment configurer et exécuter un workflow à distance à l'aide d' [Qiskit Serverless](/docs/guides/serverless)\n",
"* Comment calculer les effets du solvant implicite à l'aide d'un ordinateur quantique\n",
"\n",
"\n",
"\n",
"## Prérequis\n",
"\n",
"Nous recommandons aux utilisateurs de se familiariser avec le sujet suivant avant de suivre ce tutoriel :\n",
"\n",
"* [Diagonalisation quantique par échantillonnage](/learning/courses/quantum-diagonalization-algorithms/sqd-overview)\n",
"\n",
"\n",
"\n",
"## Arrière-plan\n",
"\n",
"Les calculs avec solvant implicite sont fréquemment utilisés en biophysique computationnelle. Ces modèles décrivent la manière dont un composé soluté interagit avec un solvant, sans modéliser directement le système solvant. On procède plutôt à une approximation dans laquelle le modèle du système de soluté est intégré dans une représentation mathématique d'un milieu diélectrique caractérisé de manière empirique. Cette approximation diélectrique interagit ensuite avec le soluté, qui est lui-même modélisé directement. Le milieu diélectrique influence les caractéristiques du système soluté, telles que l'énergie de son état fondamental, en interagissant avec son champ électronique. C'est un aspect important pour les modèles biophysiques utilisés, par exemple, dans la recherche pharmaceutique, car les composés se comportent différemment selon les environnements diélectriques. La modélisation d'un composé dans l'air (sous vide) donnera lieu à un comportement différent de celui observé lorsqu'on le modélise dans l'eau. Étant donné que les composés pharmaceutiques doivent pénétrer dans l'organisme humain, qui est lui-même composé principalement d'eau, il est utile de modéliser un composé dans une solution telle que l'eau plutôt que dans le vide. Grâce aux modèles à solvant implicite, nous pouvons obtenir ce comportement à moindre coût, même si le résultat final sera généralement moins précis que celui obtenu avec les modèles à solvant explicite, plus gourmands en ressources de calcul, qui créent des représentations directes à la fois des molécules de soluté et de solvant.\n",
"\n",
"Dans ce tutoriel, nous montrons comment un algorithme quantique, la diagonalisation quantique par échantillonnage (SQD), peut être intégré à un modèle de solvant implicite relativement peu gourmand en ressources de calcul. Dans cet exemple, nous décrivons le comportement de la méthylamine lorsqu'elle se dissout dans l'eau. Nous comparons l'algorithme quantique à une méthode classique de référence appelée CASCI et montrons que les résultats de ces calculs concordent étroitement. Nous présentons une architecture de supercalcul centrée sur la technologie quantique à petite échelle, qui transfère le post-traitement classique, très gourmand en ressources de calcul, de la partie échantillonnage quantique de la routine vers un environnement cloud au sein d' Qiskit Serverless. Le code met également en œuvre la parallélisation entre les cœurs de processeur accessibles à distance afin de réduire le temps de calcul.\n",
"\n",
"Qiskit Serverless est un cadre permettant d'exécuter des charges de travail quantiques et classiques distribuées sans avoir à gérer l'infrastructure. Il n'y a pas de provisionnement de serveurs (pas de déploiement d' EC2s s, de clusters ou de conteneurs Docker ), pas d'outils d'orchestration ( Kubernetes, Docker Swarm) et pas de surveillance ni de maintenance. Chaque tâche Serverless s'exécute dans un conteneur vierge, exécute votre code, puis s'arrête. Il n'y a pas de mémoire entre les tâches. Il vous suffit d'écrire votre code, puis de soumettre votre tâche. Dans le cadre d'une tâche sans serveur, un programme peut accéder en toute transparence aux backends d' IBM Quantum® s et y effectuer un post-traitement classique des résultats. Grâce à l' Qiskit Serverless, les utilisateurs peuvent accéder à des cœurs de processeur et à de la mémoire à distance disponibles en permanence, ce qui permet de répartir certaines charges de travail classiques sur des ressources distantes. Les utilisateurs bénéficient également d'avantages en matière de traitement parallèle des programmes, tout en évitant les problèmes courants liés à l'arrêt des appareils en cours d'exécution. Pour plus d'informations sur Qiskit Serverless, consultez sa [documentation](/docs/guides/serverless), ainsi que les ressources supplémentaires disponibles sur [GitHub](https://qiskit.github.io/qiskit-serverless/index.html).\n",
"\n",
"Ce tutoriel présente une application concrète des éléments suivants :\n",
"\n",
"* Diagonalisation quantique par échantillonnage\n",
"* Modèles de calcul client-serveur pour l'informatique quantique\n",
"\n",
"Ce tutoriel s'inspire et s'appuie sur des recherches menées à la Cleveland Clinic, telles que décrites dans [Kaliakin, Danil, et al. « Diagonalisation quantique basée sur des échantillons de solvants implicites. » The Journal of Physical Chemistry B 129.23 (2025) : 5788-5796](https://pubs.acs.org/doi/10.1021/acs.jpcb.5c01030), présentant l'ensemble du processus de calcul des solvants implicites et l'étendant à l'auto-cohérence itérative des solvants (« The Heartwood Algorithm », M. Motta, T. Pellegrini (2025), l'optimisation géométrique et la sélection automatique de la disposition des qubits. Veuillez vous reporter au [modèle de fonction SQD IEF-PCM pour Qiskit](/docs/guides/function-template-chemistry-workflow), développé conjointement par la Cleveland Clinic et l' IBM®, sur la base des recherches menées par la Cleveland Clinic, afin de bénéficier d'une interface simplifiée de type « boîte noire » permettant d'effectuer des calculs avec solvant implicite.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "55b94021",
"metadata": {},
"source": [
"\n",
"\n",
"## Exigences\n",
"\n",
"Avant de commencer ce tutoriel, assurez-vous d'avoir installé les éléments suivants :\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9456ea67",
"metadata": {},
"source": [
"* Qiskit SDK v2.0 ou une version ultérieure avec prise en charge [de la visualisation](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime v0.40 ou version ultérieure (`pip install qiskit-ibm-runtime`)\n",
"* Catalogue d' IBM s sur Qiskit `pip install qiskit_ibm_catalog`\n",
"* Qiskit : sans serveur ( IBM ) `pip install qiskit_serverless`\n",
"* Module complémentaire Qiskit : Diagonalisation quantique par échantillonnage (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",
"* Géométrique `pip install geometric`\n",
"\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "7db2e559",
"metadata": {},
"source": [
"\n",
"\n",
"## Configuration\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": [
"Créez un répertoire situé exactement à côté du programme principal du bloc-notes, nommé `source_files`. Vous placerez dans ce répertoire les fichiers d' Python s que vous souhaitez partager avec l'environnement informatique distant. Vous devez créer deux fichiers :\n",
"\n",
"* `source_files\\diagonalization_engine.py`\n",
"* `source_files\\classical_simulation.py`\n",
"\n",
"Cliquez pour développer le texte de chaque script ci-dessous, puis copiez-collez le contenu dans un fichier local en respectant ces chemins d'accès.\n",
"\n",
"\n",
" \n",
" \n",
" \n",
"\n",
" \n",
" \n",
" \n",
"\n",
"\n",
"\n",
" Pour plus d'informations, veuillez consulter le [guide sur les modèles de fonctions SQD IEF-PCM Qiskit](/docs/guides/function-template-chemistry-workflow) (élaboré conjointement par la Cleveland Clinic et IBM ) mentionné ci-dessus. Voir également la [`qiskit_addon_sqd`](/docs/addons/qiskit-addon-sqd) bibliothèque.\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_ibm_runtime import QiskitRuntimeService\n",
"from qiskit_serverless import distribute_task, get_arguments, get, save_result\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 = QiskitRuntimeService(\n",
" channel=os.environ.get(\"QISKIT_IBM_CHANNEL\"),\n",
" token=os.environ.get(\"QISKIT_IBM_TOKEN\"),\n",
" instance=os.environ.get(\"QISKIT_IBM_INSTANCE\"),\n",
")\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": [
"Nous devons partager le programme destiné à être exécuté dans l'environnement cloud et le télécharger à nouveau chaque fois que nous modifions son code 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",
"## Exemple de simulateur à petite échelle\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b0db693d",
"metadata": {},
"source": [
"Ce tutoriel n'utilise pas de simulateur à petite échelle, car son objectif est de présenter une application quantique évolutive qui dépasse le cadre de l'exploration par simulateur. Nous montrerons en revanche plus loin comment cette méthode peut être mise en œuvre à l'aide d'une méthode de comparaison classique de pointe appelée CASCI.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "431a5bd2-e6ed-471b-ad9e-c4edd27784a8",
"metadata": {},
"source": [
"\n",
"\n",
"## Exemple de matériel à grande échelle\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",
"## Étape 1 : Transposer les données classiques en un problème quantique\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "cf4bbd59",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.1: Initialiser l'objet « molécule » à l'aide de la géométrie moléculaire connue d' $\\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: Définir les effets de solvatation à l'aide du modèle de continuum polarisable (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: Optimisation géométrique à l'aide de 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: Préparez le modèle continu et l'objet de champ moyen avec les variables pertinentes\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",
"## Étape 2 : Optimiser le problème pour son exécution sur du matériel quantique\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "52506835",
"metadata": {},
"source": [
"Pour plus d'informations sur les fonctions d'aide présentées ici, consultez le tutoriel [intitulé « Diagonalisation quantique par échantillonnage d'un hamiltonien chimique](/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",
"## Étapes 3 et 4 : Exécuter à l'aide de Qiskit et effectuer le post-traitement à l'aide d' Qiskit Serverless\n",
"\n",
"Ici, nous combinons l'étape 3 (Exécution) et l'étape 4 (Post-traitement) car le contexte d'application du modèle de solvant implicite nécessite une boucle de rétroaction itérative qui exécute plusieurs cycles d'exécution et de post-traitement afin d'améliorer le résultat final du calcul.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8715bb6b",
"metadata": {},
"source": [
"\n",
"\n",
"### 3.1 Calculer l'énergie de Hartree-Fock restreinte\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: Définir une énergie de référence classique avec 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",
"## Configurer les paramètres de l'application\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 os\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",
"## Etapes suivantes\n",
"\n",
"\n",
" Si ce travail vous a paru intéressant, les ressources suivantes pourraient vous intéresser :\n",
"\n",
" * [L](https://github.com/qiskit-community/qrmi) 'interface QRMI (Quantum Resource Management Interface) — intègre les ressources quantiques et classiques dans des gestionnaires de charges de travail HPC tels que Slurm, élargissant ainsi le modèle de calcul distribué présenté dans ce tutoriel\n",
" * [Guide d' Qiskit Serverless](/docs/guides/serverless)\n",
" * [Qiskit Serverless GitHub](https://qiskit.github.io/qiskit-serverless/index.html)\n",
" * Tutoriel [sur la diagonalisation quantique par échantillonnage d'un hamiltonien chimique](/docs/tutorials/sample-based-quantum-diagonalization)\n",
" * [Présentation générale de la diagonalisation quantique par échantillonnage (SQD)](/docs/addons/qiskit-addon-sqd)\n",
" * [Déployer et exécuter un modèle de simulation de structure électronique utilisant un modèle de solvant implicite](/docs/guides/function-template-chemistry-workflow) (modèle de fonction SQD IEF-PCM Qiskit, développé conjointement par la Cleveland Clinic et l' IBM )\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"id": "a1b8767d",
"source": "© IBM Corp., 2017-2026"
}
],
"metadata": {
"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"
},
"hours": 1.5,
"qpuSeconds": 120
},
"nbformat": 4,
"nbformat_minor": 5
}