CUDA-Q Dynamics

Cavity QED

import cudaq
from cudaq import operators, Schedule, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# This example demonstrate a simulation of cavity quantum electrodynamics (interaction between light confined in a reflective cavity and atoms)

# System dimensions: atom (2-level system) and cavity (10-level system)
dimensions = {0: 2, 1: 10}

# Alias for commonly used operators
# Cavity operators
a = operators.annihilate(1)
a_dag = operators.create(1)

# Atom operators
sm = operators.annihilate(0)
sm_dag = operators.create(0)

# Defining the Hamiltonian for the system: self-energy terms and cavity-atom interaction term.
# This is the so-called Jaynes-Cummings model:
# https://en.wikipedia.org/wiki/Jaynes%E2%80%93Cummings_model
hamiltonian = 2 * np.pi * operators.number(1) + 2 * np.pi * operators.number(
    0) + 2 * np.pi * 0.25 * (sm * a_dag + sm_dag * a)

# Initial state of the system
# Atom in ground state
qubit_state = cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128)

# Cavity in a state which has 5 photons initially
cavity_state = cp.zeros((10, 10), dtype=cp.complex128)
cavity_state[5][5] = 1.0
rho0 = cudaq.State.from_data(cp.kron(qubit_state, cavity_state))

steps = np.linspace(0, 10, 201)
schedule = Schedule(steps, ["time"])

# First, evolve the system without any collapse operators (ideal).
evolution_result = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[operators.number(1), operators.number(0)],
    collapse_operators=[],
    store_intermediate_results=True,
    integrator=ScipyZvodeIntegrator())

# Then, evolve the system with a collapse operator modeling cavity decay (leaking photons)
evolution_result_decay = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[operators.number(1), operators.number(0)],
    collapse_operators=[np.sqrt(0.1) * a],
    store_intermediate_results=True,
    integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]
ideal_results = [
    get_result(0, evolution_result),
    get_result(1, evolution_result)
]
decay_results = [
    get_result(0, evolution_result_decay),
    get_result(1, evolution_result_decay)
]

fig = plt.figure(figsize=(18, 6))

