Migliora una stima SQD con l'ottimizzazione degli orbitali

La Sample-based quantum diagonalization (SQD) approssima l'energia dello stato fondamentale diagonalizzando l'Hamiltoniana in un sottospazio fisso di configurazioni elettroniche. Tale stima dipende dalla base orbitale in cui è espressa l'Hamiltoniana, e l'ottimizzazione degli orbitali (OO) sfrutta questa libertà per abbassare l'energia senza ampliare il sottospazio.

Questa guida esegue SQD su una molecola di $N_2$ e poi migliora il risultato con l'ottimizzazione degli orbitali, usando ffsim per rappresentare l'Hamiltoniana e trovare la rotazione degli orbitali che minimizza l'energia.

Esegui SQD

Costruiamo gli integrali molecolari per $N_2$ nella base degli orbitali molecolari (MO), generiamo campioni casuali uniformi ed eseguiamo SQD per ottenere un'approssimazione dello stato fondamentale.

# 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

Ottimizza gli orbitali

L'ottimizzazione degli orbitali cerca una rotazione degli orbitali che abbassi l'energia variazionale

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

dell'approssimazione SQD dello stato fondamentale $|\psi\rangle$. Una rotazione degli orbitali è specificata da una matrice unitaria $N \times N$ $\mathbf{U}$ ($N$ è il numero di orbitali spaziali), che agisce sullo stato a molti corpi tramite l'operatore

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

ffsim.optimize_orbitals restituisce la matrice $\mathbf{U}$, e applicarla alla base orbitale (tramite hamiltonian.rotated) equivale ad applicare $\mathcal{U}$ allo stato. Per i dettagli, vedi la spiegazione della rotazione degli orbitali di ffsim.

Poiché ruotare gli orbitali cambia l'Hamiltoniana vista dal sottospazio, alterniamo due passi finché l'energia non smette di migliorare:

1. Diagonalizza l'Hamiltoniana nella base corrente sull'insieme fisso di configurazioni.

2. Ottimizza gli orbitali trovando la rotazione che minimizza l'energia dello stato risultante, poi ruota gli integrali nella nuova base.

Deleghiamo il passo di rotazione degli orbitali a ffsim.optimize_orbitals, che trova la rotazione che minimizza l'energia a partire dalle matrici di densità ridotte (RDM) a uno e due corpi dello stato. Per i dettagli, vedi Sec. II A 4.

Perché qui l'ottimizzazione degli orbitali aiuta

La base degli orbitali molecolari (MO) SCF è stazionaria rispetto alle rotazioni degli orbitali per il problema CI completo. Ma SQD lavora in un piccolo sottospazio troncato (qui poche centinaia di stringhe CI su circa 19 milioni di determinanti CI completi), per il quale la base MO in generale non è ottimale, quindi ruotare gli orbitali abbassa l'energia che il sottospazio può rappresentare.

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.

Alternare diagonalizzazione e ottimizzazione degli orbitali

Manteniamo il sottospazio di diagonalizzazione fisso sulle configurazioni scoperte da SQD in precedenza, in modo che ogni iterazione isoli l'effetto della rotazione degli orbitali. Ogni iterazione:

1. Diagonalizza l'Hamiltoniana sul sottospazio fisso nella base corrente, usando il solver CI selezionato di PySCF.

2. Costruisce le RDM dello stato risultante, che è tutto ciò di cui ha bisogno ffsim.optimize_orbitals.

3. Ottimizza gli orbitali: ffsim.optimize_orbitals restituisce la rotazione che minimizza l'energia, che applichiamo agli integrali per passare alla base migliorata.

Registriamo l'energia prima di ogni passo di ottimizzazione. Poiché la base migliora a ogni iterazione, questa sequenza decresce monotonicamente verso la migliore energia ottenibile nel sottospazio fisso.

# 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

Confronta i risultati

L'ottimizzazione degli orbitali migliora la stima nel sottospazio fisso, colmando gran parte del divario con l'energia esatta pur rimanendo al di sopra di essa.

# 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