First pass of support for C++ operators

.devcontainer/devcontainer.json

@@ -241,9 +241,9 @@
 				"C_Cpp.inlayHints.referenceOperator.enabled": false,
 				"C_Cpp.doxygen.generateOnType": false,
 				"C_Cpp.default.cStandard": "c17",
-				"C_Cpp.default.cppStandard": "c++20",
+				"C_Cpp.default.cppStandard": "c++17",
 				"C_Cpp.default.configurationProvider": "ms-vscode.cmake-tools",
-				"C_Cpp.default.compileCommands": "${workspaceFolder}/build/compile_commands.json",
+				"C_Cpp.default.compileCommands": "${workspaceFolder}/build/debug/compile_commands.json",
 				"C_Cpp.default.includePath": [
 					"${workspaceFolder}/**",
 					"${CUDAQ_INSTALL_PREFIX}/**",

.github/workflows/publishing.yml

@@ -839,7 +839,7 @@ jobs:
                   fi; \
               done`
 
-          rm -rf examples applications targets && mv github-repo/docs/sphinx/examples examples && mv github-repo/docs/sphinx/applications applications && mv github-repo/docs/sphinx/targets targets
+          rm -rf examples applications targets snippets && mv github-repo/docs/sphinx/examples examples && mv github-repo/docs/sphinx/applications applications && mv github-repo/docs/sphinx/targets targets && mv github-repo/docs/sphinx/snippets/python snippets
           mv github-repo/docs/notebook_validation.py .
           GITHUB_STEP_SUMMARY=$GITHUB_STEP_SUMMARY \
           bash github-repo/scripts/validate_container.sh $backends_to_test | tee /tmp/validation.out

CMakeLists.txt

@@ -151,6 +151,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()

cmake/Modules/CMakeLists.txt

@@ -11,6 +11,7 @@ set(CONFIG_FILES
   CUDAQEmDefaultConfig.cmake
   CUDAQNloptConfig.cmake
   CUDAQSpinConfig.cmake
+  CUDAQOperatorConfig.cmake
   CUDAQConfig.cmake
   CUDAQEnsmallenConfig.cmake
   CUDAQPlatformDefaultConfig.cmake

cmake/Modules/CUDAQConfig.cmake

@@ -14,6 +14,9 @@ list(APPEND CMAKE_MODULE_PATH "${CUDAQ_CMAKE_DIR}")
 set (CUDAQSpin_DIR "${CUDAQ_CMAKE_DIR}")
 find_dependency(CUDAQSpin REQUIRED)
 
+set (CUDAQOperator_DIR "${CUDAQ_CMAKE_DIR}")
+find_dependency(CUDAQOperator REQUIRED)
+
 set (CUDAQCommon_DIR "${CUDAQ_CMAKE_DIR}")
 find_dependency(CUDAQCommon REQUIRED)
 

cmake/Modules/CUDAQEmDefaultConfig.cmake

@@ -11,6 +11,9 @@ get_filename_component(CUDAQ_EM_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" PATH)
 set (CUDAQSpin_DIR "${CUDAQ_EM_CMAKE_DIR}")
 find_dependency(CUDAQSpin REQUIRED)
 
+set (CUDAQOperator_DIR "${CUDAQ_EM_CMAKE_DIR}")
+find_dependency(CUDAQOperator REQUIRED)
+
 set (CUDAQCommon_DIR "${CUDAQ_EM_CMAKE_DIR}")
 find_dependency(CUDAQCommon REQUIRED)
 

cmake/Modules/CUDAQOperatorConfig.cmake

@@ -0,0 +1,13 @@
+# ============================================================================ #
+# 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.     #
+# ============================================================================ #
+
+get_filename_component(CUDAQ_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" PATH)
+
+if(NOT TARGET cudaq::cudaq-operator)
+  include("${CUDAQ_CMAKE_DIR}/CUDAQOperatorTargets.cmake")
+endif()

cmake/Modules/CUDAQPlatformDefaultConfig.cmake

@@ -11,6 +11,9 @@ get_filename_component(CUDAQ_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" PATH)
 set (CUDAQEmDefault_DIR "${CUDAQ_CMAKE_DIR}")
 find_dependency(CUDAQEmDefault REQUIRED)
 
+set (CUDAQOperator_DIR "${CUDAQ_CMAKE_DIR}")
+find_dependency(CUDAQOperator REQUIRED)
+
 set (CUDAQSpin_DIR "${CUDAQ_CMAKE_DIR}")
 find_dependency(CUDAQSpin REQUIRED)
 

cmake/Modules/NVQIRConfig.cmake.in

@@ -12,6 +12,7 @@ include(CMakeFindDependencyMacro)
 get_filename_component(PARENT_DIRECTORY ${NVQIR_CMAKE_DIR} DIRECTORY)
 
 find_dependency(CUDAQSpin REQUIRED HINTS "${PARENT_DIRECTORY}/cudaq")
+find_dependency(CUDAQOperator REQUIRED HINTS "${PARENT_DIRECTORY}/cudaq")
 find_dependency(CUDAQCommon REQUIRED HINTS "${PARENT_DIRECTORY}/cudaq")
 find_package(fmt QUIET)
 if (NOT fmt_FOUND)

docker/build/devdeps.ext.Dockerfile

@@ -181,7 +181,3 @@ ENV CUTENSOR_INSTALL_PREFIX="$CUTENSOR_INSTALL_PREFIX"
 ENV CUTENSOR_ROOT="$CUTENSOR_INSTALL_PREFIX"
 ENV LD_LIBRARY_PATH="$CUTENSOR_INSTALL_PREFIX/lib:$LD_LIBRARY_PATH"
 ENV CPATH="$CUTENSOR_INSTALL_PREFIX/include:$CPATH"
-
-# Active MPI support for the cuTensorNet library
-RUN cd "$CUQUANTUM_INSTALL_PREFIX/distributed_interfaces/" && source activate_mpi_cutn.sh
-ENV CUTENSORNET_COMM_LIB="$CUQUANTUM_INSTALL_PREFIX/distributed_interfaces/libcutensornet_distributed_interface_mpi.so"

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

@@ -0,0 +1,142 @@
+/*******************************************************************************
+ * 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.    *
+ ******************************************************************************/
+
+#include "cudaq/algorithms/evolve.h"
+#include "cudaq/dynamics_integrators.h"
+#include "cudaq/evolution.h"
+#include "cudaq/operators.h"
+#include "cudaq/schedule.h"
+#include "matplotlibcpp.h"
+#include <cudaq.h>
+
+namespace plt = matplotlibcpp;
+
+int main() {
+
+  cudaq::set_target_backend("dynamics");
+
+  // System dimensions:
+  // subsystem 0 (atom) has 2 levels
+  // subsystem 1 (cavity) has 10 levels
+  std::map<int, int> dimensions{{0, 2}, {1, 10}};
+
+  // For the cavity subsystem 1
+  auto a = cudaq::matrix_operator::annihilate(1);
+  auto a_dag = cudaq::matrix_operator::create(1);
+
+  // For the atom subsystem 0
+  auto sm = cudaq::matrix_operator::annihilate(0);
+  auto sm_dag = cudaq::matrix_operator::create(0);
+
+  cudaq::product_operator<cudaq::matrix_operator> atom_occ_op_t =
+      cudaq::matrix_operator::number(0);
+  cudaq::operator_sum<cudaq::matrix_operator> atom_occ_op(atom_occ_op_t);
+
+  cudaq::product_operator<cudaq::matrix_operator> cavity_occ_op_t =
+      cudaq::matrix_operator::number(1);
+  cudaq::operator_sum<cudaq::matrix_operator> cavity_occ_op(cavity_occ_op_t);
+
+  auto hamiltonian = 2 * M_PI * atom_occ_op + 2 * M_PI * cavity_occ_op +
+                     2 * M_PI * 0.25 * (sm * a_dag + sm_dag * a);
+
+  // Build the initial state
+  // For the atom, the density matrix in the ground state
+  cudaq::matrix_2 qubit_state({1.0, 0.0}, {0.0, 0.0});
+
+  // For the cavity, 1 10x10 matrix with a single photon number state at |5>
+  cudaq::matrix_2 cavity_state = cudaq::matrix_2(10, 10);
+  cavity_state[{5, 5}] = 1.0;
+
+  // Compute the tensor (kronecker) product of the atom and cavity states
+  cudaq::matrix_2 rho = qubit_state.kronecker_inplace(cavity_state);
+
+  // Flatten the matrix
+  std::vector<std::complex<double>> flat_rho;
+  for (size_t j = 0; j < rho.get_columns(); ++j) {
+    for (size_t i = 0; i < rho.get_rows(); ++i) {
+      flat_rho.push_back(rho[{i, j}]);
+    }
+  }
+
+  cudaq::state_data rho_data = flat_rho;
+
+  // Create a CUDA quantum state from a density matrix
+  auto rho0 = cudaq::state::from_data(rho_data);
+
+  // Create time steps between 0 and 10
+  const int num_steps = 201;
+  const double t0 = 0.0, t1 = 10.0;
+  std::vector<double> steps(num_steps);
+  double dt = (t1 - t0) / (num_steps - 1);
+  for (int i = 0; i < num_steps; ++i) {
+    steps[i] = t0 + i * dt;
+  }
+
+  // Create a schedule for the time steps and a label for the time parameter
+  std::vector<std::string> labels = {"time"};
+  cudaq::Schedule schedule(steps, labels);
+
+  cudaq::runge_kutta integrator;
+  integrator.dt = 0.001;
+  integrator.order = 1;
+
+  // Evolve without collapse operators
+  cudaq::evolve_result evolve_result =
+      cudaq::evolve(hamiltonian, dimensions, schedule, rho0, integrator, {},
+                    {hamiltonian}, true);
+
+  constexpr double decay_rate = 0.1;
+  auto collapse_operator = std::sqrt(decay_rate) * a;
+  // Evolve with collapse operators
+  cudaq::evolve_result evolve_result_decay =
+      cudaq::evolve(hamiltonian, dimensions, schedule, rho0, integrator,
+                    {collapse_operator}, {hamiltonian}, true);
+
+  // Lambda to extract expectation values for a given observable index
+  auto get_expectation = [](int idx,
+                            const auto &result) -> std::vector<double> {
+    std::vector<double> expectations;
+
+    auto all_exps = result.get_expectation_values().value();
+    for (const auto &exp_vals : all_exps) {
+      expectations.push_back((double)exp_vals[idx]);
+    }
+    return expectations;
+  };
+
+  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);
+
+  // Plot the results
+  plt::figure_size(1000, 600);
+
+  // Subplot 1: No decay
+  plt::subplot(1, 2, 1);
+  plt::plot(steps, ideal_result0, {{"label", "Cavity Photon Number"}});
+  plt::plot(steps, ideal_result1, {{"label", "Atom Excitation Probability"}});
+  plt::xlabel("Time");
+  plt::ylabel("Expectation value");
+  plt::legend();
+  plt::title("No decay");
+
+  // Subplot 2: With decay
+  plt::subplot(1, 2, 2);
+  plt::plot(steps, decay_result0, {{"label", "Cavity Photon Number"}});
+  plt::plot(steps, decay_result1, {{"label", "Atom Excitation Probability"}});
+  plt::xlabel("Time");
+  plt::ylabel("Expectation value");
+  plt::legend();
+  plt::title("No decay");
+
+  plt::save("cavity_qed.png");
+
+  std::cout << "Simulation complete. Plot saved to cavity_qed.png" << std::endl;
+  return 0;
+}

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

@@ -0,0 +1,138 @@
+/*******************************************************************************
+ * 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.    *
+ ******************************************************************************/
+
+#include "cudaq/algorithms/evolve.h"
+#include "cudaq/dynamics_integrators.h"
+#include "cudaq/evolution.h"
+#include "cudaq/operators.h"
+#include "cudaq/schedule.h"
+#include "matplotlibcpp.h"
+#include <cudaq.h>
+
+namespace plt = matplotlibcpp;
+
+int main() {
+
+  cudaq::set_target_backend("dynamics");
+
+  // Detuning between two qubits
+  double delta = 100 * 2 * M_PI;
+  // Static coupling between qubits
+  double J = 7 * 2 * M_PI;
+  // Spurious electromagnetic crosstalk
+  double m_12 = 0.2;
+  // Drive strength
+  double Omega = 20 * 2 * M_PI;
+
+  auto hamiltonian =
+      (delta / 2.0) * cudaq::spin_operator::z(0) +
+      J * (cudaq::spin_operator::minus(1) * cudaq::spin_operator::plus(0) +
+           cudaq::spin_operator::plus(1) * cudaq::spin_operator::minus(0)) +
+      Omega * cudaq::spin_operator::x(0) +
+      m_12 * Omega * cudaq::spin_operator::x(1);
+
+  std::map<int, int> dimensions{{0, 2}, {1, 2}};
+
+  // Build the initial state
+  cudaq::matrix_2 rho_mat({1.0, 0.0}, {0.0, 0.0});
+
+  // Flatten the matrix
+  std::vector<std::complex<double>> flat_rho;
+  for (size_t j = 0; j < rho_mat.get_columns(); ++j) {
+    for (size_t i = 0; i < rho_mat.get_rows(); ++i) {
+      flat_rho.push_back(rho_mat[{i, j}]);
+    }
+  }
+
+  cudaq::state_data rho_data = flat_rho;
+  auto rho0 = cudaq::state::from_data(rho_data);
+
+  // Two initial state vectors for the 2-qubit system (dimension 4)
+  // psi_00 corresponds to |00> and psi_10 corresponds to |10>
+  auto psi_00 = cudaq::state::from_data({1.0, 0.0, 0.0, 0.0});
+  auto psi_10 = cudaq::state::from_data({0.0, 0.0, 1.0, 0.0});
+
+  // Create a schedule of time steps
+  const int num_steps = 1001;
+  const double t0 = 0.0, t1 = 1.0;
+  std::vector<double> steps(num_steps);
+  double dt = (t1 - t0) / (num_steps - 1);
+  for (int i = 0; i < num_steps; ++i) {
+    steps[i] = t0 + i * dt;
+  }
+  std::vector<std::string> labels = {"time"};
+  cudaq::Schedule schedule(steps, labels);
+
+  cudaq::runge_kutta integrator;
+  integrator.dt = 0.001;
+  integrator.order = 1;
+
+  std::vector<cudaq::spin_operator> 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 for initial state |00>
+  auto evolution_result_00 =
+      cudaq::evolve(hamiltonian, dimensions, schedule, psi_00, observables, {},
+                    true, integrator);
+
+  // Evolution for initial state |10>
+  auto evolution_result_10 =
+      cudaq::evolve(hamiltonian, dimensions, schedule, psi_10, observables, {},
+                    true, integrator);
+
+  auto get_result = [](int idx, const auto &result) -> std::vector<double> {
+    std::vector<double> expectations;
+    auto all_exps = result.get_expectation_values().value();
+    for (const auto &exp_vals : all_exps) {
+      expectations.push_back((double)exp_vals[idx]);
+    }
+    return expectations;
+  };
+
+  // For the two evolutions, extract the six observable trajectories.
+  auto result_00_0 = get_result(0, evolution_result_00);
+  auto result_00_1 = get_result(1, evolution_result_00);
+  auto result_00_2 = get_result(2, evolution_result_00);
+  auto result_00_3 = get_result(3, evolution_result_00);
+  auto result_00_4 = get_result(4, evolution_result_00);
+  auto result_00_5 = get_result(5, evolution_result_00);
+
+  auto result_10_0 = get_result(0, evolution_result_10);
+  auto result_10_1 = get_result(1, evolution_result_10);
+  auto result_10_2 = get_result(2, evolution_result_10);
+  auto result_10_3 = get_result(3, evolution_result_10);
+  auto result_10_4 = get_result(4, evolution_result_10);
+  auto result_10_5 = get_result(5, evolution_result_10);
+
+  // Plot the results
+  plt::figure_size(1000, 600);
+
+  // Subplot 1: Expectation value <Z> for qubit 1.
+  plt::subplot(1, 2, 1);
+  plt::plot(steps, result_00_5, {{"label", "$|\\psi_0\\rangle=|00\\rangle$"}});
+  plt::plot(steps, result_10_5, {{"label", "$|\\psi_0\\rangle=|10\\rangle$"}});
+  plt::xlabel("Time");
+  plt::ylabel("$\\langle Z_2 \\rangle$");
+  plt::legend();
+
+  // Subplot 2: Expectation value <Y> for qubit 1.
+  plt::subplot(1, 2, 2);
+  plt::plot(steps, result_00_4, {{"label", "$|\\psi_0\\rangle=|00\\rangle$"}});
+  plt::plot(steps, result_10_4, {{"label", "$|\\psi_0\\rangle=|10\\rangle$"}});
+  plt::xlabel("Time");
+  plt::ylabel("$\\langle Y_2 \\rangle$");
+  plt::legend();
+
+  plt::save("cross_resonance.png");
+
+  std::cout << "Simulation complete. Plot saved to cross_resonance.png"
+            << std::endl;
+  return 0;
+}

docs/sphinx/using/backends/sims/tnsims.rst

@@ -70,11 +70,12 @@ Use the following commands to enable distribution across multiple GPUs (adjust t
         mpiexec -np 2 ./program.x
 
 .. note::
+  MPI parallelization on the :code:`tensornet` backend requires CUDA-Q's MPI support. 
+  Please refer to the instructions on how to :ref:`enable MPI parallelization <distributed-computing-with-mpi>` within CUDA-Q.  
+  CUDA-Q containers are shipped with a pre-built MPI plugin; hence no additional setup is needed.  
 
-  If the `CUTENSORNET_COMM_LIB` environment variable is not set, MPI parallelization on the :code:`tensornet` backend may fail.
-  If you are using a CUDA-Q container, this variable is pre-configured and no additional setup is needed. If you are customizing your installation or have built CUDA-Q from source, please follow the instructions for `activating the distributed interface <https://docs.nvidia.com/cuda/cuquantum/latest/getting-started/index.html#from-nvidia-devzone>`__ for the `cuTensorNet` library. This requires 
-  :ref:`installing CUDA development dependencies <additional-cuda-tools>`, and setting the `CUTENSORNET_COMM_LIB`
-  environment variable to the newly built `libcutensornet_distributed_interface_mpi.so` library.
+.. note::  
+  If the `CUTENSORNET_COMM_LIB` environment variable is set following the activation procedure described in the `cuTensorNet documentation <https://docs.nvidia.com/cuda/cuquantum/latest/getting-started/index.html#from-nvidia-devzone>`__, the cuTensorNet MPI plugin will take precedence over the builtin support from CUDA-Q.
 
 Specific aspects of the simulation can be configured by setting the following of environment variables: