Hamiltonianos para la química cuántica
Comencemos con una breve descripción del papel que desempeñan los hamiltonianos en VQE.
Descripción general: El hamiltoniano en VQE
La Dra. Victoria Lipinska nos guía a través de los hamiltonianos y explica cómo se mapean para su uso en computación cuántica.
Referencias
Los siguientes artículos se mencionan en el video anterior.
- Quantum Algorithms for Fermionic Simulations, Ortiz, et al.
- Simulating Chemistry using Quantum Computers, 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.
- Quantum computational chemistry, McArdle, et al.
- The Bravyi-Kitaev transformation for quantum computation of electronic structure, Seeley, et al., McArdle, et al.
Preparar hamiltonianos para la química cuántica
Un buen primer paso al aplicar la computación cuántica a un problema de química es definir un hamiltoniano para el sistema considerado. Aquí nos limitamos a hamiltonianos de química cuántica, ya que estos requieren un mapeo especial para sistemas de fermiones idénticos.
Como profesional que trabaja en química cuántica, probablemente ya tengas tu software preferido para modelar moléculas, con el que se puede generar un hamiltoniano que describa tu sistema de interés. Aquí usamos código que se basa exclusivamente en PySCF, numpy y Qiskit. Sin embargo, el proceso de preparación del hamiltoniano también se puede trasladar a soluciones prediseñadas. La única diferencia entre este enfoque y otro software radica en pequeñas diferencias de sintaxis; algunas de ellas se tratan en la subsección "Software de terceros" para facilitar la integración de flujos de trabajo existentes.
La generación de un hamiltoniano de química cuántica para su uso en QPUs de IBM Quantum® comprende los siguientes pasos:
- Define tu molécula (geometría, espín, espacio activo, etc.)
- Genera el hamiltoniano fermiónico (operadores de creación y aniquilación)
- Mapea el hamiltoniano fermiónico a un operador bosónico (en este contexto, mediante operadores de Pauli)
- Si usas software de terceros: maneja posibles diferencias de sintaxis entre el software generador y Qiskit
El hamiltoniano fermiónico se escribe en forma de operadores fermiónicos y tiene en cuenta, en particular, que los electrones son fermiones indistinguibles. Esto significa que obedecen una estadística completamente diferente a la de qubits distinguibles y bosónicos. Por lo tanto, el mapeo es necesario.
Quienes ya estén familiarizados con estos procesos probablemente puedan saltarse esta sección. Objetivo:
El objetivo final es obtener un hamiltoniano de la siguiente 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')]
O
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])
Comenzamos importando algunos paquetes:
import numpy as np
from pyscf import ao2mo, gto, mcscf, scf
- Define tu molécula
Aquí establecemos las propiedades de la molécula considerada. En este ejemplo hemos elegido hidrógeno diatómico (ya que los hamiltonianos resultantes son lo suficientemente cortos como para visualizarlos). El framework de simulación basado en Python para química (PySCF) ofrece una amplia colección de módulos de estructura electrónica que, entre otras cosas, pueden utilizarse para generar hamiltonianos moleculares para cálculos cuánticos. La guía de inicio rápido de PySCF es un excelente recurso para una descripción completa de todas las variables y funciones. Aquí solo damos una breve descripción general, ya que esto probablemente sea familiar para muchos de ustedes. Para más información, visita PySCF. En resumen:
distance puede usarse para moléculas diatómicas, o simplemente proporciona coordenadas cartesianas para cada átomo. Las distancias están en ángstroms.
gto genera orbitales de tipo gaussiano.
basis se refiere a las funciones utilizadas para modelar orbitales moleculares. Aquí 'sto-6g' es una base mínima común, nombrada por el ajuste de orbitales de tipo Slater con 6 gaussianos primitivos.
spin un valor entero que indica el número de electrones desapareados (corresponde a ). Ten en cuenta que algún software utiliza en su lugar la multiplicidad ().
charge la carga de la molécula.
symmetry - el grupo de simetría puntual de la molécula, especificado como cadena de texto o detectado automáticamente estableciendo "symmetry = True". "Dooh" es el grupo de simetría apropiado para moléculas diatómicas con dos átomos de la misma especie.
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>
Ten en cuenta que se puede describir la energía total (que incluye también la energía de repulsión nuclear así como la energía electrónica), la energía orbital electrónica total o la energía de un subconjunto de orbitales electrónicos (donde el subconjunto complementario está congelado). En el caso especial de , las diferentes energías se enumeran a continuación; la energía total menos la energía de repulsión nuclear da efectivamente 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(),
)
0.7199689944489797 -1.8455976628764188 -1.125628668427439 -1.8455976628764188
active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)
- Genera el hamiltoniano fermiónico
scf se refiere a una amplia gama de métodos de campo autoconsistente (Self-Consistent Field).
rhf como en mf = scf.RHF(mol): mf es un solucionador que utiliza el cálculo Hartree-Fock restringido. El kernel 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 multi-configuración.
ao2mo es una transformación de orbitales atómicos a orbitales moleculares.
También usamos 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]
Necesitamos un hamiltoniano, que a menudo se divide en la energía del núcleo electrónico (ecore, no participa en la minimización), operadores de un electrón (h1e) y energías de dos electrones (h2e). Estos 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 (de creación y aniquilación) aplicables a sistemas de fermiones (indistinguibles) y sujetos a la antisimetría bajo intercambio. Esto conduce a una estadística diferente a la de un sistema distinguible o bosónico. Para realizar cálculos en QPUs de IBM Quantum, necesitamos un operador bosónico que describa la energía. El resultado de tal mapeo se escribe convencionalmente en forma de operadores de Pauli, ya que estos son tanto hermitianos como unitarios. Existen varios mapeos que se pueden usar. Uno de los más simples es la transformación de Jordan-Wigner.
- Mapeo del hamiltoniano
Cabe señalar que existen muchas herramientas para mapear un hamiltoniano químico a uno adecuado para computadores cuánticos. Aquí implementamos el mapeo Jordan-Wigner directamente, usando exclusivamente PySCF, numpy y Qiskit. A continuación comentamos consideraciones de sintaxis para otras soluciones. La función de 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 reemplazan los operadores de creación y aniquilación en el hamiltoniano fermiónico por operadores de Pauli; creators_destructors utiliza el mapeo 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
Finalmente, build_hamiltonian utiliza las funciones cholesky, identity y creators_destructors para crear el hamiltoniano final adecuado para su uso en un computador 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()
Finalmente, usamos build_hamiltonian para construir nuestro hamiltoniano de qubits a partir de operadores de Pauli usando la transformación de Jordan-Wigner. Esto también nos proporciona la precisión de la descomposición de Cholesky utilizada.
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])
Este notebook de ejemplo de moléculas muestra la configuración y hamiltonianos para varias moléculas de diferente complejidad; con pocas adaptaciones, debería permitirte estudiar la mayoría de las moléculas pequeñas.
Brevemente, dos puntos importantes a tener en cuenta al construir los operadores fermiónicos para una molécula: cuando cambia el tipo de molécula, también cambia la simetría. Igualmente, cambia el número de orbitales con diferentes simetrías, como la "A1" con simetría cilíndrica. Estos cambios ya son evidentes en la simple extensió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 mencionar que se pierde rápidamente la intuición sobre el hamiltoniano resultante. El hamiltoniano para LiH (con el mapeador Jordan-Wigner) ya consta de 276 términos.
len(build_hamiltonian(ecore, h1e, h2e))
accuracy of Cholesky decomposition 1.1102230246251565e-16
276
Cuando hay incertidumbre respecto a las simetrías, también se puede generar 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,
)
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 otras cosas, esto proporciona tanto point group symmetry = Coov como 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 te dice necesariamente cuántos orbitales quieres incluir en tu espacio activo, pero te ayuda a reconocer qué orbitales están presentes y qué simetrías tienen.
Establecer la simetría y los orbitales suele ser útil, pero también puedes simplemente especificar el número de orbitales que deseas incluir. Considera el ejemplo del etileno a continuación. Con verbose = 4 podemos imprimir las simetrías de los diferentes 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,
)
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:
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
Sin embargo, en lugar de especificar todos los orbitales por simetría, podemos simplemente escribir:
active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)
Con este enfoque, seleccionamos varios orbitales cerca del nivel de llenado (orbitales de valencia y desocupados). Aquí se seleccionaron 5 orbitales para su inclusión en el espacio activo (del 6.o al 10.o).
print(
mol.nelectron // 2 - 2,
mol.nelectron // 2 + 2,
)
6 10
- Software de terceros
Existen varios paquetes de software para química cuántica, algunos de los cuales ofrecen múltiples mapeadores y herramientas para restringir espacios activos. Los pasos descritos anteriormente son generales y se aplican también al software de terceros. Sin embargo, dicho software puede devolver hamiltonianos en un formato que Qiskit no acepta. Por ejemplo, algunos programas devuelven hamiltonianos en la siguiente forma:
H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3]
Ten en cuenta en particular que las puertas están numeradas y los operadores de identidad no se muestran. Esto contrasta con los hamiltonianos en Qiskit, donde el término [Z2 Z3] se escribiría como ZZII (los qubits 0 y 1 son afectados por el operador identidad, los qubits 2 y 3 por el operador Z, con el qubit 0 a la derecha).
Para facilitar la integración de flujos de trabajo existentes, el siguiente bloque de código convierte de una sintaxis a otra. La función convert_openfermion_to_qiskit toma como argumentos un hamiltoniano generado en OpenFermion o Tangelo (que ya ha sido mapeado a operadores de Pauli con cualquier mapeador disponible) así como 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 notebook de Python contiene código de ejemplo completo para la migración de hamiltonianos desde otros flujos de trabajo de software a Qiskit, incluyendo la conversión anterior.
Ahora tienes una colección de herramientas para obtener el hamiltoniano que necesitas para cálculos de química cuántica en computadores cuánticos de IBM®.