Покращення оцінки SQD за допомогою оптимізації орбіталей

Квантова діагоналізація на основі вибірок (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, що знаходить поворот, який мінімізує енергію, за одно- та двочастинковими матрицями редукованої густини (RDM) стану. Докладніше дивись у розд. 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. Будує RDM отриманого стану, що єдине потрібно для 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