* Adding Rydberg hamiltonian operator

* Adding unitests for Rydberg hamiltonian operator
* Making evaluate function in scalar_operator const

Signed-off-by: Sachin Pisal <spisal@nvidia.com>

runtime/cudaq/dynamics/CMakeLists.txt

@@ -11,7 +11,7 @@ set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-ctad-maybe-unsupported")
 set(INTERFACE_POSITION_INDEPENDENT_CODE ON)
 
 set(CUDAQ_OPS_SRC
-  scalar_operators.cpp elementary_operators.cpp product_operators.cpp operator_sum.cpp schedule.cpp definition.cpp helpers.cpp
+  scalar_operators.cpp elementary_operators.cpp product_operators.cpp operator_sum.cpp schedule.cpp definition.cpp helpers.cpp rydberg_hamiltonian.cpp
 )
 
 add_library(${LIBRARY_NAME} SHARED ${CUDAQ_OPS_SRC})

runtime/cudaq/dynamics/rydberg_hamiltonian.cpp

@@ -0,0 +1,55 @@
+/*******************************************************************************
+ * 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/operators.h"
+#include <sstream>
+#include <stdexcept>
+
+namespace cudaq {
+rydberg_hamiltonian::rydberg_hamiltonian(
+    const std::vector<Coordinate> &atom_sites, const scalar_operator &amplitude,
+    const scalar_operator &phase, const scalar_operator &delta_global,
+    const std::vector<int> &atom_filling,
+    const std::optional<std::pair<scalar_operator, std::vector<double>>>
+        &delta_local)
+    : atom_sites(atom_sites), amplitude(amplitude), phase(phase),
+      delta_global(delta_global), delta_local(delta_local) {
+  if (atom_filling.empty()) {
+    this->atom_filling = std::vector<int>(atom_sites.size(), 1);
+  } else if (atom_sites.size() != atom_filling.size()) {
+    throw std::invalid_argument(
+        "Size of `atom_sites` and `atom_filling` must be equal.");
+  } else {
+    this->atom_filling = atom_filling;
+  }
+
+  if (delta_local.has_value()) {
+    throw std::runtime_error(
+        "Local detuning is an experimental feature not yet supported.");
+  }
+}
+
+const std::vector<rydberg_hamiltonian::Coordinate> &
+rydberg_hamiltonian::get_atom_sites() const {
+  return atom_sites;
+}
+
+const std::vector<int> &rydberg_hamiltonian::get_atom_filling() const {
+  return atom_filling;
+}
+
+const scalar_operator &rydberg_hamiltonian::get_amplitude() const {
+  return amplitude;
+}
+
+const scalar_operator &rydberg_hamiltonian::get_phase() const { return phase; }
+
+const scalar_operator &rydberg_hamiltonian::get_delta_global() const {
+  return delta_global;
+}
+} // namespace cudaq

runtime/cudaq/dynamics/scalar_operators.cpp

@@ -40,7 +40,7 @@ scalar_operator::scalar_operator(double value) {
 }
 
 std::complex<double> scalar_operator::evaluate(
-    std::map<std::string, std::complex<double>> parameters) {
+    std::map<std::string, std::complex<double>> parameters) const {
   return generator(parameters);
 }
 

runtime/cudaq/operators.h

@@ -64,8 +64,9 @@ class operator_sum {
   ///                      degrees of freedom: `{0:2, 1:2}`.
   /// @arg `parameters` : A map of the parameter names to their concrete,
   /// complex values.
-  matrix_2 to_matrix(const std::map<int, int> &dimensions,
-                     const std::map<std::string, double> &params = {}) const;
+  matrix_2 to_matrix(
+      const std::map<int, int> &dimensions,
+      const std::map<std::string, std::complex<double>> &params = {}) const;
 
   // Arithmetic operators
   operator_sum operator+(const operator_sum &other) const;
@@ -398,7 +399,7 @@ class scalar_operator : public product_operator {
 
   /// @brief Return the scalar operator as a concrete complex value.
   std::complex<double>
-  evaluate(std::map<std::string, std::complex<double>> parameters);
+  evaluate(std::map<std::string, std::complex<double>> parameters) const;
 
   // Return the scalar operator as a 1x1 matrix. This is needed for
   // compatability with the other inherited classes.
@@ -461,4 +462,50 @@ void operator-=(scalar_operator &self, scalar_operator other);
 void operator*=(scalar_operator &self, scalar_operator other);
 void operator/=(scalar_operator &self, scalar_operator other);
 
+/// @brief Representation of a time-dependent Hamiltonian for Rydberg system
+class rydberg_hamiltonian : public operator_sum {
+public:
+  using Coordinate = std::pair<double, double>;
+
+  /// @brief Constructor.
+  /// @param atom_sites List of 2D coordinates for trap sites.
+  /// @param amplitude Time-dependant driving amplitude, Omega(t).
+  /// @param phase Time-dependant driving phase, phi(t).
+  /// @param delta_global Time-dependant driving detuning, Delta_global(t).
+  /// @param atom_filling Optional. Marks occupied trap sites (1) and empty
+  /// sites (0). Defaults to all sites occupied.
+  /// @param delta_local Optional. A tuple of Delta_local(t) and site dependant
+  /// local detuning factors.
+  rydberg_hamiltonian(
+      const std::vector<Coordinate> &atom_sites,
+      const scalar_operator &amplitude, const scalar_operator &phase,
+      const scalar_operator &delta_global,
+      const std::vector<int> &atom_filling = {},
+      const std::optional<std::pair<scalar_operator, std::vector<double>>>
+          &delta_local = std::nullopt);
+
+  /// @brief Get atom sites.
+  const std::vector<Coordinate> &get_atom_sites() const;
+
+  /// @brief Get atom filling.
+  const std::vector<int> &get_atom_filling() const;
+
+  /// @brief Get amplitude operator.
+  const scalar_operator &get_amplitude() const;
+
+  /// @brief Get phase operator.
+  const scalar_operator &get_phase() const;
+
+  /// @brief Get global detuning operator.
+  const scalar_operator &get_delta_global() const;
+
+private:
+  std::vector<Coordinate> atom_sites;
+  std::vector<int> atom_filling;
+  scalar_operator amplitude;
+  scalar_operator phase;
+  scalar_operator delta_global;
+  std::optional<std::pair<scalar_operator, std::vector<double>>> delta_local;
+};
+
 } // namespace cudaq

unittests/CMakeLists.txt

@@ -52,6 +52,7 @@ set(CUDAQ_RUNTIME_TEST_SOURCES
   dynamics/test_runge_kutta_time_stepper.cpp
   dynamics/test_runge_kutta_integrator.cpp
   dynamics/test_helpers.cpp
+  dynamics/rydberg_hamiltonian.cpp
 )
 
 # Make it so we can get function symbols

unittests/dynamics/rydberg_hamiltonian.cpp

@@ -0,0 +1,124 @@
+/*******************************************************************************
+ * 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/operators.h"
+#include <gtest/gtest.h>
+
+using namespace cudaq;
+
+TEST(RydbergHamiltonianTest, ConstructorValidInputs) {
+  // Valid atom sites
+  std::vector<rydberg_hamiltonian::Coordinate> atom_sites = {
+      {0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}};
+
+  // Valid operators
+  scalar_operator amplitude(1.0);
+  scalar_operator phase(0.0);
+  scalar_operator delta_global(-0.5);
+
+  // Valid atom filling
+  rydberg_hamiltonian hamiltonian(atom_sites, amplitude, phase, delta_global);
+
+  EXPECT_EQ(hamiltonian.get_atom_sites().size(), atom_sites.size());
+  EXPECT_EQ(hamiltonian.get_atom_filling().size(), atom_sites.size());
+  EXPECT_EQ(hamiltonian.get_amplitude().evaluate({}),
+            std::complex<double>(1.0, 0.0));
+  EXPECT_EQ(hamiltonian.get_phase().evaluate({}),
+            std::complex<double>(0.0, 0.0));
+  EXPECT_EQ(hamiltonian.get_delta_global().evaluate({}),
+            std::complex<double>(-0.5, 0.0));
+}
+
+TEST(RydbergHamiltonianTest, ConstructorWithAtomFilling) {
+  std::vector<rydberg_hamiltonian::Coordinate> atom_sites = {
+      {0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}};
+
+  // Valid operators
+  scalar_operator amplitude(1.0);
+  scalar_operator phase(0.0);
+  scalar_operator delta_global(-0.5);
+
+  // Valid atom filling
+  std::vector<int> atom_filling = {1, 0, 1};
+
+  rydberg_hamiltonian hamiltonian(atom_sites, amplitude, phase, delta_global,
+                                  atom_filling);
+
+  EXPECT_EQ(hamiltonian.get_atom_sites().size(), atom_sites.size());
+  EXPECT_EQ(hamiltonian.get_atom_filling(), atom_filling);
+}
+
+TEST(RydbergHamiltonianTest, InvalidAtomFillingSize) {
+  std::vector<rydberg_hamiltonian::Coordinate> atom_sites = {
+      {0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}};
+
+  // Valid operators
+  scalar_operator amplitude(1.0);
+  scalar_operator phase(0.0);
+  scalar_operator delta_global(-0.5);
+
+  // Invalid atom filling size
+  std::vector<int> atom_filling = {1, 0};
+
+  EXPECT_THROW(rydberg_hamiltonian(atom_sites, amplitude, phase, delta_global,
+                                   atom_filling),
+               std::invalid_argument);
+}
+
+TEST(RydbergHamiltonianTest, UnsupportedLocalDetuning) {
+  std::vector<rydberg_hamiltonian::Coordinate> atom_sites = {
+      {0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}};
+
+  // Valid operators
+  scalar_operator amplitude(1.0);
+  scalar_operator phase(0.0);
+  scalar_operator delta_global(-0.5);
+
+  // Invalid delta_local
+  auto delta_local =
+      std::make_pair(scalar_operator(0.5), std::vector<double>{0.1, 0.2, 0.3});
+
+  EXPECT_THROW(rydberg_hamiltonian(atom_sites, amplitude, phase, delta_global,
+                                   {}, delta_local),
+               std::runtime_error);
+}
+
+TEST(RydbergHamiltonianTest, Accessors) {
+  std::vector<rydberg_hamiltonian::Coordinate> atom_sites = {
+      {0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}};
+
+  // Valid operators
+  scalar_operator amplitude(1.0);
+  scalar_operator phase(0.0);
+  scalar_operator delta_global(-0.5);
+
+  rydberg_hamiltonian hamiltonian(atom_sites, amplitude, phase, delta_global);
+
+  EXPECT_EQ(hamiltonian.get_atom_sites(), atom_sites);
+  EXPECT_EQ(hamiltonian.get_amplitude().evaluate({}),
+            std::complex<double>(1.0, 0.0));
+  EXPECT_EQ(hamiltonian.get_phase().evaluate({}),
+            std::complex<double>(0.0, 0.0));
+  EXPECT_EQ(hamiltonian.get_delta_global().evaluate({}),
+            std::complex<double>(-0.5, 0.0));
+}
+
+TEST(RydbergHamiltonianTest, DefaultAtomFilling) {
+  std::vector<rydberg_hamiltonian::Coordinate> atom_sites = {
+      {0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}};
+
+  // Valid operators
+  scalar_operator amplitude(1.0);
+  scalar_operator phase(0.0);
+  scalar_operator delta_global(-0.5);
+
+  rydberg_hamiltonian hamiltonian(atom_sites, amplitude, phase, delta_global);
+
+  std::vector<int> expected_filling(atom_sites.size(), 1);
+  EXPECT_EQ(hamiltonian.get_atom_filling(), expected_filling);
+}