diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py
index ac8b432..8c46878 100644
--- a/src/qibochem/measurement/optimization.py
+++ b/src/qibochem/measurement/optimization.py
@@ -1,29 +1,34 @@
+"""
+Functions for optimising the measurement cost of obtaining the expectation value
+"""
+
+import networkx as nx
 import numpy as np
 from qibo import gates
 
 
-def check_terms_commutativity(term1, term2, qubitwise=False):
+def check_terms_commutativity(term1, term2, qubitwise):
     """
     Check if terms 1 and 2 are mutually commuting. The 'qubitwise' flag determines if the check is for general
     commutativity (False), or the stricter qubitwise commutativity.
 
     Args:
-        term1/term2: Lists of strings representing a single Pauli term. E.g. ["X0", "Z1", "Y3"]. Obtained from a Qibo
-            SymbolicTerm as ``[factor.name for factor in term.factors]``.
-        qubitwise: Determines if the check is for general commutativity, or the stricter qubitwise commutativity
+        term1/term2: Strings representing a single Pauli term. E.g. "X0 Z1 Y3". Obtained from a Qibo SymbolicTerm as
+        ``" ".join(factor.name for factor in term.factors)``.
+        qubitwise (bool): Determines if the check is for general commutativity, or the stricter qubitwise commutativity
 
     Returns:
         bool: Do terms 1 and 2 commute?
     """
     # Get a list of common qubits for each term
-    common_qubits = sorted(
-        {_term[1:] for _term in term1 if _term[0] != "I"} & {_term[1:] for _term in term2 if _term[0] != "I"}
-    )
+    common_qubits = {_term[1:] for _term in term1.split() if _term[0] != "I"} & {
+        _term[1:] for _term in term2.split() if _term[0] != "I"
+    }
     if not common_qubits:
         return True
     # Get the single Pauli operators for the common qubits for both Pauli terms
-    term1_ops = [_op for _op in term1 if _op[1:] in common_qubits]
-    term2_ops = [_op for _op in term2 if _op[1:] in common_qubits]
+    term1_ops = [_op for _op in term1.split() if _op[1:] in common_qubits]
+    term2_ops = [_op for _op in term2.split() if _op[1:] in common_qubits]
     if qubitwise:
         # Qubitwise: Compare the Pauli terms at the common qubits. Any difference => False
         return all(_op1 == _op2 for _op1, _op2 in zip(term1_ops, term2_ops))
@@ -34,6 +39,45 @@ def check_terms_commutativity(term1, term2, qubitwise=False):
     return n_noncommuting_ops % 2 == 0
 
 
+def group_commuting_terms(terms_list, qubitwise):
+    """
+    Groups the terms in terms_list into as few groups as possible, where all the terms in each group commute
+    mutually == Finding the minimum clique cover (i.e. as few cliques as possible) for the graph whereby each node
+    is a Pauli string, and an edge exists between two nodes iff they commute.
+
+    This is equivalent to the graph colouring problem of the complement graph (i.e. edge between nodes if they DO NOT
+    commute), which this function follows.
+
+    Args:
+        terms_list: List of strings that contain X/Y. The strings should follow the output from
+            ``[factor.name for factor in term.factors]``, where term is a Qibo SymbolicTerm. E.g. ["X0", "Z1"]
+        qubitwise: Determines if the check is for general commutativity, or the stricter qubitwise commutativity
+
+    Returns:
+        list: Containing groups (lists) of Pauli strings that all commute mutually
+    """
+    # Need to convert the lists in _terms_list to strings so that they can be used as graph nodes
+    terms_str = ["".join(term) for term in terms_list]
+    # print(terms_str)
+
+    G = nx.Graph()
+    # Complement graph: Add all the terms as nodes first, then add edges between nodes if they DO NOT commute
+    G.add_nodes_from(terms_str)
+    G.add_edges_from(
+        (term1, term2)
+        for _i1, term1 in enumerate(terms_str)
+        for _i2, term2 in enumerate(terms_str)
+        if _i2 > _i1 and not check_terms_commutativity(term1, term2, qubitwise)
+    )
+
+    # Solve using Greedy Colouring on NetworkX
+    term_groups = nx.coloring.greedy_color(G)
+    group_ids = set(term_groups.values())
+    sorted_groups = [[group.split() for group, group_id in term_groups.items() if group_id == _id] for _id in group_ids]
+    print(sorted_groups)
+    return sorted_groups
+
+
 def measurement_basis_rotations(hamiltonian, n_qubits, grouping=None):
     """
     Split up and sort the Hamiltonian terms to get the basis rotation gates to be applied to a quantum circuit for the
@@ -52,17 +96,16 @@ def measurement_basis_rotations(hamiltonian, n_qubits, grouping=None):
         empty lists - ``([], [])`` - if there are no Z terms present.
     """
     result = []
