Les Hamiltoniens pour la chimie quantique
Commençons par un bref aperçu du rôle que jouent les hamiltoniens dans l'EQV.
Le hamiltonien dans VQE Aperçu
Victoria Lipinska nous présente les hamiltoniens et la manière de les cartographier pour les utiliser dans l'informatique quantique.
Références
Les articles suivants sont référencés dans la vidéo ci-dessus.
- Algorithmes quantiques pour les simulations fermioniques, Ortiz, et al.
- Simulation de la chimie à l'aide d'ordinateurs quantiques, Kassal et al.
- Comparaison des transformations de Bravyi-Kitaev et de Jordan-Wigner pour la simulation quantique de la chimie quantique, Tranter, et al.
- Chimie quantique à l'ère de l'informatique quantique, Cao, et al.
- Chimie quantique computationnelle, McArdle, et al.
- La transformation de Bravyi-Kitaev pour le calcul quantique de la structure électronique, Seeley et al.
Préparer les Hamiltoniens à la chimie quantique
Une bonne première étape dans l'application de l'informatique quantique à un problème de chimie consiste à définir un hamiltonien pour le système en question. Nous limiterons ici la discussion aux hamiltoniens de la chimie quantique, car ces hamiltoniens nécessitent une cartographie spécifique aux systèmes de fermions identiques.
Si vous travaillez dans le domaine de la chimie quantique, vous disposez probablement déjà de votre logiciel favori pour la modélisation des molécules, qui peut générer un hamiltonien décrivant le système qui vous intéresse. Nous utiliserons ici un code construit uniquement sur PySCF, numpy et Qiskit. Mais le processus de préparation des hamiltoniens s'applique également aux solutions prêtes à l'emploi. La seule différence entre cette approche et d'autres logiciels sera des différences syntaxiques mineures; certaines d'entre elles sont abordées dans la sous-section "Logiciels tiers" afin de faciliter l'intégration des flux de travail existants.
La génération d'un hamiltonien de chimie quantique pour une utilisation sur IBM Quantum® QPUs implique les étapes suivantes :
- Définir votre molécule (géométrie, spin, espace actif, etc.)
- Générer l'hamiltonien fermionique (opérateurs de création et d'annihilation)
- Carte de l'hamiltonien fermionique vers un opérateur bosonique (dans ce contexte, en utilisant les opérateurs de Pauli)
- En cas d'utilisation d'un logiciel tiers : gérer les éventuelles incohérences syntaxiques entre le logiciel de génération et Qiskit
Le 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 complètement différentes de celles des qubits bosoniques distinguables. D'où le processus de cartographie.
Ceux d'entre vous qui sont déjà familiarisés avec ces processus peuvent probablement sauter cette section.
Objectif :
L'objectif final est d'obtenir un hamiltonien de la forme :
H = [(1, "XX"), (1, "YY"), (1, "ZZ")]
print(H)Output:
[(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)Output:
SparsePauliOp(['XX', 'YY', 'ZZ'],
coeffs=[1.+0.j, 1.+0.j, 1.+0.j])
Nous commencerons par importer quelques paquets :
import numpy as np
from pyscf import ao2mo, gto, mcscf, scf- Définir votre molécule
Nous spécifierons ici les attributs de la molécule d'intérêt. Dans cet exemple, nous avons choisi l'hydrogène diatomique (parce que les hamiltoniens résultants sont suffisamment courts pour être affichés).
Le site 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 de démarrage rapide PySCF 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 beaucoup d'entre vous le connaissent déjà. Pour mieux les comprendre, consultez le site PySCF. En bref :
la distance peut être utilisée pour les molécules diatomiques, ou simplement spécifier les coordonnées cartésiennes pour chaque atome. Les distances sont exprimées en Angstrom.
gto génère des orbitales de type gaussien.
la base fait référence aux fonctions utilisées pour modéliser les orbitales moléculaires. Ici, ' sto-6g ' est une base minimale commune, nommée pour l'ajustement des orbitales de type Slater à l'aide de 6 orbitales gaussiennes primitives.
spin une valeur entière indiquant le nombre d'électrons non appariés (égal à ). Certains logiciels utilisent la multiplicité à la place ( ).
charge la charge de la molécule.
symétrie - le groupe de symétrie ponctuelle de la molécule, soit spécifié par une chaîne de caractères, soit détecté automatiquement en définissant "symmetry = True". Ici, "Dooh" est le groupe de symétrie approprié pour les molécules diatomiques comportant deux espèces d'atomes identiques.
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",
)Output:
<pyscf.gto.mole.Mole at 0x7fc718f07610>
Il faut garder à l'esprit que l'on peut décrire l'énergie totale (qui comprend l'énergie de répulsion nucléaire ainsi que l'énergie électronique), l'énergie totale des orbitales électroniques ou l'énergie d'un sous-ensemble d'orbitales électroniques (le sous-ensemble complémentaire étant gelé). Dans le cas particulier de , notez les différentes énergies ci-dessous, et notez que l'énergie totale moins l'énergie de répulsion nucléaire donne en fait 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(),
)Output:
0.7199689944489797 -1.8455976628764188 -1.125628668427439 -1.8455976628764188
active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)- Générer un hamiltonien fermionique
scf fait référence à un large éventail de méthodes de champ autoconsistant.
rhf comme dans mf = scf.RHF (mol) dans mf est un solveur qui utilise le calcul Restricted Hartree Fock. Le noyau de cette énergie (E, ci-dessous) est l'énergie totale, y compris la répulsion nucléaire et les orbitales moléculaires.
mcscf est un ensemble de champs autoconsistants à configurations multiples.
ao2mo est une transformation d'orbitales atomiques en 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]Nous voulons un hamiltonien, qui est souvent séparé en énergie d'un noyau électronique (ecore, non impliqué dans la minimisation), en opérateurs à un seul électron ( h1e ) et en énergies à deux électrons ( h2e ). Celles-ci sont explicitement extraites 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 à des systèmes de fermions (indiscernables) et, par conséquent, sujets à l'antisymétrie en cas d'échange. Il en résulte une statique différente de celle qui s'appliquerait à un système distinguable ou bosonique. Pour effectuer des calculs sur IBM Quantum QPUs, nous avons besoin d'un opérateur bosonique décrivant l'énergie. Le résultat d'un tel mappage est conventionnellement écrit en termes d'opérateurs de Pauli, puisqu'ils sont à la fois hermitiens et unitaires. Il est possible d'utiliser plusieurs mappings. L'une des plus simples est la transformation de Jordan Wigner.
- Cartographie de l'hamiltonien
Il convient de noter qu'il existe de nombreux outils permettant de convertir un hamiltonien chimique en un hamiltonien adapté à un ordinateur quantique. Ici, nous implémentons la cartographie de Jordan Wigner directement en utilisant seulement PySCF, numpy et Qiskit. Nous commentons ci-dessous les considérations syntaxiques pour d'autres solutions.
La fonction de Cholesky nous permet d'obtenir une décomposition de faible rang des termes à deux électrons 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, ngLes fonctions identity et creators_destructors remplacent les opérateurs de création et d'annihilation du hamiltonien fermionique par des opérateurs de Pauli; creators_destructors utilise le mappage 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_listEnfin, build_hamiltonian utilise les fonctions cholesky, identity et creators_destructors pour créer l'hamiltonien final pouvant être exécuté 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 qubit à partir des opérateurs de Pauli en utilisant la transformation de Jordan-Wigner. Cela nous donne également la précision de la décomposition de Cholesky que nous avons utilisée.
H = build_hamiltonian(ecore, h1e, h2e)
print(H)Output:
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])
Cet exemple de cahier de molécules montre la configuration et les hamiltoniens de plusieurs molécules de complexité variable; avec quelques modifications, cela devrait vous permettre d'examiner la plupart des petites molécules.
Notons brièvement deux points importants à prendre en compte lors de la construction des opérateurs fermioniques pour une molécule. Lorsque le type de molécule change, la symétrie change également. De même, le nombre d'orbitales présentant diverses symétries, comme l'orbitale à symétrie cylindrique " A1 ", changera. Ces changements sont évidents même avec une simple extension à LiH,, comme on peut le voir 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 que l'on peut rapidement perdre l'intuition de l'hamiltonien final obtenu. L'hamiltonien pour LiH (en utilisant le cartographe de Jordan-Wigner) se compose déjà de 276 termes.
len(build_hamiltonian(ecore, h1e, h2e))Output:
accuracy of Cholesky decomposition 1.1102230246251565e-16
276
En cas de doute concernant les symétries, il est également possible de générer des informations sur la symétrie de 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,
)Output:
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, cette méthode 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 vous indique pas nécessairement le nombre d'orbitales que vous voulez inclure dans votre espace actif, mais cela vous aide à voir quelles orbitales sont présentes et leurs symétries.
Il est souvent utile de spécifier la symétrie et les orbitales, mais vous pouvez également spécifier le nombre d'orbitales que vous souhaitez inclure. Prenons le cas de l'éthène, ci-dessous. En utilisant verbose = 4, nous pouvons imprimer 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,
)Output:
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>
Nous obtenons :
nombre d'orbitales de l'irrep Ag = 4
nombre d'orbitales de l'irrep B2g = 2
nombre d'orbitales de l'irrep B3g = 1
nombre d'orbitales de l'irrep B1u = 4
nombre d'orbitales de l'irrep B2u = 1
nombre d'orbitales de l'irrep B3u = 2
Mais plutôt que de spécifier toutes les orbitales par symétrie, nous pouvons simplement écrire :
active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)Dans cette approche, nous prenons plusieurs orbitales proches du niveau de remplissage (de valence et inoccupées). Ici, 5 orbitales ont été sélectionnées pour être incluses dans l'espace actif (de la 6ème à la 10ème).
print(
mol.nelectron // 2 - 2,
mol.nelectron // 2 + 2,
)Output:
6 10
- logiciel tiers
Plusieurs logiciels ont été développés pour la chimie quantique, certains offrant plusieurs cartographies et des outils permettant de restreindre les espaces actifs. Les étapes décrites ci-dessus sont générales et s'appliquent également aux logiciels tiers. Mais cet autre logiciel peut renvoyer les hamiltoniens dans un format qui n'est pas accepté par Qiskit. Par exemple, certains logiciels renvoient des hamiltoniens de la forme :
H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3]
Notez en particulier que les portes sont numérotées et que les opérateurs d'identité ne sont pas représentés. Cela contraste avec les hamiltoniens utilisés dans Qiskit, qui écrirait le terme [Z2 Z3] comme ZZII (les qubits 0 et 1 agissant sur l'opérateur identité, les qubits 2 et 3 agissant sur l'opérateur Z, ordonnés avec le qubit 0 le plus à droite).
Pour tenir compte des flux de travail existants, le bloc de code ci-dessous permet de passer d'une syntaxe à l'autre. La fonction convert_openfermion_to_qiskit prend comme arguments un hamiltonien généré dans OpenFermion ou Tangelo (et déjà mappé sur les 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 des exemples complets de code pour migrer des hamiltoniens d'autres logiciels vers Qiskit, y compris la conversion ci-dessus.
Vous devriez maintenant disposer d'un arsenal d'outils pour obtenir l'hamiltonien dont vous avez besoin pour effectuer des calculs de chimie quantique sur les ordinateurs quantiques IBM®.