Add uccgsd ansatz (#319)

In this PR:
1- We added UCCGSD in operator_pools for algorithm such as ADAPT-VQE
2- We added UCCGSD in stateprep for algorithm such as VQE

---------

Signed-off-by: marwafar <[email protected]>
Signed-off-by: Scott Thornton <[email protected]>
Co-authored-by: Scott Thornton <[email protected]>
Co-authored-by: Scott Thornton <[email protected]>

Author: marwafar

libs/solvers/include/cudaq/solvers/operators/operator_pools/uccgsd_operator_pool.h

@@ -0,0 +1,43 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 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.    *
+ ******************************************************************************/
+
+#pragma once
+
+#include "../operator_pool.h"
+
+namespace cudaq::solvers {
+
+/// @brief UCCGSD operator pool is a class which generates a pool
+/// of UCCGSD operators for use in quantum algorithms, like Adapt-VQE.
+/// @details This class extends the operator_pool interface
+/// therefore inherits the extension_point template, allowing for
+/// runtime extensibility.
+class uccgsd : public operator_pool {
+
+public:
+  /// @brief Generate a vector of spin operators based on the provided
+  /// configuration.
+  /// @details The UCCGSD operator pool is generated with an imaginary factor
+  /// 'i' in the coefficients of the operators, which simplifies the use for
+  /// running Adapt-VQE.
+  /// @param config A heterogeneous map containing configuration parameters for
+  /// operator generation.
+  /// @return A vector of cudaq::spin_op objects representing the generated
+  /// operator pool.
+  std::vector<cudaq::spin_op>
+  generate(const heterogeneous_map &config) const override;
+
+  /// @brief Call to macro for defining the creator function for an extension
+  /// @details This function is used by the extension point mechanism to create
+  /// instances of the uccgsd class.
+  CUDAQ_EXTENSION_CREATOR_FUNCTION(operator_pool, uccgsd)
+};
+/// @brief Register the uccgsd extension type with the CUDA-Q framework
+CUDAQ_REGISTER_TYPE(uccgsd)
+
+} // namespace cudaq::solvers

libs/solvers/include/cudaq/solvers/operators/uccgsd_excitation_utils.h

@@ -0,0 +1,55 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 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.    *
+ ******************************************************************************/
+
+#pragma once
+
+#include "cudaq/spin_op.h"
+#include <utility>
+#include <vector>
+
+namespace cudaq::solvers {
+
+/// @brief Generate all unique single excitation index pairs (p, q) with p > q
+/// @param numQubits Number of qubits in the system
+/// @return Vector of (p, q) pairs representing single excitations
+std::vector<std::pair<std::size_t, std::size_t>>
+generate_uccgsd_singles(std::size_t numQubits);
+
+/// @brief Generate all unique double excitation index pairs
+/// @details For 4 indices a < b < c < d, generates all 3 unique pairings:
+///          (a,b)-(c,d), (a,c)-(b,d), (a,d)-(b,c)
+///          Each pair is normalized so that within each pair, first > second,
+///          and pairs are ordered.
+/// @param numQubits Number of qubits in the system
+/// @return Vector of ((p,q), (r,s)) pairs where p>q, r>s, representing double
+/// excitations
+std::vector<std::pair<std::pair<std::size_t, std::size_t>,
+                      std::pair<std::size_t, std::size_t>>>
+generate_uccgsd_doubles(std::size_t numQubits);
+
+/// @brief Add a UCCGSD single excitation operator to the operator pool
+/// @details Generates the operator: 0.5 * (Y_q * Z_parity * X_p - X_q *
+/// Z_parity * Y_p)
+///          where Z_parity is the product of Z operators between q and p
+/// @param ops Vector to append the operator to
+/// @param p Higher qubit index (p > q)
+/// @param q Lower qubit index (p > q)
+void addUCCGSDSingleExcitation(std::vector<cudaq::spin_op> &ops, std::size_t p,
+                               std::size_t q);
+
+/// @brief Add a UCCGSD double excitation operator to the operator pool
+/// @details Generates the generalized double excitation operator with 8 terms
+/// @param ops Vector to append the operator to
+/// @param p First qubit index (p > q)
+/// @param q Second qubit index (p > q)
+/// @param r Third qubit index (r > s)
+/// @param s Fourth qubit index (r > s)
+void addUCCGSDDoubleExcitation(std::vector<cudaq::spin_op> &ops, std::size_t p,
+                               std::size_t q, std::size_t r, std::size_t s);
+
+} // namespace cudaq::solvers

libs/solvers/include/cudaq/solvers/stateprep/uccgsd.h

@@ -0,0 +1,41 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 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.    *
+ ******************************************************************************/
+#pragma once
+
+#include "cudaq.h"
+#include <utility>
+#include <vector>
+
+namespace cudaq::solvers::stateprep {
+
+/// @brief Generate UCCGSD operator pool (Python-style unique singles/doubles)
+/// and extract Pauli words and coefficients
+/// @param norbitals Number of spin orbitals (qubits)
+/// @param only_singles If true, only single excitations
+/// @param only_doubles If true, only double excitations
+/// @return Pair of lists: [Pauli words grouped by excitation], [coefficients
+/// grouped by excitation]
+std::pair<std::vector<std::vector<cudaq::pauli_word>>,
+          std::vector<std::vector<double>>>
+get_uccgsd_pauli_lists(std::size_t norbitals, bool only_singles = false,
+                       bool only_doubles = false);
+
+/// \pure_device_kernel
+///
+/// @brief Apply UCCGSD ansatz to a qubit register using grouped Pauli words and
+/// coefficients
+/// @param qubits Qubit register
+/// @param thetas Vector of rotation angles (one per excitation group)
+/// @param pauliWordsList Pauli words grouped by excitation
+/// @param coefficientsList Coefficients grouped by excitation
+__qpu__ void
+uccgsd(cudaq::qview<> qubits, const std::vector<double> &thetas,
+       const std::vector<std::vector<cudaq::pauli_word>> &pauliWordsList,
+       const std::vector<std::vector<double>> &coefficientsList);
+
+} // namespace cudaq::solvers::stateprep

libs/solvers/lib/CMakeLists.txt

@@ -23,6 +23,8 @@ add_library(cudaq-solvers SHARED
   operators/operator_pools/spin_complement_gsd.cpp
   operators/operator_pools/uccsd_operator_pool.cpp
   operators/operator_pools/qaoa_operator_pool.cpp
+  operators/operator_pools/uccgsd_operator_pool.cpp
+  operators/uccgsd_excitation_utils.cpp
   version.cpp
 )
 

libs/solvers/lib/operators/operator_pools/uccgsd_operator_pool.cpp

@@ -0,0 +1,42 @@
+/*******************************************************************************
+ * Copyright (c) 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/solvers/operators/operator_pools/uccgsd_operator_pool.h"
+#include "cudaq/solvers/operators/uccgsd_excitation_utils.h"
+
+using namespace cudaqx;
+
+namespace cudaq::solvers {
+
+std::vector<cudaq::spin_op>
+uccgsd::generate(const heterogeneous_map &config) const {
+
+  auto numOrbitals = config.get<std::size_t>({"num-orbitals", "num_orbitals"});
+  auto numQubits = 2 * numOrbitals;
+  if (numOrbitals == 0)
+    throw std::invalid_argument("num-orbitals must be > 0");
+
+  // For UCCGSD, we do not use alpha/beta/mixed excitations, but generate all
+  // singles and doubles
+  std::vector<cudaq::spin_op> ops;
+
+  // Generate all single excitations using shared utility
+  auto singles = generate_uccgsd_singles(numQubits);
+  for (const auto &[p, q] : singles) {
+    addUCCGSDSingleExcitation(ops, p, q);
+  }
+
+  // Generate all unique unordered double excitations using shared utility
+  auto doubles = generate_uccgsd_doubles(numQubits);
+  for (const auto &[pq, rs] : doubles) {
+    addUCCGSDDoubleExcitation(ops, pq.first, pq.second, rs.first, rs.second);
+  }
+
+  return ops;
+}
+
+} // namespace cudaq::solvers

libs/solvers/lib/operators/uccgsd_excitation_utils.cpp

@@ -0,0 +1,123 @@
+/*******************************************************************************
+ * Copyright (c) 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/solvers/operators/uccgsd_excitation_utils.h"
+#include <algorithm>
+#include <array>
+#include <set>
+
+namespace cudaq::solvers {
+
+std::vector<std::pair<std::size_t, std::size_t>>
+generate_uccgsd_singles(std::size_t numQubits) {
+  std::vector<std::pair<std::size_t, std::size_t>> singles;
+  singles.reserve(numQubits * (numQubits - 1) / 2);
+  for (std::size_t p = 1; p < numQubits; ++p)
+    for (std::size_t q = 0; q < p; ++q)
+      singles.emplace_back(p, q);
+  return singles;
+}
+
+std::vector<std::pair<std::pair<std::size_t, std::size_t>,
+                      std::pair<std::size_t, std::size_t>>>
+generate_uccgsd_doubles(std::size_t numQubits) {
+  std::set<std::pair<std::pair<std::size_t, std::size_t>,
+                     std::pair<std::size_t, std::size_t>>>
+      doubles;
+
+  // Iterate over all combinations of 4 distinct qubits
+  for (std::size_t a = 0; a < numQubits; ++a)
+    for (std::size_t b = a + 1; b < numQubits; ++b)
+      for (std::size_t c = b + 1; c < numQubits; ++c)
+        for (std::size_t d = c + 1; d < numQubits; ++d) {
+          std::array<std::size_t, 4> arr = {a, b, c, d};
+
+          // Generate all 3 unique pairings of the 4 indices
+          std::vector<std::pair<std::pair<std::size_t, std::size_t>,
+                                std::pair<std::size_t, std::size_t>>>
+              pairings = {{{arr[0], arr[1]}, {arr[2], arr[3]}},
+                          {{arr[0], arr[2]}, {arr[1], arr[3]}},
+                          {{arr[0], arr[3]}, {arr[1], arr[2]}}};
+
+          // Normalize and deduplicate each pairing
+          for (auto &pairing : pairings) {
+            auto p1 = pairing.first, p2 = pairing.second;
+
+            // Ensure within each pair: first > second
+            if (p1.first < p1.second)
+              std::swap(p1.first, p1.second);
+            if (p2.first < p2.second)
+              std::swap(p2.first, p2.second);
+
+            // Order the two pairs
+            auto sorted_pairing = std::minmax(p1, p2);
+            doubles.insert({sorted_pairing.first, sorted_pairing.second});
+          }
+        }
+
+  return std::vector<std::pair<std::pair<std::size_t, std::size_t>,
+                               std::pair<std::size_t, std::size_t>>>(
+      doubles.begin(), doubles.end());
+}
+
+void addUCCGSDSingleExcitation(std::vector<cudaq::spin_op> &ops, std::size_t p,
+                               std::size_t q) {
+  if (p > q) {
+    // Compute parity string (Z operators between q and p)
+    cudaq::spin_op_term parity;
+    for (std::size_t i = q + 1; i < p; ++i)
+      parity *= cudaq::spin::z(i);
+
+    std::complex<double> c = {0.5, 0.0};
+
+    // Single excitation: Y_q * Z_parity * X_p - X_q * Z_parity * Y_p
+    ops.emplace_back(c * cudaq::spin::y(q) * parity * cudaq::spin::x(p) -
+                     c * cudaq::spin::x(q) * parity * cudaq::spin::y(p));
+  }
+}
+
+void addUCCGSDDoubleExcitation(std::vector<cudaq::spin_op> &ops, std::size_t p,
+                               std::size_t q, std::size_t r, std::size_t s) {
+  if (p > q && r > s) {
+    // Compute parity strings
+    cudaq::spin_op_term parity_a, parity_b;
+    for (std::size_t i = q + 1; i < p; ++i)
+      parity_a *= cudaq::spin::z(i);
+    for (std::size_t i = s + 1; i < r; ++i)
+      parity_b *= cudaq::spin::z(i);
+
+    std::complex<double> c = {0.125, 0.0};
+
+    // Build the 8-term double excitation operator
+    cudaq::spin_op temp_op;
+
+    // Positive terms
+    temp_op = c * cudaq::spin::y(s) * parity_b * cudaq::spin::x(r) *
+              cudaq::spin::x(q) * parity_a * cudaq::spin::x(p);
+    temp_op += c * cudaq::spin::x(s) * parity_b * cudaq::spin::y(r) *
+               cudaq::spin::x(q) * parity_a * cudaq::spin::x(p);
+    temp_op += c * cudaq::spin::y(s) * parity_b * cudaq::spin::y(r) *
+               cudaq::spin::y(q) * parity_a * cudaq::spin::x(p);
+    temp_op += c * cudaq::spin::y(s) * parity_b * cudaq::spin::y(r) *
+               cudaq::spin::x(q) * parity_a * cudaq::spin::y(p);
+
+    // Negative terms
+    temp_op -= c * cudaq::spin::x(s) * parity_b * cudaq::spin::x(r) *
+               cudaq::spin::y(q) * parity_a * cudaq::spin::x(p);
+    temp_op -= c * cudaq::spin::x(s) * parity_b * cudaq::spin::x(r) *
+               cudaq::spin::x(q) * parity_a * cudaq::spin::y(p);
+    temp_op -= c * cudaq::spin::x(s) * parity_b * cudaq::spin::y(r) *
+               cudaq::spin::y(q) * parity_a * cudaq::spin::y(p);
+    temp_op -= c * cudaq::spin::y(s) * parity_b * cudaq::spin::x(r) *
+               cudaq::spin::y(q) * parity_a * cudaq::spin::y(p);
+
+    ops.emplace_back(temp_op);
+  }
+}
+
+} // namespace cudaq::solvers

libs/solvers/lib/stateprep/CMakeLists.txt

@@ -9,4 +9,5 @@
 cudaqx_add_device_code(cudaq-solvers
   SOURCES
     uccsd.cpp
+    uccgsd.cpp
 )