Sfrutta trim SQD per una migliore selezione delle configurazioni

Il comportamento predefinito di diagonalize_fermionic_hamiltonian è diagonalizzare ogni batch sottocampionato e conservare il risultato migliore. Trim SQD invece tratta i batch come candidati da vagliare: ciascuno viene diagonalizzato, ridotto (trimmed) alle configurazioni con il peso maggiore, e i sopravvissuti di tutti i batch vengono uniti e diagonalizzati insieme. Quella diagonalizzazione unita è quella la cui energia viene riportata dall'iterazione, e il suo output ridotto alimenta l'iterazione successiva.

L'implementazione di trim SQD si basa su una ricerca congiunta di Cleveland Clinic Foundation, RIKEN e IBM®. Se usi questa funzionalità nella tua ricerca, cita anche l'articolo. Puoi trovare la citazione nel file CITATION.bib, con la chiave di citazione merz2026crossing.

Per maggiori dettagli sul codice SQD usato in questo esempio, consulta il tutorial sull'hamiltoniana chimica.

Configura la molecola

# 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

Esegui il ciclo

Un'iterazione di trim SQD fa quanto segue:

1. Ciclo di vaglio. Suddividi le stringhe di bit in batch. Diagonalizza ogni batch e conserva la frazione trim_ratio delle sue stringhe CI, per settore di spin, ordinate per peso.

2. Ciclo di unione. Diagonalizza una sola volta i sopravvissuti uniti. Questo è il risultato che l'iterazione riporta, e le sue stringhe ridotte alimentano l'iterazione successiva.

Il sci_solver del ciclo esegue entrambi i cicli e, per impostazione predefinita, usa il solver PySCF di questo pacchetto.

Abilita trim SQD passando TrimPolicy come segue:

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

Verifica se il vaglio sta aiutando

Il vaglio ha il massimo da offrire quando le configurazioni campionate superano di gran lunga ciò che si può diagonalizzare in una sola volta, così che ordinare un pool più grande trova davvero configurazioni che un singolo batch si sarebbe perso. Ha il minimo da offrire quando un sottospazio della dimensione accessibile cattura già la maggior parte della funzione d'onda.

Il vaglio dipende dal fatto che i batch differiscano l'uno dall'altro. TrimPolicy gestisce questo aspetto modificando il comportamento di diagonalize_fermionic_hamiltonian in modo da estrarre un unico pool di samples_per_batch * num_batches stringhe di bit e partizionarlo fra i batch, così nessuna stringa di bit compare in più di un batch.

Limita il sottospazio

Un sottospazio qui è generato dal prodotto cartesiano di un array di stringhe CI di spin alfa e di uno di spin beta.

Il trimming ordina le configurazioni, ma ciò che viene riportato sono due liste di stringhe per spin. Conservare le 150 configurazioni di peso più alto di un sottospazio coinvolge tipicamente circa 150 stringhe alfa distinte e 150 stringhe beta distinte, il cui prodotto genera circa 22.500 configurazioni, non 150. Quindi il sottospazio passato all'iterazione successiva è in genere molto più grande dell'insieme che è stato ordinato e cresce quadraticamente in ciò che conservi.

Il trimming e il riporto sono limitati separatamente, perché avvengono in punti diversi: il trimming riduce ogni batch del ciclo di vaglio, mentre il riporto è ciò che passa all'iterazione successiva.

Quattro argomenti lo limitano ed è utile impostarne almeno uno:

• trim_ratio è la frazione delle stringhe di ogni batch, per settore di spin, che sopravvive al vaglio. Valori più bassi riducono il sottospazio unito ma scartano configurazioni che l'iterazione successiva deve riscoprire tramite campionamento.

• max_strings_per_trim limita le stringhe che ogni batch fornisce, per settore, indipendentemente dal rapporto. Utile con symmetrize_spin, che unisce entrambi i settori di spin in un'unica lista e quindi raddoppia circa ciò che ogni batch offre. Questo limite è per batch: il sottospazio unito può contenere fino a num_batches volte tante stringhe, meno ciò che i batch avevano in comune, che in pratica non è molto.

• carryover_threshold, carryover_ratio e max_carryover limitano le stringhe che alimentano l'iterazione successiva, applicati al sottospazio unito nel suo insieme anziché per batch. carryover_threshold è la soglia di ampiezza che usa StandardPolicy, con lo stesso valore predefinito; gli altri due limitano il numero, cosa che la sola soglia non fa.

• max_dim, su diagonalize_fermionic_hamiltonian stessa, limita ogni settore di spin di ogni sottospazio che il ciclo diagonalizza, compreso quello unito. Agisce come rete di sicurezza e si applica a ciò che la policy restituisce, dopo che il riporto è stato unito a nuovi campioni.

L'esecuzione sopra imposta insieme max_strings_per_trim=250 e max_dim=250, così che né ciò che ogni batch fornisce né il sottospazio costruito a partire da esso possano sfuggire al controllo.

Usa un eigensolver esterno

Puoi usare trim SQD con un eigensolver esterno passando l'argomento 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,
)

Lo snippet precedente richiede che il pacchetto sbd-eigensolver sia installato da PyPI.

Distribuisci il lavoro

Con TrimPolicy, il ciclo chiama sci_solver due volte per iterazione: una con tutti i batch e una con il sottospazio unito. Entrambe le chiamate passano attraverso la normale interfaccia sci_solver, che il supporto multi-processo di questo pacchetto invoca collettivamente su ogni rank MPI.

Sono potenzialmente disponibili due livelli di parallelismo e sono indipendenti:

• All'interno di una diagonalizzazione. Un solver come SBD può distribuire una singola diagonalizzazione tra i rank.

• Tra i batch. Un solver potrebbe, in linea di principio, diagonalizzare i batch in modo concorrente.

Dato un insieme di batch, il solver determina come delegare il lavoro ai suoi processi worker.