Scala i flussi di lavoro chimici SQD con Fulqrum

Fulqrum è un pacchetto software che calcola autovalori e autovettori di Hamiltoniane proiettate in un sottospazio. Costruisce una rappresentazione portabile dell'Hamiltoniana proiettata che puoi passare a un solver di autovalori di tua scelta, inclusi i comuni solver già pronti. Questi solver possono offrire un'alternativa molto più veloce al solver PySCF predefinito che l'addon SQD usa per la diagonalizzazione. Per informazioni su come installare e usare Fulqrum, consulta la documentazione di Fulqrum.

Il modulo addon_sqd_support di Fulqrum costruisce un solver che diagonalize_fermionic_hamiltonian accetta direttamente, quindi passi gli stessi integrali a entrambi e Fulqrum gestisce la proiezione e la diagonalizzazione.

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

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

_needed = {"fulqrum": "fulqrum", "numpy": "numpy", "primme": "primme", "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.cc
import pyscf.mcscf
from fulqrum.addon_sqd_support import sci_solver_from_integrals
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian

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

# Build a solver that diagonalizes with Fulqrum. It converts the integrals to a
# Fulqrum operator once, then projects onto each subspace the SQD loop requests.
sci_solver = sci_solver_from_integrals(hcore, eri, tol=1e-12)

# Run SQD
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,
    sci_solver=sci_solver,
    symmetrize_spin=True,
    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

print(f"Exact energy: {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")

Exact energy: -109.04667177808024
Estimated energy: -109.03402667559365

Usare un eigensolver diverso

Poiché Fulqrum espone l'Hamiltoniana proiettata come un scipy.sparse.linalg.LinearOperator standard, puoi diagonalizzarla con qualsiasi eigensolver di tua scelta anziché essere vincolato a scipy.sparse.linalg.eigsh. Questo rende semplice sostituirlo con un solver più adatto alle dimensioni del tuo problema o al tuo hardware.

Passa l'eigensolver che desideri come argomento eigensolver. Come esempio, questa guida mostra PRIMME, una libreria di eigensolver ad alte prestazioni. La sua interfaccia Python, primme.eigsh, segue una convenzione di chiamata simile a scipy.sparse.linalg.eigsh, quindi puoi usarla come sostituto diretto.

import primme

# Build a solver that diagonalizes with PRIMME instead of the default,
# scipy.sparse.linalg.eigsh
sci_solver = sci_solver_from_integrals(
    hcore, eri, eigensolver=primme.eigsh, tol=1e-12
)

# Run SQD using the PRIMME-based 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,
    max_iterations=5,
    sci_solver=sci_solver,
    symmetrize_spin=True,
    seed=rng,
)

print(f"Exact energy: {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")

Exact energy: -109.04667177808024
Estimated energy: -109.03937429196168