Dynamics Simulation¶

CUDA-Q enables the design, simulation and execution of quantum dynamics via the evolve API. Specifically, this API allows us to solve the time evolution of quantum systems or models. In the simulation mode, CUDA-Q provides the dynamics backend target, which is based on the cuQuantum library, optimized for performance and scale on NVIDIA GPU.

Quick Start¶

In the example below, we demonstrate a simple time evolution simulation workflow comprising of the following steps:

Define a quantum system model

A quantum system model is defined by a Hamiltonian. For example, a superconducting transmon qubit can be modeled by the following Hamiltonian

\[H = \frac{\omega_z}{2} \sigma_z + \omega_x \cos(\omega_d t)\sigma_x,\]

where \(\sigma_z\) and \(\sigma_x\) are Pauli Z and X operators, respectively.

Using CUDA-Q operator, the above time-dependent Hamiltonian can be set up as follows.

Python

omega_z = 6.5
omega_x = 4.0
omega_d = 0.5

import numpy as np
from cudaq import spin, ScalarOperator

# 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_d * t)) * spin.x(
    0)

In particular, ScalarOperator provides an easy way to model arbitrary time-dependent control signals. Details about CUDA-Q operator, including builtin operators that it supports can be found here.

C++

  // Parameters
  double omega_z = 6.5;
  double omega_x = 4.0;
  double omega_d = 0.5;

  // Qubit Hamiltonian
  auto hamiltonian = spin_op(0.5 * omega_z * spin_op::z(0));

  // Time dependent modulation
  auto mod_func =
      [omega_d](const parameter_map &params) -> std::complex<double> {
    auto it = params.find("t");
    if (it != params.end()) {
      double t = it->second.real();
      const auto result = std::cos(omega_d * t);
      return result;
    }
    throw std::runtime_error("Cannot find the time parameter.");
  };

  hamiltonian += mod_func * spin_op::x(0) * omega_x;

Details about CUDA-Q operator, including builtin operators that it supports can be found here.

Setup the evolution simulation

The below code snippet shows how to simulate the time-evolution of the above system with cudaq.evolve.

Python

import cudaq
import cupy as cp
from cudaq.dynamics import Schedule

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

# Dimensions of sub-systems: a single 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, t_final, n_steps)
schedule = Schedule(steps, ["t"])

# Run the simulation.
evolution_result = cudaq.evolve(hamiltonian,
                                dimensions,
                                schedule,
                                rho0,
                                observables=[spin.x(0),
                                             spin.y(0),
                                             spin.z(0)],
                                collapse_operators=[],
                                store_intermediate_results=True)

C++

  double t_final = 1.0;
  int n_steps = 100;

  // Define dimensions of subsystem (single two-level system)
  std::map<int, int> dimensions = {{0, 2}};

  // Initial state (ground state density matrix)
  auto rho0 = std::vector<std::complex<double>>{1.0, 0.0};

  // Schedule of time steps
  std::vector<double> steps(n_steps);
  for (int i = 0; i < n_steps; i++)
    steps[i] = i * t_final / (n_steps - 1);

  schedule schedule(steps, {"t"});

  // Numerical integrator
  // Here we choose a Runge-`Kutta` method for time evolution.
  cudaq::integrators::runge_kutta integrator(4, 0.01);

  // Observables to track
  auto observables = {cudaq::spin_op::x(0), cudaq::spin_op::y(0),
                      cudaq::spin_op::z(0)};

  // Run simulation
  // We evolve the system under the defined Hamiltonian. No collapsed operators
  // are provided (closed system evolution). The evolution returns expectation
  // values for all defined observables at each time step.
  auto evolution_result = cudaq::evolve(hamiltonian, dimensions, schedule, rho0,
                                        integrator, {}, observables, true);

Specifically, we need to set up the simulation by providing:

The system model in terms of a Hamiltonian as well as any decoherence terms, so-called collapse_operators.
The dimensionality of component systems in the model. CUDA-Q evolve allows users to model arbitrary multi-level systems, such as photonic Fock space.
The initial quantum state.
The time schedule, aka time steps, of the evolution.
Any ‘observable’ operator that we want to measure the expectation value with respect to the evolving state.

