Hamiltonianos para Química Quântica

Vamos começar com uma breve visão geral do papel que os Hamiltonianos desempenham no VQE.

Visão Geral do Hamiltoniano no VQE

A Dra. Victoria Lipinska nos guia pelos Hamiltonianos e explica como mapeá-los para uso em computação quântica.

Referências

Os artigos a seguir são referenciados no vídeo acima.

Preparando Hamiltonianos para Química Quântica

Um bom primeiro passo para aplicar computação quântica a um problema de química é definir um Hamiltoniano para o sistema de interesse. Aqui, vamos restringir a discussão aos Hamiltonianos de química quântica, já que esses Hamiltonianos exigem um mapeamento específico para sistemas de férmions idênticos.

Se você trabalha com química quântica, provavelmente já tem seu software favorito para modelar moléculas, que consegue gerar um Hamiltoniano descrevendo o sistema de interesse. Aqui, vamos usar código construído apenas com PySCF, numpy e Qiskit. Mas o processo de preparação do Hamiltoniano também se aplica a soluções pré-empacotadas. A única diferença entre essa abordagem e outros softwares serão pequenas diferenças de sintaxe; algumas delas são tratadas na subseção "Software de terceiros" para facilitar a integração com fluxos de trabalho existentes.

Gerar um Hamiltoniano de química quântica para uso em QPUs do IBM Quantum® envolve os seguintes passos:

Definir sua molécula (geometria, spin, espaço ativo etc.)
Gerar o Hamiltoniano fermiônico (operadores de criação e aniquilação)
Mapear o Hamiltoniano fermiônico para um operador bosônico (neste contexto, usando operadores de Pauli)
Se usar software de terceiros: tratar eventuais incompatibilidades de sintaxe entre o software gerador e o Qiskit

O Hamiltoniano fermiônico é escrito em termos de operadores fermiônicos e, em particular, leva em conta que elétrons são férmions indistinguíveis. Isso significa que eles obedecem a estatísticas completamente diferentes das de qubits distinguíveis e bosônicos. Daí a necessidade do processo de mapeamento.

Quem já está familiarizado com esses processos provavelmente pode pular esta seção. Objetivo:

O objetivo final é obter um Hamiltoniano na 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])

Vamos começar importando alguns pacotes:

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

Definir sua molécula

Aqui vamos especificar os atributos da molécula de interesse. Neste exemplo, escolhemos o hidrogênio diatômico (porque os Hamiltonianos resultantes são curtos o suficiente para serem exibidos). O Python-based Simulations of Chemistry Framework (PySCF) possui uma ampla coleção de módulos de estrutura eletrônica que podem ser usados para, entre outras coisas, gerar Hamiltonianos moleculares adequados para computação quântica. O guia PySCF Quickstart é um excelente recurso para uma descrição completa de todas as variáveis e funcionalidades. Vamos dar apenas uma visão bem superficial, pois isso já será familiar para muitos de vocês. Para entender melhor, visite o PySCF. Resumidamente:

distance pode ser usada para moléculas diatômicas, ou simplesmente especifique coordenadas cartesianas para cada átomo. As distâncias estão em unidades de Angstrom.

gto gera orbitais do tipo gaussiano.

basis se refere às funções usadas para modelar orbitais moleculares. Aqui, 'sto-6g' é uma base mínima comum, nomeada por ajustar Orbitais do Tipo Slater usando 6 orbitais gaussianos primitivos.

spin um valor inteiro que indica o número de elétrons desemparelhados (igual a $2S$ ). Note que alguns softwares usam multiplicidade ( $2S+1$ ) em vez disso.

charge a carga da molécula.

symmetry — o grupo de simetria pontual da molécula, especificado com uma string ou detectado automaticamente definindo "symmetry = True". Aqui "Dooh" é o grupo de simetria adequado para moléculas diatômicas com dois átomos da mesma espécie.

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>

Tenha em mente que se pode descrever a energia total (que inclui a energia de repulsão nuclear e a eletrônica), a energia orbital eletrônica total, ou a energia de algum subconjunto de orbitais eletrônicos (com o subconjunto complementar congelado). No caso específico do $\text{H}_2$ , observe as diferentes energias abaixo, e note que a energia total menos a energia de repulsão nuclear de fato resulta na energia eletrô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)

Gerar o Hamiltoniano fermiônico

scf se refere a uma ampla gama de métodos de campo autoconsistente.

rhf como em mf = scf.RHF(mol), em que mf é um solver que usa o cálculo de Hartree-Fock Restrito. O kernel disso (E, abaixo) é a energia total, incluindo repulsão nuclear e orbitais moleculares.

mcscf é um pacote de campos autoconsistentes de múltiplas configurações.

ao2mo é uma transformação de orbitais atômicos para orbitais moleculares.

Também usamos as seguintes variáveis:

ncas: número de orbitais no espaço ativo completo

nelecas: número de elétrons no espaço ativo 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 um Hamiltoniano, e ele é frequentemente separado em energia de um núcleo eletrônico (ecore, não envolvido na minimização), operadores de um elétron (h1e) e energias de dois elétrons (h2e). Estes são extraídos explicitamente abaixo nas duas últimas linhas.

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

