Add Suzuki-Trotter circuit primitive

libs/solvers/include/cudaq/solvers/hamiltonian_simulation/trotter.h

@@ -0,0 +1,64 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2026 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.    *
+ ******************************************************************************/
+#pragma once
+
+#include "cudaq.h"
+#include "cudaq/qis/pauli_word.h"
+#include "cudaq/spin_op.h"
+
+#include <cstddef>
+#include <vector>
+
+namespace cudaq::solvers {
+
+inline constexpr int first_order_trotter = 1;
+inline constexpr int second_order_trotter = 2;
+inline constexpr int fourth_order_trotter = 4;
+
+inline constexpr double forest_ruth_w1 = 1.3512071919596578;
+inline constexpr double forest_ruth_w0 = -1.7024143839193153;
+
+/// @brief Flattened Hamiltonian representation used by Trotter routines.
+/// @details Identity terms are stored separately. For H = c I + H',
+/// apply_trotter() implements the product-formula circuit for H' and omits the
+/// phase exp(-i c t). That phase is globally unobservable for a single
+/// unconditioned state evolution, but it can matter for controlled evolution,
+/// overlaps, phase estimation, Krylov/QEL moments, and other interference-based
+/// algorithms. Callers in those contexts must account for identity_coefficient.
+struct trotter_terms {
+  std::vector<double> coefficients;
+  std::vector<cudaq::pauli_word> words;
+  double identity_coefficient = 0.0;
+  std::size_t num_qubits = 0;
+};
+
+/// @brief Extract real non-identity Pauli terms from a spin operator.
+/// @details The returned identity_coefficient is not consumed by
+/// apply_trotter(). It is returned so higher-level algorithms can preserve the
+/// omitted exp(-i * identity_coefficient * time) phase when that phase becomes
+/// observable as a relative phase.
+/// @throws std::invalid_argument for coefficients with non-negligible imaginary
+/// parts.
+trotter_terms make_trotter_terms(const cudaq::spin_op &hamiltonian,
+                                 double coefficient_tolerance = 1e-12);
+
+/// @brief Apply Suzuki-Trotter evolution to a live quantum register.
+/// @details This QPU-facing primitive consumes flattened Pauli terms rather
+/// than cudaq::spin_op so that kernels only receive device-lowerable data.
+/// Identity terms are intentionally omitted. Thus, if the original Hamiltonian
+/// was H = c I + H', this applies a product-formula approximation to
+/// exp(-i H' t), not the full exp(-i H t). The omitted exp(-i c t) phase is
+/// harmless for ordinary expectation values of a single unconditioned evolved
+/// state, but must be tracked or reintroduced by controlled/interference-based
+/// algorithms.
+__qpu__ void apply_trotter(const std::vector<double> &coefficients,
+                           const std::vector<cudaq::pauli_word> &words,
+                           double time, std::size_t steps, int order,
+                           cudaq::qview<> qubits);
+
+} // namespace cudaq::solvers

libs/solvers/lib/CMakeLists.txt

@@ -31,6 +31,11 @@ add_library(cudaq-solvers SHARED
   version.cpp
 )
 
+cudaqx_add_device_code(cudaq-solvers
+  SOURCES
+    hamiltonian_simulation/trotter.cpp
+)
+
 add_subdirectory(adapt)
 add_subdirectory(optimizers)
 add_subdirectory(stateprep)

libs/solvers/lib/hamiltonian_simulation/trotter.cpp

