diff --git a/src/qibo/abstractions/abstract_gates.py b/src/qibo/abstractions/abstract_gates.py
index 3bbd349ba1..def87f7379 100644
--- a/src/qibo/abstractions/abstract_gates.py
+++ b/src/qibo/abstractions/abstract_gates.py
@@ -432,9 +432,7 @@ def parameters(self, x):
         # ``circuit.set_parameters`` method works properly.
         # pylint: disable=E1101
         if isinstance(self, BaseBackendGate):
-            self._matrix = None
-            self._native_op_matrix = None
-            self._custom_op_matrix = None
+            self._reset_unitary()
             for devgate in self.device_gates:
                 devgate.parameters = x
 
@@ -519,6 +517,13 @@ def _construct_unitary(self): # pragma: no cover
         """Constructs the gate's unitary matrix."""
         return raise_error(NotImplementedError)
 
+    def _reset_unitary(self):
+        """Resets the gate matrices back to ``None``.
+
+        Useful when the gate matrix need to be recalculated.
+        """
+        self._matrix = None
+
     @property
     @abstractmethod
     def cache(self): # pragma: no cover
diff --git a/src/qibo/abstractions/circuit.py b/src/qibo/abstractions/circuit.py
index 29a4588fa4..94e035b3e3 100644
--- a/src/qibo/abstractions/circuit.py
+++ b/src/qibo/abstractions/circuit.py
@@ -18,13 +18,13 @@ class _ParametrizedGates(list):
     total number of parameters.
     """
 
-    def __init__(self):
-        super(_ParametrizedGates, self).__init__(self)
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
         self.set = set()
         self.nparams = 0
 
-    def append(self, gate: gates.ParametrizedGate):
-        super(_ParametrizedGates, self).append(gate)
+    def append(self, gate):
+        super().append(gate)
         self.set.add(gate)
         self.nparams += gate.nparams
 
@@ -94,8 +94,6 @@ def __init__(self, nqubits):
         self.measurement_gate = None
         self.measurement_gate_result = None
 
-        self.fusion_groups = []
-
         self._final_state = None
         self.density_matrix = False
         self.repeated_execution = False
@@ -180,6 +178,16 @@ def on_qubits(self, *q):
         for gate in self.queue:
             yield gate._on_qubits(*q)
 
+    def _shallow_copy(self):
+        """Helper method for :meth:`qibo.abstractions.circuit.AbstractCircuit.copy`
+        and :meth:`qibo.core.circuit.Circuit.fuse`."""
+        new_circuit = self.__class__(**self.init_kwargs)
+        new_circuit.parametrized_gates = _ParametrizedGates(self.parametrized_gates)
+        new_circuit.trainable_gates = _ParametrizedGates(self.trainable_gates)
+        new_circuit.measurement_gate = self.measurement_gate
+        new_circuit.measurement_tuples = dict(self.measurement_tuples)
+        return new_circuit
+
     def copy(self, deep: bool = False):
         """Creates a copy of the current ``circuit`` as a new ``Circuit`` model.
 
@@ -192,9 +200,14 @@ def copy(self, deep: bool = False):
             The copied circuit object.
         """
         import copy
-        new_circuit = self.__class__(**self.init_kwargs)
         if deep:
+            new_circuit = self.__class__(**self.init_kwargs)
             for gate in self.queue:
+                if isinstance(gate, gates.FusedGate): # pragma: no cover
+                    # impractical case
+                    raise_error(NotImplementedError, "Cannot create deep copy "
+                                                     "of fused circuit.")
+
                 new_gate = copy.copy(gate)
                 new_circuit.queue.append(new_gate)
                 if isinstance(gate, gates.ParametrizedGate):
@@ -202,17 +215,10 @@ def copy(self, deep: bool = False):
                     if gate.trainable:
                         new_circuit.trainable_gates.append(new_gate)
             new_circuit.measurement_gate = copy.copy(self.measurement_gate)
-            if self.fusion_groups: # pragma: no cover
-                # impractical case
-                raise_error(NotImplementedError, "Cannot create deep copy of fused "
-                                                 "circuit.")
+            new_circuit.measurement_tuples = dict(self.measurement_tuples)
         else:
+            new_circuit = self._shallow_copy()
             new_circuit.queue = copy.copy(self.queue)
-            new_circuit.parametrized_gates = list(self.parametrized_gates)
-            new_circuit.trainable_gates = list(self.trainable_gates)
-            new_circuit.measurement_gate = self.measurement_gate
-            new_circuit.fusion_groups = list(self.fusion_groups)
-        new_circuit.measurement_tuples = dict(self.measurement_tuples)
         return new_circuit
 
     def invert(self):
@@ -517,9 +523,6 @@ def _set_parameters_list(self, parameters: List, n: int):
                                     "the circuit contains {} parametrized gates."
                                     "".format(n, len(self.trainable_gates)))
 
-        for fusion_group in self.fusion_groups:
-            fusion_group.update()
-
     def set_parameters(self, parameters):
         """Updates the parameters of the circuit's parametrized gates.
 
@@ -569,9 +572,6 @@ def set_parameters(self, parameters):
         elif isinstance(parameters, self.param_tensor_types):
             self._set_parameters_list(parameters, int(parameters.shape[0]))
         elif isinstance(parameters, dict):
