NVIDIA
diff --git a/‎docs/sphinx/api/languages/python_api.rst‎
Lines changed: 1 addition & 1 deletion b/‎docs/sphinx/api/languages/python_api.rst‎
Lines changed: 1 addition & 1 deletion
diff --git a/‎docs/sphinx/examples/cpp/dynamics/heisenberg_model_mpi.cpp‎
Lines changed: 157 additions & 0 deletions b/‎docs/sphinx/examples/cpp/dynamics/heisenberg_model_mpi.cpp‎
Lines changed: 157 additions & 0 deletions
diff --git a/‎docs/sphinx/examples/python/dynamics/cavity_qed.py‎
Lines changed: 1 addition & 1 deletion b/‎docs/sphinx/examples/python/dynamics/cavity_qed.py‎
Lines changed: 1 addition & 1 deletion
diff --git a/‎docs/sphinx/examples/python/dynamics/cross_resonance.py‎
Lines changed: 1 addition & 5 deletions b/‎docs/sphinx/examples/python/dynamics/cross_resonance.py‎
Lines changed: 1 addition & 5 deletions
diff --git a/‎python/CMakeLists.txt‎
Lines changed: 4 additions & 0 deletions b/‎python/CMakeLists.txt‎
Lines changed: 4 additions & 0 deletions
diff --git a/‎python/cudaq/__init__.py‎
Lines changed: 2 additions & 0 deletions b/‎python/cudaq/__init__.py‎
Lines changed: 2 additions & 0 deletions
@@ -150,7 +150,7 @@ Dynamics
 .. autoclass:: cudaq.dynamics.integrator.BaseIntegrator
 
 .. autoclass:: cudaq.dynamics.helpers.InitialState
-.. autofunction:: cudaq.dynamics.cudm_state.to_cupy_array
+.. autoclass:: cudaq.InitialStateType
 
 Operators
 =============================
 