@@ -0,0 +1,91 @@
+/*******************************************************************************
+ * Copyright (c) 2026 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/solvers/hamiltonian_simulation/trotter.h"
+
+#include <cmath>
+#include <stdexcept>
+
+namespace cudaq::solvers {
+
+trotter_terms make_trotter_terms(const cudaq::spin_op &hamiltonian,
+                                 double coefficient_tolerance) {
+  if (coefficient_tolerance < 0.0)
+    throw std::invalid_argument(
+        "trotter error - coefficient tolerance must be non-negative.");
+
+  trotter_terms terms;
+  terms.num_qubits = hamiltonian.num_qubits();
+
+  for (const auto &term : hamiltonian) {
+    const auto coefficient = term.evaluate_coefficient();
+    if (std::abs(coefficient.imag()) > coefficient_tolerance)
+      throw std::invalid_argument(
+          "trotter error - only real Hamiltonian coefficients are supported.");
+
+    const auto real_coefficient = coefficient.real();
+    if (term.is_identity()) {
+      // Preserve identity contributions for callers that need the omitted
+      // exp(-i c t) phase in controlled or interference-based algorithms.
+      terms.identity_coefficient += real_coefficient;
+      continue;
+    }
+
+    terms.coefficients.push_back(real_coefficient);
+    terms.words.push_back(term.get_pauli_word(terms.num_qubits));
+  }
+
+  return terms;
+}
+
+__qpu__ void apply_trotter(const std::vector<double> &coefficients,
+                           const std::vector<cudaq::pauli_word> &words,
+                           double time, std::size_t steps, int order,
+                           cudaq::qview<> qubits) {
+  if (steps == 0 || coefficients.size() != words.size())
+    return;
+
+  if (order != 1 && order != 2 && order != 4)
+    return;
+
+  const double dt = time / static_cast<double>(steps);
+  for (std::size_t step = 0; step < steps; ++step) {
+    if (order == 1) {
+      for (std::size_t i = 0; i < words.size(); ++i)
+        exp_pauli(-dt * coefficients[i], qubits, words[i]);
+    } else if (order == 4) {
+      for (std::size_t i = 0; i < words.size(); ++i)
+        exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i], qubits,
+                  words[i]);
+      for (std::size_t i = words.size(); i > 0; --i)
+        exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i - 1], qubits,
+                  words[i - 1]);
+
+      for (std::size_t i = 0; i < words.size(); ++i)
+        exp_pauli(-0.5 * -1.7024143839193153 * dt * coefficients[i], qubits,
+                  words[i]);
+      for (std::size_t i = words.size(); i > 0; --i)
+        exp_pauli(-0.5 * -1.7024143839193153 * dt * coefficients[i - 1], qubits,
+                  words[i - 1]);
+
+      for (std::size_t i = 0; i < words.size(); ++i)
+        exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i], qubits,
+                  words[i]);
+      for (std::size_t i = words.size(); i > 0; --i)
+        exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i - 1], qubits,
+                  words[i - 1]);
+    } else {
+      for (std::size_t i = 0; i < words.size(); ++i)
+        exp_pauli(-0.5 * dt * coefficients[i], qubits, words[i]);
+      for (std::size_t i = words.size(); i > 0; --i)
+        exp_pauli(-0.5 * dt * coefficients[i - 1], qubits, words[i - 1]);
+    }
+  }
+}
+
+} // namespace cudaq::solvers

libs/solvers/python/bindings/solvers/py_solvers.cpp

@@ -5,6 +5,8 @@
  * This source code and the accompanying materials are made available under    *
  * the terms of the Apache License 2.0 which accompanies this distribution.    *
  ******************************************************************************/
+#include <algorithm>
+#include <complex>
 #include <limits>
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
@@ -20,6 +22,7 @@
 #include "cudaq/python/PythonCppInterop.h"
 
 #include "cudaq/solvers/adapt.h"
+#include "cudaq/solvers/hamiltonian_simulation/trotter.h"
 #include "cudaq/solvers/qaoa.h"
 #include "cudaq/solvers/stateprep/ceo.h"
 #include "cudaq/solvers/stateprep/uccgsd.h"
@@ -818,9 +821,67 @@ operator pool. Only integer and list configuration values are currently supporte
 )#");
 }
 