-            if self.fusion_groups:
-                raise_error(TypeError, "Cannot accept new parameters as dictionary "
-                                       "for fused circuits. Use list, tuple or array.")
             diff = set(parameters.keys()) - self.trainable_gates.set
             if diff:
                 raise_error(KeyError, "Dictionary contains gates {} which are "
@@ -582,6 +582,11 @@ def set_parameters(self, parameters):
         else:
             raise_error(TypeError, "Invalid type of parameters {}."
                                    "".format(type(parameters)))
+        # Reset ``FusedGate`` matrices so that they are recalculated with the
+        # updated parameters.
+        for gate in self.queue:
+            if isinstance(gate, gates.FusedGate):
+                gate._reset_unitary()
 
     def get_parameters(self, format: str = "list",
                        include_not_trainable: bool = False
diff --git a/src/qibo/abstractions/gates.py b/src/qibo/abstractions/gates.py
index 2fdc011bb7..1e58f9bbb0 100644
--- a/src/qibo/abstractions/gates.py
+++ b/src/qibo/abstractions/gates.py
@@ -1550,3 +1550,41 @@ def __init__(self, q, t1, t2, time, excited_population=0, seed=None):
             seed=seed)
         # this case can only be applied to density matrices
         self.density_matrix = True
+
+
+class FusedGate(Gate):
+    """Collection of gates that will be fused and applied as single gate during simulation.
+
+    This gate is constructed automatically by :meth:`qibo.core.circuits.Circuit.fuse`
+    and should not be used by user.
+    :class:`qibo.abstractions.gates.FusedGate` works with arbitrary number of
+    target qubits however the backend implementation
+    :class:`qibo.core.gates.FusedGate` assumes two target qubits.
+    """
+
+    def __init__(self, *q):
+        super().__init__()
+        self.name = "fused"
+        self.target_qubits = tuple(q)
+        self.init_args = list(q)
+        self.qubit_set = set(q)
+        self.gates = []
+
+    def add(self, gate):
+        if not set(gate.qubits).issubset(self.qubit_set):
+            raise_error(ValueError, "Cannot add gate that targets {} "
+                                    "in fused gate acting on {}."
+                                    "".format(gate.qubits, self.qubits))
+        if isinstance(gate, self.__class__):
+            self.gates.extend(gate.gates)
+        else:
+            self.gates.append(gate)
+
+    def __iter__(self):
+        return iter(self.gates)
+
+    def _dagger(self):
+        dagger = self.__class__(*self.init_args)
+        for gate in self.gates[::-1]:
+            dagger.add(gate.dagger())
+        return dagger
diff --git a/src/qibo/core/circuit.py b/src/qibo/core/circuit.py
index 0df42da7c8..b6f16699dd 100644
--- a/src/qibo/core/circuit.py
+++ b/src/qibo/core/circuit.py
@@ -23,7 +23,6 @@ class Circuit(circuit.AbstractCircuit):
     Args:
         nqubits (int): Total number of qubits in the circuit.
     """
-    from qibo.core import fusion
 
     def __init__(self, nqubits):
         super(Circuit, self).__init__(nqubits)
@@ -47,37 +46,13 @@ def _add_layer(self, gate):
             self._set_nqubits(gate.additional_unitary)
             self.queue.append(gate.additional_unitary)
 
-    def _fuse_copy(self):
-        """Helper method for ``circuit.fuse``.
-
-        For standard (non-distributed) circuits this creates a copy of the
-        circuit with deep-copying the parametrized gates only.
-        For distributed circuits a fully deep copy should be created.
-        """
-        import copy
-        from qibo.abstractions.abstract_gates import ParametrizedGate
-        new_circuit = self.__class__(**self.init_kwargs)
-        for gate in self.queue:
-            if isinstance(gate, ParametrizedGate):
-                if gate.trainable:
-                    new_gate = copy.copy(gate)
-                    new_circuit.queue.append(new_gate)
-                    new_circuit.parametrized_gates.append(new_gate)
-                    new_circuit.trainable_gates.append(new_gate)
-                else:
-                    new_circuit.queue.append(gate)
-                    new_circuit.parametrized_gates.append(gate)
-            else:
-                new_circuit.queue.append(gate)
-        new_circuit.measurement_gate = copy.copy(self.measurement_gate)
-        new_circuit.measurement_tuples = dict(self.measurement_tuples)
-        return new_circuit
-
     def fuse(self):
-        """Creates an equivalent ``Circuit`` with gates fused up to two-qubits.
+        """Creates an equivalent :class:`qibo.core.circuit.Circuit` with gates fused up to two-qubits.
+
+        The current fusion algorithm with create up to two-qubit fused gates.
 
         Returns:
-            The equivalent ``Circuit`` object where the gates are fused.
+            The equivalent :class:`qibo.core.circuit.Circuit` object where the gates are fused.
 
         Example:
             ::
@@ -89,18 +64,120 @@ def fuse(self):
                 c.add([gates.Y(0), gates.Y(1)])
                 # create circuit with fused gates
                 fused_c = c.fuse()
-                # now ``fused_c`` contains only one ``gates.Unitary`` gate
-                # that is equivalent to applying the five gates of the original
-                # circuit.
+                # now ``fused_c`` contains a single gate that is equivalent
+                # to applying the five gates of the original circuit.
         """
+        from qibo import gates
         from qibo.abstractions.circuit import _Queue
-        new_circuit = self._fuse_copy()
-        new_circuit.fusion_groups = self.fusion.FusionGroup.from_queue(
-            new_circuit.queue)
-        new_circuit.queue = _Queue(self.nqubits)
-        for group in new_circuit.fusion_groups:
-            for gate in group.gates:
-                new_circuit.queue.append(gate)
+        from qibo.abstractions.abstract_gates import SpecialGate
+
+        class FusedQueue(_Queue):
+            """Helper queue implementation that checks if a gate already exists
+            in queue to avoid re-appending it.
+            """
+
+            def __init__(self, *args, **kwargs):
+                super().__init__(*args, **kwargs)
+                self.set = set()
+
+            def append(self, gate):
+                """Appends a gate in queue only if it is not already in."""
+                # Use a ``set`` instead of the original ``list`` to check if
+                # the gate already exists in queue as lookup is typically
+                # more efficient for sets
+                # (although actual performance difference is probably negligible)
+                if gate not in self.set:
+                    self.set.add(gate)
+                    super().append(gate)
+
+        # new circuit queue that will hold the fused gates
+        fused_queue = FusedQueue(self.nqubits)
+        # dictionary that maps each qubit id (int) to the corresponding
+        # active ``FusedGate`` that is part of
+        fused_gates = collections.OrderedDict()
+        # use ``OrderedDict`` so that the original gate order is not changed
+        for gate in self.queue:
+            qubits = gate.qubits
+            if len(qubits) == 1:
+                # add one-qubit gates to the active ``FusedGate`` of this qubit
+                # or create a new one if it does not exist
+                q = qubits[0]
+                if q not in fused_gates:
+                    fused_gates[q] = gates.FusedGate(q)
+                fused_gates.get(q).add(gate)
+
+            elif len(qubits) == 2:
+                # fuse two-qubit gates
+                q0, q1 = tuple(sorted(qubits))
+                if (q0 in fused_gates and q1 in fused_gates and
+                    fused_gates.get(q0) == fused_gates.get(q1)):
+                    # if the target qubit pair is compatible with the active
+                    # ``FusedGate`` of both qubits then add it to the ``FusedGate``
+                    fused_gates.get(q0).add(gate)
+                else:
+                    # otherwise we need to create a new ``FusedGate`` and
+                    # update the active gates of both target qubits
+                    fgate = gates.FusedGate(q0, q1)
+                    if q0 in fused_gates:
+                        # first qubit has existing active gate
+                        ogate = fused_gates.pop(q0)
+                        if len(ogate.target_qubits) == 1:
+                            # existing active gate is one-qubit so we just add
+                            # it to the new ``FusedGate``
+                            fgate.add(ogate)
+                        else:
+                            # existing active gate is two-qubit so we need to
+                            # add it to the new queue
+                            fused_queue.append(ogate)
+                    if q1 in fused_gates:
+                        # second qubit has existing active gate
+                        ogate = fused_gates.pop(q1)
+                        if len(ogate.target_qubits) == 1:
+                            # existing active gate is one-qubit so we just add
+                            # it to the new ``FusedGate``
+                            fgate.add(ogate)
+                        else:
+                            # existing active gate is two-qubit so we need to
+                            # add it to the new queue
+                            fused_queue.append(ogate)
+                    # add the two-qubit gate to the newly created ``FusedGate``
+                    # and update the active ``FusedGate``s of both target qubits
+                    fgate.add(gate)
+                    fused_gates[q0], fused_gates[q1] = fgate, fgate
+
+            elif isinstance(gate, SpecialGate):
+                # ``SpecialGate``s act on all qubits (like a barrier) so we
+                # so we need to temporarily stop the fusion, add all active
+                # gates in the new queue and restart fusion after the barrier
+                for g in fused_gates.values():
+                    fused_queue.append(g)
+                fused_gates = collections.OrderedDict()
+                fused_queue.append(gate)
+
+            else:
+                # gate has more than two target qubits so it cannot be included
+                # in the ``FusedGate``s which support up to two qubits.
+                # Therefore we deactivate the ``FusedGate``s of all target qubits
+                for q in qubits:
+                    if q in fused_gates:
+                        fused_queue.append(fused_gates.pop(q))
+                fused_queue.append(gate)
+
+        for gate in fused_gates.values():
+            # add remaining active ``FusedGate``s in the new queue
+            fused_queue.append(gate)
+
+        queue = _Queue(self.nqubits)
+        for gate in fused_queue:
+            if isinstance(gate, gates.FusedGate) and len(gate.gates) == 1:
+                # replace ``FusedGate``s that contain only one gate by this
+                # gate for efficiency
+                gate = gate.gates[0]
+            queue.append(gate)
+
+        # create a circuit and assign the new queue
+        new_circuit = self._shallow_copy()
+        new_circuit.queue = queue
         return new_circuit
 
     def _eager_execute(self, state):
diff --git a/src/qibo/core/distcircuit.py b/src/qibo/core/distcircuit.py
index f9d4d8f4d3..84020054e5 100644
--- a/src/qibo/core/distcircuit.py
+++ b/src/qibo/core/distcircuit.py
@@ -76,14 +76,9 @@ def copy(self, deep: bool = True):
                                     "circuits because they modify gate objects.")
         return super().copy(deep)
 
