Hamiltoniany w chemii kwantowej

Zacznijmy od krótkiego przeglądu roli, jaką hamiltoniany odgrywają w VQE.

Hamiltonian w VQE - przegląd

Dr Victoria Lipinska przeprowadzi nas przez hamiltoniany i sposób ich mapowania do wykorzystania w obliczeniach kwantowych.

Bibliografia

W powyższym materiale wideo przywoływane są następujące artykuły.

Przygotowanie hamiltonianów do chemii kwantowej

Dobrym pierwszym krokiem przy zastosowaniu obliczeń kwantowych do problemu chemicznego jest zdefiniowanie hamiltonianu dla interesującego nas układu. Tutaj ograniczymy dyskusję do hamiltonianów chemii kwantowej, ponieważ te hamiltoniany wymagają pewnego mapowania specyficznego dla układów identycznych fermionów.

Jako osoba pracująca w chemii kwantowej, prawdopodobnie masz już swoje ulubione oprogramowanie do modelowania cząsteczek, które może wygenerować hamiltonian opisujący interesujący cię układ. Tutaj użyjemy kodu zbudowanego wyłącznie na PySCF, numpy i Qiskit. Jednak proces przygotowania hamiltonianu przenosi się również na gotowe rozwiązania. Jedyna różnica między tym podejściem a innymi programami polega na drobnych różnicach w składni; niektóre z nich są omawiane w podrozdziale „Oprogramowanie zewnętrzne”, aby ułatwić integrację istniejących przepływów pracy.

Generowanie hamiltonianu chemii kwantowej do użycia na QPU IBM Quantum® obejmuje następujące kroki:

Zdefiniuj swoją cząsteczkę (geometria, spin, przestrzeń aktywna itd.)
Wygeneruj fermionowy hamiltonian (operatory kreacji i anihilacji)
Zmapuj fermionowy hamiltonian na operator bozonowy (w tym kontekście, za pomocą Pauli operators)
Jeśli używasz oprogramowania zewnętrznego: obsłuż wszelkie niezgodności składniowe między oprogramowaniem generującym a Qiskit

Fermionowy hamiltonian zapisuje się przy użyciu operatorów fermionowych, a w szczególności uwzględnia fakt, że elektrony są nierozróżnialnymi fermionami. Oznacza to, że podlegają zupełnie innej statystyce niż rozróżnialne, bozonowe kubity. Stąd proces mapowania.

Ci z was, którzy są już zaznajomieni z tymi procesami, mogą prawdopodobnie pominąć tę sekcję. Cel:

Końcowym celem jest uzyskanie hamiltonianu w postaci:

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

Lub

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

Zaczniemy od zaimportowania kilku pakietów:

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

Zdefiniuj swoją cząsteczkę

Tutaj określimy atrybuty interesującej nas cząsteczki. W tym przykładzie wybraliśmy dwuatomowy wodór (ponieważ otrzymane hamiltoniany są wystarczająco krótkie, aby je wyświetlić). Python-based Simulations of Chemistry Framework (PySCF) posiada szeroki zbiór modułów struktury elektronowej, których można używać między innymi do generowania hamiltonianów molekularnych odpowiednich do obliczeń kwantowych. Przewodnik PySCF Quickstart jest doskonałym źródłem pełnego opisu wszystkich zmiennych i funkcji. Przedstawimy jedynie najbardziej pobieżny przegląd, ponieważ wielu z was będzie to już znane. Aby lepiej to zrozumieć, odwiedź PySCF. W skrócie:

distance może być używane dla cząsteczek dwuatomowych lub po prostu określ współrzędne kartezjańskie dla każdego atomu. Odległości podawane są w jednostkach angstrema.

gto generuje orbitale typu gaussowskiego.

basis odnosi się do funkcji używanych do modelowania orbitali molekularnych. Tutaj 'sto-6g' jest popularną minimalną bazą, nazwaną tak ze względu na dopasowanie Slater-Type Orbitals przy użyciu 6 prymitywnych orbitali gaussowskich.

