Skip to main content
IBM Quantum Platform

Simulazione quantistica

Note

Yukio Kawashima (30 maggio 2024)

Scarica il pdf della lezione originale. Si noti che alcuni frammenti di codice potrebbero diventare deprecati, poiché si tratta di immagini statiche.

Il tempo approssimativo della QPU per eseguire questo esperimento è di 7 secondi.

(Questo quaderno è in gran parte tratto da un quaderno di esercitazioni per gli algoritmi di Qiskit, ormai deprecato)


1. Introduzione

Come tecnica di evoluzione in tempo reale, la trotterizzazione consiste nell'applicazione successiva di una o più porte quantistiche, scelte per approssimare l'evoluzione temporale di un sistema per una fetta di tempo. A partire dall'equazione di Schrödinger, l'evoluzione temporale di un sistema inizialmente nello stato ψ(0)\vert\psi(0)\rangle assume la forma:

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

dove HH è l'hamiltoniana indipendente dal tempo che governa il sistema. Consideriamo un'hamiltoniana che può essere scritta come una somma pesata di termini di Pauli H=jajPjH=\sum_j a_j P_j, con PjP_j che rappresenta un prodotto tensoriale di termini di Pauli che agiscono su nn qubit. In particolare, questi termini di Pauli possono commutare tra loro, oppure no. Dato uno stato al tempo t=0t=0, come si può ottenere lo stato del sistema in un tempo successivo ψ(t)|\psi(t)\rangle utilizzando un computer quantistico? L'esponenziale di un operatore può essere più facilmente compreso attraverso la sua serie di Taylor:

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

Alcuni esponenziali molto semplici, come eiZe^{iZ}, possono essere implementati facilmente sui computer quantistici utilizzando un insieme compatto di porte quantistiche. La maggior parte degli hamiltoniani di interesse non avrà un solo termine, ma molti termini. Si noti cosa succede se 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+...

Quando H1H_1 e H2H_2 commutano, si ha il caso familiare (che vale anche per i numeri e le variabili aa e bb di seguito):

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

Tuttavia, quando gli operatori non si interscambiano, i termini non possono essere riorganizzati nella serie di Taylor per semplificare in questo modo. Pertanto, esprimere Hamiltoniani complicati in porte quantistiche è una sfida.

Una soluzione consiste nel considerare un tempo molto piccolo tt, tale da rendere dominante il termine del primo ordine nell'espansione di Taylor. In base a questa ipotesi:

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}

Naturalmente, potremmo aver bisogno di evolvere il nostro stato per un periodo di tempo più lungo. Questo si ottiene utilizzando molti piccoli passi nel tempo. Questo processo è chiamato trotterizzazione:

