diff --git a/QHyper/polynomial.py b/QHyper/polynomial.py
index 60befba..f0980a6 100644
--- a/QHyper/polynomial.py
+++ b/QHyper/polynomial.py
@@ -64,11 +64,16 @@ def __init__(self, terms: dict[tuple[str, ...], float] | float | int
 
         self.terms = defaultdict(float)
 
+        non_zero_found = False
         for term, coefficient in terms.items():
             if coefficient == 0:
                 continue
+            non_zero_found = True
             self.terms[tuple(sorted(term))] += coefficient
 
+        if not non_zero_found:
+            self.terms[tuple()] = 0.0
+
     @overload
     def __add__(self, other: float | int) -> 'Polynomial': ...
 
@@ -204,6 +209,8 @@ def degree(self) -> int:
         int
             The degree of the polynomial.
         """
+        if not self.terms:
+            return 0
         return max(len(term) for term in self.terms)
 
     def get_variables(self) -> set[str]:
diff --git a/QHyper/problems/community_detection.py b/QHyper/problems/community_detection.py
index 29e5e4c..9247d94 100644
--- a/QHyper/problems/community_detection.py
+++ b/QHyper/problems/community_detection.py
@@ -22,34 +22,35 @@ class Network:
     resolution: float = 1.0
     weight: str | None = "weight"
     community: list | None = None
-    full_modularity_matrix: np.ndarray = field(init=False)
+    full_modularity_matrix: np.ndarray | None = None
     generalized_modularity_matrix: np.ndarray = field(init=False)
 
     def __post_init__(self) -> None:
         if not self.community:
             self.community = [*range(self.graph.number_of_nodes())]
-        (
-            self.full_modularity_matrix,
-            self.generalized_modularity_matrix,
-        ) = self.calculate_modularity_matrix()
+        if self.full_modularity_matrix is None:
+            self.full_modularity_matrix = self.calculate_full_modularity_matrix()
+        self.generalized_modularity_matrix = (
+            self.calculate_generalized_modularity_matrix()
+        )
 
-    def calculate_modularity_matrix(self) -> np.ndarray:
+    def calculate_full_modularity_matrix(self) -> np.ndarray:
         adj_matrix: np.ndarray = nx.to_numpy_array(self.graph, weight=self.weight)
         in_degree_matrix: np.ndarray = adj_matrix.sum(axis=1)
         out_degree_matrix: np.ndarray = adj_matrix.sum(axis=0)
         m: int = np.sum(adj_matrix)
-        full_modularity_matrix = (
+        return (
             adj_matrix
             - self.resolution * np.outer(in_degree_matrix, out_degree_matrix) / m
         )
 
-        B_bis = full_modularity_matrix[self.community, :]
+    def calculate_generalized_modularity_matrix(self) -> np.ndarray:
+        B_bis = self.full_modularity_matrix[self.community, :]
         B_community = B_bis[:, self.community]
         B_i = np.sum(B_community, axis=1)
         B_j = np.sum(B_community.T, axis=1)
         delta = np.eye(len(self.community), dtype=np.int32)
-        B_g = 0.5*( B_community + B_community.T ) - 0.5 * delta * (B_i + B_j)
-        return full_modularity_matrix, B_g
+        return 0.5 * (B_community + B_community.T) - 0.5 * delta * (B_i + B_j)
 
 
 class KarateClubNetwork(Network):
@@ -131,9 +132,9 @@ def __init__(
             self._set_objective_function()
             self._set_one_hot_constraints(communities)
         else:
-            self.variables: tuple[
-                sympy.Symbol
-            ] = self._get_discrete_variable_representation()
+            self.variables: tuple[sympy.Symbol] = (
+                self._get_discrete_variable_representation()
+            )
             self._set_objective_function()
 
     def _get_discrete_variable_representation(
@@ -164,6 +165,12 @@ def _set_objective_function(self) -> None:
 
         equation = {key: -1 * val for key, val in equation.items()}
 
+        nonzero_terms = sum(1 for v in equation.values() if not np.isclose(v, 0.0))
+        if nonzero_terms == 0:
+            raise ValueError(
+                f"The objective function is a zero polynomial - all terms in the generalized modularity matrix are 0. Try different resolution (current: {self.resolution})"
+            )
+
         self.objective_function = Polynomial(equation)
 
     def _encode_discrete_to_one_hot(
diff --git a/QHyper/solvers/quantum_annealing/dwave/advantage.py b/QHyper/solvers/quantum_annealing/dwave/advantage.py
index 481ba1d..19c05de 100644
--- a/QHyper/solvers/quantum_annealing/dwave/advantage.py
+++ b/QHyper/solvers/quantum_annealing/dwave/advantage.py
@@ -44,6 +44,8 @@ class Advantage(Solver):
 
     problem: Problem
     penalty_weights: list[float] | None = None
+    version: str | None = None
+    region: str | None = None
     num_reads: int = 1
     chain_strength: float | None = None
     token: str | None = None
@@ -51,6 +53,8 @@ class Advantage(Solver):
     def __init__(self,
                  problem: Problem,
                  penalty_weights: list[float] | None = None,
+                 version: str | None = None,
+                 region: str | None = None,
                  num_reads: int = 1,
                  chain_strength: float | None = None,
                  use_clique_embedding: bool = False,
@@ -61,12 +65,18 @@ def __init__(self,
         self.num_reads = num_reads
         self.chain_strength = chain_strength
         self.use_clique_embedding = use_clique_embedding
-        self.sampler = DWaveSampler(
+        if (self.version and not self.region) or (self.region and not self.version):
+            raise ValueError("Both 'version' and 'region' must be specified together.")
+        if self.version and self.region:
+            self.sampler = DWaveSampler(
+            solver=self.version, region=self.region,
             token=token or DWAVE_API_TOKEN, **config)
+        else:
+            self.sampler = DWaveSampler(token=token or DWAVE_API_TOKEN, **config)
         self.token = token
 
         if use_clique_embedding:
-            args = self.weigths if self.weigths else []
+            args = self.penalty_weights if self.penalty_weights else []
             qubo = Converter.create_qubo(self.problem, args)
             qubo_terms, offset = convert_qubo_keys(qubo)
             bqm = BinaryQuadraticModel.from_qubo(qubo_terms, offset=offset)
diff --git a/requirements/prod.txt b/requirements/prod.txt
index 8feee23..ad58b1b 100644
--- a/requirements/prod.txt
+++ b/requirements/prod.txt
@@ -6,4 +6,5 @@ dwave-system
 gurobipy
 types-tqdm
 pandas
+autoray==0.6.11
 git+https://github.com/wfcommons/wfcommons.git@main