Hamiltonianos para la química cuántica

Comencemos con un breve repaso del papel que desempeñan los hamiltonianos en la VQE.

El hamiltoniano en VQE Descripción general

La Dra. Victoria Lipinska nos habla de los hamiltonianos y de cómo mapearlos para su uso en computación cuántica.

Referencias

En el vídeo se hace referencia a los siguientes artículos.

Algoritmos cuánticos para simulaciones fermiónicas, Ortiz, et al.
Simulación de la química mediante ordenadores cuánticos, Kassal et al.
A Comparison of the Bravyi-Kitaev and Jordan-Wigner Transformations for the Quantum Simulation of Quantum Chemistry, Tranter, et al.
Quantum Chemistry in the Age of Quantum Computing, Cao, et al.
Química cuántica computacional, McArdle, et al.
La transformación Bravyi-Kitaev para el cálculo cuántico de la estructura electrónica, Seeley, et al., McArdle, et al.

Preparando a los habitantes de Hamilton para la química cuántica

Un buen primer paso para aplicar la computación cuántica a un problema químico es definir un Hamiltoniano para el sistema de interés. Aquí, restringiremos la discusión a los Hamiltonianos de la química cuántica, ya que esos Hamiltonianos requieren algún mapeo específico para sistemas de fermiones idénticos.

Como alguien que trabaja en química cuántica, probablemente ya tiene su software favorito para modelar moléculas, que puede generar un Hamiltoniano que describa su sistema de interés. Aquí, utilizaremos código construido únicamente en PySCF, numpy, y Qiskit. Pero el proceso de preparación hamiltoniana se traslada también a las soluciones preenvasadas. La única diferencia entre este enfoque y otro software serán pequeñas diferencias de sintaxis; algunas de ellas se abordan en la subsección "Software de terceros" para facilitar la integración de los flujos de trabajo existentes.

La generación de un Hamiltoniano de química cuántica para su uso en IBM Quantum® QPUs implica los siguientes pasos:

Defina su molécula (geometría, espín, espacio activo, etc.)
Generar el Hamiltoniano fermiónico (operadores de creación y aniquilación)
Mapa del Hamiltoniano fermiónico a un operador bosónico (en este contexto, utilizando operadores de Pauli)
Si utiliza software de terceros: Solucione cualquier desajuste de sintaxis entre el software generador y Qiskit

El hamiltoniano fermiónico se escribe en términos de operadores fermiónicos y, en particular, tiene en cuenta que los electrones son fermiones indistinguibles. Eso significa que obedecen a estadísticas completamente diferentes de las de los qubits bosónicos distinguibles. De ahí el proceso de mapeo.

Quienes ya estén familiarizados con estos procesos pueden saltarse esta sección.

Objetivo:

El objetivo final es obtener un Hamiltoniano de la forma:

H = [(1, "XX"), (1, "YY"), (1, "ZZ")]
print(H)

Output:

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

Output:

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

Empezaremos importando algunos paquetes:

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

Defina su molécula

Aquí especificaremos atributos de la molécula de interés. En este ejemplo, hemos elegido hidrógeno diatómico (porque los hamiltonianos resultantes son lo suficientemente cortos como para mostrarlos).

El marco Python -based Simulations of Chemistry Framework ( PySCF ) dispone de una amplia colección de módulos de estructura electrónica que pueden utilizarse, entre otras cosas, para generar hamiltonianos moleculares adecuados para la computación cuántica. La guía de inicio rápido PySCF es un recurso excelente para obtener una descripción completa de todas las variables y funcionalidades. Sólo daremos una visión muy somera, puesto que a muchos de ustedes ya les resultará familiar. Para entenderlos mejor, visite PySCF. Brevemente:

distancia puede utilizarse para moléculas diatómicas, o simplemente especificar coordenadas cartesianas para cada átomo. Las distancias están en unidades de Angstrom.

gto genera orbitales de tipo gaussiano.

base se refiere a las funciones utilizadas para modelar los orbitales moleculares. Aquí ' sto-6g ' es una base mínima común, llamada así para ajustar orbitales tipo Slater usando 6 orbitales Gaussianos primitivos.

spin un valor entero que indica el número de electrones no apareados (igual a $2S$ ). Tenga en cuenta que algunos programas utilizan la multiplicidad en su lugar ( $2S+1$ ).

carga la carga de la molécula.

simetría - el grupo de simetría puntual de la molécula, especificado con una cadena o detectado automáticamente estableciendo "simetría = True". Aquí "Dooh" es el grupo de simetría apropiado para moléculas diatómicas con dos de las mismas especies de átomos.

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>

Hay que tener en cuenta que se puede describir la energía total (que incluye la energía de repulsión nuclear además de la electrónica), la energía total de los orbitales electrónicos o la energía de algún subconjunto de orbitales electrónicos (con el subconjunto complementario congelado). En el caso concreto de $\text{H}_2$ , obsérvense las diferentes energías a continuación, y nótese que la energía total menos la energía de repulsión nuclear da, de hecho, la energía electrónica:

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)

Generar el Hamiltoniano fermiónico

scf hace referencia a una amplia gama de métodos de campo autoconsistente.

rhf como en mf = scf.RHF (mol) en mf es un solver que utiliza el cálculo Restricted Hartree Fock. El núcleo de esto (E, abajo) es la energía total, incluyendo la repulsión nuclear y los orbitales moleculares.

mcscf es un paquete de campos autoconsistentes multiconfiguración.

ao2mo es una transformación de orbitales atómicos a orbitales moleculares.

También utilizamos las siguientes variables:

ncas : número de orbitales en el espacio activo completo

nelecas : número de electrones en el espacio activo 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]