plt.subplot(1, 2, 1)
plt.plot(steps, ideal_results[0])
plt.plot(steps, ideal_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Atom Excitation Probability"))
plt.title("No decay")

plt.subplot(1, 2, 2)
plt.plot(steps, decay_results[0])
plt.plot(steps, decay_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Cavity Photon Number", "Atom Excitation Probability"))
plt.title("With decay")

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig('cavity_qed.png', dpi=fig.dpi)

Cross Resonance

import cudaq
from cudaq import spin, Schedule, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# This example simulates cross-resonance interactions between superconducting qubits.
# Cross-resonance interaction is key to implementing two-qubit conditional gates, e.g., the CNOT gate.
# Ref: A simple all-microwave entangling gate for fixed-frequency superconducting qubits (Physical Review Letters 107, 080502)
# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Device parameters
# Detuning between two qubits
delta = 100 * 2 * np.pi
# Static coupling between qubits
J = 7 * 2 * np.pi
# spurious electromagnetic `crosstalk` due to stray electromagnetic coupling in the device circuit and package
# see (Physical Review Letters 107, 080502)
m_12 = 0.2
# Drive strength
Omega = 20 * 2 * np.pi

# Qubit Hamiltonian (in the rotating frame w.r.t. the target qubit)
hamiltonian = delta / 2 * spin.z(0) + J * (
    spin.minus(1) * spin.plus(0) +
    spin.plus(1) * spin.minus(0)) + Omega * spin.x(0) + m_12 * Omega * spin.x(1)

# Dimensions of sub-system
dimensions = {0: 2, 1: 2}

# Initial state of the system (ground state).
rho0 = cudaq.State.from_data(
    cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128))

# Two initial states: |00> and |10>.
# We show the 'conditional' evolution when controlled qubit is in |1> state.
psi_00 = cudaq.State.from_data(
    cp.array([1.0, 0.0, 0.0, 0.0], dtype=cp.complex128))
psi_10 = cudaq.State.from_data(
    cp.array([0.0, 0.0, 1.0, 0.0], dtype=cp.complex128))

# Schedule of time steps.
steps = np.linspace(0.0, 1.0, 1001)
schedule = Schedule(steps, ["time"])

# Run the simulation.
# Control bit = 0
evolution_result_00 = cudaq.evolve(hamiltonian,
                                   dimensions,
                                   schedule,
                                   psi_00,
                                   observables=[
                                       spin.x(0),
                                       spin.y(0),
                                       spin.z(0),
                                       spin.x(1),
                                       spin.y(1),
                                       spin.z(1)
                                   ],
                                   collapse_operators=[],
                                   store_intermediate_results=True,
                                   integrator=ScipyZvodeIntegrator())

# Control bit = 1
evolution_result_10 = cudaq.evolve(hamiltonian,
                                   dimensions,
                                   schedule,
                                   psi_10,
                                   observables=[
                                       spin.x(0),
                                       spin.y(0),
                                       spin.z(0),
                                       spin.x(1),
                                       spin.y(1),
                                       spin.z(1)
                                   ],
                                   collapse_operators=[],
                                   store_intermediate_results=True,
                                   integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]
results_00 = [
    get_result(0, evolution_result_00),
    get_result(1, evolution_result_00),
    get_result(2, evolution_result_00),
    get_result(3, evolution_result_00),
    get_result(4, evolution_result_00),
    get_result(5, evolution_result_00)
]
results_10 = [
    get_result(0, evolution_result_10),
    get_result(1, evolution_result_10),
    get_result(2, evolution_result_10),
    get_result(3, evolution_result_10),
    get_result(4, evolution_result_10),
    get_result(5, evolution_result_10)
]

# The changes in recession frequencies of the target qubit when control qubit is in |1> state allow us to implement two-qubit conditional gates.
fig = plt.figure(figsize=(18, 6))

plt.subplot(1, 2, 1)
plt.plot(steps, results_00[5])
plt.plot(steps, results_10[5])
plt.ylabel(r"$\langle Z_2\rangle$")
plt.xlabel("Time")
plt.legend((r"$|\psi_0\rangle=|00\rangle$", r"$|\psi_0\rangle=|10\rangle$"))

plt.subplot(1, 2, 2)
plt.plot(steps, results_00[4])
plt.plot(steps, results_10[4])
plt.ylabel(r"$\langle Y_2\rangle$")
plt.xlabel("Time")
plt.legend((r"$|\psi_0\rangle=|00\rangle$", r"$|\psi_0\rangle=|10\rangle$"))

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig('cross_resonance.png', dpi=fig.dpi)

Gate Calibration

import cudaq
from cudaq import operators, Schedule, ScalarOperator, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# This example demonstrates the use of the dynamics simulator and optimizer to optimize a pulse.
# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Sample device parameters
# Assuming a simple transmon device Hamiltonian in rotating frame.
detuning = 0.0  # Detuning of the drive; assuming resonant drive
anharmonicity = -340.0  # Anharmonicity
sigma = 0.01  # sigma of the Gaussian pulse
cutoff = 4.0 * sigma  # total length of drive pulse

# Dimensions of sub-system
# We model `transmon` as a 3-level system to account for leakage.
dimensions = {0: 3}

# Initial state of the system (ground state).
psi0 = cudaq.State.from_data(cp.array([1.0, 0.0, 0.0], dtype=cp.complex128))


def gaussian(t):
    """
    Gaussian shape with cutoff. Starts at t = 0, amplitude normalized to one
    """
    val = (np.exp(-((t-cutoff/2)/sigma)**2/2)-np.exp(-(cutoff/sigma)**2/8)) \
           / (1-np.exp(-(cutoff/sigma)**2/8))
    return val


def dgaussian(t):
    """
    Derivative of Gaussian. Starts at t = 0, amplitude normalized to one
    """
    return -(t - cutoff / 2) / sigma * np.exp(-(
        (t - cutoff / 2) / sigma)**2 / 2 + 0.5)


# Schedule of time steps.
steps = np.linspace(0.0, cutoff, 201)
schedule = Schedule(steps, ["t"])

# We optimize for a X(pi/2) rotation
target_state = np.array([1.0 / np.sqrt(2), -1j / np.sqrt(2), 0.0],
                        dtype=cp.complex128)


# Optimize the amplitude of the drive pulse (DRAG - Derivative Removal by Adiabatic Gate)
def cost_function(amps):
    amplitude = 100 * amps[0]
    drag_amp = 100 * amps[1]
    # Qubit Hamiltonian
    hamiltonian = detuning * operators.number(0) + (
        anharmonicity / 2) * operators.create(0) * operators.create(
            0) * operators.annihilate(0) * operators.annihilate(0)
    # Drive term
    hamiltonian += amplitude * ScalarOperator(gaussian) * (
        operators.create(0) + operators.annihilate(0))

    # Drag term (leakage reduction)
    hamiltonian += 1j * drag_amp * ScalarOperator(dgaussian) * (
        operators.annihilate(0) - operators.create(0))

    # We optimize for a X(pi/2) rotation
    evolution_result = cudaq.evolve(hamiltonian,
                                    dimensions,
                                    schedule,
                                    psi0,
                                    observables=[],
                                    collapse_operators=[],
                                    store_intermediate_results=False,
                                    integrator=ScipyZvodeIntegrator())
    final_state = evolution_result.final_state()

    overlap = np.abs(final_state.overlap(target_state))
    print(
        f"Gaussian amplitude = {amplitude}, derivative amplitude = {drag_amp}, Overlap: {overlap}"
    )
    return 1.0 - overlap


# Specify the optimizer
optimizer = cudaq.optimizers.NelderMead()
optimal_error, optimal_parameters = optimizer.optimize(dimensions=2,
                                                       function=cost_function)

print("optimal overlap =", 1.0 - optimal_error)
print("optimal parameters =", optimal_parameters)

Heisenberg Model

import cudaq
from cudaq import operators, spin, Schedule, ScipyZvodeIntegrator

import numpy as np
import cupy as cp
import matplotlib.pyplot as plt
import os

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# In this example, we solve the Quantum Heisenberg model (https://en.wikipedia.org/wiki/Quantum_Heisenberg_model),
# which exhibits the so-called quantum quench effect.
# e.g., see `Quantum quenches in the anisotropic spin-1/2 Heisenberg chain: different approaches to many-body dynamics far from equilibrium`
# (New J. Phys. 12 055017)
# Number of spins
N = 9
dimensions = {}
for i in range(N):
    dimensions[i] = 2

# Initial state: alternating spin up and down
spin_state = ''
for i in range(N):
    spin_state += str(int(i % 2))

# Observable is the staggered magnetization operator
staggered_magnetization_op = operators.zero()
for i in range(N):
    if i % 2 == 0:
        staggered_magnetization_op += spin.z(i)
    else:
        staggered_magnetization_op -= spin.z(i)

staggered_magnetization_op /= N

observe_results = []
for g in [0.0, 0.25, 4.0]:
    # Heisenberg model spin coupling strength
    Jx = 1.0
    Jy = 1.0
    Jz = g

    # Construct the Hamiltonian
    H = operators.zero()

    for i in range(N - 1):
        H += Jx * spin.x(i) * spin.x(i + 1)
        H += Jy * spin.y(i) * spin.y(i + 1)
        H += Jz * spin.z(i) * spin.z(i + 1)

    steps = np.linspace(0.0, 5, 1000)
    schedule = Schedule(steps, ["time"])

    # Prepare the initial state vector
    psi0_ = cp.zeros(2**N, dtype=cp.complex128)
    psi0_[int(spin_state, 2)] = 1.0
    psi0 = cudaq.State.from_data(psi0_)

    # Run the simulation
    evolution_result = cudaq.evolve(H,
                                    dimensions,
                                    schedule,
                                    psi0,
                                    observables=[staggered_magnetization_op],
                                    collapse_operators=[],
                                    store_intermediate_results=True,
                                    integrator=ScipyZvodeIntegrator())

    exp_val = [
        exp_vals[0].expectation()
        for exp_vals in evolution_result.expectation_values()
    ]

    observe_results.append((g, exp_val))

# Plot the results
fig = plt.figure(figsize=(12, 6))
for g, exp_val in observe_results:
    plt.plot(steps, exp_val, label=f'$ g = {g}$')
plt.legend(fontsize=16)
plt.ylabel("Staggered Magnetization")
plt.xlabel("Time")
abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig("heisenberg_model.png", dpi=fig.dpi)

Landau Zener

import cudaq
from cudaq import spin, operators, ScalarOperator, Schedule, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# This example simulates the so-called Landau–Zener transition: given a time-dependent Hamiltonian such that the energy separation
# of the two states is a linear function of time, an analytical formula exists to calculate
# the probability of finding the system in the excited state after the transition.

# References:
# - https://en.wikipedia.org/wiki/Landau%E2%80%93Zener_formula
# - `The Landau-Zener formula made simple`, `Eric P Glasbrenner and Wolfgang P Schleich 2023 J. Phys. B: At. Mol. Opt. Phys. 56 104001`
# - QuTiP notebook: https://github.com/qutip/qutip-notebooks/blob/master/examples/landau-zener.ipynb

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Define some shorthand operators
sx = spin.x(0)
sz = spin.z(0)
sm = operators.annihilate(0)
sm_dag = operators.create(0)

# Dimensions of sub-system. We only have a single degree of freedom of dimension 2 (two-level system).
dimensions = {0: 2}

# Landau–Zener Hamiltonian:
# `[[-alpha*t, g], [g, alpha*t]] = g * pauli_x - alpha * t * pauli_z`
g = 2 * np.pi
# Analytical equation:
# `P(0) = exp(-pi * g ^ 2/ alpha)`
# The target ground state probability that we want to achieve
target_p0 = 0.75
# Compute `alpha` parameter:
alpha = (-np.pi * g**2) / np.log(target_p0)

# Hamiltonian
hamiltonian = g * sx - alpha * ScalarOperator(lambda t: t) * sz

# Initial state of the system (ground state)
psi0 = cudaq.State.from_data(cp.array([1.0, 0.0], dtype=cp.complex128))

# Schedule of time steps (simulating a long time range)
steps = np.linspace(-2.0, 2.0, 5000)
schedule = Schedule(steps, ["t"])

# Run the simulation.
evolution_result = cudaq.evolve(hamiltonian,
                                dimensions,
                                schedule,
                                psi0,
                                observables=[operators.number(0)],
                                collapse_operators=[],
                                store_intermediate_results=True,
                                integrator=ScipyZvodeIntegrator())

prob1 = [
    exp_vals[0].expectation()
    for exp_vals in evolution_result.expectation_values()
]

prob0 = [1 - val for val in prob1]
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(steps, prob1, 'b', steps, prob0, 'r')
ax.plot(steps, (1.0 - target_p0) * np.ones(np.shape(steps)), 'k')
ax.plot(steps, target_p0 * np.ones(np.shape(steps)), 'm')
ax.set_xlabel("Time")
ax.set_ylabel("Occupation probability")
ax.set_title("Landau-Zener transition")
ax.legend(("Excited state", "Ground state", "LZ formula (Excited state)",
           "LZ formula (Ground state)"),
          loc=0)

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
plt.savefig('landau_zener.png', dpi=fig.dpi)

Pulse

import cudaq
from cudaq import spin, operators, ScalarOperator, Schedule, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# This example simulates time evolution of a qubit (`transmon`) being driven by a pulse.
# The pulse is a modulated signal with a Gaussian envelop.
# The simulation is performed in the 'lab' frame.

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Device parameters
# Strength of the Rabi-rate in GHz.
rabi_rate = 0.1

# Frequency of the qubit transition in GHz.
omega = 5.0 * 2 * np.pi

# Define Gaussian envelope function to approximately implement a `rx(pi/2)` gate.
amplitude = 1. / 2.0  # Pi/2 rotation
sigma = 1.0 / rabi_rate / amplitude
pulse_duration = 6 * sigma


def gaussian(t, duration, amplitude, sigma):
    # Gaussian envelope function
    return amplitude * np.exp(-0.5 * (t - duration / 2)**2 / (sigma)**2)


def signal(t):
    # Modulated signal
    return np.cos(omega * t) * gaussian(t, pulse_duration, amplitude, sigma)


# Qubit Hamiltonian
hamiltonian = omega * spin.z(0) / 2
# Add modulated driving term to the Hamiltonian
hamiltonian += np.pi * rabi_rate * ScalarOperator(signal) * spin.x(0)

# Dimensions of sub-system. We only have a single degree of freedom of dimension 2 (two-level system).
dimensions = {0: 2}

# Initial state of the system (ground state).
psi0 = cudaq.State.from_data(cp.array([1.0, 0.0], dtype=cp.complex128))

# Schedule of time steps.
# Since this is a lab-frame simulation, the time step must be small to accurately capture the modulated signal.
dt = 1 / omega / 20
n_steps = int(np.ceil(pulse_duration / dt)) + 1
steps = np.linspace(0, pulse_duration, n_steps)
schedule = Schedule(steps, ["t"])

# Run the simulation.
# First, we run the simulation without any collapse operators (no decoherence).
evolution_result = cudaq.evolve(hamiltonian,
                                dimensions,
                                schedule,
                                psi0,
                                observables=[operators.number(0)],
                                collapse_operators=[],
                                store_intermediate_results=True,
                                integrator=ScipyZvodeIntegrator())

pop1 = [
    exp_vals[0].expectation()
    for exp_vals in evolution_result.expectation_values()
]
pop0 = [1.0 - x for x in pop1]
fig = plt.figure(figsize=(6, 16))
envelop = [gaussian(t, pulse_duration, amplitude, sigma) for t in steps]

plt.subplot(3, 1, 1)
plt.plot(steps, envelop)
plt.ylabel("Amplitude")
plt.xlabel("Time")
plt.title("Envelope")

modulated_signal = [
    np.cos(omega * t) * gaussian(t, pulse_duration, amplitude, sigma)
    for t in steps
]
plt.subplot(3, 1, 2)
plt.plot(steps, modulated_signal)
plt.ylabel("Amplitude")
plt.xlabel("Time")
plt.title("Signal")

plt.subplot(3, 1, 3)
plt.plot(steps, pop0)
plt.plot(steps, pop1)
plt.ylabel("Population")
plt.xlabel("Time")
plt.legend(("Population in |0>", "Population in |1>"))
plt.title("Qubit State")

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig("pulse.png", dpi=fig.dpi)

Qubit Control

import cudaq
from cudaq import spin, ScalarOperator, Schedule, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# This example simulates time evolution of a qubit (`transmon`) being driven close to resonance in the presence of noise (decoherence).
# Thus, it exhibits Rabi oscillations.
# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Qubit Hamiltonian reference: https://qiskit-community.github.io/qiskit-dynamics/tutorials/Rabi_oscillations.html
# Device parameters
# Qubit resonant frequency
omega_z = 10.0 * 2 * np.pi
# Transverse term
omega_x = 2 * np.pi
# Harmonic driving frequency
# Note: we chose a frequency slightly different from the resonant frequency to demonstrate the off-resonance effect.
omega_drive = 0.99 * omega_z

# Qubit Hamiltonian
hamiltonian = 0.5 * omega_z * spin.z(0)
# Add modulated driving term to the Hamiltonian
hamiltonian += omega_x * ScalarOperator(
    lambda t: np.cos(omega_drive * t)) * spin.x(0)

# Dimensions of sub-system. We only have a single degree of freedom of dimension 2 (two-level system).
dimensions = {0: 2}

# Initial state of the system (ground state).
rho0 = cudaq.State.from_data(
    cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128))