Note

By default, evolve will only return the final state and expectation values. To save intermediate results (at each time step specified in the schedule), the store_intermediate_results flag must be set to True.

Retrieve and plot the results

Once the simulation is complete, we can retrieve the final state and the expectation values as well as intermediate values at each time step (with store_intermediate_results=True).

For example, we can plot the Pauli expectation value for the above simulation as follows.

Python

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

import matplotlib.pyplot as plt

plt.plot(steps, get_result(0, evolution_result))
plt.plot(steps, get_result(1, evolution_result))
plt.plot(steps, get_result(2, evolution_result))
plt.ylabel("Expectation value")
plt.xlabel("Time")
plt.legend(("Sigma-X", "Sigma-Y", "Sigma-Z"))

In particular, for each time step, evolve captures an array of expectation values, one for each observable. Hence, we convert them into sequences for plotting purposes.

C++

  // Extract and print results
  for (size_t i = 0; i < steps.size(); i++) {
    double ex = evolution_result.expectation_values()[0][i].expectation();
    double ey = evolution_result.expectation_values()[1][i].expectation();
    double ez = evolution_result.expectation_values()[2][i].expectation();
    std::cout << steps[i] << " " << ex << " " << ey << " " << ez << "\n";
  }

Examples that illustrate how to use the dynamics target are available in the CUDA-Q repository.

Operator¶

CUDA-Q provides builtin definitions for commonly-used operators, such as the ladder operators (\(a\) and \(a^\dagger\)) of a harmonic oscillator, the Pauli spin operators for a two-level system, etc.

Here is a list of those operators.

Builtin Operators¶
Name	Description
`identity`	Identity operator
`zero`	Zero or null operator
`annihilate`	Bosonic annihilation operator (\(a\))
`create`	Bosonic creation operator (\(a^\dagger\))
`number`	Number operator of a bosonic mode (equivalent to \(a^\dagger a\))
`parity`	Parity operator of a bosonic mode (defined as \(e^{i\pi a^\dagger a}\))
`displace`	Displacement operator of complex amplitude \(\alpha\) (`displacement`). It is defined as \(e^{\alpha a^\dagger - \alpha^* a}\).
`squeeze`	Squeezing operator of complex squeezing amplitude \(z\) (`squeezing`). It is defined as \(\exp(\frac{1}{2}(z^*a^2 - z a^{\dagger 2}))\).
`position`	Position operator (equivalent to \((a^\dagger + a)/2\))
`momentum`	Momentum operator (equivalent to \(i(a^\dagger - a)/2\))
`spin.x`	Pauli \(\sigma_x\) operator
`spin.y`	Pauli \(\sigma_y\) operator
`spin.z`	Pauli \(\sigma_z\) operator
`spin.plus`	Pauli raising (\(\sigma_+\)) operator
`spin.minus`	Pauli lowering (\(\sigma_-\)) operator

As an example, let’s look at the Jaynes-Cummings model, which describes the interaction between a two-level atom and a light (Boson) field.

Mathematically, the Hamiltonian can be expressed as

\[H = \omega_c a^\dagger a + \omega_a \frac{\sigma_z}{2} + \frac{\Omega}{2}(a\sigma_+ + a^\dagger \sigma_-).\]

This Hamiltonian can be converted to CUDA-Q Operator representation with

Python

hamiltonian = omega_c * operators.create(1) * operators.annihilate(1) \
                + (omega_a / 2) * spin.z(0) \
                + (Omega / 2) * (operators.annihilate(1) * spin.plus(0) + operators.create(1) * spin.minus(0))

C++

  // Jaynes-Cummings Hamiltonian
  auto jc_hamiltonian =
      omega_c * boson_op::create(1) * boson_op::annihilate(1) +
      (omega_a / 2.0) * spin_op::z(0) +
      (Omega / 2.0) * (boson_op::annihilate(1) * spin_op::plus(0) +
                       boson_op::create(1) * spin_op::minus(0));

