Générer des circuits d'évolution temporelle de Krylov (SKQD)

La diagonalisation quantique de Krylov basée sur des échantillons ( SKQD ) est une variante de la diagonalisation par sous-espaces axée sur la mécanique quantique. Plutôt que d’échantillonner une seule approche variationnelle (comme dans la diagonalisation quantique par échantillonnage simple ( SQD )) ou un ensemble de circuits d’évolution temporelle aléatoires (comme dans l’approche « SqDRIFT »), elle échantillonne une base de Krylov : la famille d’états

obtenu en faisant évoluer un état de référence $\lvert \psi_0 \rangle$ par multiples croissants d'un pas de temps fixe $\Delta t$. Les chaînes de bits échantillonnées à partir de ces circuits $D$ peuplent le sous-espace de Krylov dans lequel l'hamiltonien est ensuite diagonalisé de manière classique.

Ce guide reproduit les étapes de construction du circuit de l'algorithme SKQD pour le modèle d'Anderson à impureté unique (SIAM), qui décrit une impureté magnétique couplée à un bain sans interaction, conformément à la publication sur le SKQD.

1. Construire l'hamiltonien SIAM sous la forme d'un opérateur fermionique

Le modèle d'Anderson à impureté unique place une orbitale d'impureté en interaction (indexée par l' $0$ ) en contact avec un bain sans interaction constitué de sites d' $L$ s disposés en chaîne. Son hamiltonien se décompose en une partie à un corps $H_1$ et une partie à deux corps $H_2$,

avec le terme à un seul corps représentant l'énergie locale de l'impureté (le potentiel chimique $\mu$ ), l'hybridation impureté-bain $V$ et le saut de l'impureté vers le bain de ses voisins les plus proches $t$,

et la partie à deux corps qui fait porter l' $U$ de la répulsion de Coulomb sur place uniquement à l'impureté,

Ce guide s'applique au point de symétrie particule-trou $\mu = -U/2$. Ces deux parties correspondent directement aux intégrales à un corps ( h1e sauts dans le bain, hybridation de l'impureté, énergie locale de l'impureté) et à la seule intégrale à deux corps h2e (le terme coulombien de l'impureté). Commencez par les assembler dans la position de référence, où le modèle est défini naturellement.

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. Passer à la base de l'élan

La partie « bain » du SIAM est une chaîne de fermions libres; elle est donc diagonalizable par un changement de base vers la base des impulsions. C'est une étape décisive pour SKQD. Comme le bain est invariant par translation, l'état fondamental n'est pas clairsemé dans la base des positions, mais il l'est dans la base des impulsions, où l'hamiltonien du bain est diagonal. La rotation orbitale d'une particule unique qui opère cette transformation relève de l'algèbre linéaire propre à ce problème. Elle diagonalise la matrice de « bath hopping » tout en laissant intacte l'orbitale de l'impureté, puis repositionne l'impureté sur un site central :

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

Contrairement à un workflow variationnel, cette méthode effectue une rotation des intégrales hamiltoniennes dans la base des impulsions, de sorte que l'ensemble du circuit (hamiltonien et état confondus) s'inscrit dans une base où l'état fondamental est clairsemé et où les chaînes de bits échantillonnées sont donc informatives. La rotation des intégrales à un et à deux corps correspond à un changement de base standard des tenseurs de structure électronique :

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

