Perbaiki anggaran SQD dengan pengoptimuman orbital

Sample-based quantum diagonalization (SQD) menganggar tenaga keadaan asas dengan mempepenjurukan Hamiltonian dalam subruang tetap konfigurasi elektronik. Anggaran itu bergantung pada asas orbital yang digunakan untuk menyatakan Hamiltonian, dan pengoptimuman orbital (OO) memanfaatkan kebebasan ini untuk merendahkan tenaga tanpa membesarkan subruang.

Panduan ini menjalankan SQD pada molekul $N_2$ dan kemudian memperbaiki hasilnya dengan pengoptimuman orbital, menggunakan ffsim untuk mewakili Hamiltonian dan mencari putaran orbital yang meminimumkan tenaga.

Jalankan SQD

Kami membina kamiran molekul untuk $N_2$ dalam asas orbital molekul (MO), menjana sampel rawak seragam, dan menjalankan SQD untuk memperoleh penghampiran keadaan asas.

# 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

Optimumkan orbital

Pengoptimuman orbital mencari putaran orbital yang merendahkan tenaga variasi

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

bagi penghampiran keadaan asas SQD $|\psi\rangle$. Putaran orbital ditentukan oleh matriks uniter $N \times N$ $\mathbf{U}$ ($N$ ialah bilangan orbital ruang), yang bertindak pada keadaan banyak jasad melalui operator

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

ffsim.optimize_orbitals memulangkan matriks $\mathbf{U}$, dan menggunakannya pada asas orbital (melalui hamiltonian.rotated) setara dengan menggunakan $\mathcal{U}$ pada keadaan. Lihat penjelasan putaran orbital ffsim untuk butiran.

Oleh sebab memutar orbital mengubah Hamiltonian yang dilihat oleh subruang, kami berselang-seli dua langkah sehingga tenaga berhenti bertambah baik:

1. Pepenjuruan Hamiltonian dalam asas semasa ke atas set konfigurasi yang tetap.

2. Optimumkan orbital dengan mencari putaran yang meminimumkan tenaga bagi keadaan yang terhasil, kemudian putarkan kamiran ke dalam asas baharu.

Kami menyerahkan langkah putaran orbital kepada ffsim.optimize_orbitals, yang mencari putaran peminimum tenaga daripada matriks ketumpatan terturun (RDM) satu dan dua jasad bagi keadaan itu. Lihat Sek. II A 4 untuk butiran.

Mengapa pengoptimuman orbital membantu di sini

Asas orbital molekul (MO) SCF bersifat pegun terhadap putaran orbital bagi masalah CI penuh. Tetapi SQD bekerja dalam subruang terpangkas yang kecil (di sini beberapa ratus rentetan CI daripada kira-kira 19 juta penentu CI penuh), yang mana asas MO secara amnya tidak optimum, jadi memutar orbital merendahkan tenaga yang boleh diwakili oleh subruang.

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.

Berselang-seli antara pepenjuruan dan pengoptimuman orbital

Kami mengekalkan subruang pepenjuruan tetap pada konfigurasi yang ditemui oleh SQD di atas, supaya setiap lelaran mengasingkan kesan memutar orbital. Setiap lelaran:

1. Mempepenjurukan Hamiltonian ke atas subruang tetap dalam asas semasa, menggunakan penyelesai CI terpilih PySCF.

2. Membina RDM bagi keadaan yang terhasil, iaitu semua yang diperlukan oleh ffsim.optimize_orbitals.

3. Mengoptimumkan orbital: ffsim.optimize_orbitals memulangkan putaran peminimum tenaga, yang kami gunakan pada kamiran untuk beralih ke asas yang diperbaiki.

Kami merekodkan tenaga sebelum setiap langkah pengoptimuman. Oleh sebab asas bertambah baik pada setiap lelaran, jujukan ini menurun secara monoton menuju tenaga terbaik yang boleh dicapai dalam 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 keputusan

Pengoptimuman orbital memperbaiki anggaran subruang tetap, merapatkan sebahagian besar jurang kepada tenaga tepat sambil kekal 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