Diagonalizzazione quantistica di un Hamiltoniano chimico basata su campioni

Stima di utilizzo: meno di un minuto su un processore Heron r2 (NOTA: questa è solo una stima. Il tempo di esecuzione potrebbe variare)

Risultati di apprendimento

Dopo aver seguito questo tutorial, gli utenti dovrebbero aver compreso:

• Come utilizzare l 'add-on SQD Qiskit per approssimare l'energia dello stato fondamentale di un sistema molecolare utilizzando stringhe di bit campionate da un'unità di elaborazione quantistica (QPU).
• Come utilizzare ffsim per costruire un circuito Jastrow a cluster unitario locale (LUCJ) per la simulazione di chimica quantistica.

Prerequisiti

Si consiglia agli utenti di familiarizzarsi con i seguenti argomenti prima di seguire questo tutorial:

• Chimica quantistica e seconda quantizzazione
• Utilizzo della primitiva Sampler per campionare da circuiti quantistici

Sfondo

In questo tutorial mostriamo come eseguire la post-elaborazione di campioni quantistici soggetti a rumore per approssimare lo stato fondamentale della molecola di azoto $\text{N}_2$ alla lunghezza di legame di equilibrio, utilizzando l 'add-on SQD di Qiskit per implementare l 'algoritmo di diagonalizzazione quantistica basata su campioni (SQD). Maggiori dettagli sul software sono disponibili nella documentazione corrispondente, che include anche un semplice esempio per iniziare.

Questo tutorial è consigliato agli utenti che hanno familiarità con la chimica quantistica: in particolare, è richiesta una certa dimestichezza con il calcolo delle energie dello stato fondamentale di una molecola. Per una guida dettagliata al flusso di lavoro, consulta il corso sull'algoritmo di diagonalizzazione quantistica.

L'SQD è una tecnica che consente di determinare gli autovalori e gli autovettori di operatori quantistici, come l'hamiltoniano di un sistema quantistico, ricorrendo congiuntamente al calcolo quantistico e al calcolo classico distribuito. Il calcolo distribuito classico viene utilizzato per elaborare i campioni ottenuti da un processore quantistico e per proiettare e diagonalizzare un hamiltoniano di riferimento in un sottospazio da essi generato. Un flusso di lavoro basato su SQD prevede le seguenti fasi:

1. Scegliere un'ipotesi di circuito e applicarla su un computer quantistico a uno stato di riferimento (in questo caso, lo stato Hartree-Fock ).
2. Esempi di bitstings dallo stato quantico risultante.
3. Eseguire la procedura di recupero della configurazione auto-consistente sulle stringhe di bit per ottenere l'approssimazione dello stato fondamentale.

È noto che la SQD funziona bene quando l'autostato bersaglio è rado: la funzione d'onda è supportata in un insieme di stati base $\mathcal{S} = \{|x\rangle \}$ la cui dimensione non aumenta esponenzialmente con la dimensione del problema.

Chimica quantistica

L'hamiltoniana di un sistema molecolare può essere scritta come

dove $h_{pr}$ e $h_{prqs}$ sono numeri complessi chiamati integrali molecolari che possono essere calcolati dalle specifiche della molecola utilizzando un programma informatico. In questa esercitazione, calcoliamo gli integrali utilizzando il pacchetto software PySCF pacchetto software.

Per i dettagli su come viene ricavata l'hamiltoniana molecolare, consultare un testo di chimica quantistica (per esempio, Modern Quantum Chemistry di Szabo e Ostlund). Per una spiegazione di alto livello di come i problemi di chimica quantistica vengono mappati sui computer quantistici, si veda la lezione Mapping Problems to Qubits della Qiskit Global Summer School 2024.

Approccio del cluster unitario locale di Jastrow (LUCJ)

SQD richiede un ansatz del circuito quantistico da cui estrarre i campioni. In questo tutorial utilizzeremo l'approccio del cluster unitario locale di Jastrow (LUCJ) per la sua combinazione di fondamenti fisici e facilità di implementazione hardware. Useremo ffsim per costruire il circuito di ansatz.

L'approccio LUCJ si adatta alle QPU con connettività limitata dei qubit. Gli orbitali di spin vengono mappati sui qubit in modo tale che l'ansatz non richieda l'utilizzo di porte SWAP. IBM® L'hardware presenta una topologia dei qubit a reticolo esagonale denso; in tal caso, possiamo adottare uno schema "a zig-zag", illustrato di seguito. In questo schema, gli orbitali con lo stesso spin sono associati a qubit con una topologia lineare (cerchi rossi e blu), mentre tra gli orbitali con spin diverso è presente una connessione ogni quattro orbitali spaziali, facilitata da un qubit ancilla (cerchi viola).