-    def _fuse_copy(self):
-        return self.copy(deep=True)
-
     def fuse(self):
-        if self.queues.queues:
-            raise_error(RuntimeError, "Cannot fuse distributed circuit after "
-                                      "its first execution.")
-        return super().fuse()
+        raise_error(NotImplementedError, "Fusion is not implemented for "
+                                         "distributed circuits.")
 
     def with_noise(self, noise_map, measurement_noise=None):
         raise_error(NotImplementedError, "Distributed circuit does not support "
diff --git a/src/qibo/core/fusion.py b/src/qibo/core/fusion.py
deleted file mode 100644
index bd9d774aed..0000000000
--- a/src/qibo/core/fusion.py
+++ /dev/null
@@ -1,356 +0,0 @@
-import functools
-import operator
-from qibo import K, gates
-from qibo.config import raise_error
-from qibo.abstractions.abstract_gates import ParametrizedGate
-from typing import List, Optional, Set, Tuple
-
-
-class FusionGroup:
-    """Group of one-qubit and two-qubit gates that act in two specific gates.
-
-    These gates can be fused into a single two-qubit gate represented by a
-    general 4x4 matrix.
-
-    Args:
-        qubit0 (int): Id of the first qubit that the ``FusionGroup`` act.
-        qubit1 (int): Id of the first qubit that the ``FusionGroup`` act.
-        gates0 (list): List of lists of one-qubit gates to be applied to ``qubit0``.
-            One qubit gates are split in groups according to when the two qubit gates are applied (see example).
-            Has ``len(gates0) = len(two_qubit_gates) + 1``.
-        gates1 (list): Same as ``gates0`` but for ``qubit1``.
-        two_qubit_gates (list): List of two qubit gates acting on ``qubit0`` and ``qubit1``.
-
-    Example:
-        ::
-
-            queue = [gates.H(0), gates.H(1), gates.CNOT(0, 1), gates.X(0), gates.X(1)]
-            group = fusion.FusionGroup.from_queue(queue)
-            # results to the following values for the attributes:
-            group.qubit0 = 0
-            group.qubit1 = 1
-            group.gates0 = [[gates.H(0)], [gates.X(0)]]
-            group.gates1 = [[gates.H(1)], [gates.X(1)]]
-            group.two_qubit_gates = [gates.CNOT(0, 1)]
-    """
-
-    # ``FusionGroup`` cannot start with these gates because it is more
-    # efficient to apply them on their own
-    _efficient_gates = {"CNOT", "CZ", "SWAP", "CU1"}
-
-    def __init__(self):
-        self.qubit0 = None
-        self.qubit1 = None
-        self.gates0 = [[]]
-        self.gates1 = [[]]
-        self.two_qubit_gates = []
-
-        self.completed = False
-        self.special_gate = None
-        self._fused_gates = None
-        if K.op is not None:
-            self.K = K.np
-        else:
-            self.K = K
-
-    @property
-    def qubits(self) -> Set[int]:
-        """Set of ids of the two qubits that the ``FusionGroup`` acts on."""
-        if self.qubit0 is None:
-            return {}
-        if self.qubit1 is None:
-            return {self.qubit0}
-        return {self.qubit0, self.qubit1}
-
-    @property
-    def gates(self) -> Tuple["Gate"]:
-        """Tuple with fused gates.
-
-        These gates have equivalent action with all the original gates that
-        were added in the ``FusionGroup``.
-        """
-        if self._fused_gates is None:
-            self._fused_gates = self.calculate()
-        return self._fused_gates
-
-    def update(self) -> Tuple["Gate"]:
-        """Recalculates fused gates.
-
-        This is used automatically by circuit objects in order to repeat the
-        calculation of fused gates after the parameters of original gates have
-        been changed using ``circuit.set_parameters``.
-        It assumes that the parameters of the gate objects contained in the
-        current ``FusionGroup`` have already been updated.
-        """
-        updated_gates = self.calculate()
-        for gate, new_gate in zip(self.gates, updated_gates):
-            if isinstance(gate, ParametrizedGate):
-                gate.parameters = new_gate.parameters
-        return self._fused_gates
-
-    def first_gate(self, i: int) -> Optional["Gate"]:
-        """First one-qubit gate of the group."""
-        if i < 0 or i > 1:
-            raise_error(ValueError, f"Invalid integer {i} given in FusionGroup.first_gate.")
-        gates = self.gates0 if i == 0 else self.gates1
-        for group in gates:
-            if group:
-                return group[0]
-        return None
-
-    def is_efficient(self, gate: "Gate") -> bool:
-        """Checks if given two-qubit ``gate`` is efficient.
-
-        Efficient gates are not fused if they are in the start or in the end.
-        """
-        return gate.__class__.__name__ in self._efficient_gates
-
-    @classmethod
-    def from_queue(cls, queue: List["Gate"]) -> List["FusionGroup"]:
-        """Fuses a queue of gates by combining up to two-qubit gates.
-
-        Args:
-            queue (list): List of gates.
-
-        Returns:
-            List of ``FusionGroup`` objects that correspond to the fused gates.
-        """
-        group_queue = []
-        remaining_queue = list(queue)
-        while remaining_queue:
-            gates = iter(remaining_queue)
-            gate = next(gates)
-            new_group = cls()
-            new_group.add(gate)
-            remaining_queue = []
-            for gate in gates:
-                if new_group.completed:
-                    remaining_queue.append(gate)
-                    break
-
-                if new_group.can_add(gate):
-                    commutes = True
-                    for blocking_gate in remaining_queue:
-                        commutes = commutes and gate.commutes(blocking_gate)
-                        if not commutes:
-                            break
-
-                    if commutes:
-                        new_group.add(gate)
-                    else:
-                        remaining_queue.append(gate)
-                else:
-                    remaining_queue.append(gate)
-
-            new_group.completed = True
-            remaining_queue.extend(gates)
-            group_queue.append(new_group)
-
-        return group_queue
-
-    def can_add(self, gate: "Gate") -> bool:
-        """Checks if ``gate`` can be added in the ``FusionGroup``."""
-        if self.completed:
-            return False
-        qubits = self.qubits
-        if not qubits:
-            return True
-
-        if len(gate.qubits) == 1:
-            return (len(qubits) == 1 or gate.qubits[0] in qubits)
-        if len(gate.qubits) == 2:
-            targets = set(gate.qubits)
-            if targets == qubits:
-                return True
-            if (self.qubit1 is None and self.qubit0 in targets and
-                not self.is_efficient(gate)):
-                return True
-        return False
-
-    def add(self, gate: "Gate"):
-        """Adds a gate in the group.
-
-        Raises:
-            ValueError: If the gate cannot be added in the group (eg. because
-                it acts on different qubits).
-            RuntimeError: If the group is completed.
-        """
-        if self.completed:
-            raise_error(RuntimeError, "Cannot add gates to completed FusionGroup.")
-
-        if not gate.qubits:
-            self._add_special_gate(gate)
-        elif len(gate.qubits) == 1:
-            self._add_one_qubit_gate(gate)
-        elif len(gate.qubits) == 2:
-            self._add_two_qubit_gate(gate)
-        else:
-            raise_error(ValueError, "Cannot add gate acting on {} qubits in fusion "
-                                    "group.".format(len(gate.qubits)))
-
-    def _add_special_gate(self, gate: "Gate"):
-        """Adds ``CallbackGate`` or ``Flatten`` on ``FusionGroup``."""
-        if self.qubits or self.special_gate is not None:
-            raise_error(ValueError, "Cannot add special gate on fusion group.")
-        self.special_gate = gate
-        self.completed = True
-
-    def _add_one_qubit_gate(self, gate: "Gate"):
-        """Adds one-qubit gate to ``FusionGroup``."""
-        qubit = gate.qubits[0]
-        if self.qubit0 is None or self.qubit0 == qubit:
-            self.qubit0 = qubit
-            self.gates0[-1].append(gate)
-        elif self.qubit1 is None or self.qubit1 == qubit:
-            self.qubit1 = qubit
-            self.gates1[-1].append(gate)
-        else:
-            raise_error(ValueError, "Cannot add gate on qubit {} in fusion group of "
-                                    "qubits {} and {}."
-                                    "".format(qubit, self.qubit0, self.qubit1))
-
-    def _add_two_qubit_gate(self, gate: "Gate"):
-        """Adds two-qubit gate to ``FusionGroup``."""
-        qubit0, qubit1 = gate.qubits
-        if self.qubit0 is None:
-            self.qubit0, self.qubit1 = qubit0, qubit1
-            self.two_qubit_gates.append(gate)
-            if self.is_efficient(gate):
-                self.completed = True
-
-        # case not used by the current scheme
-        elif self.qubit1 is None: # pragma: no cover
-            raise_error(NotImplementedError)
-
-            if self.is_efficient(gate):
-                raise_error(ValueError, "It is not efficient to add {} in FusionGroup "
-                                        "with only one qubit set.".format(gate))
-
-            if self.qubit0 == qubit0:
-                self.qubit1 = qubit1
-                self.two_qubit_gates.append(gate)
-            elif self.qubit0 == qubit1:
-                self.qubit1 = qubit0
-                self.two_qubit_gates.append(gate)
-            else:
-                raise_error(ValueError, "Cannot add gate on qubits {} and {} in "
-                                        "fusion group of qubit {}."
-                                        "".format(qubit0, qubit1, self.qubit0))
-
-        else:
-            if self.qubits != {qubit0, qubit1}:
-                raise_error(ValueError, "Cannot add gate on qubits {} and {} in "
-                                        "fusion group of qubits {} and {}."
-                                        "".format(qubit0, qubit1, self.qubit0, self.qubit1))
-            self.two_qubit_gates.append(gate)
-
-        self.gates0.append([])
-        self.gates1.append([])
-
-    def _one_qubit_matrix(self, gate0: "Gate", gate1: "Gate"):
-        """Calculates Kroneker product of two one-qubit gates.
-
-        Args:
-            gate0: Gate that acts on ``self.qubit0``.
-            gate1: Gate that acts on ``self.qubit1``.
-
-        Returns:
-            4x4 matrix that corresponds to the Kronecker product of the 2x2
-            gate matrices.
-        """
-        if K.op is not None:
-            return self.K.kron(gate0.matrix, gate1.matrix)
-        else:
-            matrix = self.K.tensordot(gate0.matrix, gate1.matrix, axes=0)
-            matrix = self.K.transpose(matrix, [0, 2, 1, 3])
-            return self.K.reshape(matrix, (4, 4))
-
-    def _two_qubit_matrix(self, gate: "Gate"):
-        """Calculates the 4x4 unitary matrix of a two-qubit gate.
-
-        Args:
-            gate: Two-qubit gate acting on ``(self.qubit0, self.qubit1)``.
-
-        Returns:
-            4x4 unitary matrix corresponding to the gate.
-        """
-        matrix = gate.matrix
-        if gate.qubits == (self.qubit1, self.qubit0):
-            matrix = self.K.reshape(matrix, 4 * (2,))
-            matrix = self.K.transpose(matrix, [1, 0, 3, 2])
-            matrix = self.K.reshape(matrix, (4, 4))
-        else:
-            assert gate.qubits == (self.qubit0, self.qubit1)
-        return matrix
-
-    def calculate(self):
-        """Calculates fused gate."""
-        if not self.completed:
-            raise_error(RuntimeError, "Cannot calculate fused gates for incomplete "
-                                      "FusionGroup.")
-        # Case 1: Special gate
-        if self.special_gate is not None:
-            assert not self.gates0[0] and not self.gates1[0]
-            assert not self.two_qubit_gates
-            return (self.special_gate,)
-
-        # Case 2: Two-qubit gates only (no one-qubit gates)
-        if self.first_gate(0) is None and self.first_gate(1) is None:
-            assert self.two_qubit_gates
-            # Case 2a: One two-qubit gate only
-            if len(self.two_qubit_gates) == 1:
-                return (self.two_qubit_gates[0],)
-
-            # Case 2b: Two or more two-qubit gates
-            fused_matrix = self._two_qubit_matrix(self.two_qubit_gates[0])
-            for gate in self.two_qubit_gates[1:]:
-                matrix = self._two_qubit_matrix(gate)
-                fused_matrix = self.K.matmul(matrix, fused_matrix)
-            return (gates.Unitary(fused_matrix, self.qubit0, self.qubit1),)
-
-        # Case 3: One-qubit gates exist
-        if not self.gates0[-1] and not self.gates1[-1]:
-            self.gates0.pop()
-            self.gates1.pop()
-
-        # Fuse one-qubit gates
-        ident0 = gates.I(self.qubit0)
-        gates0 = (functools.reduce(operator.matmul, reversed(gates), ident0)
-                  if len(gates) != 1 else gates[0] for gates in self.gates0)
-        if self.qubit1 is not None:
-            ident1 = gates.I(self.qubit1)
-            gates1 = (functools.reduce(operator.matmul, reversed(gates), ident1)
-                      if len(gates) != 1 else gates[0] for gates in self.gates1)
-
-        # Case 3a: One-qubit gates only (no two-qubit gates)
-        if not self.two_qubit_gates:
-            gates0 = list(gates0)[::-1]
-            fused_gate0 = (functools.reduce(operator.matmul, gates0, ident0)
-                           if len(gates0) != 1 else gates0[0])
-            if self.qubit1 is None:
-                return (fused_gate0,)
-
-            gates1 = list(gates1)[::-1]
-            fused_gate1 = (functools.reduce(operator.matmul, gates1, ident1)
-                           if len(gates1) != 1 else gates1[0])
-            return (fused_gate0, fused_gate1)
-
-        # Case 3b: One-qubit and two-qubit gates exist
-        fused_matrix = self._one_qubit_matrix(next(gates0), next(gates1))
-        for g0, g1, g2 in zip(gates0, gates1, self.two_qubit_gates):
-            matrix = self._one_qubit_matrix(g0, g1)
-            matrix2 = self._two_qubit_matrix(g2)
-            fused_matrix = self.K.matmul(self.K.matmul(matrix, matrix2),
-                                         fused_matrix)
-
-        if len(self.two_qubit_gates) == len(self.gates0):
-            g2 = self.two_qubit_gates[-1]
-            if self.is_efficient(g2):
-                fused_gate = gates.Unitary(fused_matrix, self.qubit0, self.qubit1)
-                return (fused_gate, g2)
-
-            matrix2 = self._two_qubit_matrix(g2)
-            fused_matrix = self.K.matmul(matrix2, fused_matrix)
-
-        fused_gate = gates.Unitary(fused_matrix, self.qubit0, self.qubit1)
-        return (fused_gate,)
diff --git a/src/qibo/core/gates.py b/src/qibo/core/gates.py
index 2efe419c03..84d75b3dc3 100644
--- a/src/qibo/core/gates.py
+++ b/src/qibo/core/gates.py
@@ -44,6 +44,13 @@ def _control_unitary(unitary):
         part2 = K.concatenate([zeros, unitary], axis=0)
         return K.concatenate([part1, part2], axis=1)
 
+    def _reset_unitary(self):
+        super()._reset_unitary()
+        self._native_op_matrix = None
+        self._custom_op_matrix = None
+        for gate in self.device_gates:
+            gate._reset_unitary()
+
     @property
     def cache(self):
         if self._cache is None:
@@ -1072,3 +1079,53 @@ def _density_matrix_call(self, state):
             self._target_qubits = original_targets
             return K.reshape(state, shape)
         return K._state_vector_call(self, state)
+
+
+class FusedGate(MatrixGate, abstract_gates.FusedGate):
+
+    def __init__(self, *q):
+        BackendGate.__init__(self)
+        abstract_gates.FusedGate.__init__(self, *q)
+        if self.gate_op:
+            # Custom kernels currently support up to two target qubits
+            if len(self.target_qubits) == 1:
+                self.gate_op = K.op.apply_gate
+            elif len(self.target_qubits) == 2:
+                self.gate_op = K.op.apply_two_qubit_gate
+            else:
+                raise_error(NotImplementedError, "Fused gates can target up to two qubits.")
+
+    def _construct_unitary(self):
+        """Constructs a single unitary by multiplying the matrices of the gates that are fused.
+
+        This matrix is used to perform a single update in the state during
+        simulation instead of applying the fused gates one by one.
+
+        Note that this method assumes maximum two target qubits and should be
+        updated if the fusion algorithm is extended to gates of higher rank.
+        """
+        matrix = K.qnp.eye(2 ** len(self.target_qubits))
+        for gate in self.gates:
+            # transfer gate matrix to numpy as it is more efficient for
+            # small tensor calculations
+            gmatrix = K.to_numpy(gate.matrix)
+            if len(gate.qubits) < len(self.target_qubits):
+                # fuse one-qubit gate (2x2 matrix) to a two-qubit ``FusedGate``
+                # Kronecker product with identity is needed to make the
+                # original matrix 4x4.
+                if gate.qubits[0] == self.target_qubits[0]:
+                    # gate target qubit is the first ```FusedGate`` target
+                    gmatrix = K.qnp.kron(gmatrix, K.qnp.eye(2))
+                else:
+                    # gate target qubit is the second ``FusedGate`` target
+                    gmatrix = K.qnp.kron(K.qnp.eye(2), gmatrix)
+            elif gate.qubits != self.target_qubits:
+                # fuse two-qubit gate (4x4 matrix) for which the target qubits
+                # are in opposite order compared to the ``FusedGate`` and
+                # the corresponding matrix has to be transposed before fusion
+                gmatrix = K.qnp.reshape(gmatrix, 4 * (2,))
+                gmatrix = K.qnp.transpose(gmatrix, [1, 0, 3, 2])
+                gmatrix = K.qnp.reshape(gmatrix, (4, 4))
+            # fuse the individual gate matrix to the total ``FusedGate`` matrix
+            matrix = gmatrix @ matrix
+        return K.cast(matrix)
diff --git a/src/qibo/tests/test_abstract_circuit.py b/src/qibo/tests/test_abstract_circuit.py
index 9aeb905f41..dfa3ba1ef7 100644
--- a/src/qibo/tests/test_abstract_circuit.py
+++ b/src/qibo/tests/test_abstract_circuit.py
@@ -470,9 +470,6 @@ def test_circuit_set_parameters_errors():
         c.set_parameters({0.3568})
     with pytest.raises(ValueError):
         c.queue[2].parameters = [0.1234, 0.4321, 0.156]
-    c.fusion_groups = ["test"]
-    with pytest.raises(TypeError):
-        c.set_parameters({gates.RX(0, theta=1.0): 0.568})
 
 
 def test_circuit_draw():
diff --git a/src/qibo/tests/test_abstract_gates.py b/src/qibo/tests/test_abstract_gates.py
index ecbd431587..68cdab90bc 100644
--- a/src/qibo/tests/test_abstract_gates.py
+++ b/src/qibo/tests/test_abstract_gates.py
@@ -375,3 +375,11 @@ def test_special_gate():
     assert not gate.commutes(gates.H(0))
     with pytest.raises(NotImplementedError):
         gate._on_qubits(1)
+
+
+def test_fused_gate():
+    gate = gates.FusedGate(0, 1)
+    gate.add(gates.H(0))
+    gate.add(gates.CNOT(0, 1))
+    with pytest.raises(ValueError):
+        gate.add(gates.CZ(1, 2))
diff --git a/src/qibo/tests/test_core_circuit.py b/src/qibo/tests/test_core_circuit.py
index 7fb51d5339..6a12591163 100644
--- a/src/qibo/tests/test_core_circuit.py
+++ b/src/qibo/tests/test_core_circuit.py
@@ -31,8 +31,7 @@ def test_circuit_add_layer(backend, nqubits, accelerators):
     for gate in c.queue:
         assert isinstance(gate, gates.Unitary)
 
-# TODO: Test `_fuse_copy`
-# TODO: Test `fuse`
+# :meth:`qibo.core.circuit.Circuit.fuse` is tested in `test_core_fusion.py`
 
 def test_eager_execute(backend, accelerators):
     c = Circuit(4, accelerators)
diff --git a/src/qibo/tests/test_core_circuit_features.py b/src/qibo/tests/test_core_circuit_features.py
index c6dc85ab73..30251523f9 100644
--- a/src/qibo/tests/test_core_circuit_features.py
+++ b/src/qibo/tests/test_core_circuit_features.py
@@ -78,7 +78,11 @@ def test_inverse_circuit_execution(backend, accelerators, fuse):
     c.add(gates.fSim(0, 2, theta=0.1, phi=0.3))
     c.add(gates.CU2(0, 1, phi=0.1, lam=0.1))
     if fuse:
-        c = c.fuse()
+        if accelerators:
+            with pytest.raises(NotImplementedError):
+                c = c.fuse()
+        else:
+            c = c.fuse()
     invc = c.invert()
     target_state = np.ones(2 ** 4) / 4
     final_state = invc(c(np.copy(target_state)))
diff --git a/src/qibo/tests/test_core_circuit_parametrized.py b/src/qibo/tests/test_core_circuit_parametrized.py
index da1e77b64c..e026fb9814 100644
--- a/src/qibo/tests/test_core_circuit_parametrized.py
+++ b/src/qibo/tests/test_core_circuit_parametrized.py
@@ -187,10 +187,10 @@ def test_set_parameters_with_double_variationallayer(backend, nqubits, trainable
 
 
 @pytest.mark.parametrize("trainable", [True, False])
-def test_set_parameters_with_gate_fusion(backend, trainable, accelerators):
+def test_set_parameters_with_gate_fusion(backend, trainable):
     """Check updating parameters of fused circuit."""
     params = np.random.random(9)
-    c = Circuit(5, accelerators)
+    c = Circuit(5)
     c.add(gates.RX(0, theta=params[0], trainable=trainable))
     c.add(gates.RY(1, theta=params[1]))
     c.add(gates.CZ(0, 1))
@@ -202,7 +202,9 @@ def test_set_parameters_with_gate_fusion(backend, trainable, accelerators):
     c.add(gates.RZ(1, theta=params[8]))
 
     fused_c = c.fuse()
-    K.assert_allclose(c(), fused_c())
+    final_state = fused_c()
+    target_state = c()
+    K.assert_allclose(final_state, target_state)
 
     if trainable:
         new_params = np.random.random(9)
diff --git a/src/qibo/tests/test_core_distcircuit.py b/src/qibo/tests/test_core_distcircuit.py
index e32c334eac..e235f1929a 100644
--- a/src/qibo/tests/test_core_distcircuit.py
+++ b/src/qibo/tests/test_core_distcircuit.py
@@ -56,9 +56,8 @@ def test_distributed_circuit_various_errors(backend):
 def test_distributed_circuit_fusion(backend, accelerators):
     c = DistributedCircuit(4, accelerators)
     c.add((gates.H(i) for i in range(4)))
-    final_state = c()
-    with pytest.raises(RuntimeError):
-        fused_c = c.fuse()
+    with pytest.raises(NotImplementedError):
+        c.fuse()
 
 
 def test_distributed_circuit_set_gates(backend):
diff --git a/src/qibo/tests/test_core_fusion.py b/src/qibo/tests/test_core_fusion.py
index 35b779594a..17e703af1e 100644
--- a/src/qibo/tests/test_core_fusion.py
+++ b/src/qibo/tests/test_core_fusion.py
@@ -2,97 +2,58 @@
 import numpy as np
 import pytest
 from qibo import gates, K
-from qibo.core import fusion
 from qibo.models import Circuit
 
 
-def test_fusion_group_from_queue_single_group():
-    """Check fusion that creates a single ``FusionGroup``."""
+def test_single_fusion_gate():
+    """Check circuit fusion that creates a single ``FusedGate``."""
     queue = [gates.H(0), gates.X(1), gates.CZ(0, 1)]
-    fused_groups = fusion.FusionGroup.from_queue(queue)
-    assert len(fused_groups) == 1
-    group = fused_groups[0]
-    assert group.gates0 == [[queue[0]], []]
-    assert group.gates1 == [[queue[1]], []]
-    assert group.two_qubit_gates == [queue[2]]
+    c = Circuit(2)
+    c.add(queue)
+    c = c.fuse()
+    assert len(c.queue) == 1
+    gate = c.queue[0]
+    for gate, target in zip(gate.gates, queue):
+        assert gate == target
 
 
-def test_fusion_group_from_queue_two_groups():
-    """Check fusion that creates two ``FusionGroup``s."""
+def test_two_fusion_gate():
+    """Check fusion that creates two ``FusedGate``s."""
     queue = [gates.X(0), gates.H(1),
              gates.RX(2, theta=0.1234).controlled_by(1),
              gates.H(2), gates.Y(1),
              gates.H(0)]
-    fused_groups = fusion.FusionGroup.from_queue(queue)
-    assert len(fused_groups) == 2
-    group1, group2 = fused_groups
-    assert group1.gates0 == [[queue[0], queue[5]]]
-    assert group1.gates1 == [[queue[1]]]
-    assert group1.two_qubit_gates == []
-    assert group2.gates0 == [[], [queue[4]]]
-    assert group2.gates1 == [[], [queue[3]]]
-    assert group2.two_qubit_gates == [queue[2]]
-
-
-def test_fusion_group_first_gate():
-    group = fusion.FusionGroup()
-    with pytest.raises(ValueError):
-        group.first_gate(2)
-
-
-def test_fusion_group_add():
-    group = fusion.FusionGroup()
-    group.add(gates.Flatten(np.ones(4)))
-    assert not group.can_add(gates.H(0))
-    with pytest.raises(RuntimeError):
-        group.add(gates.H(0))
-
-    group = fusion.FusionGroup()
-    group.can_add(gates.H(0))
-    group.add(gates.H(0))
-    group.can_add(gates.RX(0, theta=0.1234).controlled_by(1))
-    with pytest.raises(ValueError):
-        group.add(gates.Flatten(np.ones(4)))
-
-    group = fusion.FusionGroup()
-    group.add(gates.RX(0, theta=0.1234).controlled_by(1))
-    with pytest.raises(ValueError):
-        group.add(gates.H(2))
-
-    group = fusion.FusionGroup()
-    group.add(gates.RX(0, theta=0.1234).controlled_by(1))
-    with pytest.raises(ValueError):
-        group.add(gates.CZ(1, 2))
-
-    group = fusion.FusionGroup()
-    with pytest.raises(ValueError):
-        group.add(gates.TOFFOLI(0, 1, 2))
-
-
-def test_fusion_group_calculate(backend):
+    c = Circuit(3)
+    c.add(queue)
+    c = c.fuse()
+    assert len(c.queue) == 2
+    gate1, gate2 = c.queue
+    if len(gate1.gates) > len(gate2.gates): # pragma: no cover
+        # disabling coverage as this may not always happen
+        gate1, gate2 = gate2, gate1
+    assert gate1.gates == [queue[0], queue[-1]]
+    assert gate2.gates == queue[1:-1]
+
+
+def test_fusedgate_matrix_calculation(backend):
     queue = [gates.H(0), gates.H(1), gates.CNOT(0, 1),
              gates.X(0), gates.X(1)]
-    group = fusion.FusionGroup.from_queue(queue)
-    assert len(group) == 1
-    group = group[0]
-
-    assert len(group.gates) == 1
-    gate = group.gates[0]
+    circuit = Circuit(2)
+    circuit.add(queue)
+    circuit = circuit.fuse()
+    assert len(circuit.queue) == 1
+    fused_gate = circuit.queue[0]
 
     x = np.array([[0, 1], [1, 0]])
     h = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
     cnot = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1],
                      [0, 0, 1, 0]])
     target_matrix = np.kron(x, x) @ cnot @ np.kron(h, h)
