Hamiltoniány pro kvantovou chemii
Začněme stručným přehledem role, kterou hamiltoniány hrají ve VQE.
Přehled hamiltoniánu ve VQE
Dr. Victoria Lipinska nás provede hamiltoniány a tím, jak je mapovat pro použití v kvantových výpočtech.
Reference
Ve výše uvedeném videu se odkazuje na následující články.
- Quantum Algorithms for Fermionic Simulations, Ortiz a kol.
- Simulating Chemistry using Quantum Computers, Kassal a kol.
- A Comparison of the Bravyi–Kitaev and Jordan–Wigner Transformations for the Quantum Simulation of Quantum Chemistry, Tranter a kol.
- Quantum Chemistry in the Age of Quantum Computing, Cao a kol.
- Quantum computational chemistry, McArdle a kol.
- The Bravyi-Kitaev transformation for quantum computation of electronic structure, Seeley a kol., McArdle a kol.
Příprava hamiltoniánů pro kvantovou chemii
Dobrým prvním krokem při aplikaci kvantových výpočtů na problém chemie je definování hamiltoniánu pro zkoumaný systém. Zde omezíme diskusi na hamiltoniány kvantové chemie, protože tyto hamiltoniány vyžadují určité mapování specifické pro systémy identických fermionů.
Jako někdo pracující v kvantové chemii už pravděpodobně máš svůj oblíbený software pro modelování molekul, který dokáže vygenerovat hamiltonián popisující tvůj zkoumaný systém. Zde budeme používat kód postavený výhradně na PySCF, numpy a Qiskitu. Proces přípravy hamiltoniánu se však přenáší i na předpřipravená řešení. Jediný rozdíl mezi tímto přístupem a jiným softwarem budou drobné syntaktické odlišnosti; některé z nich jsou řešeny v podsekci „Software třetích stran" pro usnadnění integrace stávajících pracovních postupů.
Generování hamiltoniánu kvantové chemie pro použití na QPU IBM Quantum® zahrnuje následující kroky:
- Definuj svou molekulu (geometrie, spin, aktivní prostor atd.)
- Vygeneruj fermionový hamiltonián (kreační a anihilační operátory)
- Namapuj fermionový hamiltonián na bosonový operátor (v tomto kontextu pomocí Pauliho operátorů)
- Pokud používáš software třetích stran: Vyřeš případné syntaktické neshody mezi generujícím softwarem a Qiskitem
Fermionový hamiltonián je zapsán pomocí fermionových operátorů a zejména bere v úvahu, že elektrony jsou nerozlišitelné fermiony. To znamená, že se řídí zcela odlišnou statistikou než rozlišitelné, bosonové qubity. Odtud ten proces mapování.
Ti z vás, kteří jsou již s těmito procesy obeznámeni, mohou tuto sekci pravděpodobně přeskočit. Cíl:
Konečným cílem je získat hamiltonián ve tvaru:
# 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')]
Nebo
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])
Začneme importem několika balíčků:
import numpy as np
from pyscf import ao2mo, gto, mcscf, scf
- Definuj svou molekulu
Zde specifikujeme atributy zkoumané molekuly. V tomto příkladu jsme zvolili dvouatomový vodík (protože výsledné hamiltoniány jsou dostatečně krátké, aby se daly zobrazit). Python-based Simulations of Chemistry Framework (PySCF) má širokou kolekci modulů pro elektronovou strukturu, které lze mimo jiné použít ke generování molekulárních hamiltoniánů vhodných pro kvantové výpočty. Průvodce PySCF Quickstart je vynikajícím zdrojem pro úplný popis všech proměnných a funkcí. Uvedeme pouze ten nejzákladnější přehled, protože toto již bude mnoha z vás známé. Pro lepší pochopení navštiv PySCF. Stručně:
distance lze použít pro dvouatomové molekuly, nebo jednoduše zadej kartézské souřadnice pro každý atom. Vzdálenosti jsou v jednotkách Angstrom.
gto generuje orbitaly gaussovského typu.
basis odkazuje na funkce použité pro modelování molekulárních orbitalů. Zde 'sto-6g' je běžná minimální báze, pojmenovaná podle aproximace Slaterových orbitalů pomocí 6 primitivních Gaussových orbitalů.
spin celočíselná hodnota udávající počet nepárových elektronů (rovná se ). Všimni si, že některý software používá místo toho multiplicitu ().
charge náboj molekuly.
symmetry – grupa bodové symetrie molekuly, buď určená řetězcem, nebo automaticky detekovaná nastavením „symmetry = True". Zde „Dooh" je vhodná grupa symetrie pro dvouatomové molekuly se dvěma atomy stejného druhu.
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>
Měj na paměti, že lze popisovat celkovou energii (která zahrnuje jadernou odpudivou energii i elektronovou), celkovou energii elektronových orbitalů nebo energii nějaké podmnožiny elektronových orbitalů (s komplementární podmnožinou zamrznutou). V konkrétním případě si povšimni různých energií níže a všimni si, že celková energie bez jaderné odpudivé energie skutečně dává elektronovou energii:
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)
- Vygenerování fermionického hamiltoniánu
scf označuje širokou škálu metod self-consistent field (samosouhlasného pole).
rhf (jako ve mf = scf.RHF(mol) u mf) je solver, který používá výpočet Restricted Hartree Fock. Jádro (kernel) této metody (níže E) je celková energie včetně jaderného odpuzování a molekulárních orbitalů.
mcscf je balíček pro multi-configuration self-consistent fields (vícekonfigurační samosouhlasné pole).
ao2mo je transformace z atomových orbitalů na molekulární orbitaly.
Dále používáme tyto proměnné:
ncas: počet orbitalů v úplném aktivním prostoru
nelecas: počet elektronů v úplném aktivním prostoru
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]
Chceme hamiltonián a ten se často rozděluje na energii elektronového jádra (ecore, která nevstupuje do minimalizace), jednoelektronové operátory (h1e) a dvouelektronové energie (h2e). Ty jsou explicitně extrahovány na posledních dvou řádcích níže.
h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)
Tyto hamiltoniány jsou nyní fermionické (kreační a anihilační) operátory, aplikovatelné na systémy (nerozlišitelných) fermionů, a proto podléhají antisymetrii při záměně. Z toho plyne jiná statistika, než jaká by platila pro rozlišitelný nebo bosonický systém. Abychom mohli spouštět výpočty na IBM Quantum QPU, potřebujeme bosonický operátor popisující energii. Výsledek takového mapování se konvenčně zapisuje pomocí Pauliho operátorů, protože jsou zároveň hermitovské i unitární. Existuje několik mapování, která lze použít. Jedním z nejjednodušších je Jordanova-Wignerova transformace.
- Mapování hamiltoniánu
Je třeba poznamenat, že existuje mnoho nástrojů pro mapování chemického hamiltoniánu na hamiltonián vhodný pro běh na kvantovém počítači. Zde implementujeme Jordanovo-Wignerovo mapování přímo pouze pomocí PySCF, numpy a Qiskitu. Níže komentujeme syntaktické úvahy pro jiná řešení. Funkce Cholesky nám pomáhá získat nízkohodnostní rozklad (low-rank decomposition) dvouelektronových členů v hamiltoniánu.
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
Funkce identity a creators_destructors nahrazují kreační a anihilační operátory ve fermionickém hamiltoniánu Pauliho operátory; creators_destructors používá Jordanovo-Wignerovo mapování.
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
Nakonec build_hamiltonian používá funkce cholesky, identity a creators_destructors k vytvoření finálního hamiltoniánu vhodného pro běh na kvantovém počítači.
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()
Nakonec použijeme build_hamiltonian ke konstrukci našeho qubitového hamiltoniánu z Pauliho operátorů pomocí Jordanovy-Wignerovy transformace. Zároveň tím získáme přesnost Choleského rozkladu, který jsme použili.
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])
Tento ukázkový notebook s molekulami ukazuje nastavení a hamiltoniány pro několik molekul různé složitosti; s malými úpravami by ti měl umožnit prozkoumat většinu malých molekul.
Krátce si ještě všimněme dvou důležitých bodů, které je třeba zvážit při sestavování fermionových operátorů pro molekulu. Se změnou typu molekuly se mění symetrie. Podobně se mění i počet orbitalů s různými symetriemi, jako je válcově symetrická „A1". Tyto změny jsou patrné už při jednoduchém rozšíření na LiH, jak je vidět zde:
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)
Za zmínku také stojí, že člověk může rychle ztratit intuici pro výsledný hamiltonián. Hamiltonián pro LiH (s použitím Jordan-Wignerova mapperu) už obsahuje 276 členů.
len(build_hamiltonian(ecore, h1e, h2e))
accuracy of Cholesky decomposition 1.1102230246251565e-16
276
Pokud si nejsi jistý symetriemi, můžeš také vygenerovat nějaké informace o symetrii molekuly nastavením symmetry = True a 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>
Kromě jiných užitečných informací to vrací jak point group symmetry = Coov, tak i počet orbitalů v každé ireducibilní reprezentaci.
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
To ti nemusí nutně říct, kolik orbitalů chceš zahrnout do aktivního prostoru, ale pomůže ti to vidět, jaké orbitaly jsou přítomny a jaké mají symetrie.
Specifikace symetrie a orbitalů je často užitečná, ale můžeš také zadat počet orbitalů, které chceš zahrnout. Zvažme případ ethenu níže. Použitím verbose = 4 můžeme vytisknout symetrie jednotlivých orbitalů:
# 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>
Dostaneme:
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
Místo toho, abychom vyjmenovávali všechny orbitaly pomocí symetrie, můžeme jednoduše napsat:
active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)
V tomto přístupu bereš několik orbitalů v okolí úrovně zaplnění (valenční a neobsazené). Zde bylo pro zahrnutí do aktivního prostoru vybráno 5 orbitalů (šestý až desátý).
print(
mol.nelectron // 2 - 2,
mol.nelectron // 2 + 2,
)
6 10
- Software třetích stran
Pro kvantovou chemii bylo vyvinuto několik softwarových balíčků, z nichž některé nabízejí více mapperů a nástrojů pro omezování aktivních prostorů. Výše popsané kroky jsou obecné a platí i pro software třetích stran. Tento jiný software však může vracet hamiltoniány ve formátu, který Qiskit nepřijímá. Například některý software vrací hamiltoniány ve tvaru:
H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3]
Všimni si zejména, že hradla jsou očíslována a operátory identity se nezobrazují. To je v kontrastu s hamiltoniány používanými v Qiskitu, které by člen [Z2 Z3] zapsaly jako ZZII (na qubity 0 a 1 působí operátor identity, na qubity 2 a 3 působí operátor Z, přičemž qubit 0 je úplně vpravo).
Aby byly zohledněny jakékoli tvé stávající pracovní postupy, blok kódu níže převádí z jedné syntaxe do druhé. Funkce convert_openfermion_to_qiskit přijímá jako argumenty hamiltonián vygenerovaný v OpenFermion nebo Tangelo (a již namapovaný na Pauliho operátory pomocí libovolného dostupného mapperu) a počet qubitů potřebných pro danou molekulu.
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)
Dále tento Python notebook obsahuje kompletní ukázkový kód pro přenášení hamiltoniánů z jiných softwarových pracovních postupů do Qiskitu, včetně výše uvedeného převodu.
Nyní bys měl mít k dispozici arzenál nástrojů pro získání hamiltoniánu, který potřebuješ k provádění kvantově-chemických výpočtů na kvantových počítačích IBM®.