In the above code snippet, we map the cavity light field to degree index 1 and the two-level atom to degree index 0. The description of composite quantum system dynamics is independent from the Hilbert space of the system components. The latter is specified by the dimension map that is provided to the cudaq.evolve call.

Time-Dependent Dynamics¶

In the previous examples of operator construction, we assumed that the systems under consideration were described by time-independent Hamiltonian. However, we may want to simulate systems whose Hamiltonian operators have explicit time dependence.

CUDA-Q provides multiple ways to construct time-dependent operators.

Time-dependent coefficient

CUDA-Q ScalarOperator can be used to wrap a Python/C++ function that returns the coefficient value at a specific time.

As an example, we will look at a time-dependent Hamiltonian of the form \(H = H_0 + f(t)H_1\), where \(f(t)\) is the time-dependent driving strength given as \(cos(\omega t)\).

The following code sets up the problem

Python

# Define the static (drift) and control terms
H0 = spin.z(0)
H1 = spin.x(0)
H = H0 + ScalarOperator(lambda t: np.cos(omega * t)) * H1

C++

  // Hamiltonian with driving frequency
  double omega = M_PI;
  auto H0 = spin_op::z(0);
  auto H1 = spin_op::x(0);
  auto func = [omega](double t) { return std::cos(omega * t); };
  auto mod_func = [omega](double t) -> double { return std::cos(omega * t); };
  auto driven_hamiltonian = H0 + mod_func * H1;

Time-dependent operator

We can also construct a time-dependent operator from a function that returns a complex matrix representing the time dynamics of that operator.

As an example, let’s looks at the displacement operator. It can be defined as follows:

Python

import numpy
import scipy
from cudaq import operators, NumericType
from numpy.typing import NDArray


def displacement_matrix(
        dimension: int,
        displacement: NumericType) -> NDArray[numpy.complexfloating]:
    """
    Returns the displacement operator matrix.
    Args:
        displacement: Amplitude of the displacement operator.
            See also https://en.wikipedia.org/wiki/Displacement_operator.
    """
    displacement = complex(displacement)
    term1 = displacement * operators.create(0).to_matrix({0: dimension})
    term2 = numpy.conjugate(displacement) * operators.annihilate(0).to_matrix(
        {0: dimension})
    return scipy.linalg.expm(term1 - term2)


# The second argument here indicates the defined operator
# acts on a single degree of freedom, which can have any dimension.
# An argument [2], for example, would indicate that it can only
# act on a single degree of freedom with dimension two.
operators.define("displace", [0], displacement_matrix)


def displacement(degree: int) -> operators.MatrixOperatorElement:
    """
    Instantiates a displacement operator acting on the given degree of freedom.
    """
    return operators.instantiate("displace", [degree])

C++

  auto displacement_matrix =
      [](const std::vector<int> &dimensions,
         const cudaq::parameter_map &parameters) -> cudaq::complex_matrix {
    // Returns the displacement operator matrix.
    //  Args:
    //   - displacement: Amplitude of the displacement operator.
    // See also https://en.wikipedia.org/wiki/Displacement_operator.
    std::size_t dimension = dimensions[0];
    auto entry = parameters.find("displacement");
    if (entry == parameters.end())
      throw std::runtime_error("missing value for parameter 'displacement'");
    auto displacement_amplitude = entry->second;
    auto create = cudaq::complex_matrix(dimension, dimension);
    auto annihilate = cudaq::complex_matrix(dimension, dimension);
    for (std::size_t i = 0; i + 1 < dimension; i++) {
      create[{i + 1, i}] = std::sqrt(static_cast<double>(i + 1));
      annihilate[{i, i + 1}] = std::sqrt(static_cast<double>(i + 1));
    }
    auto term1 = displacement_amplitude * create;
    auto term2 = std::conj(displacement_amplitude) * annihilate;
    return (term1 - term2).exponential();
  };

  cudaq::matrix_handler::define("displace_op", {-1}, displacement_matrix);

  // Instantiate a displacement operator acting on the given degree of freedom.
  auto displacement = [](int degree) {
    return cudaq::matrix_handler::instantiate("displace_op", {degree});
  };