# Schedule of time steps.
t_final = np.pi / omega_x
dt = 2.0 * np.pi / omega_drive / 100
n_steps = int(np.ceil(t_final / dt)) + 1
steps = np.linspace(0, t_final, n_steps)
schedule = Schedule(steps, ["t"])

# Run the simulation.
# First, we run the simulation without any collapse operators (no decoherence).
evolution_result = cudaq.evolve(hamiltonian,
                                dimensions,
                                schedule,
                                rho0,
                                observables=[spin.x(0),
                                             spin.y(0),
                                             spin.z(0)],
                                collapse_operators=[],
                                store_intermediate_results=True,
                                integrator=ScipyZvodeIntegrator())

# Now, run the simulation with qubit decoherence
gamma_sm = 4.0
gamma_sz = 1.0
evolution_result_decay = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[spin.x(0), spin.y(0), spin.z(0)],
    collapse_operators=[
        np.sqrt(gamma_sm) * spin.plus(0),
        np.sqrt(gamma_sz) * spin.z(0)
    ],
    store_intermediate_results=True,
    integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]
ideal_results = [
    get_result(0, evolution_result),
    get_result(1, evolution_result),
    get_result(2, evolution_result)
]
decoherence_results = [
    get_result(0, evolution_result_decay),
    get_result(1, evolution_result_decay),
    get_result(2, evolution_result_decay)
]