-    K.assert_allclose(gate.matrix, target_matrix)
-
-    group = fusion.FusionGroup()
-    with pytest.raises(RuntimeError):
-        group.calculate()
+    K.assert_allclose(fused_gate.matrix, target_matrix)
 
 
-def test_fuse_circuit_two_qubit_only(backend):
-    """Check gate fusion in circuit with two-qubit gates only."""
+def test_fuse_circuit_two_qubit_gates(backend):
+    """Check circuit fusion in circuit with two-qubit gates only."""
     c = Circuit(2)
     c.add(gates.CNOT(0, 1))
     c.add(gates.RX(0, theta=0.1234).controlled_by(1))
@@ -103,6 +64,21 @@ def test_fuse_circuit_two_qubit_only(backend):
     K.assert_allclose(fused_c(), c())
 
 
+def test_fuse_circuit_three_qubit_gate(backend):
+    """Check circuit fusion in circuit with three-qubit gate."""
+    c = Circuit(4)
+    c.add((gates.H(i) for i in range(4)))
+    c.add(gates.CZ(0, 1))
+    c.add(gates.CZ(2, 3))
+    c.add(gates.TOFFOLI(0, 1, 2))
+    c.add(gates.SWAP(1, 2))
+    c.add((gates.H(i) for i in range(4)))
+    c.add(gates.CNOT(0, 1))
+    c.add(gates.CZ(2, 3))
+    fused_c = c.fuse()
+    K.assert_allclose(fused_c(), c(), atol=1e-12)
+
+
 @pytest.mark.parametrize("nqubits", [4, 5, 10, 11])
 @pytest.mark.parametrize("nlayers", [1, 2])
 def test_variational_layer_fusion(backend, nqubits, nlayers):