Queremos un Hamiltoniano, y esto a menudo se separa en energía de un núcleo electrónico (ecore, no involucrado en la minimización), operadores de un solo electrón ( h1e ), y energías de dos electrones ( h2e ). A continuación se extraen explícitamente en las dos últimas líneas.

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

Estos hamiltonianos son actualmente operadores fermiónicos (creación y aniquilación), aplicables a sistemas de fermiones (indistinguibles), y correspondientemente sujetos a antisimetría bajo intercambio. Esto da lugar a una estática diferente de la que se aplicaría a un sistema distinguible o bosónico. Para realizar cálculos en IBM Quantum QPUs, necesitamos un operador bosónico que describa la energía. El resultado de dicho mapeo se escribe convencionalmente en términos de operadores de Pauli, ya que ambos son hermitianos y unitarios. Se pueden utilizar varias correspondencias. Una de las más sencillas es la transformación de Jordan Wigner.

Trazado del Hamiltoniano

Hay que tener en cuenta que existen muchas herramientas para convertir un Hamiltoniano químico en uno adecuado para ser ejecutado en un ordenador cuántico. Aquí, implementamos el mapeo Jordan Wigner directamente usando sólo PySCF, numpy, y Qiskit. A continuación comentamos las consideraciones de sintaxis para otras soluciones.

La función Cholesky nos ayuda a obtener una descomposición de bajo rango de los términos de dos electrones en el 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

Las funciones identity y creators_destructors sustituyen los operadores de creación y aniquilación en el Hamiltoniano fermiónico por operadores de Pauli; creators_destructors utiliza el mapeo de 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

Por último, build_hamiltonian utiliza las funciones cholesky, identity, y creators_destructors para crear el Hamiltoniano final adecuado para su ejecución en un ordenador cuántico.

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

Por último, utilizamos build_hamiltonian para construir nuestro Hamiltoniano qubit a partir de operadores Pauli utilizando la transformación Jordan-Wigner. Esto también nos da la precisión de la descomposición Cholesky que hemos utilizado.

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

Este cuaderno de moléculas de ejemplo muestra la configuración y los Hamiltonianos para varias moléculas de complejidad variable; con una pequeña modificación, esto debería permitirle examinar la mayoría de las moléculas pequeñas.

Señalemos brevemente dos puntos importantes a tener en cuenta al construir los operadores fermiónicos de una molécula. Al cambiar el tipo de molécula, cambiará la simetría. Del mismo modo, cambiará el número de orbitales con distintas simetrías, como el de simetría cilíndrica " A1 ". Estos cambios son evidentes incluso con la simple ampliación a LiH,, como se ve aquí:

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)

También vale la pena señalar que uno puede perder rápidamente la intuición para el Hamiltoniano final resultante. El Hamiltoniano para LiH (utilizando el mapeador Jordan-Wigner) ya consta de 276 términos.

len(build_hamiltonian(ecore, h1e, h2e))

Output:

accuracy of Cholesky decomposition  1.1102230246251565e-16

En caso de duda con respecto a las simetrías, también se puede generar alguna información de simetría para la molécula estableciendo symmetry = True y 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>

Entre otra información útil, devuelve point group symmetry = Coov y también el número de orbitales en cada representación irreducible.

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

Esto no le indica necesariamente cuántos orbitales quiere incluir en su espacio activo, pero le ayuda a ver qué orbitales están presentes y sus simetrías.

Especificar la simetría y los orbitales suele ser útil, pero también puede especificar el número de orbitales que desea incluir. Consideremos el caso del eteno, a continuación. Utilizando verbose = 4, podemos imprimir las simetrías de los distintos orbitales:

# 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>

Obtenemos:

número de orbitales de irrep Ag = 4

número de orbitales de irrep B2g = 2

número de orbitales de irrep B3g = 1

número de orbitales de irrep B1u = 4

número de orbitales de irrep B2u = 1

número de orbitales de irrep B3u = 2

Pero en lugar de especificar todos los orbitales por simetría, podemos simplemente escribir:

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

En este enfoque, tomamos varios orbitales cercanos al nivel de llenado (valencia y desocupados). Aquí se han seleccionado 5 orbitales para incluirlos en el espacio activo (del 6º al 10º).

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

Output:

6 10

software de terceros

Existen varios paquetes de software desarrollados para la química cuántica, algunos de los cuales ofrecen múltiples mapeadores y herramientas para restringir los espacios activos. Los pasos descritos anteriormente son generales y se aplican también al software de terceros. Pero este otro software puede devolver Hamiltonianos en un formato que no es aceptado por Qiskit. Por ejemplo, algunos programas informáticos devuelven Hamiltonianos de la forma:

H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3]

Obsérvese, en particular, que las puertas están numeradas y que no se muestran los operadores de identidad. Esto contrasta con los Hamiltonianos utilizados en Qiskit, que escribirían el término [Z2 Z3] como ZZII (los qubits 0 y 1 actuando sobre ellos el operador identidad, los qubits 2 y 3 actuando sobre ellos el operador Z, ordenados con el qubit 0 más a la derecha).

Para acomodar cualquier flujo de trabajo existente, el bloque de código siguiente convierte de una sintaxis a la otra. La función convert_openfermion_to_qiskit toma como argumentos un Hamiltoniano generado en OpenFermion o Tangelo (y ya mapeado en operadores Pauli utilizando cualquier mapeador disponible), y el número de qubits necesarios para la molécula.

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)

Además, este cuaderno Python contiene código de ejemplo completo para migrar hamiltonianos de otros flujos de trabajo de software a Qiskit, incluida la conversión anterior.

Ahora debería tener un arsenal de herramientas para obtener el Hamiltoniano que necesita para realizar cálculos de química cuántica en IBM® quantum computers.

¿Le ha resultado útil esta página?

Informe de un error, de una errata o solicite contenido en GitHub.