Hamiltoniani per la chimica quantistica
Cominciamo con una breve panoramica del ruolo degli hamiltoniani nella VQE.
L'Hamiltoniano in VQE Panoramica
La dottoressa Victoria Lipinska ci illustra gli hamiltoniani e come mapparli per utilizzarli nell'informatica quantistica.
Riferimenti
I seguenti articoli sono citati nel video di cui sopra.
- Algoritmi quantistici per simulazioni fermioniche, Ortiz, et al.
- Simulazione della chimica con i computer quantistici, Kassal et al.
- Confronto tra le trasformazioni di Bravyi-Kitaev e Jordan-Wigner per la simulazione della chimica quantistica, Tranter, et al.
- Quantum Chemistry in the Age of Quantum Computing, Cao, et al.
- Chimica quantistica computazionale, McArdle, et al.
- La trasformazione di Bravyi-Kitaev per il calcolo quantistico della struttura elettronica, Seeley et al.
Preparare gli hamiltoniani per la chimica quantistica
Un buon primo passo per applicare l'informatica quantistica a un problema di chimica è la definizione di un'hamiltoniana per il sistema di interesse. In questa sede, limiteremo la discussione agli hamiltoniani della chimica quantistica, poiché questi ultimi richiedono una mappatura specifica per i sistemi di fermioni identici.
Se lavorate nel campo della chimica quantistica, probabilmente avete già il vostro software preferito per la modellazione delle molecole, in grado di generare un'hamiltoniana che descrive il vostro sistema di interesse. In questa sede, utilizzeremo codice costruito esclusivamente su PySCF, numpy e Qiskit. Ma il processo di preparazione hamiltoniana si trasferisce anche alle soluzioni preconfezionate. L'unica differenza tra questo approccio e altri software sarà costituita da piccole differenze di sintassi; alcune di queste sono trattate nella sottosezione "Software di terze parti" per facilitare l'integrazione dei flussi di lavoro esistenti.
La generazione di un hamiltoniano di chimica quantistica da utilizzare sulle QPU IBM Quantum® comporta le seguenti fasi:
- Definire la molecola (geometria, spin, spazio attivo e così via)
- Generare l'hamiltoniana fermionica (operatori di creazione e annichilazione)
- Mappa dall'hamiltoniano fermionico a un operatore bosonico (in questo contesto, utilizzando gli operatori di Pauli)
- Se si utilizza un software di terze parti: gestire eventuali discrepanze di sintassi tra il software di generazione e Qiskit
L'hamiltoniana fermionica è scritta in termini di operatori fermionici e, in particolare, tiene conto del fatto che gli elettroni sono fermioni indistinguibili. Ciò significa che obbediscono a statistiche completamente diverse da quelle dei qubit bosonici distinguibili. Da qui il processo di mappatura.
Chi ha già familiarità con questi processi può saltare questa sezione.
L’obiettivo:
L'obiettivo finale è ottenere un'hamiltoniana della forma:
H = [(1, "XX"), (1, "YY"), (1, "ZZ")]
print(H)Output:
[(1, 'XX'), (1, 'YY'), (1, 'ZZ')]
Oppure
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])
Inizieremo importando alcuni pacchetti:
import numpy as np
from pyscf import ao2mo, gto, mcscf, scf- Definire la propria molecola
Qui si specificano gli attributi della molecola di interesse. In questo esempio, abbiamo scelto l'idrogeno biatomico (perché gli hamiltoniani risultanti sono abbastanza brevi da poter essere visualizzati).
Il framework Python - Simulations of Chemistry ( PySCF ) dispone di un'ampia collezione di moduli di struttura elettronica che possono essere utilizzati, tra l'altro, per generare hamiltoniani molecolari adatti al calcolo quantistico. La guida rapida di PySCF è una risorsa eccellente per una descrizione completa di tutte le variabili e le funzionalità. Ne daremo solo una panoramica sommaria, poiché molti di voi ne saranno già a conoscenza. Per capirli meglio, visitate il sito PySCF. Brevemente:
la distanza può essere utilizzata per le molecole biatomiche oppure è sufficiente specificare le coordinate cartesiane per ciascun atomo. Le distanze sono in unità di Angstrom.
gto genera orbitali di tipo gaussiano.
il termine base si riferisce alle funzioni utilizzate per modellare gli orbitali molecolari. Qui ' sto-6g ' è una base minima comune, chiamata così per adattare gli orbitali di tipo Slater utilizzando 6 orbitali gaussiani primitivi.
spin un valore intero che indica il numero di elettroni spaiati (uguale a ). Si noti che alcuni software utilizzano invece la molteplicità ( ).
carica la carica della molecola.
simmetria - il gruppo di simmetria puntiforme della molecola, specificato con una stringa o rilevato automaticamente impostando "simmetria = Vero". Qui "Dooh" è il gruppo di simmetria appropriato per le molecole biatomiche con due atomi della stessa specie.
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>
Si tenga presente che si può descrivere l'energia totale (che comprende l'energia di repulsione nucleare e quella elettronica), l'energia totale degli orbitali elettronici o l'energia di un sottoinsieme di orbitali elettronici (con il sottoinsieme complementare congelato). Nel caso specifico di , si notino le diverse energie di seguito riportate e si noti che l'energia totale meno l'energia di repulsione nucleare produce di fatto l'energia elettronica:
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)- Generare l'hamiltoniana fermionica
scf si riferisce a un'ampia gamma di metodi di campo autoconsistenti.
rhf come in mf = scf.RHF (mol) in mf è un solutore che utilizza il calcolo Restricted Hartree Fock. Il nocciolo di questo (E, sotto) è l'energia totale, che include la repulsione nucleare e gli orbitali molecolari.
mcscf è un pacchetto di campi autoconsistenti multi-configurazione.
ao2mo è una trasformazione da orbitali atomici a orbitali molecolari.
Utilizziamo anche le seguenti variabili:
ncas : numero di orbitali nello spazio attivo completo
nelecas : numero di elettroni nello spazio attivo completo
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]Vogliamo un'hamiltoniana, che viene spesso separata in energia di un nucleo elettronico (ecore, non coinvolto nella minimizzazione), operatori a singolo elettrone ( h1e ) ed energie a due elettroni ( h2e ). Questi sono esplicitamente estratti nelle ultime due righe.
h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)Questi hamiltoniani sono attualmente operatori fermionici (di creazione e annichilazione), applicabili a sistemi di fermioni (indistinguibili) e, di conseguenza, soggetti ad antisimmetria per scambio. Ciò comporta una statica diversa da quella che si applicherebbe a un sistema distinguibile o bosonico. Per eseguire i calcoli sulle QPU di IBM Quantum, è necessario un operatore bosonico che descriva l'energia. Il risultato di una tale mappatura viene convenzionalmente scritto in termini di operatori di Pauli, poiché sono entrambi ermitiani e unitari. Si possono utilizzare diverse mappature. Una delle più semplici è la trasformazione di Jordan Wigner.
- Mappatura dell'hamiltoniano
Va notato che sono disponibili molti strumenti per mappare un'hamiltoniana chimica in una adatta all'esecuzione su un computer quantistico. Qui implementiamo la mappatura di Jordan Wigner direttamente utilizzando solo PySCF, numpy e Qiskit. Commentiamo di seguito le considerazioni sulla sintassi per altre soluzioni.
La funzione Cholesky ci aiuta a ottenere una decomposizione a basso rango dei termini a due elettroni dell'hamiltoniana.
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, ngLe funzioni identity e creators_destructors sostituiscono gli operatori di creazione e annichilazione nell'Hamiltoniano fermionico con operatori di Pauli; creators_destructors utilizza la mappatura di 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_listInfine, build_hamiltonian utilizza le funzioni cholesky, identity e creators_destructors per creare l'hamiltoniana finale adatta all'esecuzione su un computer quantistico.
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()Infine, utilizziamo build_hamiltonian per costruire l'hamiltoniana del qubit a partire dagli operatori di Pauli, utilizzando la trasformazione di Jordan-Wigner. Questo ci dà anche la precisione della decomposizione di Cholesky che abbiamo utilizzato.
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])
Questo quaderno di molecole di esempio mostra l'impostazione e gli hamiltoniani di diverse molecole di varia complessità; con qualche piccola modifica, questo dovrebbe consentire di esaminare la maggior parte delle piccole molecole.
Vediamo brevemente due punti importanti da considerare quando si costruiscono gli operatori fermionici per una molecola. Al variare del tipo di molecola, cambierà la simmetria. Allo stesso modo, cambierà il numero di orbitali con diverse simmetrie, come quello a simmetria cilindrica " A1 ". Questi cambiamenti sono evidenti anche con la semplice estensione a LiH,, come si vede qui:
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)Vale anche la pena di notare che si può perdere rapidamente l'intuizione per l'hamiltoniana finale risultante. L'hamiltoniana di LiH (utilizzando il mappatore di Jordan-Wigner) è già composta da 276 termini.
len(build_hamiltonian(ecore, h1e, h2e))Output:
accuracy of Cholesky decomposition 1.1102230246251565e-16
276
In caso di dubbi sulle simmetrie, si possono anche generare alcune informazioni sulla simmetria della molecola impostando symmetry = True e 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>
Tra le altre informazioni utili, restituisce sia point group symmetry = Coov sia il numero di orbitali in ogni rappresentazione irriducibile.
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
Questo non indica necessariamente quanti orbitali si vogliono includere nello spazio attivo, ma aiuta a vedere quali orbitali sono presenti e le loro simmetrie.
Spesso è utile specificare la simmetria e gli orbitali, ma è anche possibile specificare il numero di orbitali che si desidera includere. Consideriamo il caso dell'etene, qui di seguito. Utilizzando verbose = 4, possiamo stampare le simmetrie dei vari orbitali:
# 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>
Otteniamo:
numero di orbitali dell'irrep Ag = 4
num. orbitali dell'irrep B2g = 2
num. orbitali dell'irrep B3g = 1
num. orbitali dell'irrep B1u = 4
num. orbitali dell'irrep B2u = 1
num. orbitali dell'irrep B3u = 2
Ma piuttosto che specificare tutti gli orbitali per simmetria, possiamo semplicemente scrivere:
active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)In questo approccio, prendiamo diversi orbitali vicini al livello di riempimento (valenza e non occupati). Qui sono stati selezionati 5 orbitali da includere nello spazio attivo (dal 6° al 10°).
print(
mol.nelectron // 2 - 2,
mol.nelectron // 2 + 2,
)Output:
6 10
- Software di terze parti
Esistono diversi pacchetti software sviluppati per la chimica quantistica, alcuni dei quali offrono più mappatori e strumenti per limitare gli spazi attivi. I passi descritti sopra sono generali e si applicano anche al software di terze parti. Ma questo altro software potrebbe restituire gli hamiltoniani in un formato non accettato da Qiskit. Ad esempio, alcuni software restituiscono hamiltoniani della forma:
H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3]
Si noti in particolare che le porte sono numerate e che gli operatori di identità non sono mostrati. Questo è in contrasto con gli hamiltoniani utilizzati in Qiskit, che scrivono il termine [Z2 Z3] come ZZII (i qubit 0 e 1 sono agiti dall'operatore identità, i qubit 2 e 3 dall'operatore Z, ordinati con il qubit 0 più a destra).
Per adattarsi a qualsiasi flusso di lavoro esistente, il blocco di codice sottostante converte da una sintassi all'altra. La funzione convert_openfermion_to_qiskit prende come argomenti un'hamiltoniana generata in OpenFermion o Tangelo (e già mappata sugli operatori di Pauli usando qualsiasi mappatore disponibile) e il numero di qubit necessari per la molecola.
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)Inoltre, questo notebook Python contiene un codice di esempio completo per la migrazione di hamiltoniani da altri flussi di lavoro software a Qiskit, compresa la conversione di cui sopra.
Ora dovreste disporre di un arsenale di strumenti per ottenere l'hamiltoniana necessaria per eseguire calcoli di chimica quantistica sui computer quantistici IBM®.