Tirer parti du trim SQD pour améliorer la sélection des configurations

Le comportement par défaut de diagonalize_fermionic_hamiltonian consiste à diagonaliser chaque lot sous-échantillonné et à conserver le meilleur résultat. Trim SQD traite plutôt les lots comme des candidats à filtrer : chacun est diagonalisé, puis élagué pour ne garder que les configurations portant le plus grand poids, et les survivantes de tous les lots sont fusionnées et diagonalisées ensemble. C'est cette diagonalisation fusionnée dont l'énergie est rapportée par l'itération, et sa propre sortie élaguée amorce l'itération suivante.

L'implémentation de trim SQD repose sur une recherche conjointe de la Cleveland Clinic Foundation, de RIKEN et d'IBM®. Si tu utilises cette fonctionnalité dans tes recherches, cite aussi l'article. Tu trouveras la citation dans le fichier CITATION.bib, avec la clé de citation merz2026crossing.

Pour plus de détails sur le code SQD utilisé dans cet exemple, consulte le tutoriel sur le hamiltonien chimique.

Configurer la molécule

# 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

Exécuter la boucle

Une itération de trim SQD effectue les opérations suivantes :

1. Tour de filtrage. Répartis les chaînes de bits en lots. Diagonalise chaque lot et conserve la fraction trim_ratio de ses chaînes CI, par secteur de spin, classées par poids.

2. Tour fusionné. Diagonalise une fois les survivantes fusionnées. C'est le résultat que rapporte l'itération, et ses propres chaînes élaguées amorcent l'itération suivante.

Le sci_solver de la boucle effectue les deux tours et utilise par défaut le solveur PySCF de ce package.

Active trim SQD en passant le TrimPolicy comme suit :

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

Vérifier si le filtrage est utile

Le filtrage est le plus utile lorsque les configurations échantillonnées sont beaucoup plus nombreuses que ce qui peut être diagonalisé en une fois, de sorte que le classement d'un plus grand pool trouve réellement des configurations qu'un lot unique aurait manquées. Il est le moins utile lorsqu'un sous-espace de taille abordable capture déjà l'essentiel de la fonction d'onde.

Le filtrage dépend du fait que les lots diffèrent les uns des autres. TrimPolicy y veille en modifiant le comportement de diagonalize_fermionic_hamiltonian pour tirer un seul pool de samples_per_batch * num_batches chaînes de bits et le répartir entre les lots, de sorte qu'aucune chaîne de bits n'apparaisse dans plus d'un lot.

Borner le sous-espace

Un sous-espace est ici engendré par le produit cartésien d'un tableau de chaînes CI de spin alpha et d'un tableau de chaînes CI de spin bêta.

L'élagage classe les configurations, mais ce qui est reporté est constitué de deux listes de chaînes par spin. Conserver les 150 configurations de plus grand poids d'un sous-espace fait généralement intervenir près de 150 chaînes alpha distinctes et 150 chaînes bêta distinctes, dont le produit engendre environ 22 500 configurations, et non 150. Le sous-espace transmis à l'itération suivante est donc généralement bien plus grand que l'ensemble qui a été classé, et il croît de façon quadratique avec ce que tu conserves.

L'élagage et le report sont bornés séparément, car ils interviennent à des moments différents : l'élagage réduit chaque lot du tour de filtrage, tandis que le report est ce qui passe à l'itération suivante.

Quatre arguments le bornent, et il est utile d'en définir au moins un :

• trim_ratio est la fraction des chaînes de chaque lot, par secteur de spin, qui survivent au filtrage. Des valeurs plus basses réduisent le sous-espace fusionné, mais écartent des configurations que l'itération suivante devra redécouvrir par échantillonnage.

• max_strings_per_trim plafonne les chaînes que chaque lot apporte, par secteur, quel que soit le ratio. Utile avec symmetrize_spin, qui fusionne les deux secteurs de spin en une seule liste et double donc à peu près ce que chaque lot offre. Cette borne s'applique par lot : le sous-espace fusionné peut contenir jusqu'à num_batches fois plus de chaînes, moins ce que les lots avaient en commun, ce qui en pratique n'est pas beaucoup.

• carryover_threshold, carryover_ratio et max_carryover bornent les chaînes qui amorcent l'itération suivante, appliquées au sous-espace fusionné dans son ensemble plutôt que par lot. carryover_threshold est le seuil d'amplitude qu'utilise StandardPolicy, avec la même valeur par défaut ; les deux autres bornent le nombre, ce que le seuil seul ne fait pas.

• max_dim, sur diagonalize_fermionic_hamiltonian lui-même, borne chaque secteur de spin de chaque sous-espace que la boucle diagonalise, y compris le sous-espace fusionné. Il joue le rôle de filet de sécurité et s'applique à tout ce que renvoie la politique, après la fusion du report avec de nouveaux échantillons.

L'exécution ci-dessus définit ensemble max_strings_per_trim=250 et max_dim=250, de sorte que ni ce que chaque lot apporte ni le sous-espace construit à partir de cet apport ne peuvent s'emballer.

Utiliser un eigensolver externe

Tu peux utiliser trim SQD avec un eigensolver externe en passant l'argument 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,
)

L'extrait précédent nécessite l'installation du package sbd-eigensolver depuis PyPI.

Distribuer le travail

Avec TrimPolicy, la boucle appelle sci_solver deux fois par itération : une fois avec tous les lots, et une fois avec le sous-espace fusionné. Les deux appels passent par l'interface ordinaire de sci_solver, que la prise en charge multiprocessus de ce package invoque collectivement sur chaque rang MPI.

Deux niveaux de parallélisme sont potentiellement disponibles, et ils sont indépendants :

• Au sein d'une diagonalisation. Un solveur tel que SBD peut répartir une seule diagonalisation entre plusieurs rangs.

• Entre les lots. Un solveur pourrait, en principe, diagonaliser les lots de façon concurrente.

Étant donné un ensemble de lots, le solveur détermine comment déléguer le travail à ses processus de travail.