Omezení dimenze podprostoru

V tomto návodu ukážeme vliv dimenze podprostoru na techniku samokonzistentní obnovy konfigurací.

A priori nevíme, jaká je správná dimenze podprostoru pro dosažení cílové úrovně přesnosti. Víme však, že zvětšování dimenze podprostoru zvyšuje přesnost metody. Můžeme proto zkoumat přesnost předpovědí jako funkci dimenze podprostoru. Zadej molekulu a její vlastnosti.

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

_needed = {"matplotlib": "matplotlib", "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 warnings

import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
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

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

Vygeneruj několik náhodných bitových řetězců jako zástupce vzorků z QPU.

import numpy as np
from qiskit_addon_sqd.counts import generate_bit_array_uniform

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

Zavolej SQD se zvětšujícími se velikostmi dávek.

from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian

list_samples_per_batch = [50, 200, 400, 600]

# SQD options
max_iterations = 5

# Eigenstate solver options
num_batches = 10
max_davidson_cycles = 200

energies = []
subspace_dimensions = []

for samples_per_batch in list_samples_per_batch:
    result = diagonalize_fermionic_hamiltonian(
        hcore,
        eri,
        bit_array,
        samples_per_batch=samples_per_batch,
        norb=num_orbitals,
        nelec=(num_elec_a, num_elec_b),
        num_batches=num_batches,
        max_iterations=max_iterations,
        symmetrize_spin=True,
        seed=rng,
    )
    energies.append(result.energy)
    subspace_dimensions.append(np.prod(result.sci_state.amplitudes.shape))

Tento graf ukazuje, že zvětšování dimenze podprostoru vede k přesnějším výsledkům.

import matplotlib.pyplot as plt

# Data for energies plot
x1 = subspace_dimensions
y1 = np.array(energies) + nuclear_repulsion_energy

fig, axs = plt.subplots(1, 1, figsize=(12, 6))

# Plot energies
axs.plot(x1, y1, marker=".", markersize=20, label="Estimated")
axs.set_xticks(x1)
axs.set_xticklabels(x1)
axs.axhline(y=exact_energy, color="red", linestyle="--", label="Exact")
axs.set_title("Approximated Ground State Energy vs subspace dimension")
axs.set_xlabel("Subspace dimension")
axs.set_ylabel("Energy (Ha)")
axs.legend()

plt.tight_layout()
plt.show()

Nastavení explicitní hranice pomocí max_dim

Ve výše uvedeném průchodu byla dimenze podprostoru řízena jen nepřímo. samples_per_batch určuje, kolik bitových řetězců se vybere do každé dávky, ale dimenze výsledného podprostoru je dána tím, kolik různých řetězců spinu alfa a spinu beta tyto bitové řetězce obsahují, a množina konfigurací dostupných pro podvýběr roste z jedné iterace obnovy konfigurací k další. Proto jsou dimenze podprostoru na vodorovné ose výše nepravidelná čísla, která jsme nevybírali.

Argument max_dim funkce diagonalize_fermionic_hamiltonian omezuje dimenzi přímo. Omezuje počet CI řetězců zachovaných v každém spinovém sektoru, a protože podprostor je generován kartézským součinem obou sektorů, jeho dimenze je součinem dimenzí obou sektorů:

• max_dim=d omezuje oba sektory na d řetězců, takže dimenze podprostoru je nejvýše d ** 2.

• max_dim=(d_a, d_b) omezuje sektory spinu alfa a spinu beta zvlášť, takže dimenze podprostoru je nejvýše d_a * d_b. Pokud je symmetrize_spin=True, musí být obě meze stejné.

• max_dim=None, výchozí hodnota, žádné omezení neukládá.

Pokud dávka obsahuje více řetězců, než limit povoluje, nejméně důležité se zahodí: nejprve se zachovají řetězce předané přes include_configurations, poté řetězce přenesené z předchozí iterace (viz carryover_threshold) a zbývající místa se zaplní vzorkovanými řetězci v sestupném pořadí podle toho, jak často byly vzorkovány.

Protože max_dim omezuje velikost největší diagonalizace, kterou bude řešič vlastních stavů muset provést, je to pohodlný způsob, jak udržet řešič v rámci rozpočtu času běhu a paměti: zvol samples_per_batch velkoryse a nech max_dim rozhodnout, jak velký se podprostor smí stát.

Níže studii zopakujeme, tentokrát s pevným samples_per_batch a rostoucím max_dim.

list_max_dim = [100, 200, 400, 800]

# Use a generous number of samples per batch, so that max_dim is what limits
# the dimension of the subspace
samples_per_batch = 600

energies_max_dim = []
subspace_dimensions_max_dim = []

for max_dim in list_max_dim:
    result = diagonalize_fermionic_hamiltonian(
        hcore,
        eri,
        bit_array,
        samples_per_batch=samples_per_batch,
        norb=num_orbitals,
        nelec=(num_elec_a, num_elec_b),
        num_batches=num_batches,
        max_iterations=max_iterations,
        symmetrize_spin=True,
        max_dim=max_dim,
        seed=rng,
    )
    energies_max_dim.append(result.energy)
    subspace_dimensions_max_dim.append(
        int(np.prod(result.sci_state.amplitudes.shape))
    )

print(f"max_dim:            {list_max_dim}")
print(f"subspace dimension: {subspace_dimensions_max_dim}")

max_dim:            [100, 200, 400, 800]
subspace dimension: [10000, 40000, 160000, 640000]

Dimenze podprostoru je nyní v každém běhu přesně max_dim ** 2: byla zvolena, nikoli objevena. Vynesení obou průchodů do jednoho grafu ukazuje, že přesnost určuje samotná dimenze podprostoru, bez ohledu na to, jaký argument vedl k jejímu dosažení. Obě křivky se ale nemusí přesně shodovat. Běhy s max_dim berou 600 bitových řetězců na dávku a poté ponechají nejčastěji vzorkované řetězce, takže při dané dimenzi mohou zachovat lepší množinu konfigurací než běh, který této dimenze dosáhl tím, že vzorkoval méně bitových řetězců.

fig, axs = plt.subplots(1, 1, figsize=(12, 6))

axs.plot(
    subspace_dimensions,
    np.array(energies) + nuclear_repulsion_energy,
    marker=".",
    markersize=20,
    label="Varying samples_per_batch",
)
axs.plot(
    subspace_dimensions_max_dim,
    np.array(energies_max_dim) + nuclear_repulsion_energy,
    marker="^",
    markersize=12,
    label="Varying max_dim",
)
axs.axhline(y=exact_energy, color="red", linestyle="--", label="Exact")
axs.set_xscale("log")
axs.ticklabel_format(axis="y", useOffset=False)
axs.set_title("Approximated Ground State Energy vs subspace dimension")
axs.set_xlabel("Subspace dimension")
axs.set_ylabel("Energy (Ha)")
axs.legend()

plt.tight_layout()
plt.show()