diff --git a/src/tequila/grouping/binary_rep.py b/src/tequila/grouping/binary_rep.py
index 6aea305d..ac15efd0 100644
--- a/src/tequila/grouping/binary_rep.py
+++ b/src/tequila/grouping/binary_rep.py
@@ -26,7 +26,7 @@ def init_from_qubit_hamiltonian(cls, hamiltonian: QubitHamiltonian, n_qubits=Non
                 del Hof.terms[()]
                 hamiltonian = QubitHamiltonian.from_openfermion(Hof)
         if n_qubits is None:
-            n_qubits = hamiltonian.n_qubits
+            n_qubits = max(hamiltonian.qubits)+1
         binary_terms = [
             BinaryPauliString(
                 p.binary(n_qubits).binary,
diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 560c6b91..5f174bb4 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -1487,7 +1487,8 @@ def make_uccsd_ansatz(self, trotter_steps: int = 1,
         factor = 1.0 / trotter_steps
         for step in range(trotter_steps):
             for idx, angle in indices.items():
-                UCCSD += self.make_excitation_gate(indices=idx, angle=factor * angle)
+                converted = [(idx[2 * i], idx[2 * i + 1]) for i in range(len(idx) // 2)]
+                UCCSD += self.make_excitation_gate(indices=converted, angle=factor * angle)
         if hasattr(initial_amplitudes,
                    "lower") and initial_amplitudes.lower() == "mp2" and parametrized and add_singles:
             # mp2 has no singles, need to initialize them here (if not parametrized initializling as 0.0 makes no sense though)
@@ -1722,7 +1723,7 @@ def rdm2(self):
 
     def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_free: bool = True,
                      get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac", use_hcb: bool = False,
-                     rdm_trafo: QubitHamiltonian = None):
+                     rdm_trafo: QubitHamiltonian = None, decompose=None):
         """
         Computes the one- and two-particle reduced density matrices (rdm1 and rdm2) given
         a unitary U. This method uses the standard ordering in physics as denoted below.
@@ -2000,7 +2001,12 @@ def _build_2bdy_operators_hcb() -> list:
         
         # Compute expected values
         if rdm_trafo is None:
-            evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
+            if decompose is not None:
+                print("MANIPULATED")
+                X = decompose(H=qops, U=U)
+                evals = simulate(X, variables=variables)
+            else:
+                evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)
         else:
             qops = [rdm_trafo.dagger()*qops[i]*rdm_trafo for i in range(len(qops))]
             evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables)