diff --git a/notebooks/tfim_2cnot_test.py b/notebooks/tfim_2cnot_test.py
new file mode 100644
index 0000000..7f13cc4
--- /dev/null
+++ b/notebooks/tfim_2cnot_test.py
@@ -0,0 +1,62 @@
+from qibo.hamiltonians import SymbolicHamiltonian
+from boostvqe.models.dbi.double_bracket_evolution_oracles import *
+from functools import reduce
+import numpy as np
+from qibo import hamiltonians
+import matplotlib.pyplot as plt
+
+n_qubits = 3
+h_coeff = 1
+hamiltonian = SymbolicHamiltonian(nqubits=n_qubits)
+
+oracle = TFIM_EvolutionOracle(h=hamiltonian, evolution_oracle_type="trotter", steps=1, B_a=0, order=2)
+
+circuit = oracle.circuit(t_duration=1.0)
+
+unitary = circuit.unitary()
+def multikron(matrix_list):
+    """Calculates Kronecker product of a list of matrices.
+
+    Args:
+        matrix_list (list): List of matrices as ``ndarray``.
+
+    Returns:
+        ndarray: Kronecker product of all matrices in ``matrix_list``.
+    """
+    return reduce(np.kron, matrix_list)
+
+from numpy.linalg import norm
+
+def our_TFIM(nqubits, h: float = 0.0, dense: bool = True, backend=None):
+    def multikron(matrix_list):
+        """Calculates Kronecker product of a list of matrices."""
+        return reduce(np.kron, matrix_list)
+
+    from qibo.backends import matrices
+
+    matrix = (
+        - multikron([matrices.X, matrices.X]) - h * multikron([matrices.Z, matrices.I])
+    )
+    terms = [hamiltonians.terms.HamiltonianTerm(matrix, i, i + 1) for i in range(nqubits - 1)]
+    terms.append(hamiltonians.terms.HamiltonianTerm(matrix, nqubits - 1, 0))
+    ham = SymbolicHamiltonian(backend=backend)
+    ham.terms = terms
+    return ham
+
+ham = our_TFIM(nqubits=n_qubits, h=h_coeff, dense=False)
+truth = ham.exp(1)
+verification_norm = []
+for step in range(1, 21):
+    oracle = TFIM_EvolutionOracle(h=hamiltonian, evolution_oracle_type="trotter", steps=step, B_a=h_coeff, order=2)
+    circuit = oracle.circuit(t_duration=1.0)
+    unitary = circuit.unitary()
+    verification_norm.append(norm(truth-unitary))
+
+
+x = np.array([i for i in range(1, 21)])
+plt.plot(x, verification_norm, 'o')
+plt.title("verification of TFIM 2 CNOT implementation")
+plt.xlabel("steps")
+plt.ylabel("norm of difference")
+
+plt.show()
\ No newline at end of file
diff --git a/src/boostvqe/models/dbi/double_bracket_evolution_oracles.py b/src/boostvqe/models/dbi/double_bracket_evolution_oracles.py
index 7740df6..0946b51 100644
--- a/src/boostvqe/models/dbi/double_bracket_evolution_oracles.py
+++ b/src/boostvqe/models/dbi/double_bracket_evolution_oracles.py
@@ -4,8 +4,9 @@
 from functools import cached_property, reduce
 from typing import Union
 
-import hyperopt
-import matplotlib.pyplot as plt
+
+#import hyperopt
+#import matplotlib.pyplot as plt
 import numpy as np
 from qibo import Circuit, gates, symbols
 from qibo.config import raise_error
@@ -307,3 +308,51 @@ def circuit(self, t_duration, steps=None, order=None):
             steps=steps,
             order=order,
         )
+
+
+@dataclass
+class TFIM_EvolutionOracle(EvolutionOracle):
+    steps: int = None
+    B_a: float = None
+    order: int = None
+
+    def circuit(self, t_duration):
+        circuit_v = Circuit(self.h.nqubits)  # Initialize the circuit with the number of qubits
+        t_duration /= self.steps
+        def routine(tmp_circuit, enum_list, routine_t):
+            for a in enum_list:
+                tmp_circuit.add(gates.CNOT(a, (a + 1) % self.h.nqubits))
+                # Time evolution under the transverse field Ising model Hamiltonian
+                # exp(-i t (X(a) + B_a * Z(a)))
+                self._time_evolution_step(tmp_circuit, a, routine_t)
+
+                # Add second CNOT(a, a+1)
+                tmp_circuit.add(gates.CNOT(a, (a + 1) % self.h.nqubits))
+        if self.order is None:
+            for _ in range(self.steps):
+                routine(circuit_v, range(self.h.nqubits), t_duration)
+        elif self.order == 1:
+            for _ in range(self.steps):
+                routine(circuit_v, range(0, self.h.nqubits, 2), t_duration)
+                routine(circuit_v, range(1, self.h.nqubits, 2), t_duration)
+
+        elif self.order == 2:
+            for _ in range(self.steps):
+                routine(circuit_v, range(1, self.h.nqubits, 2), t_duration/2)
+                routine(circuit_v, range(0, self.h.nqubits, 2), t_duration)
+                routine(circuit_v, range(1, self.h.nqubits, 2), t_duration / 2)
+        else:
+            print("order must be either 1 or 2")
+        return circuit_v
+
+    def _time_evolution_step(self, tmp_circuit: Circuit, a: int, dt: float):
+        """Apply a single Trotter step of the time evolution operator exp(-i dt (X(a) + B_a Z(a)))."""
+
+        # Time evolution for X(a)
+        tmp_circuit.add(gates.RX(a, theta=-2*dt))  # Apply exp(-i dt X(a))
+
+        # Time evolution for Z(a)
+        tmp_circuit.add(gates.RZ(a, theta=-2*dt * self.B_a))  # Apply exp(-i dt B_a Z(a))
+
+        return tmp_circuit
+