Îmbunătățește o estimare SQD cu optimizarea orbitalelor

Diagonalizarea cuantică bazată pe eșantioane (SQD) aproximează o energie a stării fundamentale prin diagonalizarea Hamiltonianului într-un subspațiu fix de configurații electronice. Această estimare depinde de baza orbitală în care este exprimat Hamiltonianul, iar optimizarea orbitalelor (OO) exploatează această libertate pentru a reduce energia fără a mări subspațiul.

Acest ghid rulează SQD pe o moleculă de $N_2$ și apoi îmbunătățește rezultatul cu optimizarea orbitalelor, folosind ffsim pentru a reprezenta Hamiltonianul și a găsi rotația orbitală care minimizează energia.

Rulează SQD

Construim integralele moleculare pentru $N_2$ în baza orbitalelor moleculare (MO), generăm eșantioane aleatorii uniforme și rulăm SQD pentru a obține o aproximare a stării fundamentale.

# 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

Optimizează orbitalele

Optimizarea orbitalelor caută o rotație orbitală care reduce energia variațională

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

a aproximării stării fundamentale SQD $|\psi\rangle$. O rotație orbitală este specificată de o matrice unitară $N \times N$ $\mathbf{U}$ ($N$ este numărul de orbitale spațiale), care acționează asupra stării cu mai multe corpuri prin operatorul

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

ffsim.optimize_orbitals returnează matricea $\mathbf{U}$, iar aplicarea ei asupra bazei orbitale (prin hamiltonian.rotated) este echivalentă cu aplicarea lui $\mathcal{U}$ asupra stării. Vezi explicația ffsim despre rotația orbitală pentru detalii.

Deoarece rotirea orbitalelor schimbă Hamiltonianul văzut de subspațiu, alternăm doi pași până când energia nu mai se îmbunătățește:

1. Diagonalizează Hamiltonianul în baza curentă peste setul fix de configurații.

2. Optimizează orbitalele găsind rotația care minimizează energia stării rezultate, apoi rotește integralele în noua bază.

Delegăm pasul de rotație orbitală către ffsim.optimize_orbitals, care găsește rotația ce minimizează energia pornind de la matricile de densitate reduse (RDM) cu unul și două corpuri ale stării. Vezi Sec. II A 4 pentru detalii.

De ce ajută aici optimizarea orbitalelor

Baza orbitalelor moleculare (MO) SCF este staționară față de rotațiile orbitale pentru problema CI completă. Dar SQD lucrează într-un subspațiu mic trunchiat (aici câteva sute de șiruri CI din aproximativ 19 milioane de determinanți CI completi), pentru care baza MO nu este în general optimă, deci rotirea orbitalelor reduce energia pe care o poate reprezenta subspațiul.

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.

Alternarea diagonalizării și a optimizării orbitalelor

Menținem subspațiul de diagonalizare fix la configurațiile descoperite de SQD mai sus, astfel încât fiecare iterație să izoleze efectul rotirii orbitalelor. Fiecare iterație:

1. Diagonalizează Hamiltonianul peste subspațiul fix în baza curentă, folosind solverul CI selectat din PySCF.

2. Construiește RDM-urile stării rezultate, ceea ce este tot ce are nevoie ffsim.optimize_orbitals.

3. Optimizează orbitalele: ffsim.optimize_orbitals returnează rotația care minimizează energia, pe care o aplicăm integralelor pentru a trece în baza îmbunătățită.

Înregistrăm energia înainte de fiecare pas de optimizare. Deoarece baza se îmbunătățește la fiecare iterație, această secvență scade monoton spre cea mai bună energie realizabilă în subspațiul fix.

# 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

Compară rezultatele

Optimizarea orbitalelor îmbunătățește estimarea în subspațiul fix, reducând mare parte din diferența față de energia exactă, rămânând totodată peste ea.

# 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