Skip to main content
IBM Quantum Platform

Simulation quantique

Note

Yukio Kawashima (30 mai 2024)

Télécharger le pdf de la conférence originale. Notez que certains extraits de code peuvent devenir obsolètes car il s'agit d'images statiques.

Le temps approximatif d'exécution de cette expérience par la QPU est de 7 secondes.

(Ce cahier est en grande partie tiré d'un cahier de tutoriel pour les algorithmes Qiskit, aujourd'hui obsolète)


1. Introduction

En tant que technique d'évolution en temps réel, la trotterisation consiste en l'application successive d'une ou plusieurs portes quantiques, choisies pour approximer l'évolution temporelle d'un système pour une tranche de temps. A partir de l'équation de Schrödinger, l'évolution temporelle d'un système initialement dans l'état ψ(0)\vert\psi(0)\rangle prend la forme :

ψ(t)=eiHtψ(0),\vert \psi(t) \rangle = e^{-i H t} \vert \psi(0) \rangle \text{,}

HH est le hamiltonien indépendant du temps qui régit le système. Nous considérons un hamiltonien qui peut être écrit comme une somme pondérée de termes de Pauli H=jajPjH=\sum_j a_j P_j, PjP_j représentant un produit tensoriel de termes de Pauli agissant sur nn qubits. En particulier, ces termes de Pauli peuvent commuter entre eux ou non. Étant donné un état à l'instant t=0t=0, comment obtenir l'état du système à un instant ultérieur ψ(t)|\psi(t)\rangle à l'aide d'un ordinateur quantique? L'exponentielle d'un opérateur peut être plus facilement comprise à travers sa série de Taylor :

eiHt=1iHt12H2t2+...e^{-i H t} = 1-iHt-\frac{1}{2}H^2t^2+...

Certaines exponentielles très basiques, comme eiZe^{iZ}, peuvent être facilement mises en œuvre sur des ordinateurs quantiques à l'aide d'un ensemble compact de portes quantiques. La plupart des hamiltoniens qui nous intéressent n'ont pas un seul terme, mais plusieurs. Notez ce qui se passe si H=H1+H2H = H_1+H_2 :

eiHt=1i(H1+H2)t12(H1+H2)2t2+...e^{-i H t} = 1-i(H_1+H_2)t-\frac{1}{2}(H_1+H_2)^2t^2+...

Lorsque H1H_1 et H2H_2 commutent, nous avons le cas familier (qui est également vrai pour les nombres, et les variables aa et bb ci-dessous) :

ei(a+b)t=eiateibte^{-i (a+b) t} = e^{-i a t}e^{-i b t}

Mais lorsque les opérateurs ne commutent pas, les termes ne peuvent pas être réarrangés dans la série de Taylor pour simplifier de cette manière. L'expression d'hamiltoniens compliqués dans des portes quantiques constitue donc un défi.

Une solution consiste à considérer un temps très petit tt, de sorte que le terme du premier ordre dans l'expansion de Taylor domine. Dans cette hypothèse :

ei(H1+H2)t1i(H1+H2)t(1iH1t)(1iH2t)eiH1teiH2te^{-i (H_1+H_2) t} \approx 1-i(H_1+H_2)t \approx (1-i H_1 t)(1-i H_2 t) \approx e^{-i H_1 t}e^{-i H_2 t}

Bien entendu, il se peut que nous devions faire évoluer notre état pendant une période plus longue. Pour ce faire, il est nécessaire de procéder à de nombreux petits pas dans le temps. Ce processus s'appelle la trotterisation :

ψ(t)(jeiajPjt/r)rψ(0),\vert \psi(t) \rangle \approx \left(\prod_j e^{-i a_j P_j t/r} \right)^r \vert\psi(0) \rangle \text{,}

Ici, t/rt/r est la tranche de temps (étape de l'évolution) que nous choisissons. Il en résulte une porte à appliquer rr fois. Un pas de temps plus petit conduit à une approximation plus précise. Toutefois, cela conduit également à des circuits plus profonds qui, dans la pratique, entraînent une plus grande accumulation d'erreurs (un problème non négligeable pour les dispositifs quantiques à court terme).

Aujourd'hui, nous allons étudier l'évolution temporelle du modèle d'Ising sur des treillis linéaires de sites N=2N=2 et N=6N=6. Ces treillis sont constitués d'un réseau de spins σi\sigma_i qui n'interagissent qu'avec leurs voisins les plus proches. Ces spins peuvent avoir deux orientations : \uparrow et \downarrow, qui correspondent à une magnétisation de +1+1 et 1-1 respectivement.

H=Ji=0N2ZiZi+1hi=0N1Xi,H = - J \sum_{i=0}^{N-2} Z_i Z_{i+1} - h \sum_{i=0}^{N-1} X_i \text{,}

JJ décrit l'énergie d'interaction et hh l'amplitude d'un champ externe (dans la direction x ci-dessus, mais nous modifierons cette valeur). Écrivons cette expression en utilisant les matrices de Pauli, et en considérant que le champ externe a un angle α\alpha par rapport à la direction transversale,

H=Ji=0N2ZiZi+1hi=0N1(sinαZi+cosαXi).H = -J \sum_{i=0}^{N-2} Z_i Z_{i+1} -h \sum_{i=0}^{N-1} (\sin\alpha Z_i + \cos\alpha X_i) \text{.}

Ce hamiltonien est utile car il nous permet d'étudier facilement les effets d'un champ extérieur. Dans la base de calcul, le système sera codé comme suit :

État quantique
Représentation du spin
0000\lvert 0 0 0 0 \rangle\uparrow\uparrow\uparrow\uparrow
1000\lvert 1 0 0 0 \rangle\downarrow\uparrow\uparrow\uparrow
\ldots\ldots
1111\lvert 1 1 1 1 \rangle\downarrow\downarrow\downarrow\downarrow

Nous commencerons par étudier l'évolution temporelle d'un tel système quantique. Plus précisément, nous visualiserons l'évolution dans le temps de certaines propriétés du système, comme l'aimantation.

1.1 Exigences

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

  • Qiskit SDK v1.2 ou plus récent ( pip install qiskit )
  • Qiskit Runtime v0.30 ou plus tard ( pip install qiskit-ibm-runtime )
  • Numpy v1.24.1 ou plus récent < 2 ( pip install numpy )

1.2 Importer les bibliothèques

Notez que certaines bibliothèques qui pourraient être utiles ( MatrixExponential, QDrift) sont incluses même si elles ne sont pas utilisées dans ce cahier. Vous pouvez les essayer si vous avez le temps!

# Check the version of Qiskit
import qiskit

qiskit.__version__

Output:

'2.0.2'
# Import the qiskit library
import numpy as np
import matplotlib.pylab as plt
import warnings

from qiskit import QuantumCircuit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.primitives import StatevectorEstimator
from qiskit.quantum_info import Statevector, SparsePauliOp
from qiskit.synthesis import (
    SuzukiTrotter,
    LieTrotter,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2

warnings.filterwarnings("ignore")

2. Cartographier votre problème

2.1 Définition de l'hamiltonien d'Ising à champ transversal

Nous considérons ici le modèle d'Ising à champ transverse 1-D.

Tout d'abord, nous allons créer une fonction qui prend en compte les paramètres du système NN, JJ, hh et α\alpha, et qui renvoie notre hamiltonien sous la forme d'un SparsePauliOp. A SparsePauliOp est une représentation éparse d'un opérateur en termes de termes de Pauli pondérés.

def get_hamiltonian(nqubits, J, h, alpha):
    # List of Hamiltonian terms as 3-tuples containing
    # (1) the Pauli string,
    # (2) the qubit indices corresponding to the Pauli string,
    # (3) the coefficient.
    ZZ_tuples = [("ZZ", [i, i + 1], -J) for i in range(0, nqubits - 1)]
    Z_tuples = [("Z", [i], -h * np.sin(alpha)) for i in range(0, nqubits)]
    X_tuples = [("X", [i], -h * np.cos(alpha)) for i in range(0, nqubits)]

    # We create the Hamiltonian as a SparsePauliOp, via the method
    # `from_sparse_list`, and multiply by the interaction term.
    hamiltonian = SparsePauliOp.from_sparse_list(
        [*ZZ_tuples, *Z_tuples, *X_tuples], num_qubits=nqubits
    )
    return hamiltonian.simplify()

Définir l'hamiltonien

Le système que nous considérons maintenant a une taille de N=6N=6, J=0.2J=0.2, h=1.2h=1.2 et α=π8.0\alpha=\frac{\pi}{8.0} à titre d'exemple.

n_qubits = 6

hamiltonian = get_hamiltonian(nqubits=n_qubits, J=0.2, h=1.2, alpha=np.pi / 8.0)
hamiltonian

Output:

SparsePauliOp(['IIIIZZ', 'IIIZZI', 'IIZZII', 'IZZIII', 'ZZIIII', 'IIIIIZ', 'IIIIZI', 'IIIZII', 'IIZIII', 'IZIIII', 'ZIIIII', 'IIIIIX', 'IIIIXI', 'IIIXII', 'IIXIII', 'IXIIII', 'XIIIII'],
              coeffs=[-0.2       +0.j, -0.2       +0.j, -0.2       +0.j, -0.2       +0.j,
 -0.2       +0.j, -0.45922012+0.j, -0.45922012+0.j, -0.45922012+0.j,
 -0.45922012+0.j, -0.45922012+0.j, -0.45922012+0.j, -1.10865544+0.j,
 -1.10865544+0.j, -1.10865544+0.j, -1.10865544+0.j, -1.10865544+0.j,
 -1.10865544+0.j])

2.2 Définir les paramètres de la simulation d'évolution temporelle

Nous examinerons ici trois techniques différentes de trotterisation :

  • Lie-Trotter (première commande)
  • suzuki-Trotter du second ordre
  • suzuki-Trotter du quatrième ordre

Ces deux derniers seront utilisés dans l'exercice et dans l'annexe.

num_timesteps = 60
evolution_time = 30.0
dt = evolution_time / num_timesteps
product_formula_lt = LieTrotter()
product_formula_st2 = SuzukiTrotter(order=2)
product_formula_st4 = SuzukiTrotter(order=4)

2.3 Préparez le circuit quantique 1 (état initial)

Créer un état initial. Nous commencerons par une configuration de spin de \uparrow\uparrow\downarrow\downarrow\uparrow\uparrow.

initial_circuit = QuantumCircuit(n_qubits)
initial_circuit.prepare_state("001100")
# Change reps and see the difference when you decompose the circuit
initial_circuit.decompose(reps=1).draw("mpl")

Output:

Output of the previous code cell

2.4 Préparer le circuit quantique 2 (circuit unique pour l'évolution temporelle)

Nous construisons ici un circuit pour un seul pas de temps en utilisant Lie-Trotter.

La formule du produit de Lie (premier ordre) est implémentée dans la classe LieTrotter classe. Une formule du premier ordre consiste en l'approximation mentionnée dans l'introduction, où l'exponentielle matricielle d'une somme est approximée par un produit d'exponentielles matricielles :

eH1+H2eH1eH2e^{H_1+H_2} \approx e^{H_1} e^{H_2}

Comme indiqué précédemment, les circuits très profonds conduisent à l'accumulation d'erreurs et posent des problèmes aux ordinateurs quantiques modernes. Les portes à deux qubits ayant des taux d'erreur plus élevés que les portes à un qubit, la profondeur du circuit à deux qubits présente un intérêt particulier. Ce qui compte vraiment, c'est la profondeur du circuit à deux qubits après la transpilation (puisque c'est le circuit que l'ordinateur quantique exécute réellement). Mais prenons l'habitude de compter les opérations pour ce circuit, même maintenant en utilisant le simulateur.

single_step_evolution_gates_lt = PauliEvolutionGate(
    hamiltonian, dt, synthesis=product_formula_lt
)
single_step_evolution_lt = QuantumCircuit(n_qubits)
single_step_evolution_lt.append(
    single_step_evolution_gates_lt, single_step_evolution_lt.qubits
)

print(
    f"""
Trotter step with Lie-Trotter
-----------------------------
Depth: {single_step_evolution_lt.decompose(reps=3).depth()}
Gate count: {len(single_step_evolution_lt.decompose(reps=3))}
Nonlocal gate count: {single_step_evolution_lt.decompose(reps=3).num_nonlocal_gates()}
Gate breakdown: {", ".join([f"{k.upper()}: {v}" for k, v in single_step_evolution_lt.decompose(reps=3).count_ops().items()])}
"""
)
single_step_evolution_lt.decompose(reps=3).draw("mpl", fold=-1)

Output:


Trotter step with Lie-Trotter
-----------------------------
Depth: 17
Gate count: 27
Nonlocal gate count: 10
Gate breakdown: U3: 12, CX: 10, U1: 5

Output of the previous code cell

2.5 Régler les opérateurs pour mesurer

Définissons un opérateur de magnétisation iZi/N\sum_i \langle Z_i \rangle / N et un opérateur de corrélation de spin moyen iZiZi+1/(N1)\sum_i \langle Z_i Z_{i+1} \rangle/ (N - 1).

magnetization = (
    SparsePauliOp.from_sparse_list(
        [("Z", [i], 1.0) for i in range(0, n_qubits)], num_qubits=n_qubits
    )
    / n_qubits
)
correlation = SparsePauliOp.from_sparse_list(
    [("ZZ", [i, i + 1], 1.0) for i in range(0, n_qubits - 1)], num_qubits=n_qubits
) / (n_qubits - 1)
print("magnetization : ", magnetization)
print("correlation : ", correlation)

Output:

magnetization :  SparsePauliOp(['IIIIIZ', 'IIIIZI', 'IIIZII', 'IIZIII', 'IZIIII', 'ZIIIII'],
              coeffs=[0.16666667+0.j, 0.16666667+0.j, 0.16666667+0.j, 0.16666667+0.j,
 0.16666667+0.j, 0.16666667+0.j])
correlation :  SparsePauliOp(['IIIIZZ', 'IIIZZI', 'IIZZII', 'IZZIII', 'ZZIIII'],
              coeffs=[0.2+0.j, 0.2+0.j, 0.2+0.j, 0.2+0.j, 0.2+0.j])

2.6 Effectuer une simulation d'évolution temporelle

Nous surveillerons l'énergie (valeur d'espérance de l'hamiltonien), l'aimantation (valeur d'espérance de l'opérateur d'aimantation) et la corrélation moyenne de spin (valeur d'espérance de l'opérateur de corrélation moyenne de spin). La primitive StatevectorEstimator ( EstimatorV2 ) de Qiskit estime les valeurs attendues des observables, ψO^ψ\langle\psi\vert\hat{O}\vert\psi\rangle.

# Initiate the circuit
evolved_state = QuantumCircuit(initial_circuit.num_qubits)
# Start from the initial spin configuration
evolved_state.append(initial_circuit, evolved_state.qubits)
# Initiate Estimator (V2)
estimator = StatevectorEstimator()
# Set number of shots
shots = 10000
# Translate the precision required from the number of shots
precision = np.sqrt(1 / shots)
energy_list = []
mag_list = []
corr_list = []
# Estimate expectation values for t=0.0
job = estimator.run(
    [(evolved_state, [hamiltonian, magnetization, correlation])], precision=precision
)
# Get estimated expectation values
evs = job.result()[0].data.evs
energy_list.append(evs[0])
mag_list.append(evs[1])
corr_list.append(evs[2])
# Start time evolution
for n in range(num_timesteps):
    # Expand the circuit to describe delta-t
    evolved_state.append(single_step_evolution_gates_lt, evolved_state.qubits)
    # Estimate expectation values at delta-t
    job = estimator.run(
        [(evolved_state, [hamiltonian, magnetization, correlation])],
        precision=precision,
    )
    # Retrieve results (expectation values)
    evs = job.result()[0].data.evs
    energy_list.append(evs[0])
    mag_list.append(evs[1])
    corr_list.append(evs[2])
# Transform the list of expectation values (at each time step) to arrays
energy_array = np.array(energy_list)
mag_array = np.array(mag_list)
corr_array = np.array(corr_list)

2.7 Tracer l'évolution temporelle des observables

Nous traçons les valeurs d'attente que nous avons mesurées en fonction du temps.

fig, axes = plt.subplots(3, sharex=True)
times = np.linspace(0, evolution_time, num_timesteps + 1)  # includes initial state
axes[0].plot(
    times,
    energy_array,
    label="First order",
    marker="x",
    c="darkmagenta",
    ls="-",
    lw=0.8,
)
axes[1].plot(
    times, mag_array, label="First order", marker="x", c="darkmagenta", ls="-", lw=0.8
)
axes[2].plot(
    times, corr_array, label="First order", marker="x", c="darkmagenta", ls="-", lw=0.8
)
axes[0].set_ylabel("Energy")
axes[1].set_ylabel("Magnetization")
axes[2].set_ylabel("Mean spin correlation")
axes[2].set_xlabel("Time")
fig.suptitle("Observable evolution")

Output:

Text(0.5, 0.98, 'Observable evolution')
Output of the previous code cell

3. Exercice 1. Effectuer une simulation à l'aide de la méthode Suzuki-Trotter du second ordre

Essayons maintenant d'effectuer une simulation avec Suzuki-Trotter du second ordre en suivant l'exemple de Lie-Trotter présenté ci-dessus.

Le Suzuki-Trotter du second ordre peut être utilisé dans Qiskit au moyen de la classe SuzukiTrotter. En utilisant cette formule, on obtient une décomposition de second ordre :

eH1+H2eH1/2eH2eH1/2e^{H_1+H_2} \approx e^{H_1/2}e^{H_2}e^{H_1/2}

3.1 Construire un circuit pour un seul pas de temps

Utilisez product_formula_st2 ( SuzukiTrotter(order=2 )) et construisez un circuit pour un seul pas de temps en utilisant Suzuki-Trotter du second ordre. Comptez également le nombre de portes et la profondeur du circuit et comparez avec Lie-Trotter.

# Modify the line below (Use PauliEvolutionGate)
single_step_evolution_gates_st2 = PauliEvolutionGate(
    hamiltonian, dt, synthesis=product_formula_st2
)
single_step_evolution_st2 = QuantumCircuit(n_qubits)
single_step_evolution_st2.append(
    single_step_evolution_gates_st2, single_step_evolution_st2.qubits
)
# Let us print some stats
print(
    f"""
Trotter step with second-order Suzuki-Trotter
-----------------------------
Depth: {single_step_evolution_st2.decompose(reps=3).depth()}
Gate count: {len(single_step_evolution_st2.decompose(reps=3))}
Nonlocal gate count: {single_step_evolution_st2.decompose(reps=3).num_nonlocal_gates()}
Gate breakdown: {", ".join([f"{k.upper()}: {v}" for k, v in single_step_evolution_st2.decompose(reps=3).count_ops().items()])}
"""
)
single_step_evolution_st2.decompose(reps=2).draw("mpl", fold=-1)

Output:


Trotter step with second-order Suzuki-Trotter
-----------------------------
Depth: 34
Gate count: 53
Nonlocal gate count: 20
Gate breakdown: U3: 23, CX: 20, U1: 10

Output of the previous code cell

3.2 Effectuer une simulation d'évolution temporelle

Effectuer l'évolution temporelle à l'aide de Suzuki-Trotter du deuxième ordre.

# Initiate the circuit
evolved_state = QuantumCircuit(initial_circuit.num_qubits)
# Start from the initial spin configuration
evolved_state.append(initial_circuit, evolved_state.qubits)
# Initiate Estimator (V2)
estimator = StatevectorEstimator()
# Set number of shots
shots = 10000
# Translate the precision required from the number of shots
precision = np.sqrt(1 / shots)
energy_list_st2 = []
mag_list_st2 = []
corr_list_st2 = []
# Estimate expectation values for t=0.0
job = estimator.run(
    [(evolved_state, [hamiltonian, magnetization, correlation])], precision=precision
)
# Get estimated expectation values
evs = job.result()[0].data.evs
energy_list_st2.append(evs[0])
mag_list_st2.append(evs[1])
corr_list_st2.append(evs[2])
# Start time evolution
for n in range(num_timesteps):
    # Expand the circuit to describe delta-t
    evolved_state.append(single_step_evolution_gates_st2, evolved_state.qubits)
    # Estimate expectation values at delta-t
    job = estimator.run(
        [(evolved_state, [hamiltonian, magnetization, correlation])],
        precision=precision,
    )
    # Retrieve results (expectation values)
    evs = job.result()[0].data.evs
    energy_list_st2.append(evs[0])
    mag_list_st2.append(evs[1])
    corr_list_st2.append(evs[2])
# Transform the list of expectation values (at each time step) to arrays
energy_array_st2 = np.array(energy_list_st2)
mag_array_st2 = np.array(mag_list_st2)
corr_array_st2 = np.array(corr_list_st2)

3.3 Tracer les résultats de Suzuki-Trotter du second ordre

axes[0].plot(
    times,
    energy_array_st2,
    label="Second Order",
    marker="x",
    c="limegreen",
    ls="-",
    lw=0.8,
)
axes[1].plot(
    times,
    mag_array_st2,
    label="Second Order",
    marker="x",
    c="limegreen",
    ls="-",
    lw=0.8,
)
axes[2].plot(
    times,
    corr_array_st2,
    label="Second Order",
    marker="x",
    c="limegreen",
    ls="-",
    lw=0.8,
)

# Replace the legend
# legend.remove()
legend = fig.legend(
    *axes[0].get_legend_handles_labels(),
    bbox_to_anchor=(1.0, 0.5),
    loc="center left",
    framealpha=0.5,
)
fig

Output:

Output of the previous code cell

3.4 Comparer avec les résultats exacts

Les données ci-dessous sont les résultats exacts précalculés par l'ordinateur classique.

exact_times = np.array(
    [
        0.0,
        0.3,
        0.6,
        0.8999999999999999,
        1.2,
        1.5,
        1.7999999999999998,
        2.1,
        2.4,
        2.6999999999999997,
        3.0,
        3.3,
        3.5999999999999996,
        3.9,
        4.2,
        4.5,
        4.8,
        5.1,
        5.3999999999999995,
        5.7,
        6.0,
        6.3,
        6.6,
        6.8999999999999995,
        7.199999999999999,
        7.5,
        7.8,
        8.1,
        8.4,
        8.7,
        9.0,
        9.299999999999999,
        9.6,
        9.9,
        10.2,
        10.5,
        10.799999999999999,
        11.1,
        11.4,
        11.7,
        12.0,
        12.299999999999999,
        12.6,
        12.9,
        13.2,
        13.5,
        13.799999999999999,
        14.1,
        14.399999999999999,
        14.7,
        15.0,
        15.299999999999999,
        15.6,
        15.899999999999999,
        16.2,
        16.5,
        16.8,
        17.099999999999998,
        17.4,
        17.7,
        18.0,
        18.3,
        18.599999999999998,
        18.9,
        19.2,
        19.5,
        19.8,
        20.099999999999998,
        20.4,
        20.7,
        21.0,
        21.3,
        21.599999999999998,
        21.9,
        22.2,
        22.5,
        22.8,
        23.099999999999998,
        23.4,
        23.7,
        24.0,
        24.3,
        24.599999999999998,
        24.9,
        25.2,
        25.5,
        25.8,
        26.099999999999998,
        26.4,
        26.7,
        27.0,
        27.3,
        27.599999999999998,
        27.9,
        28.2,
        28.5,
        28.799999999999997,
        29.099999999999998,
        29.4,
        29.7,
        30.0,
    ]
)
exact_energy = np.array(
    [
        -1.1184402376762155,
        -1.1184402376762157,
        -1.1184402376762157,
        -1.1184402376762148,
        -1.1184402376762153,
        -1.1184402376762155,
        -1.1184402376762148,
        -1.118440237676216,
        -1.118440237676216,
        -1.1184402376762166,
        -1.1184402376762148,
        -1.118440237676216,
        -1.1184402376762153,
        -1.1184402376762148,
        -1.118440237676217,
        -1.118440237676215,
        -1.1184402376762161,
        -1.1184402376762157,
        -1.118440237676217,
        -1.1184402376762161,
        -1.1184402376762137,
        -1.1184402376762161,
        -1.1184402376762161,
        -1.118440237676218,
        -1.1184402376762155,
        -1.1184402376762166,
        -1.1184402376762155,
        -1.1184402376762137,
        -1.1184402376762186,
        -1.1184402376762215,
        -1.1184402376762148,
        -1.118440237676216,
        -1.1184402376762166,
        -1.1184402376762148,
        -1.1184402376762121,
        -1.1184402376762166,
        -1.1184402376762181,
        -1.1184402376762137,
        -1.1184402376762148,
        -1.1184402376762193,
        -1.1184402376762108,
        -1.1184402376762144,
        -1.118440237676217,
        -1.1184402376762197,
        -1.1184402376762153,
        -1.1184402376762161,
        -1.1184402376762184,
        -1.1184402376762126,
        -1.118440237676214,
        -1.118440237676214,
        -1.1184402376762161,
        -1.118440237676212,
        -1.1184402376762164,
        -1.118440237676217,
        -1.1184402376762121,
        -1.1184402376762157,
        -1.1184402376762212,
        -1.1184402376762217,
        -1.1184402376762206,
        -1.118440237676222,
        -1.1184402376762166,
        -1.118440237676212,
        -1.1184402376762137,
        -1.11844023767622,
        -1.1184402376762206,
        -1.118440237676219,
        -1.1184402376762153,
        -1.1184402376762164,
        -1.118440237676209,
        -1.1184402376762144,
        -1.1184402376762161,
        -1.118440237676216,
        -1.1184402376762173,
        -1.118440237676214,
        -1.1184402376762093,
        -1.1184402376762184,
        -1.1184402376762126,
        -1.118440237676213,
        -1.1184402376762195,
        -1.1184402376762095,
        -1.1184402376762075,
        -1.1184402376762197,
        -1.1184402376762141,
        -1.1184402376762146,
        -1.1184402376762184,
        -1.118440237676218,
        -1.1184402376762224,
        -1.118440237676219,
        -1.118440237676218,
        -1.1184402376762206,
        -1.1184402376762168,
        -1.118440237676221,
        -1.118440237676218,
        -1.1184402376762148,
        -1.1184402376762106,
        -1.1184402376762173,
        -1.118440237676216,
        -1.118440237676216,
        -1.1184402376762113,
        -1.1184402376762275,
        -1.1184402376762195,
    ]
)
exact_magnetization = np.array(
    [
        0.3333333333333333,
        0.26316769633415005,
        0.0912947227110664,
        -0.09317712543141576,
        -0.20391854332115245,
        -0.19318196655046493,
        -0.06411527074401464,
        0.12558269854206197,
        0.28252754464640606,
        0.3264196194042506,
        0.2361586169847769,
        0.060894367906122224,
        -0.10842387093076275,
        -0.18636359582538073,
        -0.1338364343947887,
        0.020284606520827753,
        0.19151142743926025,
        0.2905341647678381,
        0.2723014646745304,
        0.15147481733047252,
        -0.008179102877790292,
        -0.1242999208732406,
        -0.1372529247781061,
        -0.04083616185958952,
        0.11066094926716476,
        0.23140661570567636,
        0.2587109403786205,
        0.1868237670027325,
        0.061201779383143744,
        -0.051391248969654205,
        -0.09843899603365061,
        -0.061297056158849166,
        0.04199010081939773,
        0.15861461430963147,
        0.22336830674799552,
        0.20179555623336537,
        0.11407111438609417,
        0.01609419104778282,
        -0.04239611796730001,
        -0.04249123521065924,
        0.008850291714888112,
        0.08780898151558082,
        0.1561486776507056,
        0.17627348772811832,
        0.13870676179652253,
        0.07205869195282538,
        0.018300003064909465,
        0.0001095640839572417,
        0.015157929316037586,
        0.05077755280969454,
        0.09245534457650838,
        0.12206907551110702,
        0.12284950557969157,
        0.09570215398601932,
        0.06294378255078983,
        0.045503313813986014,
        0.043389819499542556,
        0.046725117769796744,
        0.054956411358382404,
        0.0713814528253614,
        0.08743689703248492,
        0.08951216359166674,
        0.07878386475305985,
        0.06955669116405788,
        0.06639892435963689,
        0.05890378761746903,
        0.04541796525844558,
        0.0414221088331947,
        0.05499634106912299,
        0.07409418836014572,
        0.08371859070160165,
        0.08211623987959302,
        0.07615055161378328,
        0.06702584458783024,
        0.051891407742740085,
        0.038049378383635625,
        0.03825614149768043,
        0.054183218463525695,
        0.0753534475741016,
        0.08853147112587295,
        0.08767917178542013,
        0.07709383184439536,
        0.06308595032042386,
        0.0498812359204284,
        0.04299040064096167,
        0.04769159891460652,
        0.06483569572288776,
        0.08698035745435016,
        0.10047391641776235,
        0.09747255683203637,
        0.08098863187287358,
        0.05959496723987331,
        0.04383882265040485,
        0.04232138798062125,
        0.05720514169944535,
        0.08201306299870219,
        0.10274898262000469,
        0.10707552455080133,
        0.09210856128265357,
        0.06379922105742579,
        0.03624325103307953,
    ]
)
exact_correlation = np.array(
    [
        0.2,
        0.1247704225763532,
        0.01943938494098705,
        0.03854917181332821,
        0.11196616231067426,
        0.0906546700356683,
        0.01629373561896267,
        0.011352652889791095,
        0.0636185676540077,
        0.09543834437789013,
        0.10058518161011307,
        0.11829217731417431,
        0.1397812224038133,
        0.12316460402216707,
        0.08541383059335775,
        0.06144846844403662,
        0.020246372880505827,
        -0.02693683090021662,
        0.003919250903281282,
        0.1117419430168554,
        0.19676155181256794,
        0.18594408880783336,
        0.1002673802566004,
        0.03821525827438024,
        0.04485205090247377,
        0.05348102743040269,
        0.03160026140008638,
        0.033437649060464834,
        0.10486939975320728,
        0.20249469538955758,
        0.19735507621013149,
        0.0553097261765083,
        -0.04889114490131667,
        0.011685690974970964,
        0.11705971535823065,
        0.11681165998194759,
        0.06637091239560744,
        0.10936684225958895,
        0.20225454101061405,
        0.16284420833341812,
        -0.0025823294931362067,
        -0.0763416631752919,
        0.02985268630418397,
        0.15234468006771007,
        0.14606385406970995,
        0.0935341856492092,
        0.12325421854361143,
        0.17130422930386324,
        0.10383730044042278,
        -0.031333159406547614,
        -0.05241572078596815,
        0.07722509925347705,
        0.17642188574256007,
        0.12765340239966838,
        0.06309968945093776,
        0.11574687130499339,
        0.16978282647206913,
        0.0736143632571229,
        -0.05356602733119409,
        -0.0009649396796768892,
        0.15921620111869142,
        0.17760366431811037,
        0.04736297330213485,
        0.012122870263181897,
        0.13268065586830521,
        0.1728473023503636,
        0.03999259331072221,
        -0.036997053070222885,
        0.06951528580242439,
        0.1769169993516561,
        0.12290448295710298,
        0.012897784654866427,
        0.02859435620982225,
        0.12895847695150875,
        0.13629536955485938,
        0.05394621059822597,
        0.02298040588184324,
        0.07036499900317271,
        0.11706448623132719,
        0.10435285842074606,
        0.055721236329964965,
        0.04676334743672697,
        0.08417924910022263,
        0.10611161955304965,
        0.089304171047322,
        0.06098589533081194,
        0.06314519797488709,
        0.09431492621892917,
        0.09667836915967139,
        0.0651298357290882,
        0.05176966009147416,
        0.06727229484222669,
        0.08871788283607947,
        0.09907054249093444,
        0.09785167773502176,
        0.09277216140054353,
        0.07520999642062785,
        0.05894392248382922,
        0.07236135251622376,
        0.08608284185200156,
        0.07282922961856123,
    ]
)
axes[0].plot(exact_times, exact_energy, c="k", ls=":", label="Exact")
axes[1].plot(exact_times, exact_magnetization, c="k", ls=":", label="Exact")
axes[2].plot(exact_times, exact_correlation, c="k", ls=":", label="Exact")
# Replace the legend
legend.remove()
# Select the labels of only the first axis
legend = fig.legend(
    *axes[0].get_legend_handles_labels(),
    bbox_to_anchor=(1.0, 0.5),
    loc="center left",
    framealpha=0.5,
)
fig.tight_layout()
fig

Output:

Output of the previous code cell

4. Exécution sur le matériel quantique

Nous exécutons ensuite la simulation de l'évolution temporelle sur le matériel quantique. Nous travaillerons sur un problème plus petit, de la taille du treillis N=2. Nous faisons varier le paramètre α\alpha et observons la différence dans la dynamique de la fonction d'onde.

4.1 Étape 1. Mapper les entrées classiques à un problème quantique

Choisir la configuration initiale de la simulation :

n_qubits_2 = 2
dt_2 = 1.6
product_formula = LieTrotter(reps=1)

Réglez ensuite le circuit initial :

La configuration initiale de l'essorage sera "bas-haut"

# We prepare an initial state ↓↑ (10).
# Note that Statevector and SparsePauliOp interpret the qubits from right to left
initial_circuit_2 = QuantumCircuit(n_qubits_2)
initial_circuit_2.prepare_state("10")
# Change reps and see the difference when you decompose the circuit
initial_circuit_2.decompose(reps=1).draw("mpl")

Output:

Output of the previous code cell

Calculez maintenant la valeur de référence à l'aide d'un simulateur de vecteur d'état idéal.

bar_width = 0.1
# initial_state = Statevector.from_label("10")
final_time = 1.6
eps = 1e-5

# We create the list of angles in radians, with a small epsilon
# the exactly longitudinal field, which would present no dynamics at all
alphas = np.linspace(-np.pi / 2 + eps, np.pi / 2 - eps, 5)

for i, alpha in enumerate(alphas):
    evolved_state_2 = QuantumCircuit(initial_circuit_2.num_qubits)
    evolved_state_2.append(initial_circuit_2, evolved_state_2.qubits)
    hamiltonian_2 = get_hamiltonian(nqubits=2, J=0.2, h=1.0, alpha=alpha)
    single_step_evolution_gates_2 = PauliEvolutionGate(
        hamiltonian_2, dt_2, synthesis=product_formula
    )
    evolved_state_2.append(single_step_evolution_gates_2, evolved_state_2.qubits)
    evolved_state_2 = Statevector(evolved_state_2)
    # Dictionary of probabilities
    amplitudes_dict = evolved_state_2.probabilities_dict()
    labels = list(amplitudes_dict.keys())
    values = list(amplitudes_dict.values())
    # Convert angle to degrees
    alpha_str = f"$\\alpha={int(np.round(alpha * 180 / np.pi))}^\\circ$"
    plt.bar(np.arange(4) + i * bar_width, values, bar_width, label=alpha_str, alpha=0.7)

plt.xticks(np.arange(4) + 2 * bar_width, labels)
plt.xlabel("Measurement")
plt.ylabel("Probability")
plt.suptitle(
    f"Measurement probabilities at $t={final_time}$, for various field angles $\\alpha$\n"
    f"Initial state: 10, Linear lattice of size $L=2$"
)
plt.legend()

Output:

<matplotlib.legend.Legend at 0x11c816590>
Output of the previous code cell

Nous avons préparé un système initialement avec une séquence de spins \downarrow\uparrow, qui correspond à ψ(0)=10\vert\psi(0)\rangle = \vert10\rangle. Après l'avoir laissé évoluer pour t=1.6t=1.6 sous un champ transversal ( α=0\alpha=0^\circ ), nous sommes presque assurés de mesurer \uparrow\downarrow, c'est-à-dire d'avoir un échange de spin. (Les étiquettes sont interprétées de droite à gauche). Si le champ est longitudinal ( α=±90\alpha=\pm90^\circ ), nous n'aurons pas d'évolution, nous mesurerons donc le système tel qu'il a été préparé initialement, \downarrow\uparrow. Avec des angles intermédiaires, à α=±45\alpha=\pm45^\circ, nous pourrons mesurer toutes les combinaisons avec des probabilités différentes, l'échange de spin étant le plus probable avec une probabilité de 67%.

Construire un circuit pour une expérience HW

circuit_list = []
for i, alpha in enumerate(alphas):
    evolved_state_2 = QuantumCircuit(initial_circuit_2.num_qubits)
    evolved_state_2.append(initial_circuit_2, evolved_state_2.qubits)
    hamiltonian_2 = get_hamiltonian(nqubits=2, J=0.2, h=1.0, alpha=alpha)
    single_step_evolution_gates_2 = PauliEvolutionGate(
        hamiltonian_2, dt_2, synthesis=product_formula
    )
    evolved_state_2.append(single_step_evolution_gates_2, evolved_state_2.qubits)
    evolved_state_2.measure_all()
    circuit_list.append(evolved_state_2)

4.2 Étape 2. Optimiser pour le matériel cible

Nous spécifions un backend.

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
backend.name

Output:

'ibm_strasbourg'

Ensuite, nous transposons le circuit pour le backend sélectionné.

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
circuit_isa = pm.run(circuit_list)

Vérifier le circuit.

circuit_isa[1].draw("mpl", idle_wires=False)

Output:

Output of the previous code cell

4.3 Étape 3. Exécuter à l'aide des primitives « IBM Quantum »

La primitive Sampler ( V2 ) de Qiskit fournit le nombre de chaînes de bits mesurées.

sampler = SamplerV2(mode=backend)
job = sampler.run(circuit_isa)
job_id = job.job_id()
print("job id:", job_id)

Output:

job id: d13pswfmya70008ek070

Sauvegarder les résultats

results = job.result()

4.4 Étape 4. Résultats du post-traitement

Construisez l'histogramme des chaînes de bits, ce qui correspond à l'analyse de la fonction d'onde, et comparez-les aux valeurs idéales indiquées ci-dessus.

list_temp = ["00", "01", "10", "11"]

for i, alpha in enumerate(alphas):
    # Dictionary of probabilities
    amplitudes_dict = results[i].data.meas.get_counts()
    values = []
    for str_temp in list_temp:
        values.append(
            amplitudes_dict[str_temp] / 4096.0
        )  # divided by default number of shots
    # Convert angle to degrees
    alpha_str = f"$\\alpha={int(np.round(alpha * 180 / np.pi))}^\\circ$"
    plt.bar(np.arange(4) + i * bar_width, values, bar_width, label=alpha_str, alpha=0.7)

plt.xticks(np.arange(4) + 2 * bar_width, labels)
plt.xlabel("Measurement")
plt.ylabel("Probabilities")
plt.suptitle(
    f"Measurement probabilities at $t={final_time}$, for various field angles $\\alpha$\n"
    f"Initial state: 10, Linear lattice of size $L=2$"
)
plt.legend()

Output:

<matplotlib.legend.Legend at 0x11d7af990>
Output of the previous code cell

Nous présentons ici un exemple de construction d'un circuit utilisant un Suzuki-Trotter d'ordre supérieur (quatrième ordre). Essayons maintenant de construire une simulation de circuit avec un Suzuki-Trotter du quatrième ordre en suivant les exemples présentés ci-dessus.

Le Suzuki-Trotter de quatrième ordre peut être utilisé dans Qiskit au moyen de la classe SuzukiTrotter. Le quatrième ordre peut être évalué à l'aide de la relation de récurrence suivante. Notez que l'ordre de Suzuki-Trotter est noté " 2k " dans les équations suivantes.

U^ST(2k)(t)=[U^ST(2k2)(pkt)]2U^ST(2k2)((14pk)t)[U^ST(2k2)(pkt)]2\hat{U}_{ST(2k)}\left(t\right) = \left[ \hat{U}_{ST(2k-2)}\left(p_k t\right) \right]^2 \hat{U}_{ST(2k-2)}\left( (1- 4 p_k) t\right)\left[ \hat{U}_{ST(2k-2)}\left(p_k t\right) \right]^2 pk=1/(4412k1)p_k = 1 / \left(4-4^{\frac{1}{2k-1}}\right)

Construire un circuit pour un seul pas de temps

Utilisez product_formula_st4 ( SuzukiTrotter(order=4 )) et construisez un circuit pour un seul pas de temps en utilisant Suzuki-Trotter du quatrième ordre. Comptez également le nombre de portes et la profondeur du circuit et comparez avec Lie-Trotter et Suzuki-Trotter de second ordre.

# Modify the line below (Use PauliEvolutionGate)
single_step_evolution_gates_st4 = PauliEvolutionGate(
    hamiltonian, dt, synthesis=product_formula_st4
)
single_step_evolution_st4 = QuantumCircuit(n_qubits)
single_step_evolution_st4.append(
    single_step_evolution_gates_st4, single_step_evolution_st4.qubits
)
# Let us print some stats
print(
    f"""
Trotter step with second-order Suzuki-Trotter
-----------------------------
Depth: {single_step_evolution_st4.decompose(reps=3).depth()}
Gate count: {len(single_step_evolution_st4.decompose(reps=3))}
Nonlocal gate count: {single_step_evolution_st4.decompose(reps=3).num_nonlocal_gates()}
Gate breakdown: {", ".join([f"{k.upper()}: {v}" for k, v in single_step_evolution_st4.decompose(reps=3).count_ops().items()])}
"""
)
single_step_evolution_st4.decompose(reps=2).draw("mpl", fold=-1)

Output:


Trotter step with second-order Suzuki-Trotter
-----------------------------
Depth: 170
Gate count: 265
Nonlocal gate count: 100
Gate breakdown: U3: 115, CX: 100, U1: 50

Output of the previous code cell
# Check Qiskit version
import qiskit

qiskit.__version__

Output:

'2.0.2'
Cette page a-t-elle été utile ?
Signaler un bogue, une coquille ou proposer du contenu sur GitHub.