Dynamics cpp integration

CMakeLists.txt

@@ -156,6 +156,9 @@ endif()
 if(NOT CUTENSORNET_ROOT)
   SET(CUTENSORNET_ROOT "$ENV{CUQUANTUM_INSTALL_PREFIX}")
 endif()
+if(NOT CUDENSITYMAT_ROOT)
+  SET(CUDENSITYMAT_ROOT "$ENV{CUQUANTUM_INSTALL_PREFIX}")
+endif()
 if(NOT CUTENSOR_ROOT)
   SET(CUTENSOR_ROOT "$ENV{CUTENSOR_INSTALL_PREFIX}")
 endif()

docker/build/assets.Dockerfile

@@ -150,7 +150,7 @@ RUN source /cuda-quantum/scripts/configure_build.sh && \
         libname="$(basename "$binary")" && \
         # Linking cublas dynamically necessarily adds a libgcc_s dependency to the GPU-based simulators.
         # The same holds for a range of CUDA libraries, whereas libcudart.so does not have any GCC dependencies.
-        if [ "${libname#libnvqir-custatevec}" != "$libname" ] || [ "${libname#libnvqir-tensornet}" != "$libname" ]; then \
+        if [ "${libname#libnvqir-custatevec}" != "$libname" ] || [ "${libname#libnvqir-tensornet}" != "$libname" ] || [ "${libname#libnvqir-dynamics}" != "$libname" ]; then \
             echo "Skipping validation of $libname."; \
         elif [ -n "$(ldd "${binary}" 2>/dev/null | grep gcc)" ]; then \
             has_gcc_dependencies=true && \

docs/sphinx/examples/cpp/dynamics/cavity_qed.cpp