@@ -0,0 +1,157 @@
+/*******************************************************************************
+ * 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 heisenberg_model_mpi.cpp -o a.out &&
+// mpiexec -np <N> ./a.out
+// ```
+
+#include "cudaq/algorithms/evolve.h"
+#include "cudaq/algorithms/integrator.h"
+#include "cudaq/operators.h"
+#include "export_csv_helper.h"
+#include <cudaq.h>
+
+int main() {
+  cudaq::mpi::initialize();
+  std::cout << "Number of ranks = " << cudaq::mpi::num_ranks() << "\n";
+  // Set up a 15-spin chain, where each spin is a two-level system.
+  const int num_spins = 15;
+  cudaq::dimension_map dimensions;
+  for (int i = 0; i < num_spins; i++) {
+    dimensions[i] = 2; // Each spin (site) has dimension 2.
+  }
+
+  // Initial state
+  // Prepare an initial state where the spins are arranged in a staggered
+  // configuration. Even indices get the value '0' and odd indices get '1'. For
+  // example, for 15 spins: spins: 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0
+  std::string spin_state;
+  for (int i = 0; i < num_spins; i++) {
+    spin_state.push_back((i % 2 == 0) ? '0' : '1');
+  }
+
+  // Convert the binary string to an integer index
+  // In the Hilbert space of 15 spins (size 2^15 = 32768), this index
+  // corresponds to the state |0 1 0 1 0 1 0 1 0 1 0 1 0 1 0>
+  int initial_state_index = std::stoi(spin_state, nullptr, 2);
+
+  // Build the staggered magnetization operator
+  // The staggered magnetization operator is used to measure antiferromagnetic
+  // order. It is defined as a sum over all spins of the Z operator, alternating
+  // in sign. For even sites, we add `sz`; for odd sites, we subtract `sz`.
+  auto staggered_magnetization_t = cudaq::spin_op::empty();
+  for (int i = 0; i < num_spins; i++) {
+    auto sz = cudaq::spin_op::z(i);
+    if (i % 2 == 0) {
+      staggered_magnetization_t += sz;
+    } else {
+      staggered_magnetization_t -= sz;
+    }
+  }
+
+  // Normalize the number of spins so that the observable is intensive.
+  auto staggered_magnetization_op =
+      (1 / static_cast<double>(num_spins)) * staggered_magnetization_t;
+
+  // Each entry will associate a value of g (the `anisotropy` in the Z coupling)
+  // with its corresponding time-series of expectation values of the staggered
+  // magnetization.
+  std::vector<std::pair<double, std::vector<double>>> observe_results;
+
+  // Simulate the dynamics over 100 time steps spanning from time 0 to 5.
+  const int num_steps = 100;
+  std::vector<double> steps = cudaq::linspace(0.0, 5.0, num_steps);
+
+  // For three different values of g, which sets the strength of the Z-Z
+  // interaction: g = 0.0 (isotropic in the `XY` plane), 0.25, and 4.0 (strongly
+  // `anisotropy`).
+  std::vector<double> g_values = {0.0, 0.25, 4.0};
+
+  for (auto g : g_values) {
+    // Set the coupling strengths:
+    // `Jx` and `Jy` are set to 1.0 (coupling along X and Y axes), while `Jz` is
+    // set to the current g value (coupling along the Z axis).
+    double Jx = 1.0, Jy = 1.0, Jz = g;
+
+    // The Hamiltonian is built from the nearest-neighbor interactions:
+    // H = H + `Jx` * `Sx`_i * `Sx`_{i+1}
+    // H = H + `Jy` * `Sy`_i * `Sy`_{i+1}
+    // H = H + `Jz` * `Sz`_i * `Sz`_{i+1}
+    // This is a form of the `anisotropic` Heisenberg (or `XYZ`) model.
+    auto hamiltonian = cudaq::spin_op::empty();
+    for (int i = 0; i < num_spins - 1; i++) {
+      hamiltonian =
+          hamiltonian + Jx * cudaq::spin_op::x(i) * cudaq::spin_op::x(i + 1);
+      hamiltonian =
+          hamiltonian + Jy * cudaq::spin_op::y(i) * cudaq::spin_op::y(i + 1);
+      hamiltonian =
+          hamiltonian + Jz * cudaq::spin_op::z(i) * cudaq::spin_op::z(i + 1);
+    }
+
+    // Initial state vector
+    // For a 9-spin system, the Hilbert space dimension is 2^9 = 512.
+    // Initialize the state as a vector with all zeros except for a 1 at the
+    // index corresponding to our staggered state.
+    const int state_size = 1 << num_spins;
+    std::vector<std::complex<double>> psi0_data(state_size, {0.0, 0.0});
+    psi0_data[initial_state_index] = {1.0, 0.0};
+    auto psi0 = cudaq::state::from_data(psi0_data);
+
+    // The schedule is built using the time steps array.
+    cudaq::schedule schedule(steps);
+
+    // Use a Runge-`Kutta` integrator (4`th` order) with a small time step `dt`
+    // = 0.001.
+    cudaq::integrators::runge_kutta integrator(4, 0.001);
+
+    // Evolve the initial state psi0 under the Hamiltonian, using the specified
+    // schedule and integrator. No collapse operators are included (closed
+    // system evolution). Measure the expectation value of the staggered
+    // magnetization operator at each time step.
+    auto evolve_result =
+        cudaq::evolve(hamiltonian, dimensions, schedule, psi0, integrator, {},
+                      {staggered_magnetization_op}, true);
+
+    // 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.expectation_values.value();
+      for (auto exp_vals : all_exps) {
+        expectations.push_back((double)exp_vals[idx]);
+      }
+      return expectations;
+    };
+
+    observe_results.push_back({g, get_expectation(0, evolve_result)});
+  }
+
+  if (observe_results.size() != 3) {
+    std::cerr << "Unexpected number of g values" << std::endl;
+    return 1;
+  }
+
+  if (cudaq::mpi::rank() == 0) {
+    // Only save on the first rank to prevent race condition.
+    // The `CSV` file "`heisenberg`_model.`csv`" will contain column with:
+    //    - The time steps
+    //    - The expectation values of the staggered magnetization for each g
+    //    value (labeled g_0, g_0.25, g_4).
+    export_csv("heisenberg_model_mpi_result.csv", "time", steps, "g_0",
+               observe_results[0].second, "g_0.25", observe_results[1].second,
+               "g_4", observe_results[2].second);
+
+    std::cout << "Simulation complete. The results are saved in "
+                 "heisenberg_model_mpi_result.csv file."
+              << std::endl;
+  }
+  cudaq::mpi::finalize();
+  return 0;
+}
@@ -35,7 +35,7 @@
 # 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))
+rho0 = cudaq.State.from_data(cp.kron(cavity_state, qubit_state))
 
 steps = np.linspace(0, 10, 201)
 schedule = Schedule(steps, ["time"])
 
@@ -30,16 +30,12 @@
 # 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))
+    cp.array([0.0, 1.0, 0.0, 0.0], dtype=cp.complex128))
 
 # Schedule of time steps.
 steps = np.linspace(0.0, 1.0, 1001)
 
@@ -78,3 +78,7 @@ if(CUDAQ_BUILD_TESTS)
 endif()
 
 add_subdirectory(runtime/interop)
+
+if (CUDENSITYMAT_ROOT AND CUDA_FOUND)
+  add_subdirectory(runtime/cudaq/dynamics)
+endif()
@@ -97,6 +97,8 @@
 from .dynamics.evolution import evolve, evolve_async
 from .dynamics.integrators import *
 
+InitialStateType = cudaq_runtime.InitialStateType
+
 # Optimizers + Gradients
 optimizers = cudaq_runtime.optimizers
 gradients = cudaq_runtime.gradients