Améliorer une estimation SQD avec l'optimisation d'orbitales

La diagonalisation quantique basée sur l'échantillonnage (SQD) approxime une énergie de l'état fondamental en diagonalisant l'hamiltonien dans un sous-espace fixe de configurations électroniques. Cette estimation dépend de la base d'orbitales dans laquelle l'hamiltonien est exprimé, et l'optimisation d'orbitales (OO) exploite cette liberté pour abaisser l'énergie sans agrandir le sous-espace.

Ce guide exécute SQD sur une molécule $N_2$ puis améliore le résultat grâce à l'optimisation d'orbitales, en utilisant ffsim pour représenter l'hamiltonien et trouver la rotation d'orbitales qui minimise l'énergie.

Exécuter SQD

Nous construisons les intégrales moléculaires de $N_2$ dans la base des orbitales moléculaires (MO), générons des échantillons aléatoires uniformes et exécutons SQD pour obtenir une approximation de l'état fondamental.

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

_needed = {"ffsim": "ffsim", "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.cc
import pyscf.mcscf
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import 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
)

# Run SQD
result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=100,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=1,
    max_iterations=5,
    symmetrize_spin=True,
    seed=rng,
)

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

sqd_energy = result.energy + nuclear_repulsion_energy
print(f"Exact energy:  {exact_energy:.8f}")
print(f"SQD energy:    {sqd_energy:.8f}")

Exact energy:  -109.04667178
SQD energy:    -108.98469255

Optimiser les orbitales

L'optimisation d'orbitales recherche une rotation d'orbitales qui abaisse l'énergie variationnelle

$$E = \langle \psi | \mathcal{U}^\dagger\, H\, \mathcal{U} | \psi \rangle$$

de l'approximation de l'état fondamental SQD $|\psi\rangle$. Une rotation d'orbitales est spécifiée par une matrice unitaire $N \times N$ $\mathbf{U}$ ($N$ est le nombre d'orbitales spatiales), qui agit sur l'état à N corps via l'opérateur

$$\mathcal{U} = \exp\left[\sum_{pq, \sigma} \log(\mathbf{U})_{pq}\, a^\dagger_{p\sigma} a_{q\sigma}\right].$$

ffsim.optimize_orbitals renvoie la matrice $\mathbf{U}$, et l'appliquer à la base d'orbitales (via hamiltonian.rotated) équivaut à appliquer $\mathcal{U}$ à l'état. Consulte l' explication de la rotation d'orbitales dans ffsim pour plus de détails.

Comme la rotation des orbitales modifie l'hamiltonien vu par le sous-espace, nous alternons deux étapes jusqu'à ce que l'énergie cesse de s'améliorer :

1. Diagonaliser l'hamiltonien dans la base courante sur l'ensemble fixe de configurations.

2. Optimiser les orbitales en trouvant la rotation qui minimise l'énergie de l'état obtenu, puis faire tourner les intégrales vers la nouvelle base.

Nous déléguons l'étape de rotation d'orbitales à ffsim.optimize_orbitals, qui trouve la rotation minimisant l'énergie à partir des matrices densité réduites (RDM) à un et à deux corps de l'état. Voir Sec. II A 4 pour plus de détails.

Pourquoi l'optimisation d'orbitales aide ici

La base des orbitales moléculaires (MO) SCF est stationnaire par rapport aux rotations d'orbitales pour le problème CI complet. Mais SQD travaille dans un petit sous-espace tronqué (ici quelques centaines de chaînes CI sur environ 19 millions de déterminants CI complets), pour lequel la base MO n'est généralement pas optimale ; faire tourner les orbitales abaisse donc l'énergie que le sous-espace peut représenter.

import ffsim
from pyscf import fci

# ffsim's ``MolecularHamiltonian`` uses the same "chemist" ordering for the two-body
# tensor as PySCF's ``eri``, and stores the nuclear repulsion energy as the constant
# term so that expectation values come out as total energies.

Alterner diagonalisation et optimisation d'orbitales