@@ -0,0 +1,140 @@
+/*******************************************************************************
+ * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+// Compile and run with:
+// ```
+// nvq++ --target dynamics cavity_qed.cpp -o a.out && ./a.out
+// ```
+
+#include "cudaq/algorithms/evolve.h"
+#include "cudaq/dynamics_integrators.h"
+#include "cudaq/operators.h"
+#include "export_csv_helper.h"
+#include <cudaq.h>
+
+int main() {
+
+  // Dimension of our composite quantum system:
+  // subsystem 0 (atom) has 2 levels (ground and excited states).
+  // subsystem 1 (cavity) has 10 levels (Fock states, representing photon number
+  // states).
+  cudaq::dimension_map dimensions{{0, 2}, {1, 10}};
+
+  // For the cavity subsystem 1
+  // We create the annihilation (a) and creation (a+) operators.
+  // These operators lower and raise the photon number, respectively.
+  auto a = cudaq::boson_operator::annihilate(1);
+  auto a_dag = cudaq::boson_operator::create(1);
+
+  // For the atom subsystem 0
+  // We create the annihilation (`sm`) and creation (`sm_dag`) operators.
+  // These operators lower and raise the excitation number, respectively.
+  auto sm = cudaq::boson_operator::annihilate(0);
+  auto sm_dag = cudaq::boson_operator::create(0);
+
+  // Number operators
+  // These operators count the number of excitations.
+  // For the atom (`subsytem` 0) and the cavity (`subsystem` 1) they give the
+  // population in each subsystem.
+  auto atom_occ_op = cudaq::matrix_operator::number(0);
+  auto cavity_occ_op = cudaq::matrix_operator::number(1);
+
+  // Hamiltonian
+  // The `hamiltonian` models the dynamics of the atom-cavity (cavity QED)
+  // system It has 3 parts:
+  // 1. Cavity energy: 2 * pi * photon_number_operator -> energy proportional to
+  // the number of photons.
+  // 2. Atomic energy: 2 * pi * atom_number_operator -> energy proportional to
+  // the atomic excitation.
+  // 3. Atomic-cavity interaction: 2 * pi * 0.25 * (`sm` * a_dag + `sm_dag` * a)
+  // -> represents the exchange of energy between the atom and the cavity. It is
+  // analogous to the Jaynes-Cummings model in cavity QED.
+  auto hamiltonian = (2 * M_PI * cavity_occ_op) + (2 * M_PI * atom_occ_op) +
+                     (2 * M_PI * 0.25 * (sm * a_dag + sm_dag * a));
+
+  // Build the initial state
+  // Atom (sub-system 0) in ground state.
+  // Cavity (sub-system 1) has 5 photons (Fock space).
+  // The overall Hilbert space is 2 * 10 = 20.
+  const int num_photons = 5;
+  std::vector<std::complex<double>> initial_state_vec(20, 0.0);
+  // The index is chosen such that the atom is in the ground state while the
+  // cavity is in the Fock state with 5 photons.
+  initial_state_vec[dimensions[0] * num_photons] = 1;
+
+  // Define a time evolution schedule
+  // We define a time grid from 0 to 10 (in arbitrary time units) with 201 time
+  // steps. This schedule is used by the integrator to simulate the dynamics.
+  const int num_steps = 201;
+  std::vector<double> steps = cudaq::linspace(0.0, 10.0, num_steps);
+  cudaq::schedule schedule(steps);
+
+  // Create a CUDA quantum state
+  // The initial state is converted into a quantum state object for evolution.
+  auto rho0 = cudaq::state::from_data(initial_state_vec);
+
+  // Numerical integrator
+  // Here we choose a Runge-`Kutta` method for time evolution.
+  // `dt` defines the time step for the numerical integration, and order 4
+  // indicates a 4`th` order method.
+  cudaq::integrators::runge_kutta integrator(4, 0.01);
+
+  // Evolve without collapse operators
+  // This evolution is ideal (closed system) dynamics governed solely by the
+  // Hamiltonian. The expectation values of the observables (cavity photon
+  // number and atom excitation probability) are recorded.
+  cudaq::evolve_result evolve_result =
+      cudaq::evolve(hamiltonian, dimensions, schedule, rho0, integrator, {},
+                    {cavity_occ_op, atom_occ_op}, true);
+
+  // Adding dissipation
+  // To simulate a realistic scenario, we introduce decay (dissipation).
+  // Here, the collapse operator represents photon loss from the cavity.
+  constexpr double decay_rate = 0.1;
+  auto collapse_operator = std::sqrt(decay_rate) * a;
+  // Evolve with the collapse operator to incorporate the effect of decay.
+  cudaq::evolve_result evolve_result_decay =
+      cudaq::evolve(hamiltonian, dimensions, schedule, rho0, integrator,
+                    {collapse_operator}, {cavity_occ_op, atom_occ_op}, true);
+
+  // Lambda to extract expectation values for a given observable index
+  // Here, index 0 corresponds to the cavity photon number and index 1
+  // corresponds to the atom excitation probability.
+  auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
+    std::vector<double> expectations;
+
+    auto all_exps = result.get_expectation_values().value();
+    for (auto exp_vals : all_exps) {
+      expectations.push_back((double)exp_vals[idx]);
+    }
+    return expectations;
+  };
+
+  // Retrieve expectation values from both the ideal and decaying `evolutions`.
+  auto ideal_result0 = get_expectation(0, evolve_result);
+  auto ideal_result1 = get_expectation(1, evolve_result);
+  auto decay_result0 = get_expectation(0, evolve_result_decay);
+  auto decay_result1 = get_expectation(1, evolve_result_decay);
+
+  // Export the results to `CSV` files
+  // "cavity_`qed`_ideal_result.`csv`" contains the ideal evolution results.
+  // "cavity_`qed`_decay_result.`csv`" contains the evolution results with
+  // cavity decay.
+  export_csv("cavity_qed_ideal_result.csv", "time", steps,
+             "cavity_photon_number", ideal_result0,
+             "atom_excitation_probability", ideal_result1);
+  export_csv("cavity_qed_decay_result.csv", "time", steps,
+             "cavity_photon_number", decay_result0,
+             "atom_excitation_probability", decay_result1);
+
+  std::cout << "Simulation complete. The results are saved in "
+               "cavity_qed_ideal_result.csv "
+               "and cavity_qed_decay_result.csv files."
+            << std::endl;
+  return 0;
+}

docs/sphinx/examples/cpp/dynamics/cross_resonance.cpp

@@ -0,0 +1,164 @@
+/*******************************************************************************
+ * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+// Compile and run with:
+// ```
+// nvq++ --target dynamics cavity_qed.cpp -o a.out && ./a.out
+// ```
+
+#include "cudaq/algorithms/evolve.h"
+#include "cudaq/dynamics_integrators.h"
+#include "cudaq/operators.h"
+#include "export_csv_helper.h"
+#include <cudaq.h>
+
+int main() {
+
+  // `delta` represents the detuning between the two qubits.
+  // In physical terms, detuning is the energy difference (or frequency offset)
+  // between qubit levels. Detuning term (in angular frequency units).
+  double delta = 100 * 2 * M_PI;
+  // `J` is the static coupling strength between the two qubits.
+  // This terms facilitates energy exchange between the qubits, effectively
+  // coupling their dynamics.
+  double J = 7 * 2 * M_PI;
+  // `m_12` models spurious electromagnetic `crosstalk`.
+  // `Crosstalk` is an unwanted interaction , here represented as a fraction of
+  // the drive strength applied to the second qubit.
+  double m_12 = 0.2;
+  // `Omega` is the drive strength applied to the qubits.
+  // A driving field can induce transitions between qubit states.
+  double Omega = 20 * 2 * M_PI;
+
+  // For a spin-1/2 system, the raising operator S^+ and lowering operator S^-
+  // are defined as: S^+ = 0.5 * (X + `iY`) and S^- = 0.5 * (X - `iY`) These
+  // operators allow transitions between the spin states (|0> and |1>).
+  auto spin_plus = [](int degree) {
+    return 0.5 *
+           (cudaq::spin_operator::x(degree) +
+            std::complex<double>(0.0, 1.0) * cudaq::spin_operator::y(degree));
+  };
+
+  auto spin_minus = [](int degree) {
+    return 0.5 *
+           (cudaq::spin_operator::x(degree) -
+            std::complex<double>(0.0, 1.0) * cudaq::spin_operator::y(degree));
+  };
+
+  // The Hamiltonian describes the energy and dynamics of our 2-qubit system.
+  // It consist of several parts:
+  // 1. Detuning term for qubit 0: (delta / 2) * Z. This sets the energy
+  // splitting for qubit 0.
+  // 2. Exchange interaction: J * (S^-_1 * S^+_0 + S^+_1 * S^-_0). This couples
+  // the two qubits, enabling excitation transfer.
+  // 3. Drive on qubit 0: Omega * X. A control field that drives transition in
+  // qubit 0.
+  // 4. `Crosstalk` drive on qubit 1: m_12 * Omega * X. A reduces drive on qubit
+  // 1 due to electromagnetic `crosstalk`.
+  auto hamiltonian =
+      (delta / 2.0) * cudaq::spin_operator::z(0) +
+      J * (spin_minus(1) * spin_plus(0) + spin_plus(1) * spin_minus(0)) +
+      Omega * cudaq::spin_operator::x(0) +
+      m_12 * Omega * cudaq::spin_operator::x(1);
+
+  // Each qubit is a 2-level system (dimension 2).
+  // The composite system (two qubits) has a total Hilbert space dimension of 2
+  // * 2 = 4.
+  cudaq::dimension_map dimensions{{0, 2}, {1, 2}};
+
+  // Build the initial state
+  // psi_00 represents the state |00> (both qubits in the ground state).
+  // psi_10 represents the state |10> (first qubit excited, second qubit in the
+  // ground state).
+  std::vector<std::complex<double>> psi00_data = {
+      {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
+  std::vector<std::complex<double>> psi10_data = {
+      {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
+
+  // Two initial state vectors for the 2-qubit system (dimension 4)
+  auto psi_00 = cudaq::state::from_data(psi00_data);
+  auto psi_10 = cudaq::state::from_data(psi10_data);
+
+  // Create a list of time steps for the simulation.
+  // Here we use 1001 points linearly spaced between time 0 and 1.
+  // This schedule will be used to integrate the time evolution of the system.
+  const int num_steps = 1001;
+  std::vector<double> steps = cudaq::linspace(0.0, 1.0, num_steps);
+  cudaq::schedule schedule(steps);
+
+  // Use Runge-`Kutta` integrator (4`th` order) to solve the time-dependent
+  // evolution. `dt` is the integration time step, and `order` sets the accuracy
+  // of the integrator method.
+  cudaq::integrators::runge_kutta integrator(4, 0.0001);
+
+  // The observables are the spin components along the x, y, and z directions
+  // for both qubits. These observables will be measured during the evolution.
+  auto observables = {cudaq::spin_operator::x(0), cudaq::spin_operator::y(0),
+                      cudaq::spin_operator::z(0), cudaq::spin_operator::x(1),
+                      cudaq::spin_operator::y(1), cudaq::spin_operator::z(1)};
+
+  // Evolution with 2 initial states
+  // We evolve the system under the defined Hamiltonian for both initial states
+  // simultaneously. No collapsed operators are provided (closed system
+  // evolution). The evolution returns expectation values for all defined
+  // observables at each time step.
+  auto evolution_results =
+      cudaq::evolve(hamiltonian, dimensions, schedule, {psi_00, psi_10},
+                    integrator, {}, observables, true);
+
+  // Retrieve the evolution result corresponding to each initial state.
+  auto &evolution_result_00 = evolution_results[0];
+  auto &evolution_result_10 = evolution_results[1];
+
+  // Lambda to extract expectation values for a given observable index
+  auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
+    std::vector<double> expectations;
+
+    auto all_exps = result.get_expectation_values().value();
+    for (auto exp_vals : all_exps) {
+      expectations.push_back((double)exp_vals[idx]);
+    }
+    return expectations;
+  };
+
+  // For the two `evolutions`, extract the six observable trajectories.
+  // For the |00> initial state, we extract the expectation trajectories for
+  // each observable.
+  auto result_00_0 = get_expectation(0, evolution_result_00);
+  auto result_00_1 = get_expectation(1, evolution_result_00);
+  auto result_00_2 = get_expectation(2, evolution_result_00);
+  auto result_00_3 = get_expectation(3, evolution_result_00);
+  auto result_00_4 = get_expectation(4, evolution_result_00);
+  auto result_00_5 = get_expectation(5, evolution_result_00);
+
+  // Similarly, for the |10> initial state:
+  auto result_10_0 = get_expectation(0, evolution_result_10);
+  auto result_10_1 = get_expectation(1, evolution_result_10);
+  auto result_10_2 = get_expectation(2, evolution_result_10);
+  auto result_10_3 = get_expectation(3, evolution_result_10);
+  auto result_10_4 = get_expectation(4, evolution_result_10);
+  auto result_10_5 = get_expectation(5, evolution_result_10);
+
+  // Export the results to a `CSV` file
+  // Export the Z-component of qubit 1's expectation values for both initial
+  // states. The `CSV` file "cross_resonance_z.`csv`" will have time versus (Z1)
+  // data for both |00> and |10> initial conditions.
+  export_csv("cross_resonance_z.csv", "time", steps, "<Z1>_00", result_00_5,
+             "<Z1>_10", result_10_5);
+  // Export the Y-component of qubit 1's expectation values for both initial
+  // states. The `CSV` file "cross_resonance_y.`csv`" will have time versus (Y1)
+  // data.
+  export_csv("cross_resonance_y.csv", "time", steps, "<Y1>_00", result_00_4,
+             "<Y1>_10", result_10_4);
+
+  std::cout
+      << "Simulation complete. The results are saved in cross_resonance_z.csv "
+         "and cross_resonance_y.csv files."
+      << std::endl;
+  return 0;
+}