fig = plt.figure(figsize=(18, 6))

plt.subplot(1, 2, 1)
plt.plot(steps, ideal_results[0])
plt.plot(steps, ideal_results[1])
plt.plot(steps, ideal_results[2])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Sigma-X", "Sigma-Y", "Sigma-Z"))
plt.title("No decoherence")

plt.subplot(1, 2, 2)
plt.plot(steps, decoherence_results[0])
plt.plot(steps, decoherence_results[1])
plt.plot(steps, decoherence_results[2])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Sigma-X", "Sigma-Y", "Sigma-Z"))
plt.title("With decoherence")

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig('qubit_control.png', dpi=fig.dpi)

Qubit Dynamics

import cudaq
from cudaq import spin, Schedule, RungeKuttaIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Qubit Hamiltonian
hamiltonian = 2 * np.pi * 0.1 * spin.x(0)

# Dimensions of sub-system. We only have a single degree of freedom of dimension 2 (two-level system).
dimensions = {0: 2}

# Initial state of the system (ground state).
rho0 = cudaq.State.from_data(
    cp.array([[1.0, 0.0], [0.0, 0.0]], dtype=cp.complex128))

# Schedule of time steps.
steps = np.linspace(0, 10, 101)
schedule = Schedule(steps, ["time"])