The defined operator is parameterized by the displacement amplitude. To create simulate the evolution of an operator under a time dependent displacement amplitude, we can define how the amplitude changes in time:

Python

import cudaq

# Define a system consisting of a single degree of freedom (0) with dimension 3.
system_dimensions = {0: 3}
system_operator = displacement(0)

# Define the time dependency of the system operator as a schedule that linearly
# increases the displacement parameter from 0 to 1.
time_dependence = Schedule(numpy.linspace(0, 1, 100), ['displacement'])
initial_state = cudaq.State.from_data(
    numpy.ones(3, dtype=numpy.complex128) / numpy.sqrt(3))

# Simulate the evolution of the system under this time dependent operator.
cudaq.evolve(system_operator, system_dimensions, time_dependence, initial_state)

C++

    // Define a system consisting of a single degree of freedom (0) with
    // dimension 3.
    cudaq::dimension_map system_dimensions{{0, 3}};
    auto system_operator = displacement(0);

    // Define the time dependency of the system operator as a schedule that
    // linearly increases the displacement parameter from 0 to 1.
    cudaq::schedule time_dependence(cudaq::linspace(0, 1, 100),
                                    {"displacement"});
    const std::vector<std::complex<double>> state_vec(3, 1.0 / std::sqrt(3.0));
    auto initial_state = cudaq::state::from_data(state_vec);
    cudaq::integrators::runge_kutta integrator(4, 0.01);
    // Simulate the evolution of the system under this time dependent operator.
    cudaq::evolve(system_operator, system_dimensions, time_dependence,
                  initial_state, integrator);

Let’s say we want to add a squeezing term to the system operator. We can independently vary the squeezing amplitude and the displacement amplitude by instantiating a schedule with a custom function that returns the desired value for each parameter:

Python

system_operator = displacement(0) + operators.squeeze(0)


# Define a schedule such that displacement amplitude increases linearly in time
# but the squeezing amplitude decreases, that is follows the inverse schedule.
def parameter_values(time_steps):

    def compute_value(param_name, step_idx):
        match param_name:
            case 'displacement':
                return time_steps[int(step_idx)]
            case 'squeezing':
                return time_steps[-int(step_idx + 1)]
            case _:
                raise ValueError(f"value for parameter {param_name} undefined")

    return Schedule(range(len(time_steps)), system_operator.parameters.keys(),
                    compute_value)


time_dependence = parameter_values(numpy.linspace(0, 1, 100))
cudaq.evolve(system_operator, system_dimensions, time_dependence, initial_state)

C++

    auto system_operator = displacement(0) + cudaq::matrix_op::squeeze(0);

    // Define a schedule such that displacement amplitude increases linearly in
    // time but the squeezing amplitude decreases, that is follows the inverse
    // schedule. def parameter_values(time_steps):
    auto parameter_values = [](const std::vector<double> &time_steps) {
      auto compute_value = [time_steps](const std::string &param_name,
                                        const std::complex<double> &step) {
        int step_idx = (int)step.real();
        if (param_name == "displacement")
          return time_steps[step_idx];
        if (param_name == "squeezing")
          return time_steps[time_steps.size() - (step_idx + 1)];

        throw std::runtime_error("value for parameter " + param_name +
                                 " undefined");
      };

      std::vector<std::complex<double>> steps;
      for (int i = 0; i < time_steps.size(); ++i)
        steps.emplace_back(i);
      return cudaq::schedule(steps, {"displacement", "squeezing"},
                             compute_value);
    };

    auto time_dependence = parameter_values(cudaq::linspace(0, 1, 100));
    cudaq::evolve(system_operator, system_dimensions, time_dependence,
                  initial_state, integrator);

Compile and Run C++ program

C++

nvq++ --target dynamics dynamics.cpp -o dynamics && ./dynamics

Numerical Integrators¶

