Diagonalisation quantique basée sur l'échantillonnage d'un hamiltonien chimique

Estimation de l'utilisation : moins d'une minute sur un processeur Heron r2 (NOTE : Il s'agit uniquement d'une estimation. Votre durée d'exécution peut varier.)

Résultats d'apprentissage

À l'issue de ce tutoriel, les utilisateurs devraient avoir compris :

• Comment utiliser le module complémentaire SQD Qiskit pour estimer l'énergie de l'état fondamental d'un système moléculaire à l'aide de chaînes de bits échantillonnées à partir d'une unité de traitement quantique (QPU).
• Comment utiliser ffsim pour construire un circuit Jastrow à grappes unitaires locales (LUCJ) destiné à la simulation en chimie quantique.

Prérequis

Nous recommandons aux utilisateurs de se familiariser avec les sujets suivants avant de suivre ce tutoriel :

• Chimie quantique et seconde quantification
• Utilisation de la primitive « Sampler » pour échantillonner des circuits quantiques

Arrière-plan

Dans ce tutoriel, nous montrons comment traiter des échantillons quantiques bruités afin d'approximer l'état fondamental de la molécule d'azote $\text{N}_2$ à la longueur de liaison d'équilibre, en utilisant l 'extension SQD de Qiskit pour mettre en œuvre l 'algorithme de diagonalisation quantique par échantillonnage (SQD). Vous trouverez plus de détails sur ce logiciel dans la documentation correspondante, qui comprend notamment un exemple simple pour vous aider à démarrer.

Ce tutoriel est recommandé aux utilisateurs familiarisés avec la chimie quantique, et plus particulièrement à ceux qui savent déterminer l'énergie de l'état fondamental d'une molécule. Pour un guide détaillé sur le déroulement du processus, consultez le cours sur l'algorithme de diagonalisation quantique.

La SQD est une technique permettant de déterminer les valeurs propres et les vecteurs propres d'opérateurs quantiques, tels que l'hamiltonien d'un système quantique, en combinant le calcul quantique et le calcul classique distribué. Le calcul distribué classique sert à traiter les échantillons obtenus à partir d'un processeur quantique, ainsi qu'à projeter et à diagonaliser un hamiltonien cible dans un sous-espace qu'ils génèrent. Un processus basé sur le SQD comprend les étapes suivantes :

1. Choisissez un ansatz de circuit et appliquez-le sur un ordinateur quantique à un état de référence (dans ce cas, l'état Hartree-Fock ).
2. Echantillons de chaînes de bits de l'état quantique résultant.
3. Exécutez la procédure de récupération de configuration auto-cohérente sur les chaînes de bits afin d'obtenir l'approximation de l'état fondamental.

On sait que la SQD fonctionne bien lorsque l'état propre cible est peu dense : la fonction d'onde est supportée par un ensemble d'états de base $\mathcal{S} = \{|x\rangle \}$ dont la taille n'augmente pas de manière exponentielle avec la taille du problème.

Chimie quantique

L'hamiltonien d'un système moléculaire peut être écrit comme suit

où $h_{pr}$ et $h_{prqs}$ sont des nombres complexes appelés intégrales moléculaires qui peuvent être calculées à partir des spécifications de la molécule à l'aide d'un programme informatique. Dans ce tutoriel, nous calculons les intégrales à l'aide du logiciel PySCF logiciel.

Pour plus de détails sur la manière dont le hamiltonien moléculaire est dérivé, consultez un manuel de chimie quantique (par exemple, Modern Quantum Chemistry de Szabo et Ostlund). Pour une explication de haut niveau de la manière dont les problèmes de chimie quantique sont transposés sur les ordinateurs quantiques, consultez la conférence Mapping Problems to Qubits de l'université d'été mondiale Qiskit 2024.

Approche du cluster unitaire local de Jastrow (LUCJ)

La méthode SQD nécessite un modèle de circuit quantique à partir duquel prélever des échantillons. Dans ce tutoriel, nous utiliserons l'approche LUCJ (Local Unitary Cluster Jastrow) en raison de sa combinaison entre justification physique et facilité de mise en œuvre matérielle. Nous utiliserons ffsim pour construire le circuit de référence.

L'approche LUCJ s'adapte aux QPU dont la connectivité des qubits est limitée. Les orbitales de spin sont mappées sur des qubits de telle sorte que l'ansatz ne nécessite pas de routage à l'aide de portes SWAP. IBM® Le matériel présente une topologie de qubits en réseau hexagonal dense; dans ce cas, nous pouvons adopter un motif en « zigzag », illustré ci-dessous. Dans ce schéma, les orbitales de même spin sont associées à des qubits selon une topologie en ligne (cercles rouges et bleus), et une connexion entre des orbitales de spin différent est présente tous les quatre orbitales spatiales, cette connexion étant assurée par un qubit auxiliaire (cercles violets).

Récupération de configuration auto-cohérente

La procédure de récupération de la configuration autoconsistante est conçue pour extraire autant de signaux que possible d'échantillons quantiques bruyants. Comme l'hamiltonien moléculaire conserve le nombre de particules et le spin Z, il est logique de choisir un ansatz de circuit qui conserve également ces symétries. Lorsqu'il est appliqué à l'état de Hartree-Fock, l'état résultant a un nombre de particules et un spin Z fixes dans un environnement sans bruit. Par conséquent, les moitiés de spin $\alpha$ et de spin $\beta$ de toute chaîne de bits échantillonnée à partir de cet état devraient avoir le même poids de Hamming que dans l'état Hartree-Fock. En raison de la présence de bruit dans les processeurs quantiques actuels, certaines chaînes de bits mesurées ne respectent pas cette propriété. Une forme simple de postsélection permettrait d'écarter ces chaînes de bits, mais c'est un gaspillage car ces chaînes de bits peuvent encore contenir un certain signal. La procédure de récupération autoconsistante tente de récupérer une partie de ce signal lors du post-traitement. La procédure est itérative et nécessite en entrée une estimation de l'occupation moyenne de chaque orbitale dans l'état fondamental, qui est d'abord calculée à partir des échantillons bruts. La procédure est exécutée en boucle et chaque itération comporte les étapes suivantes :

1. Pour chaque chaîne de bits qui ne respecte pas les symétries spécifiées, les bits sont retournés selon une procédure probabiliste conçue pour rapprocher la chaîne de bits de l'estimation actuelle des occupations orbitales moyennes, afin d'obtenir une nouvelle chaîne de bits.
2. Rassembler toutes les chaînes de bits anciennes et nouvelles qui satisfont aux symétries et sous-échantillonner des sous-ensembles d'une taille fixe, choisie à l'avance.
3. Pour chaque sous-ensemble de chaînes de bits, projetez l'hamiltonien dans le sous-espace couvert par les vecteurs de base correspondants (voir la section précédente pour une description de ces vecteurs de base) et calculez une estimation de l'état fondamental de l'hamiltonien projeté sur un ordinateur classique.
4. Mettre à jour l'estimation de l'occupation moyenne des orbitales avec l'estimation de l'état fondamental ayant l'énergie la plus basse.

Diagramme du flux de travail SQD

Le flux de travail de la SQD est décrit dans le diagramme suivant :

Exigences

Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :

• Qiskit SDK v1.0 ou plus tard, avec prise en charge de la visualisation
• Qiskit Runtime v0.22 ou plus tard (pip install qiskit-ibm-runtime)
• Module complémentaire SQD Qiskit v0.11 ou version ultérieure (pip install qiskit-addon-sqd)
• ffsim v0.0.75 ou version ultérieure (pip install ffsim)

Configuration

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

Exemple de simulateur à petite échelle

Dans ce tutoriel, nous allons déterminer une approximation de l'état fondamental d'une molécule d'azote à une distance de liaison proche de son équilibre. Nous utilisons d'abord un petit ensemble de bases d' STO-6G s afin de pouvoir simuler l'expérience et vérifier qu'elle fonctionne.

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

Tout d'abord, nous définissons la molécule et ses propriétés.

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

Avant de construire le circuit de l'ansatz LUCJ, nous effectuons d'abord un calcul CCSD dans la cellule de code suivante. Les amplitudes $t_1$ et $t_2$ de ce calcul seront utilisées pour initialiser les paramètres de l'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

Nous utilisons maintenant ffsim pour créer le circuit de référence. Comme notre molécule présente un état de Hartree-Fock à couche fermée, nous utilisons la variante à spin équilibré de l'approche UCJ, UCJOpSpinBalanced. Nous avons défini optimize=True dans la from_t_amplitudes méthode afin d'activer la double factorisation « compressée » des amplitudes de l' $t_2$ e (pour plus de détails, voir la section « The local unitary cluster Jastrow (LUCJ) ansatz » dans la documentation de ffsim).

Étant donné que l'ansatz LUCJ s'adapte à la connectivité disponible du QPU, nous devons initialiser le backend du QPU avant de créer l'ansatz. Pour l'instant, nous allons créer un backend générique doté d'une carte à couplage hexagonal fort et d'un ensemble de portes auquel l'approche LUCJ se décompose naturellement. Ensuite, nous utiliserons ffsim.qiskit.generate_lucj_pass_manager pour créer un gestionnaire de passes spécialisé dans la transcompilation de l'approche LUCJ vers le backend spécifié, conformément à la structure en « zigzag » décrite dans la section consacrée au contexte de l'approche LUCJ. Cette fonction utilise une heuristique de notation pour minimiser les erreurs associées à la configuration sélectionnée, ce qui est important si votre backend est un véritable QPU ou un simulateur doté d'un modèle de bruit. Outre le gestionnaire de passes, cette fonction renvoie également les paires de couplage alpha-bêta pouvant être mises en œuvre sur le matériel. Si toutes les paires ne peuvent pas être mises en œuvre, un avertissement s'affiche.

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

Étape 2 : Optimisation pour l'exécution sur du matériel quantique

Nous optimisons ensuite le circuit pour le matériel cible. En général, cette étape consiste à initialiser le backend matériel et un gestionnaire de passes pour ce backend. Cependant, comme l'approche LUCJ est adaptée à la connectivité matérielle, nous avons déjà effectué ces opérations à l'étape précédente. Il ne reste plus qu'à exécuter le gestionnaire de passes sur le circuit pour le transcompiler en un circuit ISA pouvant être exécuté directement sur le 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})

Étape 3 : Exécutez à l'aide d' Qiskit primitives

Après avoir optimisé le circuit en vue de son exécution sur le matériel, nous sommes prêts à l'exécuter sur le matériel cible et à collecter des échantillons afin d'estimer l'énergie de l'état fondamental. Comme nous ne disposons que d'un seul circuit, nous allons utiliser le mode d'exécution « Job » du service de calcul d' IBM Quantum pour exécuter notre circuit.

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]

Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité

Un indicateur utile pour évaluer la qualité du résultat du QPU est le nombre de configurations valides renvoyées. Une configuration valide possède le nombre correct de particules et le spin Z correct, ce qui signifie que la moitié droite de la chaîne binaire a un poids de Hamming égal au nombre d'électrons à spin up, et que la moitié gauche a un poids de Hamming égal au nombre d'électrons à spin down. La cellule suivante calcule la fraction des configurations échantillonnées qui sont valides.

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

Toutes les chaînes de bits sont valides, car nous effectuons un échantillonnage du circuit sur un simulateur sans bruit. Lorsqu'on exécute le programme sur un QPU bruyant, cette fraction sera inférieure à un, mais on espère qu'elle sera supérieure à celle à laquelle on s'attendrait si les chaînes de bits avaient été échantillonnées de manière aléatoire et uniforme, ce qui est calculé dans la cellule suivante.

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

Nous estimons maintenant l'énergie de l'état fondamental de l'hamiltonien à l'aide de la fonction diagonalize_fermionic_hamiltonian . Cette fonction exécute la procédure de récupération de la configuration autoconsistante pour affiner de manière itérative les échantillons quantiques bruyants afin d'améliorer l'estimation de l'énergie. Nous passons une fonction de rappel afin de pouvoir enregistrer les résultats intermédiaires pour une analyse ultérieure. Voir la documentation de l'API pour des explications sur les arguments de diagonalize_fermionic_hamiltonian.

Ici, nous utilisons initial_occupancies l'argument pour diagonalize_fermionic_hamiltonian spécifier la configuration Hartree-Fock comme estimation initiale pour les occupations orbitales dans l'état fondamental. Cette approche est judicieuse pour les systèmes dont l'état fondamental repose largement sur la configuration Hartree-Fock, mais elle peut ne pas convenir dans d'autres situations, même si des méthodes de calcul plus avancées peuvent permettre d'obtenir de meilleures estimations initiales dans ces cas-là. La spécification initial_occupancies permet également d'exécuter la récupération de configuration même si aucune configuration valide n'a été échantillonnée, comme cela peut être le cas lors de l'échantillonnage d'un grand circuit sur un QPU bruité. Sans cet argument, la récupération de la configuration échouerait et générerait une erreur si aucune configuration valide n'était fournie.

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

Visualisez les résultats

Le premier graphique montre que, dans cette simulation, nous sommes déjà très 1 mH proches de la réponse exacte après la première itération (on considère généralement que la précision chimique est de l'ordre de 1 kcal/mol$\approx$1.6 mH). Il s'agit toutefois d'un petit système et, comme les échantillons sont exempts de bruit, il n'est pas nécessaire de procéder à une reconstruction de la configuration. Sur un système plus important fonctionnant avec un QPU bruyant, plusieurs itérations de récupération de la configuration peuvent s'avérer nécessaires, et la précision finale peut s'en trouver réduite. En général, on peut améliorer l'énergie en autorisant davantage d'itérations de récupération de configuration ou en augmentant le nombre d'échantillons par lot.

Le deuxième graphique montre l'occupation moyenne de chaque orbite spatiale après la dernière itération. Nous pouvons constater que les électrons de spin-up et de spin-down occupent les cinq premières orbitales avec une forte probabilité dans nos solutions.

# 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:

Exemple de matériel à grande échelle

Nous allons maintenant exécuter un exemple plus complexe sur du matériel quantique réel. Nous allons ici dériver un espace actif pour la molécule d'azote à partir de la base de données « cc-pVDZ ».

Étapes 1 à 4

Nous regroupons ici toutes ces étapes au sein d'un processus unique à plus grande échelle, qui est ensuite exécuté sur du matériel quantique réel.

# ------------------------------ 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

Etapes suivantes