Verbesserung einer SQD-Schätzung durch Orbitaloptimierung

Bei der samplebasierten Quantendiagonalisierung (SQD) wird die Grundzustandsenergie durch Diagonalisierung des Hamilton-Operators in einem festen Unterraum elektronischer Konfigurationen approximiert. Diese Schätzung hängt von der Orbitalbasis ab, in der der Hamiltonoperator ausgedrückt wird, und die Orbitaloptimierung (OO) nutzt diese Freiheit, um die Energie zu senken, ohne den Unterraum zu vergrößern.

In dieser Anleitung wird SQD auf einem „ $N_2$ “-Molekül angewendet und das Ergebnis anschließend durch eine Orbitaloptimierung verbessert, wobei verwendet wird, um ffsim den Hamilton-Operator darzustellen und die energieminimierende Orbitalrotation zu ermitteln.

SQD ausführen

Wir erstellen die Molekülintegrale für $N_2$ in der Molekülorbital (MO)-Basis, erzeugen gleichverteilte Zufallsstichproben und führen eine SQD-Berechnung durch, um eine Näherung für den Grundzustand zu erhalten.

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

Die Orbitalen optimieren

Bei der Orbitaloptimierung wird nach einer Orbitalrotation gesucht, die die Variationsenergie verringert

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

der SQD-Grundzustandsnäherung $|\psi\rangle$ . Eine Orbitalrotation wird durch eine $N \times N$ e unitäre Matrix $\mathbf{U}$ ( $N$ ist die Anzahl der räumlichen Orbitale) festgelegt, die über den Operator auf den Vielteilchenzustand einwirkt

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

ffsim.optimize_orbitals gibt die Matrix $\mathbf{U}$ zurück, und deren Anwendung auf die Orbitalbasis (über hamiltonian.rotated) entspricht der Anwendung von $\mathcal{U}$ auf den Zustand. Weitere Informationen finden Sie in der Erklärung zur Orbitalrotation bei ffsim .

Da sich durch die Rotation der Orbitalen der für den Unterraum sichtbare Hamiltonoperator ändert, wechseln wir zwischen zwei Schritten hin und her, bis sich die Energie nicht mehr verbessert:

Diagonalisieren Sie den Hamiltonoperator in der aktuellen Basis über die feste Menge der Konfigurationen.
Optimieren Sie die Orbitale, indem Sie die Rotation ermitteln, die die Energie des resultierenden Zustands minimiert, und transformieren Sie anschließend die Integrale in die neue Basis.

Den Schritt der Orbitalrotation delegieren wir an ffsim.optimize_orbitals, das die energieminimierende Rotation anhand der reduzierten Dichtematrizen (RDMs) für ein und zwei Teilchen des Zustands ermittelt. Siehe Abschnitt II A 4 für weitere Einzelheiten.

Warum die Orbitaloptimierung hier hilft

Die SCF-Molekülorbitalbasis (MO) ist bezüglich Orbitalrotationen für das vollständige -CI-Problem stationär. Allerdings arbeitet SQD in einem kleinen, abgeschnittenen Unterraum (hier einige hundert CI-Strings von insgesamt etwa 19 Millionen vollständigen CI-Determinanten), für den die MO-Basis im Allgemeinen nicht optimal ist, sodass eine Rotation der Orbitale die Energie senkt, die der Unterraum darstellen kann.

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.

Alternierende Diagonalisierung und Orbitaloptimierung

Wir halten den Diagonalisierungsunterraum auf die oben durch SQD ermittelten Konfigurationen festgelegt, sodass bei jeder Iteration der Effekt der Rotation der Orbitale isoliert wird. Jede Iteration:

Diagonalisiert den Hamilton-Operator über den festen Unterraum in der aktuellen Basis unter Verwendung des PySCF's -Selected-CI-Lösers.
Erstellt die RDMs des resultierenden Zustands – mehr ist nichtffsim.optimize_orbitals erforderlich.
Optimiert die Orbitalen : ffsim.optimize_orbitals Gibt die energieminimierende Drehung zurück, die wir auf die Integrale anwenden, um in die verbesserte Basis überzugehen.

Vor jedem Optimierungsschritt erfassen wir den Energiezustand. Da sich die Basis mit jeder Iteration verbessert, nimmt diese Folge monoton in Richtung der besten erreichbaren Energie im festgelegten Unterraum ab.

# 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

Vergleiche die Ergebnisse

Die Orbitaloptimierung verbessert die Schätzung im festen Unterraum und schließt einen Großteil der Lücke zur exakten Energie, bleibt dabei jedoch über dieser.

# 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

War diese Seite hilfreich?

Melden Sie auf GitHub einen Fehler, einen Tippfehler oder fordern Sie Inhalte an.