Ulepszanie oszacowania SQD za pomocą optymalizacji orbitali

Sample-based quantum diagonalization (SQD) przybliża energię stanu podstawowego, diagonalizując hamiltonian w ustalonej podprzestrzeni konfiguracji elektronowych. To oszacowanie zależy od bazy orbitali, w której wyrażony jest hamiltonian, a optymalizacja orbitali (OO) wykorzystuje tę swobodę, by obniżyć energię bez powiększania podprzestrzeni.

Ten przewodnik uruchamia SQD dla cząsteczki $N_2$, a następnie poprawia wynik za pomocą optymalizacji orbitali, używając ffsim do reprezentacji hamiltonianu i znalezienia rotacji orbitali minimalizującej energię.

Uruchom SQD

Budujemy całki molekularne dla $N_2$ w bazie orbitali molekularnych (MO), generujemy jednorodne próbki losowe i uruchamiamy SQD, aby uzyskać przybliżenie stanu podstawowego.

# 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

Optymalizacja orbitali

Optymalizacja orbitali szuka rotacji orbitali, która obniża energię wariacyjną

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

przybliżenia stanu podstawowego SQD $|\psi\rangle$. Rotacja orbitali jest zadana przez macierz unitarną $N \times N$ $\mathbf{U}$ ($N$ to liczba orbitali przestrzennych), która działa na stan wielociałowy za pomocą operatora

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

ffsim.optimize_orbitals zwraca macierz $\mathbf{U}$, a zastosowanie jej do bazy orbitali (przez hamiltonian.rotated) jest równoważne zastosowaniu $\mathcal{U}$ do stanu. Szczegóły znajdziesz w wyjaśnieniu rotacji orbitali w ffsim.

Ponieważ rotacja orbitali zmienia hamiltonian widziany przez podprzestrzeń, na przemian wykonujemy dwa kroki, aż energia przestanie się poprawiać:

1. Zdiagonalizuj hamiltonian w bieżącej bazie na ustalonym zbiorze konfiguracji.

2. Zoptymalizuj orbitale, znajdując rotację minimalizującą energię otrzymanego stanu, a następnie obróć całki do nowej bazy.

Krok rotacji orbitali delegujemy do ffsim.optimize_orbitals, które znajduje rotację minimalizującą energię na podstawie jedno- i dwuciałowych zredukowanych macierzy gęstości (RDM) stanu. Szczegóły znajdziesz w Sec. II A 4.

Dlaczego optymalizacja orbitali tu pomaga

Baza orbitali molekularnych (MO) z SCF jest stacjonarna względem rotacji orbitali dla problemu pełnego CI. Jednak SQD działa w małej, obciętej podprzestrzeni (tutaj kilkaset ciągów CI spośród około 19 milionów wyznaczników pełnego CI), dla której baza MO zwykle nie jest optymalna, więc rotacja orbitali obniża energię, którą podprzestrzeń może reprezentować.

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.

Naprzemienna diagonalizacja i optymalizacja orbitali

Utrzymujemy podprzestrzeń diagonalizacji ustaloną na konfiguracjach odkrytych przez SQD powyżej, aby każda iteracja izolowała efekt rotacji orbitali. Każda iteracja:

1. Diagonalizuje hamiltonian w ustalonej podprzestrzeni w bieżącej bazie, używając solwera selected-CI z PySCF.

2. Buduje RDM otrzymanego stanu, co jest wszystkim, czego potrzebuje ffsim.optimize_orbitals.

3. Optymalizuje orbitale: ffsim.optimize_orbitals zwraca rotację minimalizującą energię, którą stosujemy do całek, aby przejść do ulepszonej bazy.

Zapisujemy energię przed każdym krokiem optymalizacji. Ponieważ baza poprawia się w każdej iteracji, ta sekwencja maleje monotonicznie w kierunku najlepszej energii osiągalnej w ustalonej podprzestrzeni.

# 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

Porównanie wyników

Optymalizacja orbitali poprawia oszacowanie dla ustalonej podprzestrzeni, zmniejszając znaczną część różnicy względem energii dokładnej, ale pozostając powyżej niej.

# 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