Esses Hamiltonianos são atualmente operadores fermiônicos (de criação e aniquilação), aplicáveis a sistemas de férmions (indistinguíveis) e, consequentemente, sujeitos à antissimetria sob troca. Isso resulta em estatísticas diferentes das que se aplicariam a um sistema distinguível ou bosônico. Para executar cálculos em QPUs do IBM Quantum, precisamos de um operador bosônico que descreva a energia. O resultado de tal mapeamento é convencionalmente escrito em termos de operadores de Pauli, pois são ao mesmo tempo Hermitianos e unitários. Existem vários mapeamentos que podem ser usados. Um dos mais simples é a transformação de Jordan-Wigner.

Mapeando o Hamiltoniano

Vale ressaltar que há muitas ferramentas disponíveis para mapear um Hamiltoniano químico para um adequado para execução em um computador quântico. Aqui, implementamos o mapeamento de Jordan-Wigner diretamente usando apenas PySCF, numpy e Qiskit. Comentamos abaixo considerações de sintaxe para outras soluções. A função Cholesky nos ajuda a obter uma decomposição de baixo posto dos termos de dois elétrons no 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

As funções identity e creators_destructors substituem os operadores de criação e aniquilação no Hamiltoniano fermiônico por operadores de Pauli; creators_destructors usa o mapeamento 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 fim, build_hamiltonian usa as funções cholesky, identity e creators_destructors para criar o Hamiltoniano final adequado para execução em um computador quâ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 fim, usamos build_hamiltonian para construir nosso Hamiltoniano de qubits a partir de operadores de Pauli, usando a transformação de Jordan-Wigner. Isso também nos dá a precisão da decomposição 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 moléculas de exemplo mostra a configuração e os Hamiltonianos para várias moléculas de complexidade variada; com algumas modificações, isso deve permitir que você examine a maioria das moléculas pequenas.

Vamos apenas mencionar brevemente dois pontos importantes a considerar ao construir os operadores fermiônicos para uma molécula. À medida que o tipo de molécula muda, a simetria também muda. Da mesma forma, o número de orbitais com diversas simetrias, como o "A1" com simetria cilíndrica, vai mudar. Essas mudanças ficam evidentes mesmo com a simples extensão ao LiH, como vemos aqui:

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 também mencionar que podemos rapidamente perder a intuição sobre o Hamiltoniano final resultante. O Hamiltoniano do LiH (usando o mapeador de Jordan-Wigner) já consiste em 276 termos.

len(build_hamiltonian(ecore, h1e, h2e))

accuracy of Cholesky decomposition  1.1102230246251565e-16

Em caso de dúvida sobre as simetrias, também é possível gerar algumas informações de simetria para a molécula definindo 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>

Entre outras informações úteis, isso retorna tanto point group symmetry = Coov quanto o número de orbitais em cada representação irredutível.

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

Isso não necessariamente diz quantos orbitais você quer incluir no seu espaço ativo, mas ajuda a ver quais orbitais estão presentes e suas simetrias.

Especificar simetria e orbitais costuma ser útil, mas você também pode especificar o número de orbitais que deseja incluir. Considere o caso do eteno, abaixo. Usando verbose = 4, podemos imprimir as simetrias dos vários orbitais:

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

Obtemos:

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

Mas em vez de especificar todos os orbitais por simetria, podemos simplesmente escrever:

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

Nessa abordagem, pegamos alguns orbitais próximos ao nível de preenchimento (valência e desocupados). Aqui, 5 orbitais foram selecionados para inclusão no espaço ativo (do 6º ao 10º).

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

6 10

Software de Terceiros

Existem vários pacotes de software desenvolvidos para química quântica, alguns oferecendo múltiplos mapeadores e ferramentas para restringir espaços ativos. Os passos descritos acima são gerais e se aplicam a softwares de terceiros também. Porém, esses outros softwares podem retornar Hamiltonianos em um formato que não é aceito pelo Qiskit. Por exemplo, alguns softwares retornam Hamiltonianos na forma:

H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3] Note em particular que as portas são numeradas e os operadores identidade não são mostrados. Isso contrasta com os Hamiltonianos usados no Qiskit, que escreveria o termo [Z2 Z3] como ZZII (qubits 0 e 1 sendo atuados pelo operador identidade, qubits 2 e 3 sendo atuados pelo operador Z, ordenados com o qubit 0 mais à direita).

Para acomodar quaisquer fluxos de trabalho existentes que você tenha, o bloco de código abaixo converte de uma sintaxe para a outra. A função convert_openfermion_to_qiskit recebe como argumentos um Hamiltoniano gerado no OpenFermion ou no Tangelo (já mapeado para operadores de Pauli usando qualquer mapeador disponível), e o número de qubits necessários para a 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)

Além disso, este notebook Python contém código de exemplo completo para migrar Hamiltonianos de outros fluxos de trabalho de software para o Qiskit, incluindo a conversão acima.

Você agora deve ter um arsenal de ferramentas para obter o Hamiltoniano necessário para realizar cálculos de química quântica em computadores quânticos IBM®.

Visão Geral do Hamiltoniano no VQE​

Referências​

Preparando Hamiltonianos para Química Quântica​

Visão Geral do Hamiltoniano no VQE

Referências

Preparando Hamiltonianos para Química Quântica