Hamiltoniani per la chimica quantistica

Cominciamo con una breve panoramica del ruolo che gli Hamiltoniani svolgono nel VQE.

Panoramica dell'Hamiltoniano nel VQE

La dott.ssa Victoria Lipinska ci guida attraverso gli Hamiltoniani e spiega come mapparli per l'uso nel calcolo quantistico.

Riferimenti

I seguenti articoli sono citati nel video precedente.

Preparare gli Hamiltoniani per la chimica quantistica

Un ottimo primo passo nell'applicare il calcolo quantistico a un problema di chimica è definire un Hamiltoniano per il sistema di interesse. In questa sede limiteremo la discussione agli Hamiltoniani di chimica quantistica, poiché tali Hamiltoniani richiedono una mappatura specifica per i sistemi di fermioni identici.

Se lavori nel campo della chimica quantistica, probabilmente hai già il tuo software preferito per la modellazione delle molecole, in grado di generare un Hamiltoniano che descriva il sistema di interesse. Qui utilizzeremo codice basato esclusivamente su PySCF, numpy e Qiskit. Il processo di preparazione dell'Hamiltoniano si applica tuttavia anche a soluzioni preconfezionate. L'unica differenza tra questo approccio e altri software consiste in lievi differenze sintattiche; alcune di queste vengono affrontate nella sottosezione "Software di terze parti" per facilitare l'integrazione dei flussi di lavoro esistenti.

La generazione di un Hamiltoniano di chimica quantistica per l'uso su QPU IBM Quantum® prevede i seguenti passaggi:

Definire la molecola (geometria, spin, spazio attivo e così via)
Generare l'Hamiltoniano fermionico (operatori di creazione e distruzione)
Mappare dall'Hamiltoniano fermionico a un operatore bosonico (in questo contesto, usando operatori di Pauli)
In caso di software di terze parti: gestire eventuali incompatibilità sintattiche tra il software generatore e Qiskit

L'Hamiltoniano fermionico è scritto in termini di operatori fermionici e, in particolare, tiene conto del fatto che gli elettroni sono fermioni indistinguibili. Ciò significa che obbediscono a una statistica completamente diversa da quella di qubit distinguibili e bosonici. Da qui la necessità della mappatura.

Chi ha già familiarità con questi processi può probabilmente saltare questa sezione. Obiettivo:

L'obiettivo finale è ottenere un Hamiltoniano della forma:

# 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')]

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)

SparsePauliOp(['XX', 'YY', 'ZZ'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j])

Iniziamo importando alcuni pacchetti:

import numpy as np
from pyscf import ao2mo, gto, mcscf, scf

Definire la molecola

Qui specifichiamo gli attributi della molecola di interesse. In questo esempio abbiamo scelto l'idrogeno biatomico (poiché gli Hamiltoniani risultanti sono abbastanza brevi da essere visualizzati). Il Python-based Simulations of Chemistry Framework (PySCF) dispone di un'ampia raccolta di moduli di struttura elettronica che possono essere utilizzati per, tra le altre cose, generare Hamiltoniani molecolari adatti al calcolo quantistico. La guida PySCF Quickstart è un'ottima risorsa per una descrizione completa di tutte le variabili e funzionalità. Forniremo solo una panoramica molto sintetica, dal momento che questo argomento sarà già familiare a molti di voi. Per approfondire, visita PySCF. In breve:

distance può essere usato per molecole biatomiche, oppure si possono specificare direttamente le coordinate cartesiane di ciascun atomo. Le distanze sono espresse in Angstrom.

gto genera orbitali di tipo gaussiano.

basis si riferisce alle funzioni utilizzate per modellare gli orbitali molecolari. Qui 'sto-6g' è una base minimale comune, denominata così per l'adattamento degli orbitali di tipo Slater usando 6 gaussiane primitive.

spin un valore intero che indica il numero di elettroni spaiati (uguale a $2S$ ). Nota che alcuni software usano invece la molteplicità ( $2S+1$ ).

charge la carica della molecola.

symmetry - il gruppo di simmetria puntuale della molecola, specificato tramite una stringa oppure rilevato automaticamente impostando "symmetry = True". Qui "Dooh" è il gruppo di simmetria appropriato per le molecole biatomiche omonucleari.

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>

Tieni presente che si può descrivere l'energia totale (che include l'energia di repulsione nucleare oltre a quella elettronica), l'energia orbitale elettronica totale, o l'energia di un sottoinsieme di orbitali elettronici (con il sottoinsieme complementare congelato). Nel caso specifico di $\text{H}_2$ , nota le diverse energie riportate di seguito, e osserva che l'energia totale meno l'energia di repulsione nucleare restituisce effettivamente 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(),
)

0.7199689944489797 -1.8455976628764188 -1.125628668427439 -1.8455976628764188

active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)

Generare l'Hamiltoniano fermionico

scf si riferisce a un'ampia gamma di metodi a campo autoconsistente.

rhf come in mf = scf.RHF(mol), è un solver che utilizza il calcolo Hartree-Fock ristretto. Il kernel di questo (E, di seguito) è l'energia totale, compresa la repulsione nucleare e gli orbitali molecolari.

mcscf è un pacchetto per campi autoconsistenti multiconfigurazione.

