diff --git a/.github/workflows/rules.yml b/.github/workflows/rules.yml
index 177ed86..c0f2917 100644
--- a/.github/workflows/rules.yml
+++ b/.github/workflows/rules.yml
@@ -12,7 +12,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        python-version: [3.7, 3.8, 3.9, "3.10"]
+        python-version: [3.8, 3.9, "3.10"]
     uses: qiboteam/workflows/.github/workflows/rules.yml@main
     with:
       os: ${{ matrix.os }}
diff --git a/.gitignore b/.gitignore
index 68bc17f..d903da0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -146,6 +146,7 @@ dmypy.json
 # Pyre type checker
 .pyre/
 
+
 # pytype static type analyzer
 .pytype/
 
diff --git a/setup.py b/setup.py
index aefac89..2b1184c 100644
--- a/setup.py
+++ b/setup.py
@@ -42,6 +42,8 @@ setup(
         "qibo>=0.1.10",
         "qibojit>=0.0.7",
         "quimb[tensor]>=1.4.0",
+        "cupy>=11.6.0",
+        "cuquantum-python-cu11",
     ],
     extras_require={
         "docs": [],
@@ -54,7 +56,7 @@ setup(
             "pylint>=2.16.0",
         ],
     },
-    python_requires=">=3.7.0",
+    python_requires=">=3.8.0",
     long_description=(HERE / "README.md").read_text(encoding="utf-8"),
     long_description_content_type="text/markdown",
 )
diff --git a/src/qibotn/QiboCircuitConvertor.py b/src/qibotn/QiboCircuitConvertor.py
new file mode 100644
index 0000000..c30cfb6
--- /dev/null
+++ b/src/qibotn/QiboCircuitConvertor.py
@@ -0,0 +1,111 @@
+import cupy as cp
+import numpy as np
+
+
+class QiboCircuitToEinsum:
+    """Convert a circuit to a Tensor Network (TN) representation.
+    The circuit is first processed to an intermediate form by grouping each gate
+    matrix with its corresponding qubit it is acting on to a list. It is then
+    converted it to an equivalent TN expression through the class function
+    state_vector_operands() following the Einstein summation convention in the
+    interleave format.
+
+    See document for detail of the format: https://docs.nvidia.com/cuda/cuquantum/python/api/generated/cuquantum.contract.html
+
+    The output is to be used by cuQuantum's contract() for computation of the
+    state vectors of the circuit.
+    """
+
+    def __init__(self, circuit, dtype="complex128"):
+        self.backend = cp
+        self.dtype = getattr(self.backend, dtype)
+        self.init_basis_map(self.backend, dtype)
+        self.init_intermediate_circuit(circuit)
+
+    def state_vector_operands(self):
+        input_bitstring = "0" * len(self.active_qubits)
+
+        input_operands = self._get_bitstring_tensors(input_bitstring)
+
+        (
+            mode_labels,
+            qubits_frontier,
+            next_frontier,
+        ) = self._init_mode_labels_from_qubits(self.active_qubits)
+
+        gate_mode_labels, gate_operands = self._parse_gates_to_mode_labels_operands(
+            self.gate_tensors, qubits_frontier, next_frontier
+        )
+
+        operands = input_operands + gate_operands
+        mode_labels += gate_mode_labels
+
+        out_list = []
+        for key in qubits_frontier:
+            out_list.append(qubits_frontier[key])
+
+        operand_exp_interleave = [x for y in zip(
+            operands, mode_labels) for x in y]
+        operand_exp_interleave.append(out_list)
+        return operand_exp_interleave
+
+    def _init_mode_labels_from_qubits(self, qubits):
+        n = len(qubits)
+        frontier_dict = {q: i for i, q in enumerate(qubits)}
+        mode_labels = [[i] for i in range(n)]
+        return mode_labels, frontier_dict, n
+
+    def _get_bitstring_tensors(self, bitstring):
+        return [self.basis_map[ibit] for ibit in bitstring]
+
+    def _parse_gates_to_mode_labels_operands(
+        self, gates, qubits_frontier, next_frontier
+    ):
+        mode_labels = []
+        operands = []
+
+        for tensor, gate_qubits in gates:
+            operands.append(tensor)
+            input_mode_labels = []
+            output_mode_labels = []
+            for q in gate_qubits:
+                input_mode_labels.append(qubits_frontier[q])
+                output_mode_labels.append(next_frontier)
+                qubits_frontier[q] = next_frontier
+                next_frontier += 1
+            mode_labels.append(output_mode_labels + input_mode_labels)
+        return mode_labels, operands
+
+    def op_shape_from_qubits(self, nqubits):
+        """Modify tensor to cuQuantum shape
+        (qubit_states,input_output) * qubits_involved
+        """
+        return (2, 2) * nqubits
+
+    def init_intermediate_circuit(self, circuit):
+        self.gate_tensors = []
+        gates_qubits = []
+
+        for gate in circuit.queue:
+            gate_qubits = gate.control_qubits + gate.target_qubits
+            gates_qubits.extend(gate_qubits)
+
+            # self.gate_tensors is to extract into a list the gate matrix together with the qubit id that it is acting on
+            # https://github.com/NVIDIA/cuQuantum/blob/6b6339358f859ea930907b79854b90b2db71ab92/python/cuquantum/cutensornet/_internal/circuit_parser_utils_cirq.py#L32
+            required_shape = self.op_shape_from_qubits(len(gate_qubits))
+            self.gate_tensors.append(
+                (
+                    cp.asarray(gate.matrix).reshape(required_shape),
+                    gate_qubits,
+                )
+            )
+
+        # self.active_qubits is to identify qubits with at least 1 gate acting on it in the whole circuit.
+        self.active_qubits = np.unique(gates_qubits)
+
+    def init_basis_map(self, backend, dtype):
+        asarray = backend.asarray
+        state_0 = asarray([1, 0], dtype=dtype)
+        state_1 = asarray([0, 1], dtype=dtype)
+
+        self.basis_map = {"0": state_0, "1": state_1}
diff --git a/src/qibotn/__main__.py b/src/qibotn/__main__.py
index 8ed7439..0476be5 100644
--- a/src/qibotn/__main__.py
+++ b/src/qibotn/__main__.py
@@ -1,5 +1,6 @@
 import argparse
