Pagbutihin ang estimate ng SQD gamit ang orbital optimization

Tinatantya ng Sample-based quantum diagonalization (SQD) ang ground-state energy sa pamamagitan ng pag-diagonalize ng Hamiltonian sa isang nakapirming subspace ng mga electronic configuration. Ang estimate na iyon ay nakadepende sa orbital basis kung saan ipinapahayag ang Hamiltonian, at ginagamit ng orbital optimization (OO) ang kalayaang ito para pababain ang energy nang hindi pinalalaki ang subspace.

Pinapatakbo ng gabay na ito ang SQD sa isang molecule ng $N_2$ at pagkatapos ay pinapabuti ang resulta gamit ang orbital optimization, gamit ang ffsim para irepresenta ang Hamiltonian at hanapin ang orbital rotation na nagpapababa ng energy.

Patakbuhin ang SQD

Binubuo namin ang mga molecular integral para sa $N_2$ sa molecular-orbital (MO) basis, gumagawa ng mga uniform na random sample, at pinapatakbo ang SQD para makakuha ng ground-state approximation.

# 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

I-optimize ang mga orbital

Naghahanap ang orbital optimization ng orbital rotation na nagpapababa sa variational energy

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

ng SQD ground-state approximation na $|\psi\rangle$. Ang orbital rotation ay tinutukoy ng isang $N \times N$ na unitary matrix na $\mathbf{U}$ ($N$ ang bilang ng mga spatial orbital), na umaaksyon sa many-body state sa pamamagitan ng operator na

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

Ibinabalik ng ffsim.optimize_orbitals ang matrix na $\mathbf{U}$, at ang pag-apply nito sa orbital basis (sa pamamagitan ng hamiltonian.rotated) ay katumbas ng pag-apply ng $\mathcal{U}$ sa state. Tingnan ang paliwanag ng ffsim tungkol sa orbital-rotation para sa mga detalye.

Dahil binabago ng pag-rotate ng mga orbital ang Hamiltonian na nakikita ng subspace, salitan nating ginagawa ang dalawang hakbang hanggang tumigil ang pagbuti ng energy:

1. I-diagonalize ang Hamiltonian sa kasalukuyang basis sa ibabaw ng nakapirming hanay ng mga configuration.

2. I-optimize ang mga orbital sa pamamagitan ng paghahanap ng rotation na nagpapababa sa energy ng nagreresultang state, pagkatapos ay i-rotate ang mga integral sa bagong basis.

Ipinapaubaya namin ang hakbang ng orbital-rotation sa ffsim.optimize_orbitals, na nakakahanap ng rotation na nagpapababa ng energy mula sa one- at two-body reduced density matrix (RDM) ng state. Tingnan ang Sec. II A 4 para sa mga detalye.

Bakit nakakatulong dito ang orbital optimization

Ang SCF molecular-orbital (MO) basis ay stationary kaugnay ng mga orbital rotation sa problemang full-CI. Ngunit gumagana ang SQD sa isang maliit na pinutol na subspace (dito ay ilang daang CI string mula sa humigit-kumulang 19 milyong full-CI determinant), kung saan ang MO basis ay karaniwang hindi optimal, kaya ang pag-rotate ng mga orbital ay nagpapababa sa energy na kayang irepresenta ng subspace.

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.

Salitan na diagonalization at orbital optimization

Pinapanatili naming nakapirmi ang diagonalization subspace sa mga configuration na natuklasan ng SQD sa itaas, para ihiwalay ng bawat iteration ang epekto ng pag-rotate ng mga orbital. Ang bawat iteration ay:

1. Nagdi-diagonalize ng Hamiltonian sa ibabaw ng nakapirming subspace sa kasalukuyang basis, gamit ang selected-CI solver ng PySCF.

2. Bumubuo ng mga RDM ng nagreresultang state, na siya lamang kailangan ng ffsim.optimize_orbitals.

3. Nag-o-optimize ng mga orbital: ibinabalik ng ffsim.optimize_orbitals ang rotation na nagpapababa ng energy, na ina-apply namin sa mga integral para lumipat sa pinahusay na basis.

Itinatala namin ang energy bago ang bawat hakbang ng optimization. Dahil bumubuti ang basis sa bawat iteration, ang sequence na ito ay monotonic na bumababa patungo sa pinakamahusay na energy na makakamit sa nakapirming subspace.

# 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

Ihambing ang mga resulta

Pinapabuti ng orbital optimization ang estimate sa nakapirming subspace, na nagpapaliit sa malaking bahagi ng agwat sa tamang energy habang nananatiling mas mataas dito.

# 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