# Run the simulation.
# First, we run the simulation without any collapse operators (ideal).
evolution_result = cudaq.evolve(hamiltonian,
                                dimensions,
                                schedule,
                                rho0,
                                observables=[spin.y(0), spin.z(0)],
                                collapse_operators=[],
                                store_intermediate_results=True,
                                integrator=RungeKuttaIntegrator())

# Now, run the simulation with qubit decaying due to the presence of a collapse operator.
evolution_result_decay = cudaq.evolve(
    hamiltonian,
    dimensions,
    schedule,
    rho0,
    observables=[spin.y(0), spin.z(0)],
    collapse_operators=[np.sqrt(0.05) * spin.x(0)],
    store_intermediate_results=True,
    integrator=RungeKuttaIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]
ideal_results = [
    get_result(0, evolution_result),
    get_result(1, evolution_result)
]
decay_results = [
    get_result(0, evolution_result_decay),
    get_result(1, evolution_result_decay)
]

fig = plt.figure(figsize=(18, 6))

plt.subplot(1, 2, 1)
plt.plot(steps, ideal_results[0])
plt.plot(steps, ideal_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Sigma-Y", "Sigma-Z"))
plt.title("No decay")

plt.subplot(1, 2, 2)
plt.plot(steps, decay_results[0])
plt.plot(steps, decay_results[1])
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Sigma-Y", "Sigma-Z"))
plt.title("With decay")

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig('qubit_dynamics.png', dpi=fig.dpi)

Silicon Spin Qubit

import cudaq
from cudaq import spin, operators, Schedule, ScalarOperator, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt
import matplotlib as mpl

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# This example demonstrates simulation of an electrically-driven silicon spin qubit.
# The system dynamics is taken from https://journals.aps.org/prapplied/pdf/10.1103/PhysRevApplied.19.044078

dimensions = {0: 2}
resonance_frequency = 2 * np.pi * 10  # 10 Ghz

# Run the simulation:
evolution_results = []
get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]
# Sweep the amplitude
amplitudes = np.linspace(0.0, 0.5, 20)
for amplitude in amplitudes:
    rho0 = cudaq.State.from_data(cp.array([1.0, 0.0], dtype=cp.complex128))
    t_final = 100
    dt = 0.005
    n_steps = int(np.ceil(t_final / dt)) + 1
    steps = np.linspace(0, t_final, n_steps)
    schedule = Schedule(steps, ["t"])
    # Electric dipole spin resonance (`EDSR`) Hamiltonian
    H = 0.5 * resonance_frequency * spin.z(0) + ScalarOperator(
        lambda t: 0.5 * np.sin(resonance_frequency * t) * amplitude) * spin.x(0)
    evolution_result = cudaq.evolve(H,
                                    dimensions,
                                    schedule,
                                    rho0,
                                    observables=[operators.number(0)],
                                    collapse_operators=[],
                                    store_intermediate_results=True,
                                    integrator=ScipyZvodeIntegrator())
    evolution_results.append(get_result(0, evolution_result))

