Eine SQD-Schätzung mit Orbitaloptimierung verbessern

Die Sample-based quantum diagonalization (SQD) approximiert eine Grundzustandsenergie, indem sie den Hamiltonian in einem festen Unterraum elektronischer Konfigurationen diagonalisiert. Diese Schätzung hängt von der Orbitalbasis ab, in der der Hamiltonian ausgedrückt wird, und die Orbitaloptimierung (OO) nutzt diese Freiheit, um die Energie zu senken, ohne den Unterraum zu vergrößern.

Diese Anleitung führt SQD an einem $N_2$-Molekül aus und verbessert das Ergebnis dann mit Orbitaloptimierung, wobei ffsim verwendet wird, um den Hamiltonian darzustellen und die energieminimierende Orbitalrotation zu finden.

SQD ausführen

Wir erstellen die molekularen Integrale für $N_2$ in der Molekülorbital-Basis (MO), erzeugen gleichverteilte Zufalls-Samples und führen SQD aus, um eine Grundzustandsnäherung zu erhalten.

# 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

Die Orbitale optimieren

Die Orbitaloptimierung sucht nach einer Orbitalrotation, die die variationelle Energie

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

der SQD-Grundzustandsnäherung $|\psi\rangle$ senkt. Eine Orbitalrotation wird durch eine unitäre $N \times N$-Matrix $\mathbf{U}$ angegeben ($N$ ist die Anzahl der räumlichen Orbitale), die über den Operator auf den Vielteilchenzustand wirkt

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

ffsim.optimize_orbitals gibt die Matrix $\mathbf{U}$ zurück, und ihre Anwendung auf die Orbitalbasis (über hamiltonian.rotated) ist gleichbedeutend mit der Anwendung von $\mathcal{U}$ auf den Zustand. Details findest du in der ffsim-Erklärung zur Orbitalrotation.

Da das Rotieren der Orbitale den Hamiltonian verändert, den der Unterraum sieht, wechseln wir zwei Schritte ab, bis sich die Energie nicht mehr verbessert:

1. Diagonalisiere den Hamiltonian in der aktuellen Basis über der festen Menge von Konfigurationen.

2. Optimiere die Orbitale, indem du die Rotation findest, die die Energie des resultierenden Zustands minimiert, und rotiere dann die Integrale in die neue Basis.

Den Schritt der Orbitalrotation überlassen wir ffsim.optimize_orbitals, welches die energieminimierende Rotation aus den Ein- und Zweikörper-reduzierten Dichtematrizen (RDMs) des Zustands findet. Details siehe Abschn. II A 4.

Warum Orbitaloptimierung hier hilft

Die SCF-Molekülorbital-Basis (MO) ist stationär gegenüber Orbitalrotationen für das volle CI-Problem. SQD arbeitet jedoch in einem kleinen, abgeschnittenen Unterraum (hier einige hundert CI-Strings von etwa 19 Millionen Full-CI-Determinanten), für den die MO-Basis im Allgemeinen nicht optimal ist, sodass das Rotieren der Orbitale die Energie senkt, die der Unterraum darstellen kann.

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.

Diagonalisierung und Orbitaloptimierung abwechseln

Wir halten den Diagonalisierungs-Unterraum fest bei den Konfigurationen, die SQD oben gefunden hat, sodass jede Iteration die Wirkung der Orbitalrotation isoliert. Jede Iteration:

1. Diagonalisiert den Hamiltonian über dem festen Unterraum in der aktuellen Basis mit PySCFs Selected-CI-Solver.

2. Erstellt die RDMs des resultierenden Zustands, was alles ist, was ffsim.optimize_orbitals benötigt.

3. Optimiert die Orbitale: ffsim.optimize_orbitals gibt die energieminimierende Rotation zurück, die wir auf die Integrale anwenden, um in die verbesserte Basis zu wechseln.

Wir erfassen die Energie vor jedem Optimierungsschritt. Da sich die Basis in jeder Iteration verbessert, nimmt diese Folge monoton ab und nähert sich der besten Energie, die im festen Unterraum erreichbar ist.

# 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

Die Ergebnisse vergleichen

Die Orbitaloptimierung verbessert die Schätzung im festen Unterraum und schließt einen Großteil der Lücke zur exakten Energie, bleibt dabei aber oberhalb von ihr.

# 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