Tingkatkan estimasi SQD dengan optimasi orbital

Sample-based quantum diagonalization (SQD) memperkirakan energi ground-state dengan mendiagonalisasi Hamiltonian dalam subruang tetap dari konfigurasi elektronik. Estimasi itu bergantung pada basis orbital tempat Hamiltonian diekspresikan, dan optimasi orbital (OO) memanfaatkan kebebasan ini untuk menurunkan energi tanpa memperbesar subruang.

Panduan ini menjalankan SQD pada molekul $N_2$ lalu memperbaiki hasilnya dengan optimasi orbital, memakai ffsim untuk merepresentasikan Hamiltonian dan mencari rotasi orbital yang meminimalkan energi.

Jalankan SQD

Kita bangun integral molekuler untuk $N_2$ dalam basis molecular-orbital (MO), bangkitkan sampel acak uniform, lalu jalankan SQD untuk mendapatkan aproksimasi ground-state.

# 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

Optimalkan orbitalnya

Optimasi orbital mencari rotasi orbital yang menurunkan energi variasional

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

dari aproksimasi ground-state SQD $|\psi\rangle$. Sebuah rotasi orbital ditentukan oleh matriks unitary $N \times N$ $\mathbf{U}$ ($N$ adalah jumlah orbital spasial), yang bekerja pada state banyak-benda lewat operator

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

ffsim.optimize_orbitals mengembalikan matriks $\mathbf{U}$, dan menerapkannya pada basis orbital (lewat hamiltonian.rotated) setara dengan menerapkan $\mathcal{U}$ pada state. Lihat penjelasan rotasi orbital ffsim untuk detailnya.

Karena merotasi orbital mengubah Hamiltonian yang dilihat oleh subruang, kita bergantian melakukan dua langkah sampai energi tidak membaik lagi:

1. Diagonalisasi Hamiltonian dalam basis saat ini atas himpunan konfigurasi yang tetap.

2. Optimalkan orbital dengan mencari rotasi yang meminimalkan energi dari state yang dihasilkan, lalu rotasikan integral ke basis yang baru.

Kita serahkan langkah rotasi orbital ke ffsim.optimize_orbitals, yang mencari rotasi peminimal energi dari reduced density matrix (RDM) satu dan dua benda dari state tersebut. Lihat Sec. II A 4 untuk detailnya.

Mengapa optimasi orbital membantu di sini

Basis molecular-orbital (MO) SCF bersifat stasioner terhadap rotasi orbital untuk masalah CI penuh. Tapi SQD bekerja di subruang terpotong yang kecil (di sini beberapa ratus string CI dari sekitar 19 juta determinan full-CI), yang basis MO-nya umumnya tidak optimal, sehingga merotasi orbital menurunkan energi yang bisa direpresentasikan subruang itu.

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.

Bergantian antara diagonalisasi dan optimasi orbital

Kita jaga subruang diagonalisasi tetap fixed pada konfigurasi yang ditemukan SQD di atas, supaya setiap iterasi mengisolasi pengaruh rotasi orbital. Setiap iterasi:

1. Mendiagonalisasi Hamiltonian atas subruang tetap dalam basis saat ini, memakai solver selected-CI milik PySCF.

2. Membangun RDM dari state yang dihasilkan, yang merupakan satu-satunya yang dibutuhkan ffsim.optimize_orbitals.

3. Mengoptimalkan orbital: ffsim.optimize_orbitals mengembalikan rotasi peminimal energi, yang kita terapkan ke integral untuk berpindah ke basis yang lebih baik.

Kita catat energi sebelum setiap langkah optimasi. Karena basis membaik di setiap iterasi, urutan ini menurun secara monoton menuju energi terbaik yang bisa dicapai di subruang tetap.

# 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

Bandingkan hasilnya

Optimasi orbital memperbaiki estimasi subruang tetap, menutup sebagian besar selisih dengan energi eksak sambil tetap berada di atasnya.

# 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