fig, ax = plt.subplots()
im = ax.contourf(steps, amplitudes, evolution_results)
ax.set_xlabel("Time (ns)")
ax.set_ylabel(f"Amplitude (a.u.)")
fig.suptitle(f"Excited state probability")
fig.colorbar(im)
abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
# For reference, see figure 5 in https://journals.aps.org/prapplied/pdf/10.1103/PhysRevApplied.19.044078
fig.savefig("spin_qubit_edsr.png", dpi=fig.dpi)

Tensor Callback

import cudaq
from cudaq import ElementaryOperator, spin, operators, Schedule, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# This example demonstrates the use of callback functions to define time-dependent operators.

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Consider a simple 2-level system Hamiltonian exhibits the Landau–Zener transition:
# `[[-alpha*t, g], [g, alpha*t]]
# This can be defined as a callback tensor:
g = 2.0 * np.pi
alpha = 10.0 * 2 * np.pi


def callback_tensor(t):
    return np.array([[-alpha * t, g], [g, alpha * t]], dtype=np.complex128)


# Analytical formula
lz_formula_p0 = np.exp(-np.pi * g**2 / (alpha))
lz_formula_p1 = 1.0 - lz_formula_p0

# Let's define the control term as a callback tensor that acts on 2-level systems
ElementaryOperator.define("lz_op", [2], callback_tensor)

# Hamiltonian
hamiltonian = ElementaryOperator("lz_op", [0])

# Dimensions of sub-system. We only have a single degree of freedom of dimension 2 (two-level system).
dimensions = {0: 2}

# Initial state of the system (ground state)
psi0 = cudaq.State.from_data(cp.array([1.0, 0.0], dtype=cp.complex128))

# Schedule of time steps (simulating a long time range)
steps = np.linspace(-4.0, 4.0, 10000)
schedule = Schedule(steps, ["t"])

# Run the simulation.
evolution_result = cudaq.evolve(hamiltonian,
                                dimensions,
                                schedule,
                                psi0,
                                observables=[operators.number(0)],
                                collapse_operators=[],
                                store_intermediate_results=True,
                                integrator=ScipyZvodeIntegrator())

prob1 = [
    exp_vals[0].expectation()
    for exp_vals in evolution_result.expectation_values()
]

prob0 = [1 - val for val in prob1]
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(steps, prob1, 'b', steps, prob0, 'r')
ax.plot(steps, lz_formula_p1 * np.ones(np.shape(steps)), 'k')
ax.plot(steps, lz_formula_p0 * np.ones(np.shape(steps)), 'm')
ax.set_xlabel("Time")
ax.set_ylabel("Occupation probability")
ax.set_title("Landau-Zener transition")
ax.legend(("Excited state", "Ground state", "LZ formula (Excited state)",
           "LZ formula (Ground state)"),
          loc=0)

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig('tensor_callback.png', dpi=fig.dpi)