+void bindTrotter(nb::module_ &mod) {
+  auto terms_to_tuple = [](const cudaq::spin_op &hamiltonian,
+                           double coefficient_tolerance) {
+    auto terms =
+        cudaq::solvers::make_trotter_terms(hamiltonian, coefficient_tolerance);
+    return nb::make_tuple(terms.coefficients, terms.words,
+                          terms.identity_coefficient, terms.num_qubits);
+  };
+
+  mod.def(
+      "make_trotter_terms", terms_to_tuple, nb::arg("hamiltonian"),
+      nb::arg("coefficient_tolerance") = 1e-12,
+      R"(Return flattened non-identity coefficients, Pauli words, identity coefficient, and qubit count.
+
+The identity coefficient is returned separately because apply_trotter() omits
+that global phase. For H = c I + H', apply_trotter() approximates exp(-i H' t).
+Callers using controlled evolution, overlaps, phase estimation, Krylov/QEL
+moments, or other interference-based algorithms must account for the omitted
+exp(-i c t) phase when it becomes a relative phase.)");
+
+  mod.def(
+      "make_trotter_terms",
+      [terms_to_tuple](const cudaq::spin_op_term &hamiltonian,
+                       double coefficient_tolerance) {
+        return terms_to_tuple(cudaq::spin_op(hamiltonian),
+                              coefficient_tolerance);
+      },
+      nb::arg("hamiltonian"), nb::arg("coefficient_tolerance") = 1e-12);
+
+  mod.def(
+      "apply_trotter",
+      [](const std::vector<double> &, const std::vector<cudaq::pauli_word> &,
+         double, std::size_t, int, cudaq::qview<>) {},
+      R"(Apply Suzuki-Trotter evolution inside a CUDA-Q kernel.
+
+The QPU-facing form consumes flattened terms:
+
+    apply_trotter(coefficients, pauli_words, time, steps, order, qubits)
+
+Use make_trotter_terms(H) on the host to extract coefficients and words from a
+SpinOperator before passing them to a kernel. Identity terms are not applied by
+this primitive; use the returned identity coefficient to track or reintroduce
+the omitted phase in controlled or interference-based algorithms.)");
+
+  cudaq::python::registerDeviceKernel(
+      nb::cast<std::string>(mod.attr("__name__")), "apply_trotter",
+      cudaq::python::getMangledArgsString<
+          const std::vector<double> &, const std::vector<cudaq::pauli_word> &,
+          double, std::size_t, int, cudaq::qview<>>());
+
+  cudaq::python::addDeviceKernelInterop<
+      const std::vector<double> &, const std::vector<cudaq::pauli_word> &,
+      double, std::size_t, int, cudaq::qview<>>(
+      mod, "hamiltonian_simulation", "apply_trotter",
+      "Apply Suzuki-Trotter evolution inside a CUDA-Q kernel.");
+}
+
 void bindSolvers(nb::module_ &mod) {
 
   addStatePrepKernels(mod);
+  bindTrotter(mod);
 
   auto solvers = mod; //.def_submodule("solvers");
   bindOperators(solvers);

libs/solvers/python/cudaq_solvers/__init__.py

@@ -9,10 +9,82 @@
 try:
     from ._pycudaqx_solvers_the_suffix_matters_cudaq_solvers import *
     from ._pycudaqx_solvers_the_suffix_matters_cudaq_solvers import __version__
+    from ._pycudaqx_solvers_the_suffix_matters_cudaq_solvers import make_trotter_terms as _cpp_make_trotter_terms
 except ImportError as _e:
     from ._system_dep_check import raise_if_missing_system_dep
     raise_if_missing_system_dep(_e, "cudaq-solvers")
     raise
+
+
+def _maybe_call(value):
+    return value() if callable(value) else value
+
+
+def _is_python_spin_operator(value):
+    return hasattr(value, "term_count") and hasattr(value, "__iter__")
+
+
+def _is_python_spin_term(value):
+    return hasattr(value, "evaluate_coefficient") and hasattr(
+        value, "get_pauli_word")
+
+
+def _term_coefficient(term, coefficient_tolerance):
+    coefficient = term.evaluate_coefficient()
+    if abs(coefficient.imag) > coefficient_tolerance:
+        raise ValueError(
+            "trotter error - only real Hamiltonian coefficients are supported.")
+    return float(coefficient.real)
+
+
+def _term_qubit_extent(term):
+    max_degree = _maybe_call(getattr(term, "max_degree", -1))
+    return max_degree + 1 if max_degree >= 0 else 0
+
+
+def make_trotter_terms(hamiltonian, coefficient_tolerance=1e-12):
+    """Return flattened terms for Suzuki-Trotter circuit primitives.
+
+    Returns ``(coefficients, words, identity_coefficient, num_qubits)`` where
+    ``words`` are padded Pauli strings suitable for CUDA-Q kernel arguments.
+
+    ``apply_trotter`` omits identity terms. For ``H = c I + H'``, it applies a
+    product-formula approximation to ``exp(-i H' t)`` and leaves the phase
+    ``exp(-i c t)`` to the caller. This phase cancels in ordinary expectation
+    values of one unconditioned evolved state, but it can matter for controlled
+    evolution, overlaps, phase estimation, Krylov/QEL moments, and other
+    interference-based algorithms.
+    """
+    if coefficient_tolerance < 0.0:
+        raise ValueError(
+            "trotter error - coefficient tolerance must be non-negative.")
+
+    if _is_python_spin_term(hamiltonian):
+        num_qubits = _term_qubit_extent(hamiltonian)
+        terms = [hamiltonian]
+    elif _is_python_spin_operator(hamiltonian):
+        num_qubits = int(_maybe_call(getattr(hamiltonian, "qubit_count", 0)))
+        terms = list(hamiltonian)
+    else:
+        return _cpp_make_trotter_terms(hamiltonian, coefficient_tolerance)
+
+    coefficients = []
+    words = []
+    identity_coefficient = 0.0
+
+    for term in terms:
+        coefficient = _term_coefficient(term, coefficient_tolerance)
+        term_extent = _term_qubit_extent(term)
+        num_qubits = max(num_qubits, term_extent)
+        if term.is_identity():
+            identity_coefficient += coefficient
+            continue
+        words.append(term.get_pauli_word(num_qubits))
+        coefficients.append(coefficient)
+
+    return coefficients, words, identity_coefficient, num_qubits
+
+
 try:
     from .gqe_algorithm.gqe import gqe
 except ImportError: