Hamiltonieni pentru chimia cuantică

Să începem cu o scurtă prezentare a rolului pe care Hamiltonienii îl joacă în VQE.

Prezentare generală a Hamiltonianului în VQE

Dr. Victoria Lipinska ne explică Hamiltonienii și cum să îi mapăm pentru a-i utiliza în calcul cuantic.

Referințe

Articolele de mai jos sunt menționate în videoclipul de mai sus.

Pregătirea Hamiltonienilor pentru chimia cuantică

Un prim pas bun în aplicarea calculului cuantic la o problemă de chimie este definirea unui Hamiltonian pentru sistemul de interes. Aici vom restrânge discuția la Hamiltonienii din chimia cuantică, deoarece acești Hamiltonieni necesită o mapare specifică sistemelor de fermioni identici.

Ca cineva care lucrează în chimia cuantică, probabil că ai deja software-ul tău preferat pentru modelarea moleculelor, care poate genera un Hamiltonian ce descrie sistemul de interes. Aici vom folosi cod construit exclusiv pe PySCF, numpy și Qiskit. Dar procesul de pregătire a Hamiltonianului se transferă și la soluțiile prepackaged. Singura diferență dintre această abordare și alt software vor fi diferențe minore de sintaxă; câteva dintre acestea sunt abordate în subsecțiunea „Software terț" pentru a facilita integrarea fluxurilor de lucru existente.

Generarea unui Hamiltonian de chimie cuantică pentru utilizare pe QPU-urile IBM Quantum® implică următorii pași:

Definește molecula ta (geometrie, spin, spațiu activ și altele)
Generează Hamiltonianul fermionic (operatori de creare și anihilare)
Mapează de la Hamiltonianul fermionic la un operator bosonic (în acest context, folosind operatori Pauli)
Dacă folosești software terț: gestionează orice nepotriviri de sintaxă între software-ul generator și Qiskit

Hamiltonianul fermionic este scris în termeni de operatori fermionici și, în particular, ține cont de faptul că electronii sunt fermioni indisolubili. Asta înseamnă că respectă statistici complet diferite față de qubiții distingebili, bosonici. De aceea este nevoie de procesul de mapare.

Cei dintre voi care sunt deja familiarizați cu aceste procese pot sări probabil peste această secțiune. Obiectiv:

Obiectivul final este obținerea unui Hamiltonian de 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')]

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

Vom începe prin importarea unor pachete:

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

Definește molecula ta

Aici vom specifica atributele moleculei de interes. În acest exemplu, am ales hidrogenul diatomic (deoarece Hamiltonienii rezultanți sunt suficient de scurți pentru a fi afișați). Python-based Simulations of Chemistry Framework (PySCF) dispune de o colecție vastă de module de structură electronică ce pot fi utilizate, printre altele, pentru a genera Hamiltonieni moleculari potriviți pentru calculul cuantic. Ghidul PySCF Quickstart este o resursă excelentă pentru o descriere completă a tuturor variabilelor și funcționalităților. Vom oferi doar cea mai sumară prezentare, deoarece acest lucru va fi deja familiar multora dintre voi. Pentru a înțelege mai bine, vizitează PySCF. Pe scurt:

distance poate fi utilizat pentru molecule diatomice sau pur și simplu specifică coordonatele carteziene pentru fiecare atom. Distanțele sunt exprimate în unități de Angstrom.

gto generează orbitali de tip gaussian.

basis se referă la funcțiile utilizate pentru modelarea orbitalilor moleculari. Aici sto-6g este o bază minimală comună, numită după ajustarea Orbitalilor de Tip Slater folosind 6 orbitali gaussieni primitivi.

spin o valoare întreagă care indică numărul de electroni nepereche (egal cu $2S$ ). Reține că unele programe folosesc multiplicitatea în schimb ( $2S+1$ ).

charge sarcina moleculei.

symmetry - grupul de simetrie punctuală al moleculei, specificat fie printr-un șir de caractere, fie detectat automat prin setarea „symmetry = True". Aici „Dooh" este grupul de simetrie potrivit pentru moleculele diatomice cu două specii de atomi identice.

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>

Reține că se poate descrie energia totală (care include energia de repulsie nucleară, precum și cea electronică), energia orbitală electronică totală sau energia unui subset de orbitali electronici (cu subsetul complementar înghețat). În cazul specific al $\text{H}_2$ , observă diferitele energii de mai jos și că energia totală minus energia de repulsie nucleară dă de fapt energia electronică:

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)

Generarea Hamiltonianului fermionic

scf se referă la o gamă largă de metode de câmp autoconsistent.

rhf ca în mf = scf.RHF(mol) — mf este un solver care folosește calculul Hartree-Fock restricționat. Nucleul acestuia (E, mai jos) este energia totală, incluzând repulsia nucleară și orbitalii moleculari.

mcscf este un pachet de câmpuri autoconsistente cu configurații multiple.

ao2mo este o transformare din orbitali atomici în orbitali moleculari.

Folosim și următoarele variabile:

ncas: numărul de orbitali din spațiul activ complet

nelecas: numărul de electroni din spațiul activ 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]

Dorim un Hamiltonian, care este adesea separat în energia unui miez electronic (ecore, neimplicat în minimizare), operatori cu un singur electron (h1e) și energii cu doi electroni (h2e). Acestea sunt extrase explicit mai jos în ultimele două linii.

h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)