spin wartość całkowitoliczbowa wskazująca liczbę niesparowanych elektronów (równa $2S$ ). Zauważ, że niektóre programy używają multipletowości zamiast tego ( $2S+1$ ).

charge ładunek cząsteczki.

symmetry - grupa symetrii punktowej cząsteczki, określona jako łańcuch znaków lub wykrywana automatycznie przez ustawienie "symmetry = True". Tutaj "Dooh" jest odpowiednią grupą symetrii dla cząsteczek dwuatomowych z dwoma atomami tego samego typu.

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>

Należy pamiętać, że można opisywać energię całkowitą (która obejmuje zarówno energię odpychania jądrowego, jak i elektronową), całkowitą energię orbitali elektronowych lub energię pewnego podzbioru orbitali elektronowych (z zamrożonym podzbiorem komplementarnym). W szczególnym przypadku $\text{H}_2$ zwróć uwagę na różne energie poniżej i zauważ, że energia całkowita pomniejszona o energię odpychania jądrowego daje w istocie energię elektronową:

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)

Wygeneruj fermionowy Hamiltonian

scf odnosi się do szerokiego zakresu metod pola samouzgodnionego (self-consistent field).

rhf jak w mf = scf.RHF(mol) to solver wykorzystujący obliczenia metodą Restricted Hartree Fock. Jądro tego obliczenia (E, poniżej) to energia całkowita, obejmująca odpychanie jądrowe i molecular orbitals.

mcscf to pakiet dla wielokonfiguracyjnych metod pola samouzgodnionego.

ao2mo to transformacja z orbitali atomowych do molecular orbitals.

Używamy również następujących zmiennych:

ncas: liczba orbitali w pełnej active space

nelecas: liczba elektronów w pełnej active space

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]

Chcemy uzyskać Hamiltonian, a ten jest często rozdzielany na energię rdzenia elektronowego (ecore, nieuczestniczącą w minimalizacji), operatory jednoelektronowe (h1e) oraz energie dwuelektronowe (h2e). Są one jawnie ekstrahowane poniżej w dwóch ostatnich liniach.

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

Te Hamiltoniany są obecnie operatorami fermionowymi (creation i annihilation), mającymi zastosowanie do układów (nierozróżnialnych) fermions, i odpowiednio podlegają antysymetrii względem wymiany. Powoduje to inną statystykę niż ta, która miałaby zastosowanie do układu rozróżnialnego lub bozonowego. Aby uruchomić obliczenia na QPU IBM Quantum, potrzebujemy operatora bozonowego opisującego energię. Wynik takiego mapowania jest konwencjonalnie zapisywany w kategoriach operatorów Pauli, ponieważ są one zarówno hermitowskie, jak i unitarne. Istnieje kilka mapowań, których można użyć. Jednym z najprostszych jest transformacja Jordan-Wigner.

Mapowanie Hamiltonianu

Należy zauważyć, że istnieje wiele narzędzi dostępnych do mapowania chemicznego Hamiltonianu na taki, który nadaje się do uruchamiania na komputerze kwantowym. Tutaj implementujemy mapowanie Jordan-Wigner bezpośrednio, używając tylko PySCF, numpy i Qiskit. Poniżej komentujemy rozważania składniowe dla innych rozwiązań. Funkcja Cholesky pomaga nam uzyskać niskorangową dekompozycję składników dwuelektronowych w Hamiltonianie.

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

Funkcje identity i creators_destructors zastępują creation i annihilation operators w fermionowym Hamiltonianie operatorami Pauli; creators_destructors używa mapowania 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

Na koniec, build_hamiltonian używa funkcji cholesky, identity i creators_destructors, aby utworzyć ostateczny Hamiltonian nadający się do uruchomienia na komputerze kwantowym.

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

Na koniec używamy build_hamiltonian do skonstruowania naszego kubitowego Hamiltonianu z operatorów Pauli, stosując transformację Jordan-Wigner. To daje nam również dokładność dekompozycji Cholesky, której użyliśmy.

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