Transmon Resonator

import cudaq
from cudaq import operators, spin, Schedule, ScipyZvodeIntegrator
import numpy as np
import cupy as cp
import os
import matplotlib.pyplot as plt

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# This example demonstrates a simulation of a superconducting transmon qubit coupled to a resonator (i.e., cavity).
# References:
# - "Charge-insensitive qubit design derived from the Cooper pair box", PRA 76, 042319
# - QuTiP lecture: https://github.com/jrjohansson/qutip-lectures/blob/master/Lecture-10-cQED-dispersive-regime.ipynb

# Number of cavity photons
N = 20

# System dimensions: transmon + cavity
dimensions = {0: 2, 1: N}

# See III.B of PRA 76, 042319
# System parameters
# Unit: GHz
omega_01 = 3.0 * 2 * np.pi  # transmon qubit frequency
omega_r = 2.0 * 2 * np.pi  # resonator frequency
# Dispersive shift
chi_01 = 0.025 * 2 * np.pi
chi_12 = 0.0

omega_01_prime = omega_01 + chi_01
omega_r_prime = omega_r - chi_12 / 2.0
chi = chi_01 - chi_12 / 2.0

# System Hamiltonian
hamiltonian = 0.5 * omega_01_prime * spin.z(0) + (
    omega_r_prime + chi * spin.z(0)) * operators.number(1)

# Initial state of the system
# Transmon in a superposition state
transmon_state = cp.array([1. / np.sqrt(2.), 1. / np.sqrt(2.)],
                          dtype=cp.complex128)


# Helper to create a coherent state in Fock basis truncated at `num_levels`.
# Note: There are a couple of ways of generating a coherent state,
# e.g., see https://qutip.readthedocs.io/en/v5.0.3/apidoc/functions.html#qutip.core.states.coherent
# or https://en.wikipedia.org/wiki/Coherent_state
# Here, in this example, we use a the formula: `|alpha> = D(alpha)|0>`,
# i.e., apply the displacement operator on a zero (or vacuum) state to compute the corresponding coherent state.
def coherent_state(num_levels, amplitude):
    displace_mat = operators.displace(0).to_matrix({0: num_levels},
                                                   displacement=amplitude)
    # `D(alpha)|0>` is the first column of `D(alpha)` matrix
    return cp.array(np.transpose(displace_mat)[0])


# Cavity in a coherent state
cavity_state = coherent_state(N, 2.0)
psi0 = cudaq.State.from_data(cp.kron(transmon_state, cavity_state))

steps = np.linspace(0, 250, 1000)
schedule = Schedule(steps, ["time"])

# Evolve the system
evolution_result = cudaq.evolve(hamiltonian,
                                dimensions,
                                schedule,
                                psi0,
                                observables=[
                                    operators.number(1),
                                    operators.number(0),
                                    operators.position(1),
                                    operators.position(0)
                                ],
                                collapse_operators=[],
                                store_intermediate_results=True,
                                integrator=ScipyZvodeIntegrator())

get_result = lambda idx, res: [
    exp_vals[idx].expectation() for exp_vals in res.expectation_values()
]
count_results = [
    get_result(0, evolution_result),
    get_result(1, evolution_result)
]

quadrature_results = [
    get_result(2, evolution_result),
    get_result(3, evolution_result)
]

fig = plt.figure(figsize=(18, 6))

plt.subplot(1, 2, 1)
plt.plot(steps, count_results[0])
plt.plot(steps, count_results[1])
plt.ylabel("n")
plt.xlabel("Time [ns]")
plt.legend(("Cavity Photon Number", "Transmon Excitation Probability"))
plt.title("Excitation Numbers")

plt.subplot(1, 2, 2)
plt.plot(steps, quadrature_results[0])
plt.ylabel("x")
plt.xlabel("Time [ns]")
plt.legend(("cavity"))
plt.title("Resonator Quadrature")

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
fig.savefig('transmon_resonator.png', dpi=fig.dpi)

Initial State (Multi-GPU Multi-Node)

import cudaq
from cudaq import operators, spin, Schedule, RungeKuttaIntegrator

import numpy as np
import matplotlib.pyplot as plt
import os

# On a system with multiple GPUs, `mpiexec` can be used as follows:
# `mpiexec -np <N> python3 multi_gpu.py `
cudaq.mpi.initialize()

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# Large number of spins
N = 20
dimensions = {}
for i in range(N):
    dimensions[i] = 2

# Observable is the average magnetization operator
avg_magnetization_op = operators.zero()
for i in range(N):
    avg_magnetization_op += (spin.z(i) / N)

