تحسين تقدير 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$$

لتقريب الحالة الأرضية $|\psi\rangle$ الناتج عن SQD. ويُحدَّد دوران المدارات بمصفوفة وحدوية $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، التي تجد الدوران الذي يقلل الطاقة انطلاقًا من مصفوفات الكثافة المختزلة (RDMs) أحادية وثنائية الجسم للحالة. راجع القسم 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. نبني RDMs للحالة الناتجة، وهو كل ما تحتاجه 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