Ripristino della configurazione auto-coerente

La procedura di recupero della configurazione autoconsistente è progettata per estrarre il maggior segnale possibile da campioni quantistici rumorosi. Poiché l'hamiltoniana molecolare conserva il numero di particelle e lo spin Z, ha senso scegliere un'ipotesi di circuito che conservi anche queste simmetrie. Quando viene applicato allo stato Hartree-Fock, lo stato risultante ha un numero di particelle e uno spin Z fissi nell'impostazione senza rumore. Pertanto, le metà di spin $\alpha$ e di spin $\beta$ di qualsiasi stringa di bit campionata da questo stato dovrebbero avere lo stesso peso di Hamming dello stato Hartree-Fock. A causa della presenza di rumore negli attuali processori quantistici, alcune stringhe di bit misurate violeranno questa proprietà. Una semplice forma di post-selezione consentirebbe di scartare queste stringhe di bit, ma è uno spreco perché queste stringhe di bit potrebbero contenere ancora qualche segnale. La procedura di recupero autoconsistente tenta di recuperare parte del segnale in post-elaborazione. La procedura è iterativa e richiede come input una stima delle occupazioni medie di ciascun orbitale nello stato fondamentale, che viene prima calcolata dai campioni grezzi. La procedura viene eseguita in un ciclo e ogni iterazione prevede le seguenti fasi:

1. Per ogni stringa di bit che viola le simmetrie specificate, si capovolgono i suoi bit con una procedura probabilistica volta ad avvicinare la stringa di bit alla stima corrente delle occupazioni orbitali medie, per ottenere una nuova stringa di bit.
2. Raccogliere tutte le bitstringhe vecchie e nuove che soddisfano le simmetrie e sotto-campionare sottoinsiemi di dimensioni fisse, scelte in anticipo.
3. Per ogni sottoinsieme di stringhe di bit, proiettare l'hamiltoniana nel sottospazio coperto dai vettori base corrispondenti (si veda la sezione precedente per una descrizione di questi vettori base) e calcolare una stima dello stato fondamentale dell'hamiltoniana proiettata su un computer classico.
4. Aggiornare la stima delle occupazioni orbitali medie con la stima dello stato fondamentale con l'energia più bassa.

Diagramma del flusso di lavoro SQD

Il flusso di lavoro SQD è rappresentato nel seguente diagramma:

Requisiti

Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:

• Qiskit SDK v1.0 o versioni successive, con supporto alla visualizzazione
• Qiskit Runtime v0.22 o successivamente (pip install qiskit-ibm-runtime)
• Componente aggiuntivo SQD Qiskit v0.11 o versioni successive (pip install qiskit-addon-sqd)
• ffsim v0.0.75 o versioni successive (pip install ffsim)

Configura

import math

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.primitives import StatevectorSampler
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Esempio di simulatore su piccola scala

In questo tutorial cercheremo di trovare un'approssimazione dello stato fondamentale di una molecola di azoto in prossimità della sua distanza di legame di equilibrio. Per prima cosa utilizziamo un piccolo set di basi " STO-6G " in modo da poter simulare l'esperimento e assicurarci che funzioni.

Fase 1: mappare gli input classici su un problema quantistico

Per prima cosa, identifichiamo la molecola e le sue proprietà.

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="sto-6g",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
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), norb)

# Compute exact energy using FCI
reference_energy = cas.run().e_tot

print(f"norb = {norb}")
print(f"nelec = {nelec}")

Output:

converged SCF energy = -108.464957764796
CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000
norb = 8
nelec = (5, 5)

Prima di costruire il circuito di ansatz LUCJ, eseguiamo un calcolo CCSD nella seguente cella di codice. Le ampiezze $t_1$ e $t_2$ di questo calcolo saranno utilizzate per inizializzare i parametri dell'ansatz.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

Output:

E(CCSD) = -108.5933309085008  E_corr = -0.1283731437052354