-    # Split up the Z and X/Y terms first
+    # Split up the Z and X/Y terms
     z_only_terms = [
         term for term in hamiltonian.terms if not any(factor.name[0] in ("X", "Y") for factor in term.factors)
     ]
     xy_terms = [term for term in hamiltonian.terms if term not in z_only_terms]
-    # Add the Z terms into result first
+    # Add the Z terms into result first, followed by the terms with X/Y's
     if z_only_terms:
         result.append(([gates.M(_i) for _i in range(n_qubits)], z_only_terms))
     else:
-        result.append(([], []))  # No terms with only Z's
-    # Then add the X/Y terms in
+        result.append(([], []))
     if xy_terms:
         if grouping is None:
             result += [
diff --git a/tests/test_measurement_optimisation.py b/tests/test_measurement_optimisation.py
index e3c217f..e96b650 100644
--- a/tests/test_measurement_optimisation.py
+++ b/tests/test_measurement_optimisation.py
@@ -6,7 +6,10 @@
 
 # from qibochem.driver import Molecule
 # from qibochem.measurement import expectation
-from qibochem.measurement.optimization import check_terms_commutativity
+from qibochem.measurement.optimization import (
+    check_terms_commutativity,
+    group_commuting_terms,
+)
 
 # from qibo import Circuit, gates
 # from qibo.hamiltonians import SymbolicHamiltonian
@@ -16,15 +19,31 @@
 @pytest.mark.parametrize(
     "term1,term2,qwc_expected,gc_expected",
     [
-        (["X0"], ["Z0"], False, False),
-        (["X0"], ["Z1"], True, True),
-        (["X0", "X1"], ["Y0", "Y1"], False, True),
-        (["X0", "Y1"], ["Y0", "Y1"], False, False),
+        ("X0", "Z0", False, False),
+        ("X0", "Z1", True, True),
+        ("X0 X1", "Y0 Y1", False, True),
+        ("X0 Y1", "Y0 Y1", False, False),
     ],
 )
 def test_check_terms_commutativity(term1, term2, qwc_expected, gc_expected):
     """Do two Pauli strings commute (qubit-wise or generally)?"""
     qwc_result = check_terms_commutativity(term1, term2, qubitwise=True)
     assert qwc_result == qwc_expected
-    gc_result = check_terms_commutativity(term1, term2)
+    gc_result = check_terms_commutativity(term1, term2, qubitwise=False)
     assert gc_result == gc_expected
+
+
+# @pytest.mark.parametrize(
+#     "term_list,qwc_expected,gc_expected",
+#     [
+#         ([["X0"], ["Z0"], ["Z1"], ["Z0", "Z1"]], {[["Z0"], ["Z1"], ["Z0", "Z1"]], ["X0"]}, {[["Z0"], ["Z1"], ["Z0", "Z1"]], ["X0"]}),
+#         # (["X0"], ["Z1"], True, True),
+#         # (["X0", "X1"], ["Y0", "Y1"], False, True),
+#         # (["X0", "Y1"], ["Y0", "Y1"], False, False),
+#     ],
+# )
+# def test_group_commuting_terms(term_list, qwc_expected, gc_expected):
+#     qwc_result = group_commuting_terms(term_list, qubitwise=True)
+#     assert set(qwc_result) == qwc_expected
+#     gc_result = group_commuting_terms(term_list, qubitwise=False)
+#     assert set(gc_result) == gc_expected