# Arbitrary coupling constant
g = 1.0
# Construct the Hamiltonian
H = operators.zero()
for i in range(N):
    H += 2 * np.pi * spin.x(i)
    H += 2 * np.pi * spin.y(i)
for i in range(N - 1):
    H += 2 * np.pi * g * spin.x(i) * spin.x(i + 1)
    H += 2 * np.pi * g * spin.y(i) * spin.z(i + 1)

steps = np.linspace(0.0, 1, 200)
schedule = Schedule(steps, ["time"])

# Initial state (expressed as an enum)
psi0 = cudaq.operator.InitialState.ZERO
# This can also be used to initialize a uniformly-distributed wave-function instead.
# `psi0 = cudaq.operator.InitialState.UNIFORM`

# Run the simulation
evolution_result = cudaq.evolve(H,
                                dimensions,
                                schedule,
                                psi0,
                                observables=[avg_magnetization_op],
                                collapse_operators=[],
                                store_intermediate_results=True,
                                integrator=RungeKuttaIntegrator())

exp_val = [
    exp_vals[0].expectation()
    for exp_vals in evolution_result.expectation_values()
]

if cudaq.mpi.rank() == 0:
    # Plot the results
    fig = plt.figure(figsize=(12, 6))
    plt.plot(steps, exp_val)
    plt.ylabel("Average Magnetization")
    plt.xlabel("Time")
    abspath = os.path.abspath(__file__)
    dname = os.path.dirname(abspath)
    os.chdir(dname)
    fig.savefig("spin_model.png", dpi=fig.dpi)

cudaq.mpi.finalize()

Heisenberg Model (Multi-GPU Multi-Node)

import cudaq
from cudaq import operators, spin, Schedule, RungeKuttaIntegrator

import numpy as np
import cupy as cp
import matplotlib.pyplot as plt
import os

# On a system with multiple GPUs, `mpiexec` can be used as follows:
# `mpiexec -np <N> python3 multi_gpu.py `
cudaq.mpi.initialize()

# Set the target to our dynamics simulator
cudaq.set_target("dynamics")

# In this example, we solve the Quantum Heisenberg model (https://en.wikipedia.org/wiki/Quantum_Heisenberg_model),
# which exhibits the so-called quantum quench effect.
# e.g., see `Quantum quenches in the anisotropic spin-1/2 Heisenberg chain: different approaches to many-body dynamics far from equilibrium`
# (New J. Phys. 12 055017)
# Large number of spins
N = 21
dimensions = {}
for i in range(N):
    dimensions[i] = 2

# Initial state: alternating spin up and down
spin_state = ''
for i in range(N):
    spin_state += str(int(i % 2))

# Observable is the staggered magnetization operator
staggered_magnetization_op = operators.zero()
for i in range(N):
    if i % 2 == 0:
        staggered_magnetization_op += spin.z(i)
    else:
        staggered_magnetization_op -= spin.z(i)

staggered_magnetization_op /= N

observe_results = []
for g in [0.25, 4.0]:
    # Heisenberg model spin coupling strength
    Jx = 1.0
    Jy = 1.0
    Jz = g

    # Construct the Hamiltonian
    H = operators.zero()

    for i in range(N - 1):
        H += Jx * spin.x(i) * spin.x(i + 1)
        H += Jy * spin.y(i) * spin.y(i + 1)
        H += Jz * spin.z(i) * spin.z(i + 1)

    steps = np.linspace(0.0, 1, 100)
    schedule = Schedule(steps, ["time"])

    # Prepare the initial state vector
    psi0_ = cp.zeros(2**N, dtype=cp.complex128)
    psi0_[int(spin_state, 2)] = 1.0
    psi0 = cudaq.State.from_data(psi0_)

    # Run the simulation
    evolution_result = cudaq.evolve(H,
                                    dimensions,
                                    schedule,
                                    psi0,
                                    observables=[staggered_magnetization_op],
                                    collapse_operators=[],
                                    store_intermediate_results=True,
                                    integrator=RungeKuttaIntegrator())

    exp_val = [
        exp_vals[0].expectation()
        for exp_vals in evolution_result.expectation_values()
    ]

    observe_results.append((g, exp_val))

if cudaq.mpi.rank() == 0:
    # Plot the results
    fig = plt.figure(figsize=(12, 6))
    for g, exp_val in observe_results:
        plt.plot(steps, exp_val, label=f'$ g = {g}$')
    plt.legend(fontsize=16)
    plt.ylabel("Staggered Magnetization")
    plt.xlabel("Time")
    abspath = os.path.abspath(__file__)
    dname = os.path.dirname(abspath)
    os.chdir(dname)
    fig.savefig("heisenberg_model.png", dpi=fig.dpi)

cudaq.mpi.finalize()