@@ -115,7 +91,7 @@ def test_variational_layer_fusion(backend, nqubits, nlayers):
         c.add((gates.RY(i, next(theta_iter)) for i in range(nqubits)))
         c.add((gates.CZ(i, i + 1) for i in range(0, nqubits - 1, 2)))
         c.add((gates.RY(i, next(theta_iter)) for i in range(nqubits)))
-        c.add((gates.CZ(i, i + 1) for i in range(1, nqubits - 2, 2)))
+        c.add((gates.CZ(i, i + 1) for i in range(1, nqubits - 1, 2)))
         c.add(gates.CZ(0, nqubits - 1))
 
     fused_c = c.fuse()
@@ -139,13 +115,12 @@ def test_random_circuit_fusion(backend, nqubits, ngates):
         while q0 == q1:
             q0, q1 = np.random.randint(0, nqubits, (2,))
         c.add(gate(q0, q1))
-
     fused_c = c.fuse()
     K.assert_allclose(fused_c(), c())
 
 
 def test_controlled_by_gates_fusion(backend):
-    """Check gate fusion in circuit that contains ``controlled_by`` gates."""
+    """Check circuit fusion that contains ``controlled_by`` gates."""
     c = Circuit(4)
     c.add((gates.H(i) for i in range(4)))
     c.add(gates.RX(1, theta=0.1234).controlled_by(0))
