Améliorer une estimation SQD grâce à l'optimisation orbitale

La diagonalisation quantique par échantillonnage (SQD) permet d'estimer l'énergie de l'état fondamental en diagonalisant l'hamiltonien dans un sous-espace fixe de configurations électroniques. Cette estimation dépend de la base orbitale dans laquelle l'hamiltonien est exprimé, et l'optimisation orbitale (OO) exploite cette liberté pour réduire l'énergie sans élargir le sous-espace.

Ce guide explique comment exécuter SQD sur une molécule de type « $N_2$ », puis comment améliorer le résultat grâce à une optimisation orbitale, en utilisant ffsim pour représenter l'hamiltonien et déterminer la rotation orbitale minimisant l'énergie.

Exécuter SQD

Nous construisons les intégrales moléculaires pour l’ $N_2$ e dans la base des orbitales moléculaires (MO), générons des échantillons aléatoires uniformes, puis exécutons SQD afin d’obtenir une approximation de l’état fondamental.

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,
)

Output:

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}")

Output:

Exact energy:  -109.04667178
SQD energy:    -108.98469255

Optimiser les orbitales

L'optimisation orbitale consiste à rechercher une rotation orbitale qui minimise l'énergie variationnelle

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

de l'approximation de l'état fondamental SQD $|\psi\rangle$ . Une rotation orbitale est définie par une matrice unitaire d' $N \times N$ s $\mathbf{U}$ (où $N$ est le nombre d'orbitales spatiales), qui agit sur l'état à plusieurs corps par l'intermédiaire de l'opérateur

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

ffsim.optimize_orbitals renvoie la matrice d' $\mathbf{U}$ ; l'appliquer à la base orbitale (via hamiltonian.rotated) revient à appliquer $\mathcal{U}$ à l' état. Pour plus de détails, consultez l' explication de ffsim sur la rotation orbitale .

Étant donné que la rotation des orbitales modifie l'hamiltonien associé au sous-espace, nous alternons ces deux étapes jusqu'à ce que l'énergie cesse de s'améliorer :

Diagonaliser l'hamiltonien dans la base actuelle sur l'ensemble fixe de configurations.
Optimisez les orbitales en déterminant la rotation qui minimise l'énergie de l' état résultant, puis effectuez une rotation des intégrales dans la nouvelle base.

Nous confions l'étape de rotation orbitale à ffsim.optimize_orbitals, qui détermine la rotation minimisant l'énergie à partir des matrices de densité réduites (RDM) à un et à deux corps de l'état. Voir la section Voir II A 4 pour plus de détails.

Pourquoi l'optimisation orbitale s'avère utile dans ce cas précis

La base d'orbitales moléculaires (OM) SCF est stationnaire par rapport aux rotations orbitales pour le problème -CI complet . Mais la méthode SQD opère dans un petit sous-espace tronqué (ici, quelques centaines de chaînes CI sur environ 19 millions de déterminants CI complets), pour lequel la base MO n'est généralement pas optimale; la rotation des orbitales réduit donc l'énergie que ce sous-espace peut représenter.

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.

Diagonalisation alternée et optimisation orbitale

Nous conservons le sous-espace de diagonalisation défini par les configurations mises en évidence par la méthode SQD ci-dessus, de sorte que chaque itération isole l'effet de la rotation des orbitales. À chaque itération :

Diagonalise l'hamiltonien sur le sous-espace fixe dans la base actuelle, à l'aide du solveur CI sélectionné de PySCF's.
Crée les RDM de l'état résultant, ce quiffsim.optimize_orbitals suffit.
Optimise les orbitales : ffsim.optimize_orbitals renvoie la rotation minimisant l'énergie, que nous appliquons aux intégrales afin de passer à la base améliorée.

Nous enregistrons l'énergie avant chaque étape d'optimisation. Comme la base s'améliore à chaque itération, cette suite diminue de manière monotone vers la meilleure énergie pouvant être atteinte dans le sous-espace 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())

Output:

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

Comparez les résultats

L'optimisation orbitale améliore l'estimation du sous-espace fixe, comblant ainsi une grande partie de l'écart par rapport à l'énergie exacte tout en restant au-dessus de celle-ci.

# 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}")

Output:

Exact energy:      -109.04667178
SQD energy (MO):   -108.98469255
Energy after OO:   -109.01178727

Cette page a-t-elle été utile ?

Signaler un bogue, une coquille ou proposer du contenu sur GitHub.