Hello AHS: Run your first Analog Hamiltonian Simulation
Interacting spin chain
For a canonical example of a system of many interacting particles, let us consider a ring of eight spins (each of which can be in “up” ∣↑⟩ and “down” ∣↓⟩ states). Albeit small, this model system already exhibits a handful of interesting phenomena of naturally occurring magnetic materials. In this example, we will show how to prepare a so-called anti-ferromagnetic order, where consecutive spins point in opposite directions.
![Diagram connecting 8 circle nodes that contain inversing up and down arrows.](/images/braket/latest/developerguide/images/AntiFerromagnetic.png)
Arrangement
We will use one neutral atom to stand for each spin, and the “up” and “down” spin states will be encoded in excited Rydberg state and ground state of the atoms, respectively. First, we create the 2-d arrangement. We can program the above ring of spins with the following code.
Prerequisites: You need to pip install the Braket SDKpip install
matplotlib
.
import numpy as np import matplotlib.pyplot as plt # required for plotting from braket.ahs.atom_arrangement import AtomArrangement a = 5.7e-6 # nearest-neighbor separation (in meters) register = AtomArrangement() register.add(np.array([0.5, 0.5 + 1/np.sqrt(2)]) * a) register.add(np.array([0.5 + 1/np.sqrt(2), 0.5]) * a) register.add(np.array([0.5 + 1/np.sqrt(2), - 0.5]) * a) register.add(np.array([0.5, - 0.5 - 1/np.sqrt(2)]) * a) register.add(np.array([-0.5, - 0.5 - 1/np.sqrt(2)]) * a) register.add(np.array([-0.5 - 1/np.sqrt(2), - 0.5]) * a) register.add(np.array([-0.5 - 1/np.sqrt(2), 0.5]) * a) register.add(np.array([-0.5, 0.5 + 1/np.sqrt(2)]) * a)
which we can also plot with
fig, ax = plt.subplots(1, 1, figsize=(7,7)) xs, ys = [register.coordinate_list(dim) for dim in (0, 1)] ax.plot(xs, ys, 'r.', ms=15) for idx, (x, y) in enumerate(zip(xs, ys)): ax.text(x, y, f" {idx}", fontsize=12) plt.show() # this will show the plot below in an ipython or jupyter session
![Scatter plot showing points distributed across positive and negative values on both axes.](/images/braket/latest/developerguide/images/PlotNeutralAtoms.png)
Interaction
To prepare the anti-ferromagnetic phase, we need to induce interactions between
neighboring spins. We use the van der Waals
interaction
![Hamiltonian interaction equation showing this this interaction as expressed as a sum over all spin pairs (j,k).](/images/braket/latest/developerguide/images/HInteraction.png)
Here, nj=∣↑j⟩⟨↑j∣ is an operator that takes the value of 1 only if spin j is in the “up” state, and 0 otherwise. The strength is Vj,k=C6/(dj,k)6, where C6 is the fixed coefficient, and dj,k is the Euclidean distance between spins j and k. The immediate effect of this interaction term is that any state where both spin j and spin k are “up” have elevated energy (by the amount Vj,k). By carefully designing the rest of the AHS program, this interaction will prevent neighboring spins from both being in the “up” state, an effect commonly known as "Rydberg blockade."
Driving field
At the beginning of the AHS program, all spins (by default) start in their “down” state, they are in a so-called ferromagnetic phase. Keeping an eye on our goal to prepare the anti-ferromagnetic phase, we specify a time-dependent coherent driving field that smoothly transitions the spins from this state to a many-body state where the “up” states are preferred. The corresponding Hamiltonian can be written as
![Mathematical equation depicting the calculation of a Hamiltonian drive function.](/images/braket/latest/developerguide/images/HDrive.png)
where Ω(t),ϕ(t),Δ(t) are the time-dependent, global amplitude (aka Rabi frequency
To program a smooth transition from the ferromagnetic phase to the anti-ferromagnetic phase, we specify the driving field with the following code.
from braket.timings.time_series import TimeSeries from braket.ahs.driving_field import DrivingField # smooth transition from "down" to "up" state time_max = 4e-6 # seconds time_ramp = 1e-7 # seconds omega_max = 6300000.0 # rad / sec delta_start = -5 * omega_max delta_end = 5 * omega_max omega = TimeSeries() omega.put(0.0, 0.0) omega.put(time_ramp, omega_max) omega.put(time_max - time_ramp, omega_max) omega.put(time_max, 0.0) delta = TimeSeries() delta.put(0.0, delta_start) delta.put(time_ramp, delta_start) delta.put(time_max - time_ramp, delta_end) delta.put(time_max, delta_end) phi = TimeSeries().put(0.0, 0.0).put(time_max, 0.0) drive = DrivingField( amplitude=omega, phase=phi, detuning=delta )
We can visualize the time series of the driving field with the following script.
fig, axes = plt.subplots(3, 1, figsize=(12, 7), sharex=True) ax = axes[0] time_series = drive.amplitude.time_series ax.plot(time_series.times(), time_series.values(), '.-'); ax.grid() ax.set_ylabel('Omega [rad/s]') ax = axes[1] time_series = drive.detuning.time_series ax.plot(time_series.times(), time_series.values(), '.-'); ax.grid() ax.set_ylabel('Delta [rad/s]') ax = axes[2] time_series = drive.phase.time_series # Note: time series of phase is understood as a piecewise constant function ax.step(time_series.times(), time_series.values(), '.-', where='post'); ax.set_ylabel('phi [rad]') ax.grid() ax.set_xlabel('time [s]') plt.show() # this will show the plot below in an ipython or jupyter session
![Three graphs showing phi, delta, and omega over time. The top subplot shows the growth to just above 6 rads/s where it stays for 4 seconds until it drops back to 0. The middle subplot depicts the associated linear growth of the derivative, and the bottom subplot illustrates a flat line near zero.](/images/braket/latest/developerguide/images/DrivingTimeSeries.png)
AHS program
The register, the driving field, (and the implicit van der Waals interactions) make
up the Analog Hamiltonian Simulation program ahs_program
.
from braket.ahs.analog_hamiltonian_simulation import AnalogHamiltonianSimulation ahs_program = AnalogHamiltonianSimulation( register=register, hamiltonian=drive )
Running on local simulator
Since this example is small (less than 15 spins), before running it on an AHS-compatible QPU, we can run it on the local AHS simulator which comes with the Braket SDK. Since the local simulator is available for free with the Braket SDK, this is best practice to ensure that our code can correctly execute.
Here, we can set the number of shots to a high value (say, 1 million) because the local simulator tracks the time evolution of the quantum state and draws samples from the final state; hence, increasing the number of shots, while increasing the total runtime only marginally.
from braket.devices import LocalSimulator device = LocalSimulator("braket_ahs") result_simulator = device.run( ahs_program, shots=1_000_000 ).result() # takes about 5 seconds
Analyzing simulator results
We can aggregate the shot results with the following function that infers the state of each spin (which may be “d” for “down”, “u” for “up”, or “e” for empty site), and counts how many times each configuration occurred across the shots.
from collections import Counter def get_counts(result): """Aggregate state counts from AHS shot results A count of strings (of length = # of spins) are returned, where each character denotes the state of a spin (site): e: empty site u: up state spin d: down state spin Args: result (braket.tasks.analog_hamiltonian_simulation_quantum_task_result.AnalogHamiltonianSimulationQuantumTaskResult) Returns dict: number of times each state configuration is measured """ state_counts = Counter() states = ['e', 'u', 'd'] for shot in result.measurements: pre = shot.pre_sequence post = shot.post_sequence state_idx = np.array(pre) * (1 + np.array(post)) state = "".join(map(lambda s_idx: states[s_idx], state_idx)) state_counts.update((state,)) return dict(state_counts) counts_simulator = get_counts(result_simulator) # takes about 5 seconds print(counts_simulator)
{'udududud': 330944, 'dudududu': 329576, 'dududdud': 38033, ...}
Here counts
is a dictionary that counts the number of times each state
configuration is observed across the shots. We can also visualize them with the
following code.
from collections import Counter def has_neighboring_up_states(state): if 'uu' in state: return True if state[0] == 'u' and state[-1] == 'u': return True return False def number_of_up_states(state): return Counter(state)['u'] def plot_counts(counts): non_blockaded = [] blockaded = [] for state, count in counts.items(): if not has_neighboring_up_states(state): collection = non_blockaded else: collection = blockaded collection.append((state, count, number_of_up_states(state))) blockaded.sort(key=lambda _: _[1], reverse=True) non_blockaded.sort(key=lambda _: _[1], reverse=True) for configurations, name in zip((non_blockaded, blockaded), ('no neighboring "up" states', 'some neighboring "up" states')): plt.figure(figsize=(14, 3)) plt.bar(range(len(configurations)), [item[1] for item in configurations]) plt.xticks(range(len(configurations))) plt.gca().set_xticklabels([item[0] for item in configurations], rotation=90) plt.ylabel('shots') plt.grid(axis='y') plt.title(f'{name} configurations') plt.show() plot_counts(counts_simulator)
![Bar chart showing a large number of shots with no neighboring "up" states configurations.](/images/braket/latest/developerguide/images/AHSCounts1.png)
![Bar chart showing shots of some neighboring "up" states configurations, with states descending from high to low values.](/images/braket/latest/developerguide/images/AHSCounts2.png)
From the plots, we can read the following observations the verify that we successfully prepared the anti-ferromagnetic phase.
-
Generally, non-blockaded states (where no two neighboring spins are in the “up” state) are more common than states where at least one pair of neighboring spins are both in “up” states.
-
Generally, states with more "up" excitations are favored, unless the configuration is blockaded.
-
The most common states are indeed the perfect anti-ferromagnetic states
"dudududu"
and"udududud"
. -
The second most common states are the ones where there is only 3 “up” excitations with consecutive separations of 1, 2, 2. This shows that the van der Waals interaction has an affect (albeit much smaller) on next-nearest neighbors too.
Running on QuEra’s Aquila QPU
Prerequisites: Apart from pip installing the Braket
SDK
Note
If you are using a Braket hosted notebook instance, the Braket SDK comes pre-installed with the instance.
With all dependencies installed, we can connect to the Aquila QPU.
from braket.aws import AwsDevice aquila_qpu = AwsDevice("arn:aws:braket:us-east-1::device/qpu/quera/Aquila")
To make our AHS program suitable for the QuEra machine, we need to
round all values to comply with the levels of precision allowed by the
Aquila QPU. (These requirements are governed by the device parameters
with “Resolution” in their name. We can see them by executing
aquila_qpu.properties.dict()
in a notebook. For more details of
capabilities and requirements of Aquila, see the Introduction to Aquiladiscretize
method.
discretized_ahs_program = ahs_program.discretize(aquila_qpu)
Now we can run the program (running only 100 shots for now) on the Aquila QPU.
Note
Running this program on the Aquila processor will incur a cost. The
Amazon Braket SDK includes a Cost Tracker
task = aquila_qpu.run(discretized_ahs_program, shots=100) metadata = task.metadata() task_arn = metadata['quantumTaskArn'] task_status = metadata['status'] print(f"ARN: {task_arn}") print(f"status: {task_status}")
task ARN: arn:aws:braket:us-east-1:123456789012:quantum-task/12345678-90ab-cdef-1234-567890abcdef task status: CREATED
Due to the large variance of how long a quantum task may take to run (depending on availability windows and QPU utilization), it is a good idea to note down the quantum task ARN, so we can check its status at a later time with the following code snippet.
# Optionally, in a new python session from braket.aws import AwsQuantumTask SAVED_TASK_ARN = "arn:aws:braket:us-east-1:123456789012:quantum-task/12345678-90ab-cdef-1234-567890abcdef" task = AwsQuantumTask(arn=SAVED_TASK_ARN) metadata = task.metadata() task_arn = metadata['quantumTaskArn'] task_status = metadata['status'] print(f"ARN: {task_arn}") print(f"status: {task_status}")
*[Output]* task ARN: arn:aws:braket:us-east-1:123456789012:quantum-task/12345678-90ab-cdef-1234-567890abcdef task status: COMPLETED
Once the status is COMPLETED (which can also be checked from the quantum tasks page of the
Amazon Braket console
result_aquila = task.result()
Analyzing QPU results
Using the same get_counts
functions as before, we can compute the
counts:
counts_aquila = get_counts(result_aquila) print(counts_aquila)
*[Output]* {'udududud': 24, 'dudududu': 17, 'dududdud': 3, ...}
and plot them with plot_counts
:
plot_counts(counts_aquila)
![Bar chart showing rate of shots with and without neighboring "up" states.](/images/braket/latest/developerguide/images/QPUPlotCounts.png)
Note that a small fraction of shots have empty sites (marked with “e”). This is due to a 1—2% per atom preparation imperfections of the Aquila QPU. Apart from this, the results match with the simulation within the expected statistical fluctuation due to small number of shots.
Next
Congratulations, you have now run your first AHS workload on Amazon Braket using the local AHS simulator and the Aquila QPU.
To learn more about Rydberg physics, Analog Hamiltonian Simulation and the
Aquila device, refer to our example notebooks