Ten przykładowy notatnik z molekułami przedstawia konfigurację i Hamiltonians dla kilku molekuł o różnej złożoności; po niewielkich modyfikacjach powinno to umożliwić badanie większości małych molekuł.

Zauważmy krótko dwie istotne kwestie, które należy brać pod uwagę przy konstruowaniu operatorów fermionowych dla molekuły. Gdy zmienia się typ molekuły, zmienia się również symetria. Podobnie zmienia się liczba orbitali o różnych symetriach, takich jak cylindrycznie symetryczna "A1". Zmiany te są widoczne już przy prostym rozszerzeniu do LiH, co pokazano poniżej:

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)

Warto również zauważyć, że można szybko stracić intuicję co do ostatecznego Hamiltonianu. Hamiltonian dla LiH (przy użyciu qubit mappera Jordan-Wigner) składa się już z 276 wyrazów.

len(build_hamiltonian(ecore, h1e, h2e))

accuracy of Cholesky decomposition  1.1102230246251565e-16

W razie wątpliwości co do symetrii można również wygenerować informacje o symetrii molekuły, ustawiając symmetry = True oraz 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>

Wśród innych przydatnych informacji zwraca to zarówno point group symmetry = Coov, jak i liczbę orbitali w każdej reprezentacji nieprzywiedlnej.

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

Niekoniecznie informuje to, ile orbitali chcesz uwzględnić w swojej active space, ale pomaga zobaczyć, jakie orbitale są obecne i jakie mają symetrie.

Określanie symetrii i orbitali jest często pomocne, ale można również określić liczbę orbitali, które chcesz uwzględnić. Rozważmy poniższy przypadek etenu. Korzystając z verbose = 4, możemy wypisać symetrie różnych 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>

Otrzymujemy:

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

Zamiast jednak określać wszystkie orbitale na podstawie symetrii, możemy po prostu napisać:

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

W tym podejściu bierzemy kilka orbitali w pobliżu poziomu zapełnienia (walencyjne i nieobsadzone). Tutaj do uwzględnienia w przestrzeni aktywnej wybrano 5 orbitali (od 6. do 10.).

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

6 10

Oprogramowanie firm trzecich

Istnieje kilka pakietów oprogramowania opracowanych na potrzeby chemii kwantowej, z których część oferuje wiele mapperów oraz narzędzia do ograniczania przestrzeni aktywnych. Opisane powyżej kroki mają charakter ogólny i stosują się również do oprogramowania firm trzecich. Jednak to inne oprogramowanie może zwracać Hamiltoniany w formacie nieakceptowanym przez Qiskit. Na przykład niektóre oprogramowanie zwraca Hamiltoniany w postaci:

H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3] Zwróć szczególną uwagę na to, że bramki są numerowane, a operatory tożsamości nie są pokazywane. Kontrastuje to z Hamiltonianami używanymi w Qiskit, w których człon [Z2 Z3] byłby zapisany jako ZZII (kubity 0 i 1 są poddawane działaniu operatora tożsamości, a kubity 2 i 3 operatora Z, przy czym kubit 0 znajduje się najbardziej po prawej stronie).

Aby ułatwić integrację z istniejącymi przepływami pracy, poniższy blok kodu dokonuje konwersji z jednej składni na drugą. Funkcja convert_openfermion_to_qiskit przyjmuje jako argumenty Hamiltonian wygenerowany w OpenFermion lub Tangelo (i już odwzorowany na operatory Pauli za pomocą dowolnego dostępnego mappera) oraz liczbę kubitów potrzebnych dla danej cząsteczki.

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)

Ponadto ten notatnik Pythona zawiera kompletny przykładowy kod do migracji Hamiltonianów z innych przepływów pracy oprogramowania do Qiskit, w tym powyższą konwersję.

Hamiltonian w VQE - przegląd​

Bibliografia​

Przygotowanie hamiltonianów do chemii kwantowej​

Hamiltonian w VQE - przegląd

Bibliografia

Przygotowanie hamiltonianów do chemii kwantowej