שיפור אומדן SQD באמצעות אופטימיזציית אורביטלים

Sample-based quantum diagonalization (SQD) מקרב אנרגיית מצב יסוד על ידי האלכסון ההמילטוניאן בתת-מרחב קבוע של תצורות אלקטרוניות. האומדן הזה תלוי בבסיס האורביטלים שבו ההמילטוניאן מבוטא, ואופטימיזציית אורביטלים (OO) מנצלת את החופש הזה כדי להוריד את האנרגיה בלי להגדיל את תת-המרחב.

המדריך הזה מריץ SQD על מולקולת $N_2$ ואז משפר את התוצאה באמצעות אופטימיזציית אורביטלים, תוך שימוש ב-ffsim כדי לייצג את ההמילטוניאן ולמצוא את סיבוב האורביטלים שממזער את האנרגיה.

הרצת SQD

אנחנו בונים את האינטגרלים המולקולריים של $N_2$ בבסיס האורביטלים המולקולריים (MO), יוצרים דגימות אקראיות אחידות, ומריצים SQD כדי לקבל קירוב למצב היסוד.

# 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

אופטימיזציה של האורביטלים

אופטימיזציית אורביטלים מחפשת סיבוב אורביטלים שמוריד את האנרגיה הווריאציונית

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

של קירוב מצב היסוד של SQD‏ $|\psi\rangle$. סיבוב אורביטלים מוגדר על ידי מטריצה אוניטרית $N \times N$‏ $\mathbf{U}$ (כש-$N$ הוא מספר האורביטלים המרחביים), שפועלת על המצב מרובה-הגופים באמצעות האופרטור

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

ffsim.optimize_orbitals מחזיר את המטריצה $\mathbf{U}$, והפעלתה על בסיס האורביטלים (באמצעות hamiltonian.rotated) שקולה להפעלת $\mathcal{U}$ על המצב. לפרטים ראו את ההסבר של ffsim על סיבוב אורביטלים.

מכיוון שסיבוב האורביטלים משנה את ההמילטוניאן שתת-המרחב רואה, אנחנו מחליפים בין שני צעדים עד שהאנרגיה מפסיקה להשתפר:

1. מאלכסנים את ההמילטוניאן בבסיס הנוכחי מעל הקבוצה הקבועה של התצורות.

2. מבצעים אופטימיזציה לאורביטלים על ידי מציאת הסיבוב שממזער את האנרגיה של המצב שהתקבל, ואז מסובבים את האינטגרלים לבסיס החדש.

אנחנו מאצילים את שלב סיבוב האורביטלים ל-ffsim.optimize_orbitals, שמוצא את הסיבוב שממזער את האנרגיה מתוך מטריצות הצפיפות המצומצמות החד-גופיות והדו-גופיות (RDMs) של המצב. לפרטים ראו Sec. II A 4.

למה אופטימיזציית אורביטלים עוזרת כאן

בסיס האורביטלים המולקולריים (MO) של SCF נייח ביחס לסיבובי אורביטלים עבור בעיית ה-CI המלאה. אבל SQD עובד בתת-מרחב קטוע קטן (כאן כמה מאות מחרוזות CI מתוך כ-19 מיליון דטרמיננטות של CI מלא), ועבורו בסיס ה-MO בדרך כלל אינו אופטימלי, ולכן סיבוב האורביטלים מוריד את האנרגיה שתת-המרחב יכול לייצג.

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.

מחליפים בין האלכסון לאופטימיזציית אורביטלים

אנחנו משאירים את תת-מרחב ההאלכסון קבוע על התצורות ש-SQD גילה למעלה, כך שכל איטרציה מבודדת את ההשפעה של סיבוב האורביטלים. כל איטרציה:

1. מאלכסנת את ההמילטוניאן מעל תת-המרחב הקבוע בבסיס הנוכחי, באמצעות פותר ה-CI הנבחר של PySCF.

2. בונה את ה-RDMs של המצב שהתקבל, וזה כל מה ש-ffsim.optimize_orbitals צריך.

3. מבצעת אופטימיזציה לאורביטלים: ffsim.optimize_orbitals מחזיר את הסיבוב שממזער את האנרגיה, ואנחנו מפעילים אותו על האינטגרלים כדי לעבור לבסיס המשופר.

אנחנו מתעדים את האנרגיה לפני כל צעד אופטימיזציה. מכיוון שהבסיס משתפר בכל איטרציה, הסדרה הזו יורדת באופן מונוטוני לעבר האנרגיה הטובה ביותר שאפשר להשיג בתת-המרחב הקבוע.

# 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

השוואת התוצאות

אופטימיזציית אורביטלים משפרת את אומדן תת-המרחב הקבוע, וסוגרת חלק גדול מהפער לאנרגיה המדויקת תוך כדי הישארות מעליה.

# 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