-from qibotn import qasm_quimb
+
+import qibotn.quimb
 
 
 def parser():
@@ -12,7 +13,7 @@ def parser():
 
 def main(args: argparse.Namespace):
     print("Testing for %d nqubits" % (args.nqubits))
-    qasm_quimb.eval_QI_qft(args.nqubits, args.qasm_circ, args.init_state)
+    qibotn.quimb.eval(args.nqubits, args.qasm_circ, args.init_state)
 
 
 if __name__ == "__main__":
diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
new file mode 100644
index 0000000..e6f3e8c
--- /dev/null
+++ b/src/qibotn/cutn.py
@@ -0,0 +1,8 @@
+# from qibotn import quimb as qiboquimb
+from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
+from cuquantum import contract
+
+
+def eval(qibo_circ, datatype):
+    myconvertor = QiboCircuitToEinsum(qibo_circ, dtype=datatype)
+    return contract(*myconvertor.state_vector_operands())
diff --git a/tests/test_cuquantum_cutensor_backend.py b/tests/test_cuquantum_cutensor_backend.py
new file mode 100644
index 0000000..e7f2804
--- /dev/null
+++ b/tests/test_cuquantum_cutensor_backend.py
@@ -0,0 +1,48 @@
+from timeit import default_timer as timer
+
+import config
+import numpy as np
+import pytest
+import qibo
+from qibo.models import QFT
+
+
+def qibo_qft(nqubits, swaps):
+    circ_qibo = QFT(nqubits, swaps)
+    state_vec = np.array(circ_qibo())
+    return circ_qibo, state_vec
+
+
+def time(func):
+    start = timer()
+    res = func()
+    end = timer()
+    time = end - start
+    return time, res
+
+
+@pytest.mark.gpu
+@pytest.mark.parametrize("nqubits", [1, 2, 5, 10])
+def test_eval(nqubits: int, dtype="complex128"):
+    """Evaluate QASM with cuQuantum.
+
+    Args:
+        nqubits (int): Total number of qubits in the system.
+        dtype (str): The data type for precision, 'complex64' for single,
+            'complex128' for double.
+    """
+    import qibotn.cutn
+
+    # Test qibo
+    qibo.set_backend(backend=config.qibo.backend,
+                     platform=config.qibo.platform)
+    qibo_time, (qibo_circ, result_sv) = time(
+        lambda: qibo_qft(nqubits, swaps=True))
+
+    # Test Cuquantum
+    cutn_time, result_tn = time(
+        lambda: qibotn.cutn.eval(qibo_circ, dtype).flatten())
+
+    assert 1e-2 * qibo_time < cutn_time < 1e2 * qibo_time
+    assert np.allclose(
+        result_sv, result_tn), "Resulting dense vectors do not match"