Acești Hamiltonieni sunt în prezent operatori fermionici (de creare și anihilare), aplicabili sistemelor de fermioni (indistinguibili) și, corespunzător, supuși antisimetriei la schimb. Aceasta duce la statistici diferite față de cele care s-ar aplica unui sistem distingibil sau bosonic. Pentru a rula calcule pe QPU-urile IBM Quantum, avem nevoie de un operator bosonic care să descrie energia. Rezultatul unui astfel de mapaj este scris în mod convențional în termeni de operatori Pauli, deoarece aceștia sunt atât hermitici, cât și unitari. Există mai multe mapaje pe care le poți folosi. Unul dintre cele mai simple este transformarea Jordan-Wigner.

Maparea Hamiltonianului

Trebuie menționat că există multe instrumente disponibile pentru maparea unui Hamiltonian chimic la unul potrivit pentru rularea pe un calculator cuantic. Aici, implementăm mapajul Jordan-Wigner direct, folosind doar PySCF, numpy și Qiskit. Comentăm mai jos considerațiile de sintaxă pentru alte soluții. Funcția Cholesky ne ajută să obținem o descompunere de rang scăzut a termenilor cu doi electroni din Hamiltonian.

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

Funcțiile identity și creators_destructors înlocuiesc operatorii de creare și anihilare din Hamiltonianul fermionic cu operatori Pauli; creators_destructors folosește mapajul 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

În final, build_hamiltonian folosește funcțiile cholesky, identity și creators_destructors pentru a crea Hamiltonianul final, adecvat rulării pe un calculator cuantic.

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

În final, folosim build_hamiltonian pentru a construi Hamiltonianul nostru de qubiți din operatori Pauli, prin transformarea Jordan-Wigner. Aceasta ne oferă și acuratețea descompunerii Cholesky utilizate.

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

Acest notebook cu molecule exemplu prezintă configurarea și Hamiltonienii pentru mai multe molecule de complexitate variată; cu mici modificări, acesta ar trebui să îți permită să examinezi majoritatea moleculelor mici.

Să notăm pe scurt două aspecte importante de luat în considerare la construirea operatorilor fermionici pentru o moleculă. Pe măsură ce tipul moleculei se schimbă, se va schimba și simetria. În mod similar, se va schimba și numărul de orbitali cu diverse simetrii, cum ar fi „A1" cu simetrie cilindrică. Aceste schimbări sunt evidente chiar și cu simpla extindere la LiH, după cum se vede aici:

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)

De asemenea, merită menționat că intuiția pentru Hamiltonianul rezultat final se pierde rapid. Hamiltonianul pentru LiH (folosind mapper-ul Jordan-Wigner) constă deja din 276 de termeni.

len(build_hamiltonian(ecore, h1e, h2e))

accuracy of Cholesky decomposition  1.1102230246251565e-16

Când există dubii cu privire la simetrii, se pot genera și informații despre simetria moleculei setând symmetry = True și 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>

Printre alte informații utile, aceasta returnează atât point group symmetry = Coov, cât și numărul de orbitali din fiecare reprezentare ireductibilă.

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

Acest lucru nu îți spune neapărat câți orbitali vrei să incluzi în spațiul tău activ, dar te ajută să vezi ce orbitali sunt prezenți și simetriile lor.

Specificarea simetriei și a orbitalilor este adesea utilă, dar poți specifica și numărul de orbitali pe care vrei să îi incluzi. Consideră cazul etenei, de mai jos. Folosind verbose = 4, putem afișa simetriile diferiților 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>

Obținem:

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

Dar în loc să specificăm toți orbitalii după simetrie, putem scrie pur și simplu:

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

Prin această abordare, selectăm câțiva orbitali din apropierea nivelului de umplere (valenți și neocupați). Aici, 5 orbitali au fost selectați pentru a fi incluși în spațiul activ (al 6-lea până la al 10-lea).

print(
    mol.nelectron // 2 - 2,
    mol.nelectron // 2 + 2,
)

6 10

Software terț

Există mai multe pachete software dezvoltate pentru chimie cuantică, unele oferind mai mulți mapperi și instrumente pentru restricționarea spațiilor active. Pașii descriși mai sus sunt generali și se aplică și software-ului terț. Dar acest alt software poate returna Hamiltonieni într-un format care nu este acceptat de Qiskit. De exemplu, unele software-uri returnează Hamiltonieni de forma:

H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3] Observă în special că porțile sunt numerotate, iar operatorii identitate nu sunt afișați. Aceasta este în contrast cu Hamiltonienii folosiți în Qiskit, care ar scrie termenul [Z2 Z3] ca ZZII (qubiții 0 și 1 acționați de operatorul identitate, qubiții 2 și 3 acționați de operatorul Z, ordonați cu qubit-ul 0 cel mai la dreapta).

Pentru a acomoda orice fluxuri de lucru existente, blocul de cod de mai jos convertește dintr-o sintaxă în cealaltă. Funcția convert_openfermion_to_qiskit primește ca argumente un Hamiltonian generat în OpenFermion sau Tangelo (deja mapat pe operatori Pauli folosind orice mapper disponibil) și numărul de qubiți necesari pentru moleculă.

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)

În plus, acest notebook Python conține cod de exemplu complet pentru migrarea Hamiltonienilor din alte fluxuri de lucru software în Qiskit, inclusiv conversia de mai sus.

Prezentare generală a Hamiltonianului în VQE​

Referințe​

Pregătirea Hamiltonienilor pentru chimia cuantică​

Prezentare generală a Hamiltonianului în VQE

Referințe

Pregătirea Hamiltonienilor pentru chimia cuantică