diff --git a/libs/solvers/include/cudaq/solvers/adapt.h b/libs/solvers/include/cudaq/solvers/adapt.h
index d25be1d6..eba17627 100644
--- a/libs/solvers/include/cudaq/solvers/adapt.h
+++ b/libs/solvers/include/cudaq/solvers/adapt.h
@@ -11,6 +11,7 @@
 #include "cudaq/qis/qubit_qis.h"
 #include "cudaq/solvers/operators/operator_pool.h"
 #include "cudaq/solvers/vqe.h"
+#include "cudaq/solvers/adapt/variant.h"
 
 #include <functional>
 
@@ -121,6 +122,7 @@ class adapt_impl : public cudaqx::extension_point<adapt_impl> {
 ///      It can be either "warm", or "cold". [default: "cold"]
 ///  - "tol" (double): Tolerance for optimization [default: 1e-12]
 /// @return Result of the ADAPT-VQE algorithm
+[[deprecated("cudaq::solvers::adapt_vqe is deprecated; use cudaq::solvers::adapt(..., AdaptVariant::VQE, ...) instead.")]]
 static inline adapt::result
 adapt_vqe(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
           const spin_op &H, const std::vector<spin_op> &poolList,
@@ -154,6 +156,7 @@ adapt_vqe(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
 ///      It can be either "warm", or "cold". [default: "cold"]
 ///  - "tol" (double): Tolerance for optimization [default: 1e-12]
 /// @return Result of the ADAPT-VQE algorithm
+[[deprecated("cudaq::solvers::adapt_vqe is deprecated; use cudaq::solvers::adapt(..., AdaptVariant::VQE, ...) instead.")]]
 static inline adapt::result
 adapt_vqe(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
           const spin_op &H, const std::vector<spin_op> &poolList,
@@ -188,6 +191,7 @@ adapt_vqe(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
 ///      It can be either "warm", or "cold". [default: "cold"]
 ///  - "tol" (double): Tolerance for optimization [default: 1e-12]
 /// @return Result of the ADAPT-VQE algorithm
+[[deprecated("cudaq::solvers::adapt_vqe is deprecated; use cudaq::solvers::adapt(..., AdaptVariant::VQE, ...) instead.")]]
 static inline adapt::result
 adapt_vqe(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
           const spin_op &H, const std::vector<spin_op> &poolList,
@@ -199,4 +203,37 @@ adapt_vqe(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
   return impl->run(initialState, H, poolList, optimizer, gradient, options);
 }
 
+// Generic ADAPT front-door with optimizer & gradient
+static inline adapt::result
+adapt(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
+      const spin_op &H, const std::vector<spin_op> &poolList,
+      const optim::optimizer &optimizer, AdaptVariant variant,
+      const std::string &gradient = "",
+      const heterogeneous_map options = heterogeneous_map()) {
+  auto opts = options;
+  opts.insert("adapt_variant", static_cast<int>(variant));
+  auto &platform = cudaq::get_platform();
+  auto impl =
+      adapt::adapt_impl::get(platform.is_simulator() ? "simulator" : "remote");
+  return impl->run(initialState, H, poolList, optimizer, gradient, opts);
+}
+
+// Overload without gradient
+static inline adapt::result
+adapt(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
+      const spin_op &H, const std::vector<spin_op> &poolList,
+      const optim::optimizer &optimizer, AdaptVariant variant,
+      const heterogeneous_map options = heterogeneous_map()) {
+  return adapt(initialState, H, poolList, optimizer, variant, "", options);
+}
+
+// Overload with default optimizer
+static inline adapt::result
+adapt(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
+      const spin_op &H, const std::vector<spin_op> &poolList,
+      AdaptVariant variant, const heterogeneous_map options = heterogeneous_map()) {
+  auto opt = optim::optimizer::get("cobyla");
+  return adapt(initialState, H, poolList, *opt, variant, "", options);
+}
+
 } // namespace cudaq::solvers
diff --git a/libs/solvers/include/cudaq/solvers/adapt/variant.h b/libs/solvers/include/cudaq/solvers/adapt/variant.h
new file mode 100644
index 00000000..afb5ccc7
--- /dev/null
+++ b/libs/solvers/include/cudaq/solvers/adapt/variant.h
@@ -0,0 +1,52 @@
+#pragma once
+#include <string_view>
+
+namespace cudaq::solvers::adapt {
+
+enum class AdaptVariant : int {
+  VQE  = 0,
+  QAOA = 1,
+};
+
+inline constexpr std::string_view to_string(AdaptVariant v) {
+  switch (v) {
+    case AdaptVariant::VQE:  return "vqe";
+    case AdaptVariant::QAOA: return "qaoa";
+    default:                 return "unknown";
+  }
+}
+
+// ASCII helpers (locale-independent)
+constexpr unsigned char ascii_lower(unsigned char c) {
+  return (c >= 'A' && c <= 'Z') ? static_cast<unsigned char>(c + ('a' - 'A')) : c;
+}
+constexpr bool ascii_isspace(unsigned char c) {
+  return c == ' ' || c == '\t' || c == '\n' || c == '\r' || c == '\f' || c == '\v';
+}
+
+inline constexpr std::string_view trim_ascii(std::string_view s) {
+  size_t b = 0, e = s.size();
+  while (b < e && ascii_isspace(static_cast<unsigned char>(s[b]))) ++b;
+  while (e > b && ascii_isspace(static_cast<unsigned char>(s[e - 1]))) --e;
+  return s.substr(b, e - b);
+}
+
+// ASCII-only case-insensitive equals
+inline constexpr bool iequals_ascii(std::string_view a, std::string_view b) {
+  if (a.size() != b.size()) return false;
+  for (size_t i = 0; i < a.size(); ++i) {
+    if (ascii_lower(static_cast<unsigned char>(a[i])) !=
+        ascii_lower(static_cast<unsigned char>(b[i])))
+      return false;
+  }
+  return true;
+}
+
+inline bool parse_variant(std::string_view s, AdaptVariant &out) {
+  s = trim_ascii(s);
+  if (iequals_ascii(s, "vqe"))  { out = AdaptVariant::VQE;  return true; }
+  if (iequals_ascii(s, "qaoa")) { out = AdaptVariant::QAOA; return true; }
+  return false;
+}
+
+} // namespace cudaq::solvers::adapt
diff --git a/libs/solvers/lib/adapt/adapt_simulator.cpp b/libs/solvers/lib/adapt/adapt_simulator.cpp
index 047e0e56..9d9c4614 100644
--- a/libs/solvers/lib/adapt/adapt_simulator.cpp
+++ b/libs/solvers/lib/adapt/adapt_simulator.cpp
@@ -15,6 +15,7 @@
 #include "device/prepare_state.h"
 #include "cudaq/solvers/adapt/adapt_simulator.h"
 #include "cudaq/solvers/vqe.h"
+#include "cudaq/solvers/adapt/variant.h"
 
 #include <nlohmann/json.hpp>
 
@@ -36,6 +37,26 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
   double latestEnergy = std::numeric_limits<double>::max();
   double ediff = std::numeric_limits<double>::max();
 
+  auto variant = adapt::AdaptVariant::VQE;
+  if (options.contains("adapt_variant")) {
+    // String first
+    try {
+      auto s = options.get<std::string>("adapt_variant");
+      adapt::AdaptVariant tmp;
+      if (adapt::parse_variant(s, tmp)) variant = tmp;
+    } catch (...) {
+      // Int fallback
+      try {
+        int v = options.get<int>("adapt_variant");
+        variant = static_cast<adapt::AdaptVariant>(v);
+      } catch (...) {}
+    }
+  }
+
+  if (variant == adapt::AdaptVariant::QAOA) {
+    throw std::runtime_error("ADAPT-QAOA not yet implemented.");
+  }
+
   int maxIter = options.get<int>("max_iter", 30);
   auto grad_norm_tolerance = options.get<double>("grad_norm_tolerance", 1e-5);
   auto tolNormDiff = options.get<double>("grad_norm_diff_tolerance", 1e-5);
@@ -44,6 +65,7 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
   auto mutable_options = options;
 
   auto numQubits = H.num_qubits();
+
   // Assumes each rank can see numQpus, models a distributed
   // architecture where each rank is a compute node, and each node
   // has numQpus GPUs available. Each GPU is indexed 0, 1, 2, ..
@@ -51,17 +73,11 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
   std::size_t numRanks =
       cudaq::mpi::is_initialized() ? cudaq::mpi::num_ranks() : 1;
   std::size_t rank = cudaq::mpi::is_initialized() ? cudaq::mpi::rank() : 0;
+  
   double energy = 0.0, lastNorm = std::numeric_limits<double>::max();
 
-  // poolList is split into numRanks chunks, and each chunk can be
-  // further parallelized across numQpus.
-  // Compute the [H,Oi]
   std::vector<spin_op> commutators;
   std::size_t total_elements = pool.size();
-  std::size_t elements_per_rank = total_elements / numRanks;
-  std::size_t remainder = total_elements % numRanks;
-  std::size_t start = rank * elements_per_rank + std::min(rank, remainder);
-  std::size_t end = start + elements_per_rank + (rank < remainder ? 1 : 0);
 
   // Check if operator has only imaginary coefficients
   // checking the first one is enough, we assume the pool is homogeneous
@@ -70,7 +86,7 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
       (std::abs(c.real()) <= 1e-9) && (std::abs(c.imag()) > 1e-9);
   auto coeff = (!isImaginary) ? std::complex<double>{0.0, 1.0}
                               : std::complex<double>{1.0, 0.0};
-
+  
   for (auto &op : pool) {
     auto commutator = H * op - op * H;
     commutator.canonicalize().trim();
@@ -78,22 +94,11 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
       commutators.push_back(coeff * commutator);
   }
 
-  nlohmann::json initInfo = {{"num-qpus", numQpus},
-                             {"numRanks", numRanks},
-                             {"num-pool-elements", pool.size()},
-                             {"num-elements-per-rank", end - start}};
-  if (rank == 0)
-    cudaq::info("[adapt] init info: {}", initInfo.dump(4));
-
-  // We'll need to know the local to global index map
-  std::vector<std::size_t> localToGlobalMap(end - start);
-  for (int i = 0; i < end - start; i++)
-    localToGlobalMap[i] = start + i;
-
   // Start of with the initial |psi_n>
   cudaq::state state = get_state(adapt_kernel, numQubits, initialState, thetas,
                                  coefficients, pauliWords, poolIndices);
 
+
   int step = 0;
   while (true) {
     if (options.get("verbose", false))
@@ -105,29 +110,16 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
       break;
     }
     step++;
-
     // Step 1 - compute <psi|[H,Oi]|psi> vector
     std::vector<double> gradients;
     std::vector<observe_result> results;
     double gradNorm = 0.0;
-    std::vector<async_observe_result> resultHandles;
 
     if (numQpus == 1) {
       for (std::size_t i = 0; i < commutators.size(); i++) {
         cudaq::info("Compute commutator {}", i);
         results.emplace_back(observe(prepare_state, commutators[i], state));
       }
-    } else {
-      for (std::size_t i = 0, qpuCounter = 0; i < commutators.size(); i++) {
-        if (rank == 0)
-          cudaq::info("Compute commutator {}", i);
-        if (qpuCounter % numQpus == 0)
-          qpuCounter = 0;
-        resultHandles.emplace_back(
-            observe_async(qpuCounter++, prepare_state, commutators[i], state));
-      }
-      for (auto &handle : resultHandles)
-        results.emplace_back(handle.get());
     }
 
     // Get the gradient results
@@ -141,43 +133,10 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
       norm += g * g;
     norm = std::sqrt(norm);
 
-    // All ranks have a norm, need to reduce that across all
-    if (mpi::is_initialized())
-      norm = cudaq::mpi::all_reduce(norm, std::plus<double>());
-
     auto iter = std::max_element(gradients.begin(), gradients.end());
     double maxGrad = *iter;
     auto maxOpIdx = std::distance(gradients.begin(), iter);
 
-    if (mpi::is_initialized()) {
-      std::vector<int> allMaxOpIndices(numRanks);
-      std::vector<double> allMaxGrads(numRanks);
-      // Distribute the max gradient from this rank to others
-      cudaq::mpi::all_gather(allMaxGrads, {*iter});
-      // Distribute the corresponding idx from this rank to others,
-      // make sure we map back to global indices
-      cudaq::mpi::all_gather(allMaxOpIndices,
-                             {static_cast<int>(localToGlobalMap[maxOpIdx])});
-
-      // Everyone has the indices, loop over and pick out the
-      // max from all calculations
-      std::size_t cachedIdx = 0;
-      double cachedGrad = 0.0;
-      for (std::size_t i = 0; i < allMaxGrads.size(); i++)
-        if (allMaxGrads[i] > cachedGrad) {
-          cachedGrad = allMaxGrads[i];
-          cachedIdx = allMaxOpIndices[i];
-        }
-
-      maxOpIdx = cachedIdx;
-    }
-
-    if (rank == 0) {
-      cudaq::info("[adapt] index of element with max gradient is {}", maxOpIdx);
-      cudaq::info("current norm is {} and last iteration norm is {}", norm,
-                  lastNorm);
-    }
-
     // Convergence is reached if gradient values are small
     if (norm < grad_norm_tolerance ||
         std::fabs(lastNorm - norm) < tolNormDiff || ediff < thresholdE)
@@ -259,6 +218,7 @@ simulator::run(const cudaq::qkernel<void(cudaq::qvector<> &)> &initialState,
   }
 
   return std::make_tuple(energy, thetas, chosenOps);
-}
 
-} // namespace cudaq::solvers::adapt
+  }
+
+}
\ No newline at end of file