ao2mo è una trasformazione dagli orbitali atomici agli 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 Hamiltoniano, e questo viene spesso scomposto in: energia del nucleo elettronico (ecore, non coinvolta nella minimizzazione), operatori a un elettrone (h1e) ed energie a due elettroni (h2e). Questi vengono estratti esplicitamente di seguito 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 distruzione), applicabili a sistemi di fermioni (indistinguibili) e, di conseguenza, soggetti all'antisimmetria sotto scambio. Questo comporta una statistica diversa da quella applicabile a un sistema distinguibile o bosonico. Per eseguire calcoli su QPU IBM Quantum, è necessario un operatore bosonico che descriva l'energia. Il risultato di tale mappatura è convenzionalmente scritto in termini di operatori di Pauli, poiché sono sia hermitiani che unitari. Esistono diverse mappature utilizzabili. Una delle più semplici è la trasformazione Jordan-Wigner.

Mappare l'Hamiltoniano

È opportuno notare che esistono molti strumenti disponibili per mappare un Hamiltoniano chimico in uno adatto all'esecuzione su un computer quantistico. Qui implementiamo la mappatura Jordan-Wigner direttamente usando solo PySCF, numpy e Qiskit. Di seguito commentiamo alcune considerazioni sintattiche per altre soluzioni. La funzione Cholesky ci aiuta a ottenere una decomposizione a basso rango dei termini a due elettroni nell'Hamiltoniano.

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

Le funzioni identity e creators_destructors sostituiscono gli operatori di creazione e distruzione nell'Hamiltoniano fermionico con operatori di Pauli; creators_destructors utilizza la mappatura 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

Infine, build_hamiltonian utilizza le funzioni cholesky, identity e creators_destructors per costruire l'Hamiltoniano finale adatto 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, usiamo build_hamiltonian per costruire il nostro Hamiltoniano di qubit dagli operatori di Pauli tramite la trasformazione Jordan-Wigner. Questo ci fornisce anche l'accuratezza della decomposizione di Cholesky utilizzata.

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])

Questo notebook di esempio sulle molecole mostra la configurazione e gli Hamiltoniani per diverse molecole di complessità variabile; con poche modifiche, dovrebbe consentirti di esaminare la maggior parte delle piccole molecole.

Notiamo brevemente due punti importanti da considerare quando si costruiscono gli operatori fermionici per una molecola. Al variare del tipo di molecola, cambia anche la simmetria. Analogamente, cambia il numero di orbitali con varie simmetrie, come la 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 notare che si può perdere rapidamente l'intuizione sull'Hamiltoniano finale risultante. L'Hamiltoniano per LiH (usando il mapper Jordan-Wigner) è già composto da 276 termini.

len(build_hamiltonian(ecore, h1e, h2e))

accuracy of Cholesky decomposition  1.1102230246251565e-16

In caso di dubbi riguardo alle simmetrie, è possibile 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,
)

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, questo restituisce sia point group symmetry = Coov sia il numero di orbitali in ciascuna 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 ti dice necessariamente quanti orbitali vuoi includere nel tuo spazio attivo, ma ti aiuta a vedere quali orbitali sono presenti e quali sono le loro simmetrie.

Specificare la simmetria e gli orbitali è spesso utile, ma puoi anche specificare il numero di orbitali da includere. Considera il caso dell'etene, di seguito. Usando 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,
)

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:

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

Ma invece di specificare tutti gli orbitali per simmetria, possiamo semplicemente scrivere:

active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)

Con questo approccio, selezioniamo alcuni orbitali vicini al livello di riempimento (di 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,
)

6 10

Software di terze parti

Esistono diversi pacchetti software sviluppati per la chimica quantistica, alcuni dei quali offrono più mapper e strumenti per limitare gli spazi attivi. I passaggi descritti in precedenza sono generali e si applicano anche al software di terze parti. Tuttavia, tale software potrebbe restituire Hamiltoniani in un formato non accettato da Qiskit. Ad esempio, alcuni software restituiscono Hamiltoniani nella forma:

H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3] Nota in particolare che i gate sono numerati e gli operatori identità non vengono mostrati. Questo contrasta con gli Hamiltoniani utilizzati in Qiskit, in cui il termine [Z2 Z3] sarebbe scritto come ZZII (i qubit 0 e 1 vengono agiti dall'operatore identità, i qubit 2 e 3 vengono agiti dall'operatore Z, ordinati con il qubit 0 più a destra).

Per adattarsi agli eventuali flussi di lavoro esistenti, il blocco di codice seguente converte da una sintassi all'altra. La funzione convert_openfermion_to_qiskit prende come argomenti un Hamiltoniano generato in OpenFermion o Tangelo (già mappato sugli operatori di Pauli con qualsiasi mapper 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 codice di esempio completo per migrare Hamiltoniani da altri flussi di lavoro software a Qiskit, inclusa la conversione sopra descritta.

Dovresti ora disporre di un arsenale di strumenti per ottenere l'Hamiltoniano necessario per eseguire calcoli di chimica quantistica su computer quantistici IBM®.

Panoramica dell'Hamiltoniano nel VQE​

Riferimenti​

Preparare gli Hamiltoniani per la chimica quantistica​

Panoramica dell'Hamiltoniano nel VQE

Riferimenti

Preparare gli Hamiltoniani per la chimica quantistica