Skip to content
Pasqal Documentation

Getting started with emu-mps

import matplotlib.pyplot as plt
from emu_mps import (
MPS,
MPSConfig,
MPSBackend,
BitStrings,
Fidelity,
Occupation,
)
from utils_examples import afm_sequence_from_register, square_perimeter_points
import pulser
from pulser.devices import AnalogDevice
import numpy as np

Emu-mps is a Pulser backend, designed to Emulate the dynamics of programmable arrays of neutral atoms. It works in a very intuitive way: the main method is called run which is used to run the simulation. This method requires as argument a Pulser sequence and a configuration file.

Here are the steps that we will go through in the following:

  1. using Pulser we will create first the atomic Register
  2. again with Pulser we will then generate the Pulse sequence that produce the AFM state for the created register
  3. we create the configuration object (MPSConfig) and we fill it with the needed observables and time step
  4. we instantiate the backend class (MPSBackend) and run the simulation
  5. we inspect the results obtained
Omega_max = 2 * 2 * np.pi
delta_0 = -6 * Omega_max / 2
delta_f = 1 * Omega_max / 2
t_rise = 500
t_fall = 1500
sweep_factor = 2
square_length = 3
R_interatomic = AnalogDevice.rydberg_blockade_radius(Omega_max / 2)
coords = R_interatomic * square_perimeter_points(square_length)
reg = pulser.Register.from_coordinates(coords)
reg.draw(blockade_radius=R_interatomic, draw_graph=True, draw_half_radius=True)
seq = afm_sequence_from_register(
reg, Omega_max, delta_0, delta_f, t_rise, t_fall, sweep_factor, AnalogDevice
)
seq.draw("input")

As mentioned earlier, to run a simulation with the emu-mps backend we need to provide as input a Pulse sequence - which we have just created - and a configuration object.

We still have to create the configuration for the emu-mps backend. This is done via an instantiation of the configuration class MPSConfig which contains all the observables that we wish to measure and the time step chosen for the simulation, along with various algorithm specific parameters that are explained in the documentation.

We start by setting a bigger discretization time than the default one provided ($dt=10$) and enforcing that the times at which we compute the observables are an integer multiple of $dt$. For simplicity, we will measure the observables only at the final time of the simulation.

We also fix the basis along which the measurements will be done. For the details regarding the conventions used we refer to the Pulser documentation (external).

dt = 100
eval_times = [1.0]
basis = ("r","g")

It samples at desired time steps the evolved state, returning the bitStrings in a counter.

sampling_times = 1000
bitstrings = BitStrings(evaluation_times=eval_times, num_shots=sampling_times)

The fidelity is computed as $\langle \psi_{evolved} | \phi_{given} \rangle$ where

  • $\psi_{evolved}$ is the system state at desired steps
  • $\phi_{given}$ is the state in which we want to project the state system.

In this tutorial we will compute the fidelity against the dominant of the two antiferromagnetic states $\phi_{given} = |rgrgrgrg>$

nqubits = len(seq.register.qubit_ids)
afm_string_pure = {"rgrgrgrg": 1.0}
afm_mps_state = MPS.from_state_amplitudes(
eigenstates=basis, amplitudes=afm_string_pure
)
fidelity_mps_pure = Fidelity(evaluation_times=eval_times, state=afm_mps_state)

It is computed as $\langle \psi_{evolved} |\frac{(1+Z_i)}{2}|\psi_{evolved}\rangle$ and often informally referred to as the magnetization at each atom site.

density = Occupation(
evaluation_times=[x/seq.get_duration() for x in range(0, seq.get_duration(), dt)]
)
mpsconfig = MPSConfig(
dt=dt,
observables=[
bitstrings,
fidelity_mps_pure,
density,
],
)

We instantiate the backend class MPSBackend and use its run method where we pass as argument the sequence and the backend configuration objects.

sim = MPSBackend(seq, config=mpsconfig)
results = sim.run()

In the following lines, we are going to give a brief code examples of how you can get the information from the results object

results.get_result_tags()

Below, we retrieve the bistrings computed. We observe that they were indeed computed only at a single time $ns = 3400$, and we find those that were sampled the highest number of times with their relative counting.

results.get_result_times(bitstrings)
bitstrings_final = results.get_result(bitstrings, 1.0)
max_val = max(bitstrings_final.values()) # max number of counts in the bitstring
max_string = [key for key, value in bitstrings_final.items() if value == max_val]
print(
"The most frequent bitstring is {} which was sampled {} times".format(
max_string, max_val
)
)
filtered_counts = [count for count in bitstrings_final.values() if count > 20]
filtered_bitstrings = [
bitstring for bitstring, count in bitstrings_final.items() if count > 20
]
x_labels = range(len(filtered_bitstrings))
with plt.style.context("seaborn-v0_8-darkgrid"):
fig, ax = plt.subplots()
ax.bar(x_labels, filtered_counts, color="teal", alpha=0.8)
ax.set_xlabel("Bitstrings")
ax.set_ylabel("Counts")
ax.set_title("Histogram of Bitstring Counts (Counts > 20)")
ax.set_xticks(x_labels)
ax.set_xticklabels(filtered_bitstrings, rotation="vertical")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.tight_layout()
plt.show()

Here we compute the fidelity of the system against the two different AFM state realizations defined above.

fidelity_pure = results.get_result(fidelity_mps_pure,1.0)
print(
"The fidelity computed for the system final state against the pure state |rgrgrgr> is {}.\nThe probability of the system being in that sate is equal to {} ".format(
fidelity_pure, abs(fidelity_pure) ** 2
)
)

Here, we plot the time evolution of the magnetization of the system sites, and we observe how the system slowly reaches the AFM state.

magnetization_values = np.array(list(results.occupation))
magnetization_times = results.get_result_times(density)
fig, ax = plt.subplots(figsize=(8, 4), layout="constrained")
num_time_points, positions = magnetization_values.shape
x, y = np.meshgrid(np.arange(num_time_points), np.arange(1, positions + 1))
im = plt.pcolormesh(magnetization_times, y, magnetization_values.T, shading="auto")
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Qubit")
ax.set_title("State Density")
ax.set_yticks(np.arange(1, positions + 1))
cb = fig.colorbar(im, ax=ax)