Adaptez vos flux de travail chimiques SQD à grande échelle grâce à Fulqrum

Fulqrum est un logiciel qui calcule les valeurs propres et les vecteurs propres des hamiltoniens projetés dans un sous-espace. Il construit une représentation portable de l'hamiltonien projeté que vous pouvez transmettre à un solveur de valeurs propres de votre choix, y compris les solveurs standard disponibles dans le commerce. Ces solveurs peuvent constituer une alternative bien plus rapide au solveur PySCF par défaut utilisé par l'extension SQD pour la diagonalisation. Pour plus d'informations sur l'installation et l'utilisation de Fulqrum, consultez la documentation de Fulqrum.

Pour plus de détails sur le code SQD utilisé dans cet exemple, consultez le tutoriel sur la chimie de Hamilton.

import fulqrum as fq
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
import scipy.sparse.linalg
from fulqrum.convert.integrals import integrals_to_fq_fermionic_op
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import (
    SCIResult,
    SCIState,
    diagonalize_fermionic_hamiltonian,
)

# Specify molecule properties
num_orbitals = 16
num_elec_a = num_elec_b = 5
spin_sq = 0

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

# Convert Hamiltonian to Fulqrum format
fermionic_op = integrals_to_fq_fermionic_op(
    one_body_integrals=hcore, two_body_integrals=eri
)
fulqrum_operator = fermionic_op.extended_jw_transformation()

# Initialize a Hamiltonian with a trivial subspace to be updated later
Hsub = fq.SubspaceHamiltonian(
    fulqrum_operator,
    fq.Subspace([["0" * num_orbitals], ["0" * num_orbitals]]),
)


# Define fulqrum solver
def solve_sci_batch(ci_batches, _h1, _h2, norb, nelec, tol=1e-12):
    # Ignore the _h1 and _h2 inputs because we store the Hamiltonian separately
    results = []
    for ints_a, ints_b in ci_batches:
        strs_a = [format(x, f"0{norb}b") for x in ints_a]
        strs_b = [format(x, f"0{norb}b") for x in ints_b]
        subspace = fq.Subspace([strs_a, strs_b])
        Hsub.update_subspace(subspace)
        Hsub_csr_linop = Hsub.to_csr_linearoperator_fast()
        diag = Hsub.diagonal_vector()
        v0 = np.zeros(len(subspace), dtype=Hsub.dtype)
        v0[np.argmin(diag)] = 1.0
        eigs, vecs = scipy.sparse.linalg.eigsh(
            Hsub_csr_linop, k=1, which="SA", tol=tol, v0=v0
        )
        vec = vecs[:, 0]
        occs = subspace.get_orbital_occupancies(
            probs=np.abs(vec) ** 2, norb=norb
        )
        # SQD ordering and Fulqrum ordering are transpose of each other.
        # ascontiguousarray() is needed for matching C contiguous ordering inside PySCF
        state = SCIState(
            amplitudes=np.ascontiguousarray(
                vec.reshape(len(strs_b), len(strs_a)).T
            ),
            ci_strs_a=ints_a,
            ci_strs_b=ints_b,
            norb=norb,
            nelec=nelec,
        )
        results.append(
            SCIResult(
                energy=eigs[0], sci_state=state, orbital_occupancies=occs
            )
        )
    return results


# 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=solve_sci_batch,
    symmetrize_spin=True,
    seed=rng,
)

Output:

converged SCF energy = -108.835236570775
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}")

Output:

Exact energy: -109.04667177808028
Estimated energy: -109.0340266755937

Utilisez un autre solveur d'équations propres

Comme Fulqrum expose l'hamiltonien projeté sous une forme standard scipy.sparse.linalg.LinearOperator, vous pouvez le diagonaliser à l'aide de n'importe quel solveur d'eigenvaleurs de votre choix, sans être limité à scipy.sparse.linalg.eigsh. Cela permet de remplacer facilement un solveur par un autre mieux adapté à la taille de votre problème ou à votre matériel.

À titre d'exemple, ce guide présente PRIMME, une bibliothèque de solveurs d'équations propres hautement performante. Son interface Python, primme.eigsh, suit une convention d'appel similaire à celle de scipy.sparse.linalg.eigsh, ce qui vous permet de l'utiliser comme remplacement direct.

La seule différence par rapport au solveur précédent réside dans l'appel à primme.eigsh au lieu de scipy.sparse.linalg.eigsh.

import primme


# Define a fulqrum solver that uses PRIMME to diagonalize the projected Hamiltonian
def solve_sci_batch_primme(ci_batches, _h1, _h2, norb, nelec, tol=1e-12):
    # Ignore the _h1 and _h2 inputs because we store the Hamiltonian separately
    results = []
    for ints_a, ints_b in ci_batches:
        strs_a = [format(x, f"0{norb}b") for x in ints_a]
        strs_b = [format(x, f"0{norb}b") for x in ints_b]
        subspace = fq.Subspace([strs_a, strs_b])
        Hsub.update_subspace(subspace)
        Hsub_csr_linop = Hsub.to_csr_linearoperator_fast()
        diag = Hsub.diagonal_vector()
        v0 = np.zeros(len(subspace), dtype=Hsub.dtype)
        v0[np.argmin(diag)] = 1.0
        # primme.eigsh follows the same calling convention as scipy.sparse.linalg.eigsh,
        # except that it expects the initial guess v0 as a 2D array of column vectors.
        eigs, vecs = primme.eigsh(
            Hsub_csr_linop, k=1, which="SA", tol=tol, v0=v0.reshape(-1, 1)
        )
        vec = vecs[:, 0]
        occs = subspace.get_orbital_occupancies(
            probs=np.abs(vec) ** 2, norb=norb
        )
        # SQD ordering and Fulqrum ordering are transpose of each other.
        # ascontiguousarray() is needed for matching C contiguous ordering inside PySCF
        state = SCIState(
            amplitudes=np.ascontiguousarray(
                vec.reshape(len(strs_b), len(strs_a)).T
            ),
            ci_strs_a=ints_a,
            ci_strs_b=ints_b,
            norb=norb,
            nelec=nelec,
        )
        results.append(
            SCIResult(
                energy=eigs[0], sci_state=state, orbital_occupancies=occs
            )
        )
    return results


# 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=solve_sci_batch_primme,
    symmetrize_spin=True,
    seed=rng,
)

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

Output:

Exact energy: -109.04667177808028
Estimated energy: -109.03937429196165

Cette page a-t-elle été utile ?

Signaler un bogue, une coquille ou proposer du contenu sur GitHub.