Hamiltoniens pour la chimie quantique
Commençons par un bref aperçu du rôle que jouent les hamiltoniens dans le VQE.
Aperçu de l'hamiltonien dans le VQE
Dr. Victoria Lipinska nous guide à travers les hamiltoniens et explique comment les mapper pour les utiliser en informatique quantique.
Références
Les articles suivants sont cités dans la vidéo ci-dessus.
- Quantum Algorithms for Fermionic Simulations, Ortiz, et al.
- Simulating Chemistry using Quantum Computers, Kassal et al.
- A Comparison of the Bravyi–Kitaev and Jordan–Wigner Transformations for the Quantum Simulation of Quantum Chemistry, Tranter, et al.
- Quantum Chemistry in the Age of Quantum Computing, Cao, et al.
- Quantum computational chemistry, McArdle, et al.
- The Bravyi-Kitaev transformation for quantum computation of electronic structure, Seeley, et al., McArdle, et al.
Préparer les hamiltoniens pour la chimie quantique
Une bonne première étape pour appliquer l'informatique quantique à un problème de chimie consiste à définir un hamiltonien pour le système étudié. Ici, nous limiterons la discussion aux hamiltoniens de chimie quantique, car ces hamiltoniens nécessitent un certain mapping spécifique aux systèmes de fermions identiques.
En tant que spécialiste de la chimie quantique, tu as probablement déjà ton logiciel préféré pour modéliser les molécules, capable de générer un hamiltonien décrivant ton système d'intérêt. Ici, nous utiliserons du code reposant uniquement sur PySCF, numpy et Qiskit. Mais le processus de préparation de l'hamiltonien s'applique également aux solutions prépackagées. La seule différence entre cette approche et d'autres logiciels réside dans de légères différences de syntaxe ; certaines sont abordées dans la sous-section « Logiciels tiers » pour faciliter l'intégration des workflows existants.
Générer un hamiltonien de chimie quantique pour une utilisation sur les QPUs IBM Quantum® implique les étapes suivantes :
- Définir ta molécule (géométrie, spin, espace actif, etc.)
- Générer l'hamiltonien fermionique (opérateurs de création et d'annihilation)
- Mapper depuis l'hamiltonien fermionique vers un opérateur bosonique (ici, en utilisant des opérateurs de Pauli)
- Si tu utilises un logiciel tiers : gérer les éventuelles incompatibilités de syntaxe entre le logiciel générateur et Qiskit
L'hamiltonien fermionique est écrit en termes d'opérateurs fermioniques, et tient notamment compte du fait que les électrons sont des fermions indiscernables. Cela signifie qu'ils obéissent à des statistiques totalement différentes de celles des qubits discernables et bosoniques. D'où la nécessité du mapping.
Ceux et celles qui sont déjà familiers avec ces processus peuvent probablement passer cette section. Objectif :
L'objectif final est d'obtenir un hamiltonien de la forme :
# Added by doQumentation — required packages for this notebook
!pip install -q numpy openfermion pyscf qiskit
H = [(1, "XX"), (1, "YY"), (1, "ZZ")]
print(H)
[(1, 'XX'), (1, 'YY'), (1, 'ZZ')]
Ou
from qiskit.quantum_info import SparsePauliOp
H = SparsePauliOp(["XX", "YY", "ZZ"], coeffs=[1.0 + 0.0j, 1.0 + 0.0j, 1.0 + 0.0j])
print(H)
SparsePauliOp(['XX', 'YY', 'ZZ'],
coeffs=[1.+0.j, 1.+0.j, 1.+0.j])
On commence par importer quelques packages :
import numpy as np
from pyscf import ao2mo, gto, mcscf, scf
- Définir ta molécule
Ici, nous allons spécifier les attributs de la molécule d'intérêt. Dans cet exemple, nous avons choisi l'hydrogène diatomique (car les hamiltoniens résultants sont suffisamment courts pour être affichés). Le Python-based Simulations of Chemistry Framework (PySCF) dispose d'une vaste collection de modules de structure électronique qui peuvent être utilisés, entre autres, pour générer des hamiltoniens moléculaires adaptés au calcul quantique. Le guide PySCF Quickstart est une excellente ressource pour une description complète de toutes les variables et fonctionnalités. Nous n'en donnerons qu'un aperçu très succinct, car tout cela sera déjà familier pour beaucoup d'entre vous. Pour mieux comprendre, consulte PySCF. En bref :
distance peut être utilisé pour les molécules diatomiques, ou spécifie simplement les coordonnées cartésiennes pour chaque atome. Les distances sont exprimées en Angström.
gto génère des orbitales de type gaussien.
basis désigne les fonctions utilisées pour modéliser les orbitales moléculaires. Ici, 'sto-6g' est une base minimale courante, dont le nom provient de l'ajustement d'orbitales de type Slater (Slater-Type Orbitals) à l'aide de 6 orbitales gaussiennes primitives.
spin une valeur entière indiquant le nombre d'électrons non appariés (égal à ). Note que certains logiciels utilisent la multiplicité à la place ().
charge la charge de la molécule.
symmetry - le groupe de symétrie ponctuelle de la molécule, spécifié par une chaîne de caractères ou détecté automatiquement en définissant "symmetry = True". Ici, "Dooh" est le groupe de symétrie approprié pour les molécules diatomiques avec deux atomes de la même espèce.
distance = 0.735
a = distance / 2
mol = gto.Mole()
mol.build(
verbose=0,
atom=[
["H", (0, 0, -a)],
["H", (0, 0, a)],
],
basis="sto-6g",
spin=0,
charge=0,
symmetry="Dooh",
)
<pyscf.gto.mole.Mole at 0x7fc718f07610>
Garde à l'esprit qu'on peut décrire l'énergie totale (qui inclut l'énergie de répulsion nucléaire ainsi que l'énergie électronique), l'énergie orbitale électronique totale, ou l'énergie d'un sous-ensemble d'orbitales électroniques (le sous-ensemble complémentaire étant gelé). Dans le cas spécifique de , note les différentes énergies ci-dessous, et observe que l'énergie totale diminuée de l'énergie de répulsion nucléaire donne bien l'énergie électronique :
mf = scf.RHF(mol)
mf.scf()
print(
mf.energy_nuc(),
mf.energy_elec()[0],
mf.energy_tot(),
mf.energy_tot() - mol.energy_nuc(),
)
0.7199689944489797 -1.8455976628764188 -1.125628668427439 -1.8455976628764188
active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)
- Générer l'hamiltonien fermionique
scf désigne une large gamme de méthodes de champ auto-cohérent.
rhf comme dans mf = scf.RHF(mol) : mf est un solveur qui utilise le calcul Hartree-Fock restreint. Le noyau de celui-ci (E, ci-dessous) est l'énergie totale, incluant la répulsion nucléaire et les orbitales moléculaires.
mcscf est un package de champs auto-cohérents multi-configurations.
ao2mo est une transformation des orbitales atomiques vers les orbitales moléculaires.
Nous utilisons également les variables suivantes :
ncas : nombre d'orbitales dans l'espace actif complet
nelecas : nombre d'électrons dans l'espace actif complet
E1 = mf.kernel()
mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))
mo = mx.sort_mo(active_space, base=0)
E2 = mx.kernel(mo)[:2]
On souhaite obtenir un hamiltonien, généralement décomposé en énergie d'un cœur électronique (ecore, non impliqué dans la minimisation), opérateurs mono-électroniques (h1e) et énergies bi-électroniques (h2e). Ces termes sont explicitement extraits ci-dessous dans les deux dernières lignes.
h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)
Ces hamiltoniens sont actuellement des opérateurs fermioniques (création et annihilation), applicables aux systèmes de fermions (indiscernables), et donc soumis à l'antisymétrie par échange. Cela engendre des statistiques différentes de celles qui s'appliqueraient à un système discernable ou bosonique. Pour effectuer des calculs sur les QPUs IBM Quantum, nous avons besoin d'un opérateur bosonique décrivant l'énergie. Le résultat d'un tel mapping s'écrit conventionnellement en termes d'opérateurs de Pauli, car ils sont à la fois hermitiens et unitaires. Plusieurs mappings sont possibles. L'un des plus simples est la transformation de Jordan-Wigner.
- Mapper l'hamiltonien
Il convient de noter qu'il existe de nombreux outils pour mapper un hamiltonien chimique vers un hamiltonien adapté à l'exécution sur un ordinateur quantique. Ici, nous implémentons directement le mapping de Jordan-Wigner en utilisant uniquement PySCF, numpy et Qiskit. Nous commentons ci-dessous les considérations de syntaxe pour d'autres solutions. La fonction Cholesky nous aide à obtenir une décomposition de faible rang des termes bi-électroniques de l'hamiltonien.
def cholesky(V, eps):
# see https://arxiv.org/pdf/1711.02242.pdf section B2
# see https://arxiv.org/abs/1808.02625
# see https://arxiv.org/abs/2104.08957
no = V.shape[0]
chmax, ng = 20 * no, 0
W = V.reshape(no**2, no**2)
L = np.zeros((no**2, chmax))
Dmax = np.diagonal(W).copy()
nu_max = np.argmax(Dmax)
vmax = Dmax[nu_max]
while vmax > eps:
L[:, ng] = W[:, nu_max]
if ng > 0:
L[:, ng] -= np.dot(L[:, 0:ng], (L.T)[0:ng, nu_max])
L[:, ng] /= np.sqrt(vmax)
Dmax[: no**2] -= L[: no**2, ng] ** 2
ng += 1
nu_max = np.argmax(Dmax)
vmax = Dmax[nu_max]
L = L[:, :ng].reshape((no, no, ng))
print(
"accuracy of Cholesky decomposition ",
np.abs(np.einsum("prg,qsg->prqs", L, L) - V).max(),
)
return L, ng
Les fonctions identity et creators_destructors remplacent les opérateurs de création et d'annihilation de l'hamiltonien fermionique par des opérateurs de Pauli ; creators_destructors utilise le mapping de Jordan-Wigner.
def identity(n):
return SparsePauliOp.from_list([("I" * n, 1)])
def creators_destructors(n, mapping="jordan_wigner"):
c_list = []
if mapping == "jordan_wigner":
for p in range(n):
if p == 0:
ell, r = "I" * (n - 1), ""
elif p == n - 1:
ell, r = "", "Z" * (n - 1)
else:
ell, r = "I" * (n - p - 1), "Z" * p
cp = SparsePauliOp.from_list([(ell + "X" + r, 0.5), (ell + "Y" + r, -0.5j)])
c_list.append(cp)
else:
raise ValueError("Unsupported mapping.")
d_list = [cp.adjoint() for cp in c_list]
return c_list, d_list
Enfin, build_hamiltonian utilise les fonctions cholesky, identity et creators_destructors pour construire l'hamiltonien final adapté à l'exécution sur un ordinateur quantique.
def build_hamiltonian(ecore: float, h1e: np.ndarray, h2e: np.ndarray) -> SparsePauliOp:
ncas, _ = h1e.shape
C, D = creators_destructors(2 * ncas, mapping="jordan_wigner")
Exc = []
for p in range(ncas):
Excp = [C[p] @ D[p] + C[ncas + p] @ D[ncas + p]]
for r in range(p + 1, ncas):
Excp.append(
C[p] @ D[r]
+ C[ncas + p] @ D[ncas + r]
+ C[r] @ D[p]
+ C[ncas + r] @ D[ncas + p]
)
Exc.append(Excp)
# low-rank decomposition of the Hamiltonian
Lop, ng = cholesky(h2e, 1e-6)
t1e = h1e - 0.5 * np.einsum("pxxr->pr", h2e)
H = ecore * identity(2 * ncas)
# one-body term
for p in range(ncas):
for r in range(p, ncas):
H += t1e[p, r] * Exc[p][r - p]
# two-body term
for g in range(ng):
Lg = 0 * identity(2 * ncas)
for p in range(ncas):
for r in range(p, ncas):
Lg += Lop[p, r, g] * Exc[p][r - p]
H += 0.5 * Lg @ Lg
return H.chop().simplify()
Enfin, nous utilisons build_hamiltonian pour construire notre hamiltonien de qubits à partir des opérateurs de Pauli via la transformation de Jordan-Wigner. Cela nous donne également la précision de la décomposition de Cholesky utilisée.
H = build_hamiltonian(ecore, h1e, h2e)
print(H)
accuracy of Cholesky decomposition 2.220446049250313e-16
SparsePauliOp(['IIII', 'IIIZ', 'IZII', 'IIZI', 'ZIII', 'IZIZ', 'IIZZ', 'ZIIZ', 'IZZI', 'ZZII', 'ZIZI', 'YYYY', 'XXYY', 'YYXX', 'XXXX'],
coeffs=[-0.09820182+0.j, -0.1740751 +0.j, -0.1740751 +0.j, 0.2242933 +0.j,
0.2242933 +0.j, 0.16891402+0.j, 0.1210099 +0.j, 0.16631441+0.j,
0.16631441+0.j, 0.1210099 +0.j, 0.17504456+0.j, 0.04530451+0.j,
0.04530451+0.j, 0.04530451+0.j, 0.04530451+0.j])
Ce carnet d'exemples de molécules présente la configuration et les hamiltoniens pour plusieurs molécules de complexité variable ; avec quelques modifications, il devrait te permettre d'étudier la plupart des petites molécules.
Notons brièvement deux points importants à prendre en compte lors de la construction des opérateurs fermioniques d'une molécule. À mesure que le type de molécule change, sa symétrie change également. De même, le nombre d'orbitales avec diverses symétries, comme la symétrie cylindrique "A1", varie. Ces changements sont évidents même avec l'extension simple à LiH, comme on le voit ici :
distance = 1.56
mol = gto.Mole()
mol.build(
verbose=0,
atom=[["Li", (0, 0, 0)], ["H", (0, 0, distance)]],
basis="sto-6g",
spin=0,
charge=0,
symmetry="Coov",
)
mf = scf.RHF(mol)
E1 = mf.kernel()
# %% ----------------------------------------------------------------------------------------------
mx = mcscf.CASCI(mf, ncas=5, nelecas=(1, 1))
cas_space_symmetry = {"A1": 3, "E1x": 1, "E1y": 1}
mo = mcscf.sort_mo_by_irrep(mx, mf.mo_coeff, cas_space_symmetry)
E2 = mx.kernel(mo)[:2]
h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)
Il convient également de noter qu'on perd rapidement l'intuition sur l'hamiltonien final obtenu. L'hamiltonien de LiH (utilisant le mappeur Jordan-Wigner) se compose déjà de 276 termes.
len(build_hamiltonian(ecore, h1e, h2e))
accuracy of Cholesky decomposition 1.1102230246251565e-16
276
En cas de doute sur les symétries, on peut aussi générer des informations de symétrie pour la molécule en définissant symmetry = True et verbose = 4 :
distance = 1.56
mol = gto.Mole()
mol.build(
verbose=4,
atom=[["Li", (0, 0, 0)], ["H", (0, 0, distance)]],
basis="sto-6g",
spin=0,
charge=0,
symmetry=True,
)
System: uname_result(system='Linux', node='IBM-R912JTRT', release='5.10.102.1-microsoft-standard-WSL2', version='#1 SMP Wed Mar 2 00:30:59 UTC 2022', machine='x86_64') Threads 16
Python 3.11.12 (main, May 16 2025, 02:33:32) [GCC 11.4.0]
numpy 2.3.1 scipy 1.16.0 h5py 3.14.0
Date: Mon Jun 30 12:56:55 2025
PySCF version 2.9.0
PySCF path /home/porter284/.pyenv/versions/3.11.12/lib/python3.11/site-packages/pyscf
[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 2
[INPUT] num. electrons = 4
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry True subgroup None
[INPUT] Mole.unit = angstrom
[INPUT] Symbol X Y Z unit X Y Z unit Magmom
[INPUT] 1 Li 0.000000000000 0.000000000000 0.000000000000 AA 0.000000000000 0.000000000000 0.000000000000 Bohr 0.0
[INPUT] 2 H 0.000000000000 0.000000000000 1.560000000000 AA 0.000000000000 0.000000000000 2.947972754321 Bohr 0.0
nuclear repulsion = 1.01764848253846
point group symmetry = Coov
symmetry origin: [0. 0. 0.73699319]
symmetry axis x: [1. 0. 0.]
symmetry axis y: [0. 1. 0.]
symmetry axis z: [0. 0. 1.]
num. orbitals of irrep A1 = 4
num. orbitals of irrep E1x = 1
num. orbitals of irrep E1y = 1
number of shells = 4
number of NR pGTOs = 36
number of NR cGTOs = 6
basis = sto-6g
ecp = {}
CPU time: 9.85
<pyscf.gto.mole.Mole at 0x7fc719f94850>
Parmi d'autres informations utiles, cela renvoie à la fois point group symmetry = Coov et le nombre d'orbitales dans chaque représentation irréductible.
point group symmetry = Coov
num. orbitals of irrep A1 = 4
num. orbitals of irrep E1x = 1
num. orbitals of irrep E1y = 1
number of shells = 4
Cela ne t'indique pas nécessairement combien d'orbitales tu veux inclure dans ton espace actif, mais cela te permet de voir quelles orbitales sont présentes et quelles sont leurs symétries.
Spécifier la symétrie et les orbitales est souvent utile, mais tu peux aussi simplement spécifier le nombre d'orbitales à inclure. Considère le cas de l'éthène, ci-dessous. En utilisant verbose = 4, on peut afficher les symétries des différentes orbitales :
# Replace these variables with correct distances:
a = 1
b = 1
c = 1
# Build
mol = gto.Mole()
mol.build(
verbose=4,
atom=[
["C", (0, 0, a)],
["C", (0, 0, -a)],
["H", (0, c, b)],
["H", (0, -c, b)],
["H", (0, c, -b)],
["H", (0, -c, -b)],
],
basis="sto-6g",
spin=0,
charge=0,
symmetry=True,
)
System: uname_result(system='Linux', node='IBM-R912JTRT', release='5.10.102.1-microsoft-standard-WSL2', version='#1 SMP Wed Mar 2 00:30:59 UTC 2022', machine='x86_64') Threads 16
Python 3.11.12 (main, May 16 2025, 02:33:32) [GCC 11.4.0]
numpy 2.3.1 scipy 1.16.0 h5py 3.14.0
Date: Mon Jun 30 12:57:07 2025
PySCF version 2.9.0
PySCF path /home/porter284/.pyenv/versions/3.11.12/lib/python3.11/site-packages/pyscf
[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 6
[INPUT] num. electrons = 16
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry True subgroup None
[INPUT] Mole.unit = angstrom
[INPUT] Symbol X Y Z unit X Y Z unit Magmom
[INPUT] 1 C 0.000000000000 0.000000000000 1.000000000000 AA 0.000000000000 0.000000000000 1.889726124565 Bohr 0.0
[INPUT] 2 C 0.000000000000 0.000000000000 -1.000000000000 AA 0.000000000000 0.000000000000 -1.889726124565 Bohr 0.0
[INPUT] 3 H 0.000000000000 1.000000000000 1.000000000000 AA 0.000000000000 1.889726124565 1.889726124565 Bohr 0.0
[INPUT] 4 H 0.000000000000 -1.000000000000 1.000000000000 AA 0.000000000000 -1.889726124565 1.889726124565 Bohr 0.0
[INPUT] 5 H 0.000000000000 1.000000000000 -1.000000000000 AA 0.000000000000 1.889726124565 -1.889726124565 Bohr 0.0
[INPUT] 6 H 0.000000000000 -1.000000000000 -1.000000000000 AA 0.000000000000 -1.889726124565 -1.889726124565 Bohr 0.0
nuclear repulsion = 29.3377079104231
point group symmetry = D2h
symmetry origin: [0. 0. 0.]
symmetry axis x: [0. 1. 0.]
symmetry axis y: [1. 0. 0.]
symmetry axis z: [-0. -0. -1.]
num. orbitals of irrep Ag = 4
num. orbitals of irrep B2g = 2
num. orbitals of irrep B3g = 1
num. orbitals of irrep B1u = 4
num. orbitals of irrep B2u = 1
num. orbitals of irrep B3u = 2
number of shells = 10
number of NR pGTOs = 84
number of NR cGTOs = 14
basis = sto-6g
ecp = {}
CPU time: 9.92
<pyscf.gto.mole.Mole at 0x7fc719fa9290>
On obtient :
num. orbitals of irrep Ag = 4
num. orbitals of irrep B2g = 2
num. orbitals of irrep B3g = 1
num. orbitals of irrep B1u = 4
num. orbitals of irrep B2u = 1
num. orbitals of irrep B3u = 2
Mais plutôt que de spécifier toutes les orbitales par symétrie, on peut simplement écrire :
active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)
Dans cette approche, on sélectionne plusieurs orbitales proches du niveau de remplissage (valence et inoccupées). Ici, 5 orbitales ont été choisies pour être incluses dans l'espace actif (les 6e à 10e).
print(
mol.nelectron // 2 - 2,
mol.nelectron // 2 + 2,
)
6 10
- Logiciels tiers
Il existe plusieurs packages logiciels développés pour la chimie quantique, certains offrant plusieurs mappeurs et des outils pour restreindre les espaces actifs. Les étapes décrites ci-dessus sont générales et s'appliquent également aux logiciels tiers. Mais ces autres logiciels peuvent retourner des hamiltoniens dans un format non accepté par Qiskit. Par exemple, certains logiciels retournent des hamiltoniens de la forme :
H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3]
Note en particulier que les portes sont numérotées et que les opérateurs identité ne sont pas affichés. C'est le contraire des hamiltoniens utilisés dans Qiskit, qui écriraient le terme [Z2 Z3] sous la forme ZZII (les qubits 0 et 1 étant affectés par l'opérateur identité, les qubits 2 et 3 étant affectés par l'opérateur Z, ordonnés avec le qubit 0 à l'extrême droite).
Pour s'adapter à tout workflow existant, le bloc de code ci-dessous effectue la conversion d'une syntaxe à l'autre. La fonction convert_openfermion_to_qiskit prend en arguments un hamiltonien généré dans OpenFermion ou Tangelo (déjà mappé sur des opérateurs de Pauli à l'aide de n'importe quel mappeur disponible), et le nombre de qubits nécessaires pour la molécule.
from openfermion import QubitOperator
from qiskit.quantum_info import SparsePauliOp
def convert_openfermion_to_qiskit(
openfermion_operator: QubitOperator, num_qubits: int
) -> SparsePauliOp:
terms = openfermion_operator.terms
labels = []
coefficients = []
for term, constant in terms.items():
# Default set to identity
operator = list("I" * num_qubits)
# Iterate through PauliSum and replace I with Pauli
for index, pauli in term:
operator[index] = pauli
label = "".join(operator)
labels.append(label)
coefficients.append(constant)
return SparsePauliOp(labels, coefficients)
De plus, ce carnet Python contient un exemple de code complet pour migrer des hamiltoniens d'autres workflows logiciels vers Qiskit, y compris la conversion ci-dessus.
Tu disposes maintenant d'un arsenal d'outils pour obtenir l'hamiltonien dont tu as besoin pour effectuer des calculs de chimie quantique sur les ordinateurs quantiques IBM®.