Ora utilizziamo ffsim per creare il circuito di ansatz. Poiché la nostra molecola presenta uno stato di Hartree-Fock a guscio chiuso, utilizziamo la variante a spin bilanciato dell’ansatz UCJ, UCJOpSpinBalanced. Abbiamo impostato optimize=True nel metodo from_t_amplitudes per abilitare la doppia fattorizzazione "compressa" delle ampiezze di $t_2$ (per ulteriori dettagli, consultare la documentazione di ffsim relativa all'approccio LUCJ (Local Unitary Cluster Jastrow)).

Poiché l'ansatz LUCJ si adatta alla connettività disponibile della QPU, è necessario inizializzare il backend della QPU prima di creare l'ansatz. Per ora, creeremo un backend generico con una mappa ad accoppiamento esagonale forte e un insieme di gate in cui l'ansatz LUCJ si scompone naturalmente. Quindi, useremo ffsim.qiskit.generate_lucj_pass_manager per creare un gestore di passaggi specializzato nella trasposizione dell'approccio LUCJ al backend specificato, secondo lo schema "a zig-zag" descritto nella sezione introduttiva dedicata all'approccio LUCJ. Questa funzione utilizza un algoritmo euristico di valutazione per ridurre al minimo gli errori associati al layout selezionato, il che è importante se il backend è una vera QPU o un simulatore con un modello di rumore. Oltre a restituire il gestore dei pass, questa funzione restituisce anche le coppie di accoppiamento alfa-beta che possono essere implementate sull'hardware. Se non è possibile implementare tutte le coppie, viene emesso un avviso.

import warnings

from qiskit.transpiler import CouplingMap

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
coupling_map = CouplingMap.from_heavy_hex(3)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
    # Limit the number of optimization iterations to prevent the code cell
    # from running too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it
# to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

Fase 2: Ottimizzazione per l'esecuzione su hardware quantistico

Successivamente, ottimizziamo il circuito per un hardware specifico. In genere, questa fase prevede l'inizializzazione del backend hardware e di un gestore di passaggi per quel backend. Tuttavia, poiché l'approccio LUCJ è adattato alla connettività hardware, abbiamo già eseguito queste operazioni nella fase precedente. Non resta che eseguire il gestore di passaggi sul circuito per traspilarlo in un circuito ISA che possa essere eseguito direttamente sulla QPU.

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

Output:

Gate counts: OrderedDict({'xx_plus_yy': 86, 'p': 16, 'measure': 16, 'cp': 15, 'x': 10, 'swap': 2, 'barrier': 1})

Passaggio 3: eseguire utilizzando Qiskit primitives

Dopo aver ottimizzato il circuito per l'esecuzione su hardware, siamo pronti a eseguirlo sull'hardware di destinazione e a raccogliere campioni per la stima dell'energia dello stato fondamentale. Poiché disponiamo di un solo circuito, utilizzeremo la modalità di esecuzione “ IBM Quantum ” di Compute Service per eseguire il nostro circuito.

rng = np.random.default_rng()
sampler = StatevectorSampler(seed=rng)
job = sampler.run([isa_circuit], shots=100_000)

Output:

Warning: Trying to add QuantumRegister to a QuantumCircuit having a layout

primitive_result = job.result()
pub_result = primitive_result[0]

Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato

Un parametro utile per valutare la qualità dell'output della QPU è il numero di configurazioni valide restituite. Una configurazione valida ha il numero corretto di particelle e lo spin Z, il che significa che la metà destra della stringa di bit ha un peso di Hamming pari al numero di elettroni con spin up, mentre la metà sinistra ha un peso di Hamming pari al numero di elettroni con spin down. La cella seguente calcola la frazione di configurazioni campionate che sono valide.

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )


bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")

Output:

Fraction of sampled configurations that are valid: 1.0

Tutte le stringhe di bit sono valide perché stiamo effettuando il campionamento del circuito su un simulatore privo di rumore. Se l'operazione viene eseguita su una QPU soggetta a rumore, la frazione sarà inferiore a uno, ma si spera che sia maggiore di quella che ci si aspetterebbe se le stringhe di bit fossero campionate in modo casuale e uniforme, come calcolato nella cella seguente.

expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)

Output:

Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

Ora stimiamo l'energia di stato fondamentale dell'hamiltoniana utilizzando la funzione diagonalize_fermionic_hamiltonian . Questa funzione esegue la procedura di recupero della configurazione autoconsistente per affinare iterativamente i campioni quantistici rumorosi e migliorare la stima dell'energia. Si passa una funzione di callback in modo da poter salvare i risultati intermedi per un'analisi successiva. Vedere la documentazione dell'API per le spiegazioni degli argomenti di diagonalize_fermionic_hamiltonian.

