Aprovecha trim SQD para mejorar la selección de configuraciones

El comportamiento predeterminado de diagonalize_fermionic_hamiltonian es diagonalizar cada lote submuestreado y quedarse con el mejor resultado. Trim SQD, en cambio, trata los lotes como candidatos que se deben filtrar: cada uno se diagonaliza, se recorta hasta las configuraciones con mayor peso, y los supervivientes de todos los lotes se fusionan y se diagonalizan juntos. Esa diagonalización fusionada es la que da la energía que informa la iteración, y su propia salida recortada inicializa la siguiente iteración.

La implementación de trim SQD se basa en una investigación conjunta de la Cleveland Clinic Foundation, RIKEN e IBM®. Si usas esta función en tu investigación, cita también el artículo. Puedes encontrar la cita en el archivo CITATION.bib, con la clave de cita merz2026crossing.

Para más detalles sobre el código de SQD usado en este ejemplo, consulta el tutorial del hamiltoniano químico.

Configura la molécula

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit_addon_sqd": "qiskit-addon-sqd"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import numpy as np
import pyscf
import pyscf.mcscf
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
from qiskit_addon_sqd.trim import TrimPolicy

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

# Create a seed to control randomness throughout this workflow
rng = np.random.default_rng(24)

# Generate random samples
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)

WARN: Unable to to identify input symmetry using original axes.
Different symmetry axes will be used.

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

Ejecuta el bucle

Una iteración de trim SQD hace lo siguiente:

1. Ronda de filtrado. Divide las cadenas de bits en lotes. Diagonaliza cada lote y conserva la fracción trim_ratio de sus cadenas CI, por sector de espín, ordenadas por peso.

2. Ronda fusionada. Diagonaliza una vez los supervivientes fusionados. Este es el resultado que informa la iteración, y sus propias cadenas recortadas inicializan la siguiente iteración.

El sci_solver del bucle realiza ambas rondas y, por defecto, usa el solver PySCF de este paquete.

Habilita trim SQD pasando TrimPolicy de la siguiente manera:

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    max_iterations=5,
    policy=TrimPolicy(trim_ratio=0.5, max_strings_per_trim=250),
    symmetrize_spin=True,
    max_dim=250,
    seed=rng,
)

sci_state = result.sci_state
print(
    f"Merged subspace: {len(sci_state.ci_strs_a)} x {len(sci_state.ci_strs_b)} CI strings"
)
print(f"Exact energy:     {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")

Merged subspace: 148 x 148 CI strings
Exact energy:     -109.04667177808024
Estimated energy: -109.034369147063

Comprueba si el filtrado está ayudando

El filtrado ofrece más cuando las configuraciones muestreadas superan con creces lo que se puede diagonalizar a la vez, de modo que ordenar un conjunto más grande encuentra realmente configuraciones que un solo lote habría pasado por alto. Ofrece menos cuando un subespacio del tamaño asequible ya captura la mayor parte de la función de onda.

El filtrado depende de que los lotes difieran entre sí. TrimPolicy se encarga de esto modificando el comportamiento de diagonalize_fermionic_hamiltonian para extraer un único conjunto de samples_per_batch * num_batches cadenas de bits y repartirlo entre los lotes, de modo que ninguna cadena de bits aparezca en más de uno.

Acota el subespacio

Un subespacio aquí está generado por el producto cartesiano de un arreglo de cadenas CI de espín alfa y otro de espín beta.

El recorte ordena configuraciones, pero lo que se transfiere son dos listas de cadenas por espín. Conservar las 150 configuraciones de mayor peso de un subespacio suele implicar cerca de 150 cadenas alfa distintas y 150 cadenas beta distintas, cuyo producto abarca aproximadamente 22.500 configuraciones, no 150. Por tanto, el subespacio que se pasa a la siguiente iteración es en general mucho mayor que el conjunto que se ordenó, y crece de forma cuadrática con lo que conserves.

El recorte y la transferencia se acotan por separado, porque ocurren en puntos distintos: el recorte reduce cada lote de la ronda de filtrado, mientras que la transferencia es lo que pasa a la siguiente iteración.

Cuatro argumentos lo acotan, y conviene establecer al menos uno:

• trim_ratio es la fracción de las cadenas de cada lote, por sector de espín, que sobreviven al filtrado. Los valores más bajos reducen el subespacio fusionado, pero descartan configuraciones que la siguiente iteración debe redescubrir mediante muestreo.

• max_strings_per_trim limita las cadenas que aporta cada lote, por sector, independientemente de la proporción. Es útil con symmetrize_spin, que fusiona ambos sectores de espín en una sola lista y por tanto duplica aproximadamente lo que ofrece cada lote. Este límite es por lote: el subespacio fusionado puede contener hasta num_batches veces esa cantidad de cadenas, menos lo que los lotes tuvieran en común, que en la práctica no es mucho.

• carryover_threshold, carryover_ratio y max_carryover acotan las cadenas que inicializan la siguiente iteración, y se aplican al subespacio fusionado en su conjunto y no por lote. carryover_threshold es el umbral de amplitud que usa StandardPolicy, con el mismo valor predeterminado; los otros dos acotan el número de cadenas, algo que el umbral por sí solo no hace.

• max_dim, en la propia diagonalize_fermionic_hamiltonian, acota cada sector de espín de todos los subespacios que diagonaliza el bucle, incluido el fusionado. Actúa como respaldo y se aplica a lo que devuelva la política, después de fusionar la transferencia con muestras nuevas.

La ejecución anterior establece max_strings_per_trim=250 y max_dim=250 a la vez, de modo que ni lo que aporta cada lote ni el subespacio construido a partir de ello pueden descontrolarse.

Usa un eigensolver externo

Puedes usar trim SQD con un eigensolver externo pasando el argumento sci_solver:

from sbd.sbd_solver import create_sbd_solver

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    policy=TrimPolicy(trim_ratio=0.1),
    sci_solver=create_sbd_solver(sbd_config={"method": 0, "eps": 1e-10}),
    seed=rng,
)

El fragmento anterior requiere que el paquete sbd-eigensolver esté instalado desde PyPI.

Distribuye el trabajo

Con TrimPolicy, el bucle llama a sci_solver dos veces por iteración: una con todos los lotes y otra con el subespacio fusionado. Ambas llamadas pasan por la interfaz ordinaria de sci_solver, que el soporte multiproceso de este paquete invoca de forma colectiva en cada rango de MPI.

Potencialmente hay dos niveles de paralelismo disponibles, y son independientes:

• Dentro de una diagonalización. Un solver como SBD puede repartir una sola diagonalización entre rangos.

• Entre los lotes. En principio, un solver podría diagonalizar los lotes de forma concurrente.

Dado un conjunto de lotes, el solver determina cómo delegar el trabajo a sus procesos de trabajo.