diff --git a/src/qibo/tests/test_core_gates.py b/src/qibo/tests/test_core_gates.py
index dec04ea2e6..52da58ac8f 100644
--- a/src/qibo/tests/test_core_gates.py
+++ b/src/qibo/tests/test_core_gates.py
@@ -1,4 +1,4 @@
-"""Test gates defined in `qibo/core/cgates.py` and `qibo/core/gates.py`."""
+"""Test gates defined in `qibo/core/gates.py`."""
 import pytest
 import numpy as np
 from qibo import gates, K
@@ -522,3 +522,22 @@ def test_thermal_relaxation_channel_errors(backend, t1, t2, time, excpop):
     with pytest.raises(ValueError):
         gate = gates.ThermalRelaxationChannel(
             0, t1, t2, time, excited_population=excpop)
+
+
+def test_fused_gate_init(backend):
+    gate = gates.FusedGate(0)
+    gate = gates.FusedGate(0, 1)
+    if K.op is not None:
+        with pytest.raises(NotImplementedError):
+            gate = gates.FusedGate(0, 1, 2)
+
+
+def test_fused_gate_construct_unitary(backend):
+    gate = gates.FusedGate(0, 1)
+    gate.add(gates.H(0))
+    gate.add(gates.H(1))
+    gate.add(gates.CZ(0, 1))
+    hmatrix = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
+    czmatrix = np.diag([1, 1, 1, -1])
+    target_matrix = czmatrix @ np.kron(hmatrix, hmatrix)
+    K.assert_allclose(gate.matrix, target_matrix)