Qui, utilizziamo initial_occupancies l'argomento per diagonalize_fermionic_hamiltonian specificare la configurazione di Hartree-Fock come ipotesi iniziale per le occupazioni orbitali nello stato fondamentale. Questo approccio è sensato per i sistemi in cui lo stato fondamentale ha un supporto significativo sulla configurazione di Hartree-Fock, ma potrebbe non essere appropriato in altre situazioni, anche se metodi computazionali più avanzati potrebbero fornire ipotesi iniziali migliori in questi casi. Specificando è initial_occupancies anche possibile eseguire il ripristino della configurazione anche se non sono state campionate configurazioni valide, come può accadere quando si campiona un circuito di grandi dimensioni su una QPU rumorosa. Senza questo argomento, il ripristino della configurazione non andrebbe a buon fine e genererebbe un errore se non fossero state fornite configurazioni valide.

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []


def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )


result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

Output:

Iteration 1
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Iteration 2
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Final energy: -108.59275573641656
Final energy error: 0.0032316145694579745

Visualizza i risultati

Il primo grafico mostra che in questa simulazione siamo già vicini 1 mH alla risposta esatta dopo la prima iterazione (l'accuratezza chimica è generalmente considerata pari a 1 kcal/mol$\approx$1.6 mH). Si tratta però di un sistema di piccole dimensioni e, poiché i campioni sono privi di rumore, non è necessario recuperare la configurazione. Su un sistema più grande in esecuzione su una QPU soggetta a rumore, potrebbero essere necessarie più iterazioni di recupero della configurazione e la precisione finale potrebbe risultare inferiore. In generale, è possibile migliorare l'energia consentendo un maggior numero di iterazioni di recupero della configurazione oppure aumentando il numero di campioni per lotto.

Il secondo grafico mostra l'occupazione media di ciascun orbitale spaziale dopo l'iterazione finale. Possiamo notare che sia gli elettroni con spin-up che quelli con spin-down occupano i primi cinque orbitali con alta probabilità nelle nostre soluzioni.

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

Output:

Esempio di hardware su larga scala

Ora eseguiamo un esempio più ampio su un hardware quantistico reale. In questa sede, ricaveremo uno spazio attivo per la molecola di azoto a partire dal set di basi " cc-pVDZ ".

Passaggi da 1 a 4

Qui riuniamo tutte le fasi in un unico flusso di lavoro su scala più ampia, che viene poi eseguito su hardware quantistico reale.

# ------------------------------ Step 1 ------------------------------
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="cc-pvdz",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
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), norb)

# Store reference energy from SCI calculation performed separately
reference_energy = -109.22802921665716

print(f"norb = {norb}")
print(f"nelec = {nelec}")

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)
print(f"Using backend {backend.name}")

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
    # Limit the number of optimization iterations to prevent the code cell
    # from running too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it
# to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()


# ------------------------------ Step 2 ------------------------------

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")


# ------------------------------ Step 3 ------------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_SQD"]
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]


# ------------------------------ Step 4 ------------------------------

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")
expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)
# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the
# orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the
# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []


result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

Output:

converged SCF energy = -108.929838385609
norb = 26
nelec = (5, 5)
E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450045
Using backend ibm_boston

Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].
Removing interaction (24, 24) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].
Removing interaction (20, 20) from the end.

Gate counts: OrderedDict({'sx': 7039, 'rz': 6990, 'cz': 1858, 'x': 61, 'measure': 52, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.02124
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.13889134249762
		Subspace dimension: 120409
	Subsample 1
		Energy: -109.11785470455858
		Subspace dimension: 110889
	Subsample 2
		Energy: -109.13234360554011
		Subspace dimension: 130321
Iteration 2
	Subsample 0
		Energy: -109.16392179579177
		Subspace dimension: 223729
	Subsample 1
		Energy: -109.16281938332986
		Subspace dimension: 223729
	Subsample 2
		Energy: -109.16955816711932
		Subspace dimension: 233289
Iteration 3
	Subsample 0
		Energy: -109.17905772999075
		Subspace dimension: 324900
	Subsample 1
		Energy: -109.17532445048462
		Subspace dimension: 357604
	Subsample 2
		Energy: -109.1733168689756
		Subspace dimension: 348100
Iteration 4
	Subsample 0
		Energy: -109.18437778820451
		Subspace dimension: 474721
	Subsample 1
		Energy: -109.18450164209159
		Subspace dimension: 476100
	Subsample 2
		Energy: -109.18493571190754
		Subspace dimension: 487204
Iteration 5
	Subsample 0
		Energy: -109.18616522497996
		Subspace dimension: 622521
	Subsample 1
		Energy: -109.18652868888333
		Subspace dimension: 644809
	Subsample 2
		Energy: -109.18753326484406
		Subspace dimension: 585225
Final energy: -109.18753326484406
Final energy error: 0.040495951813099396

Passi successivi