diff --git a/setup.py b/setup.py
index e781945..13285eb 100644
--- a/setup.py
+++ b/setup.py
@@ -41,7 +41,7 @@ setup(
     install_requires=[
         "qibo>=0.1.10",
         "qibojit>=0.0.7",
-        "quimb[tensor]>=1.4.0",
+        "quimb[tensor]>=1.6.0",
     ],
     extras_require={
         "docs": [],
diff --git a/src/qibotn/quimb.py b/src/qibotn/quimb.py
index 8414540..a9ac510 100644
--- a/src/qibotn/quimb.py
+++ b/src/qibotn/quimb.py
@@ -3,16 +3,25 @@ import quimb.tensor as qtn
 from qibo.models import Circuit as QiboCircuit
 
 
-def from_qibo(circuit: QiboCircuit, psi0=None):
+def from_qibo(circuit: QiboCircuit, is_mps: False, psi0=None, method='svd',
+              cutoff=1e-6, cutoff_mode='abs'):
     nqubits = circuit.nqubits
-    tncirc = qtn.Circuit(nqubits, psi0=psi0)
+    gate_opt = {}
+    if is_mps:
+        tncirc = qtn.CircuitMPS(nqubits, psi0=psi0)
+        gate_opt["method"] = method
+        gate_opt["cutoff"] = cutoff
+        gate_opt["cutoff_mode"] = cutoff_mode
+    else:
+        tncirc = qtn.Circuit(nqubits, psi0=psi0)
 
     for gate in circuit.queue:
         tncirc.apply_gate(
             gate.name,
             *gate.parameters,
             *gate.qubits,
-            parametrize=len(gate.parameters) > 0
+            parametrize=False if is_mps else (len(gate.parameters) > 0),
+            **gate_opt
         )
 
     return tncirc
@@ -24,7 +33,7 @@ def init_state_tn(nqubits, init_state_sv):
     return qtn.tensor_1d.MatrixProductState.from_dense(init_state_sv, dims)
 
 
-def eval(qasm: str, init_state, backend="numpy"):
+def eval(qasm: str, init_state, is_mps, backend="numpy"):
     """Evaluate QASM with Quimb
 
     backend (quimb): numpy, cupy, jax. Passed to ``opt_einsum``.
@@ -32,7 +41,7 @@ def eval(qasm: str, init_state, backend="numpy"):
     """
     circuit = QiboCircuit.from_qasm(qasm)
     init_state_mps = init_state_tn(circuit.nqubits, init_state)
-    circ_quimb = from_qibo(circuit, psi0=init_state_mps)
+    circ_quimb = from_qibo(circuit, is_mps, psi0=init_state_mps)
     interim = circ_quimb.psi.full_simplify(seq="DRC")
     amplitudes = interim.to_dense(backend=backend).flatten()
 
diff --git a/tests/test_qasm_quimb_backend.py b/tests/test_qasm_quimb_backend.py
index e962018..f48d48a 100644
--- a/tests/test_qasm_quimb_backend.py
+++ b/tests/test_qasm_quimb_backend.py
@@ -10,7 +10,8 @@ from qibo.models import QFT
 
 
 def create_init_state(nqubits):
-    init_state = np.random.random(2**nqubits) + 1j * np.random.random(2**nqubits)
+    init_state = np.random.random(2**nqubits) + \
+        1j * np.random.random(2**nqubits)
     init_state = init_state / np.sqrt((np.abs(init_state) ** 2).sum())
     return init_state
 
@@ -29,8 +30,9 @@ def time(func):
     return time, res
 
 
-@pytest.mark.parametrize("nqubits", [1, 2, 5, 10])
-def test_eval(nqubits: int):
+@pytest.mark.parametrize("nqubits, tolerance, is_mps",
+                         [(1, 1e-6, True), (2, 1e-6, False), (5, 1e-3, True), (10, 1e-3, False)])
+def test_eval(nqubits: int, tolerance: float, is_mps: bool):
     # hack quimb to use the correct number of processes
     # TODO: remove completely, or at least delegate to the backend
     # implementation
@@ -41,7 +43,8 @@ def test_eval(nqubits: int):
     init_state_tn = copy.deepcopy(init_state)
 
     # Test qibo
-    qibo.set_backend(backend=config.qibo.backend, platform=config.qibo.platform)
+    qibo.set_backend(backend=config.qibo.backend,
+                     platform=config.qibo.platform)
     qibo_time, (qibo_circ, result_sv) = time(
         lambda: qibo_qft(nqubits, init_state, swaps=True)
     )
@@ -52,9 +55,10 @@ def test_eval(nqubits: int):
     # Test quimb
     quimb_time, result_tn = time(
         lambda: qibotn.quimb.eval(
-            qasm_circ, init_state_tn, backend=config.quimb.backend
+            qasm_circ, init_state_tn, is_mps, backend=config.quimb.backend
         )
     )
 
     assert 1e-2 * qibo_time < quimb_time < 1e2 * qibo_time
-    assert np.allclose(result_sv, result_tn), "Resulting dense vectors do not match"
+    assert np.allclose(result_sv, result_tn,
+                       atol=tolerance), "Resulting dense vectors do not match"