Comme la rotation laisse l'orbitale d'impureté fixe (elle n'agit que sur le bain), l'interaction locale reste inchangée. h2e_momentum reste un terme composé d'un seul nombre, tout comme h2e. Le seul changement réside dans le fait que l'impureté a été déplacée vers un indice central. Cela se voit directement dans les tenseurs intégraux, qui comportent chacun un élément non nul :

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

C'est cette localité qui permet de synthétiser l'interaction à moindre coût par la suite.

Regroupez les intégrales à un corps dans un FermionOperator à l'aide du constructeur d'intégrales électroniques from_1body_tril_spin_sym(), qui les attend regroupées selon un ordre triangulaire inférieur. La partie à deux corps est un terme de type « nombre-nombre »; par conséquent, plutôt que de la traiter via le mécanisme d'intégrale à deux corps, construisez-la directement à l'aide de cre()/ann() : l'opérateur de Coulomb sur site $U \, n_{p\uparrow} n_{p\downarrow}$ pour le mode d'impureté $p$ (son mode $\uparrow$ à l'indice p, son mode $\downarrow$ à l'indice dans p + norb la configuration à spin en bloc). Les deux constructeurs transforment les orbitales spatiales norb en modes fermioniques 2 * norb de spin double.

Conservez l' $H_1$ (partie à un corps) et l' $H_2$ (partie à deux corps) en tant qu'opérateurs distincts plutôt que de les additionner immédiatement, en faisant correspondre les termes $H_1$ et $H_2$ des équations section-1. En créant chaque élément une seule fois (dans les deux bases), on peut le réutiliser par la suite. Une fois additionnés, ils donnent l'hamiltonien complet permettant la diagonalisation exacte dans l'une ou l'autre base. Plus précisément, l'opérateur « $H_2$ » (base de la quantité de mouvement) est l'opérateur à partir duquel évolue l'étape de Trotter à l'étape 4.

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

Le changement de base est une transformation de similitude unitaire; il ne modifie donc pas le spectre. Vérifiez que l'opérateur est correct en comparant sa valeur propre la plus faible dans le secteur « $(\text{norb}, \text{nelec})$ » à celle obtenue par une diagonalisation exacte. La fonction expose FermionOperator une interface SciPy ( LinearOperator soutenue par un noyau matriciel-vectoriel FCI (Full Configuration Interaction) natif), ce qui permet de la transmettre directement à scipy.sparse.linalg.eigsh().

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

L'avantage de ce changement de base réside dans la parcimonie. Pour illustrer cela, comptez le nombre de déterminants de base de calcul nécessaires pour rendre compte de 99 % du poids de l'état fondamental, puis comparez la base d'impulsion à la base de position. L'énergie ne dépend pas des bases, contrairement au nombre de déterminants :

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

Les deux opérateurs partagent la même énergie, mais l'état fondamental dans la base de position est réparti sur des milliers de déterminants, tandis que la base d'impulsion le concentre sur une vingtaine d'entre eux; c'est pourquoi l'échantillonnage SKQD fonctionne dans la base d'impulsion et non dans la base de position.

3. Préparer l'état de référence

La référence SKQD $\lvert \psi_0 \rangle$ correspond à la superposition de toutes les excitations des électrons les plus proches du niveau de Fermi dans les modes d'impulsion vides les plus proches. La publication SKQD la construit à partir d'un réseau de « XXPlusYYGates » écrits à la main. Il semble fortement corrélé (il comporte des centaines d'amplitudes non nulles dans la base de calcul), mais un réseau de rotations bimodales conservant le nombre est une opération à une seule particule (gaussienne fermionique); l'état qu'il produit est donc un déterminant de Slater unique. Cela signifie que vous pouvez le préparer à l'aide d'une seule porte PrepareSlaterDeterminant , compte tenu de la rotation à particule unique mise en œuvre par le réseau.

Étant donné que chacune d'entre elles XXPlusYYGate agit comme une rotation de Givens sur une paire de modes, cette rotation peut être construite directement dans le cadre de la théorie à particule unique (une matrice de Givens par porte) au lieu de générer manuellement les portes à deux qubits. Les grilles forment un maillage de rotations « plus proche voisin » qui s'épanouit à partir du niveau de Fermi en couches de plus en plus larges, répartissant chaque électron proche du niveau de Fermi entre les modes vides situés juste au-dessus de lui :

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

La préparation du déterminant de Slater s'effectue secteur de spin par secteur, car chaque secteur possède sa propre forme (électrons, orbitales). Une seule porte PrepareSlaterDeterminant (une occupation, les modes nocc les plus bas comblés, ainsi que la rotation d’une seule particule) prépare un secteur; cette même porte est appliquée une fois aux modes alpha et 0..norb une fois aux modes bêta norb..2*norb, puisque les deux secteurs de spin partagent ici la même occupation de référence et la même rotation :

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

Cette porte déclarative remplace le réseau XXPlusYYGate codé manuellement; le transpileur génère la synthèse par rotation de Givens.

4. Assembler les circuits de Krylov

Chaque circuit de Krylov consiste en la préparation d'un état de référence suivie d'une évolution sous l' $e^{-i k \, \Delta t \, H}$. On décompose l'hamiltonien en ses composantes à un corps et à deux corps, $H = H_1 + H_2$, puis on applique la méthode de Trotter à l'évolution obtenue à partir de cette décomposition. C'est le choix déterminant pour la profondeur du circuit. L'évolution à un corps $e^{-i t H_1}$ est une opération de fermions libres (gaussienne fermionique); il s'agit donc d'une rotation de base à une particule (une opération OrbitalRotation avec une $e^{-i t \, \mathbf{h}_1}$ unitaire) plutôt que d'un calcul nécessitant une application de la règle de Pauli-Trotter terme par terme. L'évolution à deux corps $e^{-i t H_2}$ n'agit que sur le terme d'interaction locale unique. En construisant cette étape de cette manière, le transpilateur génère un enchaînement simple de rotations de Givens ainsi qu’une seule porte d’interaction à deux qubits, soit la structure simple que la publication SKQD construit manuellement.

Nous avons déjà la partie « deux corps ». Il s'agit de l'opérateur two_body_momentum défini à l'étape 2 (le terme « sur site » seul). L'évolution à un corps $e^{-i t H_1}$ est indépendante du spin; il s'agit donc d'une seule opération par secteur, OrbitalRotation avec une $e^{-i t \, \mathbf{h}_1}$ unitaire appliquée une fois à chaque secteur de spin, ce qui correspond à la même structure parallèle que la préparation de l'état de référence.

Chaque circuit de Krylov évolue pendant un temps total de $k \, \Delta t$ s. Cela est réalisé à l'aide d'un produit de Trotter du second ordre, composé de $k$ étapes de taille $\Delta t$, de sorte que l'erreur par étape reste fixe à mesure que la dimension de Krylov augmente. À chaque étape, une rotation complète autour de l'axe central (une par secteur de rotation) est intercalée entre deux demi-étapes de l'évolution à deux corps (simplifiée, à un seul terme). La dimension de Krylov $D$ correspond au nombre de ces circuits, avec des valeurs comprises dim entre ( 0 état de référence uniquement) et $D - 1$.

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

5. Conversion en circuits de qubits

La transpilation FermionicCircuit associe les portes fermioniques à des portes de qubits. On utilise ici la transformation jordan_wigner() de Jordan-Wigner. Les portes OrbitalRotation se synthétisent en un ensemble de rotations de Givens, et le terme unique se Evolution réduit à une seule porte d'interaction à deux qubits.

Comme chaque porte fermionique dispose ici d'une synthèse par défaut, vous pouvez utiliser directement generate_preset_jw_pass_manager() les éléments prêts à l'emploi plutôt que d'assembler vous-même les étapes. Tous les arguments de mots-clés sont transmis à l’étape « qubit »; utilisez « pass » pour optimization_level=0 obtenir une représentation fidèle, mais non optimisée, de la profondeur synthétisée. La génération de l'ensemble complet des circuits de Krylov se fait alors à l'aide d'une boucle sur la dimension de Krylov. Conservez les circuits de qubits « FermionicCircuits » non transpilés à côté des circuits de qubits « qu » transpilés. Les versions fermioniques servent à la simulation exacte (vecteur d'état) à l'étape 6, tandis que les versions transpilées (avec les mesures ajoutées) correspondent à ce qu'un backend exécuterait.

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         circuit = pass_manager.run(fermionic_circuit)
...         circuit.measure_all()
...         circuits.append(circuit)
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

La dimension de Krylov se reflète directement dans la profondeur à deux qubits des circuits. La préparation de l’état de référence représente un coût fixe, et chaque puissance de Krylov supplémentaire ajoute une étape de Trotter d’ordre deux à l’opérateur d’évolution temporelle.

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

Le premier circuit ne calcule que le déterminant de référence, tandis que le dernier comporte le plus grand nombre d'étapes de Trotter (quatre dans ce cas) :

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. Chaînes de bits échantillonnées à partir des circuits de Krylov

Au niveau matériel, vous exécuteriez les circuits transpilés et effectueriez les mesures. Il s'agit de la version « silencieuse » (vecteur d'état) de cette étape, qui montre les chaînes de bits que SKQD collecterait. Simulez chaque élément non transposé FermionicCircuit avec ffsim.apply_unitary() et échantillonnez le vecteur d'état résultant avec ffsim.sample_state_vector().

En simulation, une porte PrepareSlaterDeterminant fonctionne selon le principe « validation puis rotation » (elle vérifie que l'état entrant correspond à son déterminant de référence, puis applique la rotation). Par conséquent, la donnée d'entrée correspond au déterminant d'occupation simple : les modes nocc les plus bas occupés par spin, sans rotation. La porte s'ouvre reference_rotation d'elle-même.

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

Si l'on regroupe les échantillons provenant des cinq circuits de Krylov, l'ensemble de support est minuscule : quelques centaines de chaînes binaires distinctes parmi les déterminants d' $\binom{8}{4}^2 = 4900$ s du secteur (4, 4) . Cette concentration, conséquence directe de la parcimonie de la base d'impulsion établie à l'étape 2, rend la diagonalisation classique qui suit plus facile à traiter.

En représentant graphiquement les nombres d'occurrences à l'aide de, on constate qiskit.visualization.plot_histogram() clairement la rareté des données. Quelques configurations prédominent, avec en tête le déterminant de référence.

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. Diagonaliser dans le sous-espace échantillonné

Les chaînes de bits échantillonnées constituent l'entrée de la partie classique du SKQD. L'hamiltonien est projeté sur le sous-espace engendré par les configurations échantillonnées, puis y est diagonalisé. Cette valeur est transmise au paquet qiskit-addon-sqd, qui diagonalize_fermionic_hamiltonian() exécute la boucle complète de diagonalisation quantique basée sur des échantillons. Il construit un sous-espace à partir de lots de configurations échantillonnées, diagonalise l'hamiltonien dans ce sous-espace, puis affine celui-ci de manière itérative en recourant à la récupération de configuration : il inverse les occupations pour corriger les configurations qui enfreignent la symétrie connue du nombre de particules. La restauration de la configuration est conçue pour corriger les dommages causés par le bruit matériel, qui altère les chaînes de bits échantillonnées en les plaçant dans un secteur comportant un nombre de particules erroné; l'échantillonnage étant ici exempt de bruit, il n'y a donc rien à réparer de ce côté-là.

Le solveur utilise les nombres échantillonnés sous la forme d'un BitArray, construit directement à partir du pool de counts l'étape 6. Elle se diagonalise dans la même base d'impulsion que celle utilisée pour le reste de la construction; elle reprend donc directement les intégrales h2e_momentum h1e_momentum et de l'étape 2, sans réorganisation. L'énergie récupérée correspond exactement à la valeur de référence, à bien moins d'un milli-Hartree près, à partir de quelques centaines seulement de configurations échantillonnées. L'énergie SQD constituant une borne supérieure variationnelle de l'état fondamental réel, elle est toujours égale ou supérieure à reference_energy (ce qui reste dans le cadre de l'étape 2) :

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts))
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

Visualisez le calcul : la convergence de l'énergie au fil des itérations, ainsi que l'occupation moyenne de chaque orbitale spatiale dans l'état fondamental récupéré. L'énergie par itération est la plus faible parmi tous les lots de cette itération, et l'occupation correspond à la somme des orbital_occupancies sur les deux secteurs de spin. Même en l'absence de bruit à corriger, l'énergie diminue d'itération en itération. À chaque itération, la diagonalisation s'effectue au sein de lots num_batches ne comprenant que configurations samples_per_batch , sélectionnées par sous-échantillonnage à partir de l'ensemble; le raffinement itératif résultant du report de la chaîne binaire améliore progressivement la sélection des configurations qui se retrouvent dans ces lots, ce qui permet d'augmenter l'énergie récupérée :

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

Cela permet de boucler la boucle SKQD. Les circuits basés sur la quantité de mouvement des étapes 1 à 5 produisent un échantillon clairsemé (étape 6), et la diagonalisation basée sur cet échantillon décrite ci-dessus permet de reconvertir cet échantillon en une énergie d'état fondamental.

Etapes suivantes

Ce guide a exécuté l'intégralité du pipeline SKQD sur un simulateur de vecteurs d'état sans bruit. Au niveau matériel, le code transpilé circuits serait exécuté et mesuré à la place de l'échantillonnage du vecteur d'état à l'étape 6, les comptes obtenus alimentant sans modification la diagonalisation de step-7. Consultez la documentation de Qiskit pour obtenir de l'aide sur l'exécution des circuits, et les tutoriels du module complémentaire SQD pour en savoir plus sur le post-traitement par diagonalisation en sous-espace, notamment son comportement sur des échantillons bruités, où la récupération de la configuration joue un rôle prépondérant.

Consultez le guide des relations ffsim pour comprendre pourquoi et ffsim.apply_unitary() ffsim.sample_state_vector() fonctionnent de manière native sur les circuits fermioniques de ce package, et pourquoi le secteur à nombre de particules fixe utilisé à l'étape 6 rend l'éparsité basée sur la quantité de mouvement, établie à l'étape 2, efficace en termes d'échantillonnage.