Nous gardons le sous-espace de diagonalisation fixe sur les configurations découvertes par SQD ci-dessus, afin que chaque itération isole l'effet de la rotation des orbitales. Chaque itération :

1. Diagonalise l'hamiltonien sur le sous-espace fixe dans la base courante, en utilisant le solveur CI sélectionné de PySCF.

2. Construit les RDM de l'état obtenu, ce dont ffsim.optimize_orbitals a seul besoin.

3. Optimise les orbitales : ffsim.optimize_orbitals renvoie la rotation qui minimise l'énergie, que nous appliquons aux intégrales pour passer dans la base améliorée.

Nous enregistrons l'énergie avant chaque étape d'optimisation. Comme la base s'améliore à chaque itération, cette séquence décroît de façon monotone vers la meilleure énergie atteignable dans le sous-espace fixe.

# Fix the diagonalization subspace to the configurations found by SQD.
ci_strings = (result.sci_state.ci_strs_a, result.sci_state.ci_strs_b)
nelec = (num_elec_a, num_elec_b)

# Start from the MO basis in which we ran SQD.
hamiltonian_opt = ffsim.MolecularHamiltonian(
    hcore, eri, constant=nuclear_repulsion_energy
)

num_iters = 10
for i in range(num_iters):
    # Diagonalize over the fixed subspace in the current basis.
    myci = fci.selected_ci.SelectedCI()
    myci = fci.addons.fix_spin_(myci, ss=spin_sq)
    _, amplitudes = fci.selected_ci.kernel_fixed_space(
        myci,
        hamiltonian_opt.one_body_tensor,
        hamiltonian_opt.two_body_tensor,
        num_orbitals,
        nelec,
        ci_strs=ci_strings,
    )

    # Build the RDMs and record the energy before re-optimizing the orbitals.
    dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)
    rdm = ffsim.ReducedDensityMatrix(dm1, dm2)
    energy = rdm.expectation(hamiltonian_opt).real
    print(f"Iteration {i}: energy = {energy:.8f}")

    # Rotate the Hamiltonian into the energy-minimizing basis for the next iteration.
    # optimize_orbitals returns the unitary matrix U minimizing
    # rdm.rotated(U).expectation(hamiltonian), equivalently
    # rdm.expectation(hamiltonian.rotated(U.conj().T)), so we rotate by U^dagger.
    orbital_rotation = ffsim.optimize_orbitals(rdm, hamiltonian_opt)
    hamiltonian_opt = hamiltonian_opt.rotated(orbital_rotation.T.conj())

Iteration 0: energy = -108.98452447
Iteration 1: energy = -108.99981993
Iteration 2: energy = -109.00585329
Iteration 3: energy = -109.00816569
Iteration 4: energy = -109.00936616
Iteration 5: energy = -109.01014322
Iteration 6: energy = -109.01069439
Iteration 7: energy = -109.01109308
Iteration 8: energy = -109.01138928
Iteration 9: energy = -109.01161411

Comparer les résultats

L'optimisation d'orbitales améliore l'estimation sur sous-espace fixe, comblant une grande partie de l'écart avec l'énergie exacte tout en restant au-dessus.

# Diagonalize once more in the final optimized basis to report the improved energy.
myci = fci.selected_ci.SelectedCI()
myci = fci.addons.fix_spin_(myci, ss=spin_sq)
_, amplitudes = fci.selected_ci.kernel_fixed_space(
    myci,
    hamiltonian_opt.one_body_tensor,
    hamiltonian_opt.two_body_tensor,
    num_orbitals,
    nelec,
    ci_strs=ci_strings,
)
dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)
energy_after_oo = (
    ffsim.ReducedDensityMatrix(dm1, dm2).expectation(hamiltonian_opt).real
)

print(f"Exact energy:      {exact_energy:.8f}")
print(f"SQD energy (MO):   {sqd_energy:.8f}")
print(f"Energy after OO:   {energy_after_oo:.8f}")

Exact energy:      -109.04667178
SQD energy (MO):   -108.98469255
Energy after OO:   -109.01178727