diff --git a/src/qrisp/operators/fermionic/fermionic_operator.py b/src/qrisp/operators/fermionic/fermionic_operator.py
index 44f5b700..cc5e7458 100644
--- a/src/qrisp/operators/fermionic/fermionic_operator.py
+++ b/src/qrisp/operators/fermionic/fermionic_operator.py
@@ -708,14 +708,19 @@ def get_measurement(
     #
     
 
-    def trotterization(self, t = 1, steps = 1, iter = 1):
+    def trotterization(self, forward_evolution = True):
         r"""
-        Returns a function for performing Hamiltonian simulation, i.e., approximately implementing the unitary operator $e^{itH}$ via Trotterization.
+        Returns a function for performing Hamiltonian simulation, i.e., approximately implementing the unitary operator $U(t) = e^{-itH}$ via Trotterization.
         Note that this method will always simulate the **hermitized** operator, i.e.
         
         .. math::
             
             H = (O + O^\dagger)/2
+            
+        Parameters
+        ----------
+        forward_evolution, bool, optional
+            If set to False $U(t)^\dagger = e^{itH}$ will be executed (usefull for quantum phase estimation). The default is True.
 
         Returns
         -------
@@ -801,7 +806,7 @@ def trotter_step(qarg, t, steps):
                         coeff = reduced_H.terms_dict[ferm_term]
                         pauli_hamiltonian = ferm_term.fermionic_swap(permutation).to_qubit_term()
                         pauli_term = list(pauli_hamiltonian.terms_dict.keys())[0]
-                        pauli_term.simulate(coeff*t/steps*pauli_hamiltonian.terms_dict[pauli_term], new_qarg)
+                        pauli_term.simulate(-coeff*t/steps*pauli_hamiltonian.terms_dict[pauli_term]*(-1)**int(forward_evolution), new_qarg)
                 
 
         def U(qarg, t=1, steps=1, iter=1):
@@ -812,8 +817,9 @@ def U(qarg, t=1, steps=1, iter=1):
         return U
     
     def group_up(self, denominator):
-        
         term_groups = group_up_terms(self, denominator)
+        if len(term_groups) == 0:
+            return [self]
         groups = []
         for term_group in term_groups:
             O = FermionicOperator({term : self.terms_dict[term] for term in term_group})
diff --git a/src/qrisp/operators/hamiltonian_tools.py b/src/qrisp/operators/hamiltonian_tools.py
index 983faa39..202b7fd8 100644
--- a/src/qrisp/operators/hamiltonian_tools.py
+++ b/src/qrisp/operators/hamiltonian_tools.py
@@ -178,6 +178,9 @@ def dsatur_coloring(num_vertices, adjacency_matrix):
     return colors
 
 def find_coloring(G):
+    if len(G) == 0:
+        return []
+    
     adjacency_matrix = nx.to_numpy_array(G)
     
     coloring_1 = rlf_coloring(len(G), adjacency_matrix)
@@ -202,6 +205,9 @@ def group_up_terms(H, group_denominator):
     
     coloring = find_coloring(G)
     
+    if len(coloring) == 0:
+        return []
+    
     groups = []
     for i in range(np.max(coloring)+1): groups.append([])
     
diff --git a/src/qrisp/operators/qubit/bound_qubit_operator.py b/src/qrisp/operators/qubit/bound_qubit_operator.py
index 95a82d16..b9d51ba6 100644
--- a/src/qrisp/operators/qubit/bound_qubit_operator.py
+++ b/src/qrisp/operators/qubit/bound_qubit_operator.py
@@ -646,7 +646,6 @@ def get_measurement(
         return unbound_operator.get_measurement(qarg,
                                                 precision=precision, 
                                                 backend=backend, 
-                                                shots=shots, 
                                                 compile=compile, 
                                                 compilation_kwargs=compilation_kwargs, 
                                                 subs_dic=subs_dic,
diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py
index 0ecfa99e..e8ed523c 100644
--- a/src/qrisp/operators/qubit/qubit_operator.py
+++ b/src/qrisp/operators/qubit/qubit_operator.py
@@ -778,6 +778,8 @@ def commuting_groups(self):
     
     def group_up(self, group_denominator):
         term_groups = group_up_terms(self, group_denominator)
+        if len(term_groups) == 0:
+            return [self]
         groups = []
         for term_group in term_groups:
             H = QubitOperator({term : self.terms_dict[term] for term in term_group})
@@ -1423,11 +1425,11 @@ def get_measurement(
     # Trotterization
     #
 
-    def trotterization(self, method='commuting_qw'):
+    def trotterization(self, method='commuting_qw', forward_evolution = True):
         r"""
         .. _ham_sim:
         
-        Returns a function for performing Hamiltonian simulation, i.e., approximately implementing the unitary operator $e^{itH}$ via Trotterization.
+        Returns a function for performing Hamiltonian simulation, i.e., approximately implementing the unitary operator $U(t) = e^{-itH}$ via Trotterization.
         Note that this method will always simulate the **hermitized** operator, i.e.
         
         .. math::
@@ -1441,6 +1443,8 @@ def trotterization(self, method='commuting_qw'):
             The method for grouping the QubitTerms. 
             Available are ``commuting`` (groups such that all QubitTerms mutually commute) and ``commuting_qw`` (groups such that all QubitTerms mutually commute qubit-wise).
             The default is ``commuting_qw``.
+        forward_evolution : bool, optional
+            If set to False $U(t)^\dagger = e^{itH}$ will be executed (usefull for quantum phase estimation). The default is True.
 
         Returns
         -------
@@ -1512,7 +1516,7 @@ def trotter_step(qarg, t, steps):
                             intersect_groups = diagonal_operator.group_up(lambda a, b: not a.intersect(b))
                             for intersect_group in intersect_groups:
                                 for term,coeff in intersect_group.terms_dict.items():
-                                    term.simulate(coeff*t/steps, qarg)
+                                    term.simulate(-coeff*t/steps*(-1)**int(forward_evolution), qarg)
         
         if method=='commuting':
             def trotter_step(qarg, t, steps):
@@ -1521,7 +1525,7 @@ def trotter_step(qarg, t, steps):
                         intersect_groups = diagonal_operator.group_up(lambda a, b: not a.intersect(b))
                         for intersect_group in intersect_groups:
                             for term,coeff in intersect_group.terms_dict.items():
-                                term.simulate(coeff*t/steps, qarg)
+                                term.simulate(-coeff*t/steps*(-1)**int(forward_evolution), qarg)
 
         def U(qarg, t=1, steps=1, iter=1):
             merge([qarg])
diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py
index fc8901f3..0dae6576 100644
--- a/src/qrisp/operators/qubit/qubit_term.py
+++ b/src/qrisp/operators/qubit/qubit_term.py
@@ -243,7 +243,7 @@ def simulate(self, coeff, qv, ctrl = None):
         # If no non-trivial indices are found, we perform a global phase
         # and are done.
         if len(Z_indices + projector_qubits) == 0:
-            gphase(coeff, qv[0])
+            gphase(-coeff, qv[0])
             return
         
         # If there are only projectors, the circuit is a mcp gate
@@ -260,14 +260,13 @@ def simulate(self, coeff, qv, ctrl = None):
             # Perform the mcp            
             if len(projector_qubits) == 1:
                 if len(flip_qubits) == 1:
-                    p(-coeff, projector_qubits[0])
-                    gphase(coeff, projector_qubits[0])
-                else:
                     p(coeff, projector_qubits[0])
+                    gphase(-coeff, projector_qubits[0])
+                else:
+                    p(-coeff, projector_qubits[0])
             else:
                 with env:
-                    
-                        mcp(coeff, projector_qubits, method = "balauca")
+                    mcp(-coeff, projector_qubits, method = "balauca")
                     
             return
         
@@ -400,12 +399,12 @@ def flip_anchor_qubit(qv, anchor_index, Z_indices):
                 if control_qubit_available:
                     # Use Selinger's circuit (page 5)
                     with conjugate(cx)(qv[anchor_index], hs_anc):
-                        rz(-coeff, qv[anchor_index])
+                        rz(coeff, qv[anchor_index])
                         if flip_control_phase:
                             coeff = -coeff
-                        rz(coeff, hs_anc)
+                        rz(-coeff, hs_anc)
                 else:
-                    rz(-coeff*2, qv[anchor_index])
+                    rz(coeff*2, qv[anchor_index])
                 
     
         if len(projector_indices) >= 2:
diff --git a/tests/operator_tests/test_hamiltonian.py b/tests/operator_tests/test_hamiltonian.py
index c35aaa27..925686b2 100644
--- a/tests/operator_tests/test_hamiltonian.py
+++ b/tests/operator_tests/test_hamiltonian.py
@@ -63,7 +63,7 @@ def test_trotterization():
     E0 = G.ground_state_energy()
     assert np.abs(E0-(-0.804899065613056)) < 2e-2
 
-    U = G.trotterization()
+    U = G.trotterization(forward_evolution = False)
 
     qv = QuantumVariable(2)
     x(qv) 
diff --git a/tests/operator_tests/test_qubit_hamiltonian_simulation.py b/tests/operator_tests/test_qubit_hamiltonian_simulation.py
index ddc6a298..e4371d9c 100644
--- a/tests/operator_tests/test_qubit_hamiltonian_simulation.py
+++ b/tests/operator_tests/test_qubit_hamiltonian_simulation.py
@@ -17,7 +17,7 @@
 """
 
 from qrisp import QuantumVariable, x, QPE
-from qrisp.operators import X, Y, Z, A, C, P0, P1
+from qrisp.operators import X, Y, Z, A, C, P0, P1, QubitOperator
 import numpy as np
 
 def test_qubit_hamiltonian_simulation():
@@ -29,7 +29,7 @@ def test_qubit_hamiltonian_simulation():
 
     for method in ['commuting_qw', 'commuting']:
 
-        U = H.trotterization(method=method)
+        U = H.trotterization(method=method, forward_evolution = False)
 
         # Find minimum eigenvalue of H with Hamiltonian simulation and QPE
 
@@ -38,7 +38,7 @@ def test_qubit_hamiltonian_simulation():
         x(qv[0])
         E1 = H.get_measurement(qv)
         assert abs(E1-E0)<5e-2
-
+        
         qpe_res = QPE(qv,U,precision=10,kwargs={"steps":3},iter_spec=True)
 
         results = qpe_res.get_measurement()    
@@ -53,10 +53,10 @@ def test_qubit_hamiltonian_simulation():
     def verify_trotterization(H,method):
         
         # Compute hermitean matrix
-        H_matrix = H.to_sparse_matrix().todense()
+        H_matrix = H.to_sparse_matrix(4).todense()
         H_matrix = (H_matrix + H_matrix.transpose().conjugate())/2
         # Compute unitary matrix
-        U_matrix = expm(1j*H_matrix)
+        U_matrix = expm(-1j*H_matrix)
         
         # Perform trotterization
         qv = QuantumVariable(int(np.log2(H_matrix.shape[0])))
@@ -96,7 +96,7 @@ def verify_trotterization(H,method):
                 for O3 in operator_list:
                     H = O0(0)*O1(1)*O2(2)*O3(3)
                     if H is 1:
-                        continue
+                        H = QubitOperator() + 1
                     print(H)
                     verify_trotterization(H,'commuting_qw')
                     verify_trotterization(H,'commuting')