ψ(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{,}

Qui t/rt/r è la fetta di tempo (passo evolutivo) che stiamo scegliendo. Di conseguenza, viene creato un cancello da applicare rr volte. Un passo temporale minore porta a un'approssimazione più accurata. Tuttavia, questo porta anche a circuiti più profondi che, in pratica, comportano un maggiore accumulo di errori (una preoccupazione non trascurabile per i dispositivi quantistici a breve termine).

Oggi studieremo l'evoluzione temporale del modello di Ising su reticoli lineari di siti N=2N=2 e N=6N=6. Questi reticoli sono costituiti da una serie di spin σi\sigma_i che interagiscono solo con i loro vicini più prossimi. Questi spin possono avere due orientamenti: \uparrow e \downarrow, che corrispondono a una magnetizzazione di +1+1 e 1-1 rispettivamente.

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{,}

dove JJ descrive l'energia di interazione e hh la grandezza di un campo esterno (nella direzione x di cui sopra, ma che modificheremo). Scriviamo questa espressione utilizzando le matrici di Pauli e considerando che il campo esterno ha un angolo α\alpha rispetto alla direzione trasversale,

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{.}

Questa hamiltoniana è utile perché ci permette di studiare facilmente gli effetti di un campo esterno. Nella base computazionale, il sistema sarà codificato come segue:

Stato quantistico
Rappresentazione dello 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

Inizieremo a studiare l'evoluzione temporale di un sistema quantistico di questo tipo. In particolare, visualizzeremo l'evoluzione temporale di alcune proprietà del sistema, come la magnetizzazione.

1.1 Requisiti

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

  • Qiskit SDK v1.2 o successivo ( pip install qiskit )
  • Qiskit Runtime v0.30 o più tardi ( pip install qiskit-ibm-runtime )
  • Numpy v1.24.1 o successivo < 2 ( pip install numpy )

1.2 Importa le librerie

Si noti che alcune librerie che potrebbero essere utili ( MatrixExponential, QDrift) sono incluse anche se non sono utilizzate in questo notebook. Se avete tempo, potete provarli!

# 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. Mappare il problema

2.1 Definizione dell'Hamiltoniano di Ising a campo trasversale

Consideriamo qui il modello di Ising a campo trasverso 1-D.

Per prima cosa, creeremo una funzione che accetta i parametri del sistema NN, JJ, hh e α\alpha, e restituisce la nostra hamiltoniana come SparsePauliOp. A SparsePauliOp è una rappresentazione rada di un operatore in termini di termini di Pauli ponderati.

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

Definire l'Hamiltoniano

Il sistema che consideriamo ora ha una dimensione di N=6N=6, J=0.2J=0.2, h=1.2h=1.2 e α=π8.0\alpha=\frac{\pi}{8.0} come esempio.

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 Imposta i parametri della simulazione dell'evoluzione temporale

Qui prenderemo in considerazione tre diverse tecniche di trotterizzazione:

  • Lie-Trotter (primo ordine)
  • suzuki-Trotter del secondo ordine
  • suzuki-Trotter del quarto ordine

Questi ultimi due saranno utilizzati nell'esercizio e nell'appendice.

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 Preparare il circuito quantistico 1 (stato iniziale)

Creare uno stato iniziale. Inizieremo con una configurazione di spin di \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 Preparare il circuito quantistico 2 (circuito singolo per l'evoluzione temporale)

Costruiamo qui un circuito per un singolo passo temporale utilizzando Lie-Trotter.

La formula del prodotto di Lie (del primo ordine) è implementata nella classe LieTrotter classe. Una formula del primo ordine consiste nell'approssimazione indicata nell'introduzione, in cui la matrice esponenziale di una somma è approssimata da un prodotto di matrici esponenziali:

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

Come già detto, i circuiti molto profondi portano all'accumulo di errori e causano problemi ai moderni computer quantistici. Poiché le porte a due equbit hanno tassi di errore più elevati rispetto a quelle a singolo equbit, una quantità di particolare interesse è la profondità del circuito a due equbit. Ciò che conta davvero è la profondità del circuito a due qubit dopo la transpilazione (poiché è il circuito che il computer quantistico esegue effettivamente). Ma prendiamo l'abitudine di contare le operazioni per questo circuito, anche ora utilizzando il simulatore.

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 Impostare gli operatori per misurare

Definiamo un operatore di magnetizzazione iZi/N\sum_i \langle Z_i \rangle / N e un operatore di correlazione media di spin 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 Eseguire una simulazione dell'evoluzione nel tempo

Controlleremo l'energia (valore di aspettativa dell'hamiltoniana), la magnetizzazione (valore di aspettativa dell'operatore di magnetizzazione) e la correlazione media di spin (valore di aspettativa dell'operatore di correlazione media di spin). La primitiva di Qiskit StatevectorEstimator ( EstimatorV2 ) stima i valori di aspettativa delle osservabili, ψ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 Tracciare l'evoluzione temporale delle grandezze osservabili

Tracciamo i valori di aspettativa misurati rispetto al tempo.

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. Esercizio 1. Eseguire la simulazione utilizzando il metodo Suzuki-Trotter di secondo ordine

Proviamo ora a eseguire una simulazione con Suzuki-Trotter del secondo ordine, seguendo l'esempio di Lie-Trotter illustrato in precedenza.

Il Suzuki-Trotter del secondo ordine può essere utilizzato in Qiskit per mezzo della classe SuzukiTrotter. Utilizzando questa formula, si ottiene una decomposizione del secondo ordine:

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

3.1 Costruisci un circuito per un singolo passo temporale

Utilizzare product_formula_st2 ( SuzukiTrotter(order=2 )) e costruire un circuito per un singolo passo temporale utilizzando Suzuki-Trotter del secondo ordine. Inoltre, contate il numero di porte e la profondità del circuito e confrontatelo con 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 Eseguire una simulazione dell'evoluzione nel tempo

Eseguire l'evoluzione temporale utilizzando Suzuki-Trotter del secondo ordine.

# 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 Tracciare i risultati di secondo ordine di Suzuki-Trotter

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 Confronta con risultati esatti

I dati riportati di seguito sono i risultati esatti precalcolati dal computer classico.

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. Esecuzione sull'hardware quantistico

Eseguiamo quindi la simulazione dell'evoluzione temporale sull'hardware quantistico. Lavoreremo su un problema più piccolo, di dimensione reticolare N=2. Variando il parametro α\alpha si osserva la differenza nella dinamica della funzione d'onda.

4.1 Passaggio 1. Mappare gli input classici su un problema quantistico

Scegliere l'impostazione iniziale della simulazione:

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

Quindi impostare il circuito iniziale:

La configurazione iniziale della rotazione sarà "down-up"

# 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

Ora si calcola il valore di riferimento utilizzando un simulatore di vettore di stato ideale.

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

Abbiamo preparato un sistema inizialmente con una sequenza di spin \downarrow\uparrow, che corrisponde a ψ(0)=10\vert\psi(0)\rangle = \vert10\rangle. Dopo averlo lasciato evolvere per t=1.6t=1.6 sotto un campo trasversale ( α=0\alpha=0^\circ ), abbiamo quasi la garanzia di misurare \uparrow\downarrow, cioè di avere uno scambio di spin. (Si noti che le etichette sono interpretate da destra a sinistra). Se il campo è longitudinale ( α=±90\alpha=\pm90^\circ ), non avremo alcuna evoluzione, quindi misureremo il sistema come è stato preparato inizialmente, \downarrow\uparrow. Con angoli intermedi, a α=±45\alpha=\pm45^\circ, potremo misurare tutte le combinazioni con diverse probabilità, essendo lo scambio di spin il più probabile con una probabilità del 67%.

Costruire un circuito per l'esperimento 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 Fase 2. Ottimizza per l'hardware di destinazione

Specifichiamo un backend.

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

Output:

'ibm_strasbourg'

Quindi si traspila il circuito per il backend selezionato.

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

Controllare il circuito.

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

Output:

Output of the previous code cell

4.3 Fase 3. Esegui con le primitive di IBM Quantum

La primitiva Sampler ( V2 ) di Qiskit fornisce il conteggio delle stringhe di bit misurate.

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

Output:

job id: d13pswfmya70008ek070

Salvare i risultati

results = job.result()

4.4 Fase 4. Risultati della post-elaborazione

Costruire l'istogramma delle stringhe di bit, che corrisponde all'analisi della funzione d'onda, e confrontarlo con i valori ideali mostrati sopra.

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

Qui mostriamo un esempio di costruzione di un circuito con Suzuki-Trotter di ordine superiore (quarto ordine). Proviamo ora a costruire una simulazione di circuito con Suzuki-Trotter del quarto ordine seguendo gli esempi illustrati sopra.

Il Suzuki-Trotter del quarto ordine può essere utilizzato in Qiskit per mezzo della classe SuzukiTrotter. Il quarto ordine può essere valutato utilizzando la seguente relazione di ricorsione. Si noti che l'ordine di Suzuki-Trotter è indicato come " 2k " nelle seguenti equazioni.

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)

Costruisci un circuito per un singolo passo temporale

Utilizzare product_formula_st4 ( SuzukiTrotter(order=4 )) e costruire un circuito per un singolo passo temporale utilizzando Suzuki-Trotter del quarto ordine. Inoltre, contate il numero di porte e la profondità del circuito e confrontatelo con Lie-Trotter e Suzuki-Trotter del secondo ordine.

# 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'
Questa pagina è stata utile?
Segnala un bug, un errore di battitura o richiedi contenuti su GitHub.