For Python, CUDA-Q provides a set of numerical integrators, to be used with the dynamics backend target.

Numerical Integrators¶
Name	Description
`RungeKuttaIntegrator`	Explicit 4th-order Runge-Kutta method (default integrator)
`ScipyZvodeIntegrator`	Complex-valued variable-coefficient ordinary differential equation solver (provided by SciPy)
`CUDATorchDiffEqDopri5Integrator`	Runge-Kutta of order 5 of Dormand-Prince-Shampine (provided by `torchdiffeq`)
`CUDATorchDiffEqAdaptiveHeunIntegrator`	Runge-Kutta of order 2 (provided by `torchdiffeq`)
`CUDATorchDiffEqBosh3Integrator`	Runge-Kutta of order 3 of Bogacki-Shampine (provided by `torchdiffeq`)
`CUDATorchDiffEqDopri8Integrator`	Runge-Kutta of order 8 of Dormand-Prince-Shampine (provided by `torchdiffeq`)
`CUDATorchDiffEqEulerIntegrator`	Euler method (provided by `torchdiffeq`)
`CUDATorchDiffEqExplicitAdamsIntegrator`	Explicit Adams-Bashforth method (provided by `torchdiffeq`)
`CUDATorchDiffEqImplicitAdamsIntegrator`	Implicit Adams-Bashforth-Moulton method (provided by `torchdiffeq`)
`CUDATorchDiffEqMidpointIntegrator`	Midpoint method (provided by `torchdiffeq`)
`CUDATorchDiffEqRK4Integrator`	Fourth-order Runge-Kutta with 3/8 rule (provided by `torchdiffeq`)

Note

To use Torch-based integrators, users need to install torchdiffeq (e.g., with pip install torchdiffeq). This is an optional dependency of CUDA-Q, thus will not be installed by default.

Note

If you are using CUDA 12.8 on Blackwell, you may need to install nightly torch.

See Blackwell Torch Dependencies for more information.

Warning

Torch-based integrators require a CUDA-enabled Torch installation. Depending on your platform (e.g., aarch64), the default Torch pip package may not have CUDA support.

The below command can be used to verify your installation:

python3 -c "import torch; print(torch.version.cuda)"

If the output is a ‘None’ string, it indicates that your Torch installation does not support CUDA. In this case, you need to install a CUDA-enabled Torch package via other mechanisms, e.g., building Torch from source or using their Docker images.

For C++, CUDA-Q provides Runge-Kutta integrator, to be used with the dynamics backend target.

Numerical Integrators¶
Name	Description
`runge_kutta`	1st-order (Euler method), 2nd-order (Midpoint method), and 4th-order (classical Runge-Kutta method).

Multi-GPU Multi-Node Execution¶

CUDA-Q dynamics target supports parallel execution on multiple GPUs. To enable parallel execution, the application must initialize MPI as follows.

Python

cudaq.mpi.initialize()

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

# Initial state (expressed as an enum)
psi0 = cudaq.dynamics.InitialState.ZERO

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

cudaq.mpi.finalize()

mpiexec -np <N> python3 program.py

where N is the number of processes.

By initializing the MPI execution environment (via cudaq.mpi.initialize()) in the application code and invoking it via an MPI launcher, we have activated the multi-node multi-GPU feature of the dynamics target. Specifically, it will detect the number of processes (GPUs) and distribute the computation across all available GPUs.

Note

The number of MPI processes must be a power of 2, one GPU per process.

Note

Not all integrators are capable of handling distributed state. Errors will be raised if parallel execution is activated but the selected integrator does not support distributed state.

Warning

As of cuQuantum version 24.11, there are a couple of known limitations for parallel execution:

Computing the expectation value of a mixed quantum state is not supported. Thus, collapse_operators are not supported if expectation calculation is required.
Some combinations of quantum states and quantum many-body operators are not supported. Errors will be raised in those cases.

Examples¶

The Dynamics Examples section of the docs contains a number of excellent dynamics examples demonstrating how to simulate basic physics models, specific qubit modalities, and utilize multi-GPU multi-Node capabilities.