diff --git a/scratchpad/qtensor_MPS/constants.py b/scratchpad/qtensor_MPS/constants.py
index e3f54288..415ae0a5 100644
--- a/scratchpad/qtensor_MPS/constants.py
+++ b/scratchpad/qtensor_MPS/constants.py
@@ -1,8 +1,7 @@
 import numpy as np
 
-hmatrix = (
-    1 / np.sqrt(2) * np.array([[1.0, 1.0], [1.0, -1.0]], dtype=np.complex64)
-)
+pauli_matrices = ["X", "Z", "I", "Y"]
+hmatrix = 1 / np.sqrt(2) * np.array([[1.0, 1.0], [1.0, -1.0]], dtype=np.complex64)
 imatrix = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.complex64)
 xmatrix = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex64)
 ymatrix = np.array([[0.0, -1j], [1j, 0.0]], dtype=np.complex64)
@@ -28,10 +27,10 @@
 )
 swap_matrix = np.reshape(_swap_matrix, newshape=(2, 2, 2, 2))
 
-_hbar = 1.0545718*10e-34
-_sigma_z =  _hbar * 0.5 * zmatrix
-_sigma_x_pos =  _hbar * 0.5 * (xmatrix + 1j*ymatrix)
-_sigma_x_neg =  _hbar * 0.5 * (xmatrix - 1j*ymatrix)
+_hbar = 1.0545718 * 10e-34
+_sigma_z = _hbar * 0.5 * zmatrix
+_sigma_x_pos = _hbar * 0.5 * (xmatrix + 1j * ymatrix)
+_sigma_x_neg = _hbar * 0.5 * (xmatrix - 1j * ymatrix)
 
 _sigma_z_sigma_z_gate_matrix = np.kron(_sigma_z, _sigma_z)
 sigma_z_sigma_z_gate_matrix = _sigma_z_sigma_z_gate_matrix
@@ -40,10 +39,14 @@
 # TODO: What does it look if hamiltonian is e^(I) / e^(0)
 
 _sigma_x_pos_sigma_x_neg_gate_matrix = np.kron(_sigma_x_pos, _sigma_x_neg)
-sigma_x_pos_sigma_x_neg_gate_matrix = np.reshape(_sigma_x_pos_sigma_x_neg_gate_matrix, newshape=(2,2,2,2))
+sigma_x_pos_sigma_x_neg_gate_matrix = np.reshape(
+    _sigma_x_pos_sigma_x_neg_gate_matrix, newshape=(2, 2, 2, 2)
+)
 
 _sigma_x_neg_sigma_x_pos_gate_matrix = np.kron(_sigma_x_neg, _sigma_x_pos)
-sigma_x_neg_sigma_x_pos_gate_matrix = np.reshape(_sigma_x_neg_sigma_x_pos_gate_matrix, newshape=(2,2,2,2))
+sigma_x_neg_sigma_x_pos_gate_matrix = np.reshape(
+    _sigma_x_neg_sigma_x_pos_gate_matrix, newshape=(2, 2, 2, 2)
+)
 # _ising_hamiltonian_matrix = sigma_z_sigma_z_gate_matrix + 0.5 * sigma_x_pos_sigma_x_neg_gate_matrix + 0.5 * sigma_x_neg_sigma_x_pos_gate_matrix
 _ising_hamiltonian_matrix = sigma_z_sigma_z_gate_matrix
-ising_hamiltonian_matrix = np.reshape(_ising_hamiltonian_matrix, newshape=(2,2,2,2))
\ No newline at end of file
+ising_hamiltonian_matrix = np.reshape(_ising_hamiltonian_matrix, newshape=(2, 2, 2, 2))
diff --git a/scratchpad/qtensor_MPS/gates.py b/scratchpad/qtensor_MPS/gates.py
index d65d507a..0377904f 100644
--- a/scratchpad/qtensor_MPS/gates.py
+++ b/scratchpad/qtensor_MPS/gates.py
@@ -6,6 +6,11 @@
 from scipy.linalg import qr
 
 
+def get_gate(gc):
+    gates = {"X": xgate, "I": igate, "Y": ygate, "Z": zgate}
+    return gates.get(gc, lambda: None)().tensor
+
+
 def xgate() -> tn.Node:
     return tn.Node(deepcopy(xmatrix), name="xgate")
 
@@ -15,7 +20,11 @@ def igate() -> tn.Node:
 
 
 def zgate() -> tn.Node:
-    return tn.Node(deepcopy(zmatrix), name="xgate")
+    return tn.Node(deepcopy(zmatrix), name="zgate")
+
+
+def ygate() -> tn.Node:
+    return tn.Node(deepcopy(ymatrix), name="ygate")
 
 
 def cnot() -> tn.Node:
@@ -29,7 +38,7 @@ def hgate() -> tn.Node:
 def sigmaRx(t) -> tn.Node:
     _hbar = 1.0545718 * 10e-34
     # theta =  -1j * t * 0.5 * _hbar * 0.5
-    theta = 1j * 0.1
+    theta = 0.1 * t
 
     gate_matrix = np.array(
         [
@@ -45,13 +54,13 @@ def sigmaRzz(t) -> tn.Node:
     _hbar = 1.0545718 * 10e-34
     # theta =  -1j * t * 0.5 * _hbar * 0.5 * _hbar * 0.5
     # J = 0.1Hz
-    theta = -1j * 0.7
+    theta = -1j * 0.7 * t
 
     gate_matrix = np.array(
         [
             [np.exp(theta / 2), 0.0, 0.0, 0.0],
-            [0.0, np.exp(theta / 2), 0.0, 0.0],
-            [0.0, 0.0, np.exp(theta / 2), 0.0],
+            [0.0, np.exp(-1 * theta / 2), 0.0, 0.0],
+            [0.0, 0.0, np.exp(-1 * theta / 2), 0.0],
             [0.0, 0.0, 0.0, np.exp(theta / 2)],
         ]
     )
diff --git a/scratchpad/qtensor_MPS/mpo.py b/scratchpad/qtensor_MPS/mpo.py
index b70bda3c..57e5cbcf 100644
--- a/scratchpad/qtensor_MPS/mpo.py
+++ b/scratchpad/qtensor_MPS/mpo.py
@@ -1,7 +1,9 @@
 import tensornetwork as tn
 import numpy as np
-from gates import igate
+from gates import igate, get_gate, xgate
+from constants import pauli_matrices, xmatrix
 from typing import List, Union, Text, Optional, Any, Type
+from copy import deepcopy
 
 
 class MPOLayer:
@@ -62,62 +64,87 @@ def get_mpo_nodes(self, original) -> list[tn.Node]:
     def get_mpo_node(self, index, original) -> list[tn.Node]:
         return self.get_mpo_nodes(original)[index]
 
-    def construct_mpo(self, gate_function) -> "MPO":
-        # IIZ
-        gate_function = gate_function.reshape([self.physical_dim] * 2 * self.N)
-        to_split = tn.Node(
-            gate_function,
-            axis_names=["u" + str(i) for i in range(self.N)]
-            + ["d" + str(i) for i in range(self.N)],
-        )
+    def construct_mpo(self, pauli_string) -> "MPO":
+        # check if all elements are pauli matrices
+        paulis = list(pauli_string)
+
+        if any(pstring not in pauli_matrices for pstring in paulis):
+            print("Error: Not all are pauli string")
+            # Comments: change gatefunction - paulistring
+            pauli_string = pauli_string.reshape([self.physical_dim] * 2 * self.N)
+            to_split = tn.Node(
+                pauli_string,
+                axis_names=["u" + str(i) for i in range(self.N)]
+                + ["d" + str(i) for i in range(self.N)],
+            )
 
-        nodes = []
+            nodes = []
 
-        for i in range(self.N - 1):
-            left_edges = []
-            right_edges = []
+            for i in range(self.N - 1):
+                left_edges = []
+                right_edges = []
 
-            for edge in to_split.get_all_dangling():
-                if edge.name == "u" + str(i) or edge.name == "d" + str(i):
-                    left_edges.append(edge)
-                else:
-                    right_edges.append(edge)
-
-            if nodes:
-                for edge in nodes[-1].get_all_nondangling():
-                    if to_split in edge.get_nodes():
+                for edge in to_split.get_all_dangling():
+                    if edge.name == "u" + str(i) or edge.name == "d" + str(i):
                         left_edges.append(edge)
+                    else:
+                        right_edges.append(edge)
+
+                if nodes:
+                    for edge in nodes[-1].get_all_nondangling():
+                        if to_split in edge.get_nodes():
+                            left_edges.append(edge)
+
+                left, right, _ = tn.split_node(
+                    to_split,
+                    left_edges,
+                    right_edges,
+                    left_name=self.tensor_name + str(i),
+                    # max_singular_values=1,
+                )
+                nodes.append(left)
+                to_split = right
+            to_split.name = self.tensor_name + str(self.N - 1)
+            nodes.append(to_split)
+
+            self._nodes = nodes
 
-            left, right, _ = tn.split_node(
-                to_split,
-                left_edges,
-                right_edges,
-                left_name=self.tensor_name + str(i),
-                max_singular_values=1,
-            )
-            nodes.append(left)
-            to_split = right
-        to_split.name = self.tensor_name + str(self.N - 1)
-        nodes.append(to_split)
-
-        self._nodes = nodes
-
-    def add_single_qubit_gate(self, gate_matrix, idx):
-        if idx in range(1, self.N):
-            gate_matrix = gate_matrix.reshape(
-                self.physical_dim, 1, 1, self.physical_dim
-            )
         else:
-            gate_matrix = gate_matrix.reshape(self.physical_dim, 1, self.physical_dim)
+            if len(pauli_string) != self.N:
+                # pad pauli string
+                pauli_string += "I" * (self.N - len(pauli_string))
 
-        gate_node = tn.Node(
-            gate_matrix,
-            name=self.tensor_name + str(idx + 1),
-        )
+            nodes = self.get_mpo_nodes(False)
+            paulis = list(pauli_string)
 
-        self._nodes[idx] = gate_node
+            for i, (node, pauli) in enumerate(zip(nodes, paulis)):
+                pauli_gate = get_gate(pauli)
 
-    def two_qubit_svd(self, new_node, operating_qubits):
+                if i == 0 or i == self.N - 1:
+                    pauli_gate = pauli_gate.reshape(
+                        self.physical_dim, 1, self.physical_dim
+                    )
+                else:
+                    pauli_gate = pauli_gate.reshape(
+                        self.physical_dim, 1, 1, self.physical_dim
+                    )
+
+                node.set_tensor(pauli_gate)
+
+            self._nodes = nodes
+
+    def add_single_qubit_gate(self, gate, idx):
+        node = self._nodes[idx]
+        lst = list(node.get_all_dangling())
+        mps_index_edge = lst[0]
+        gate_edge = gate[1]
+        temp_node = tn.connect(mps_index_edge, gate_edge)
+        new_node = tn.contract(temp_node, name=self._nodes[idx].name)
+        self._nodes[idx] = new_node
+
+    def two_qubit_svd(
+        self, new_node, operating_qubits, left_gate_edge, right_gate_edge
+    ):
         left_connected_edge = None
         right_connected_edge = None
 
@@ -137,21 +164,13 @@ def two_qubit_svd(self, new_node, operating_qubits):
         left_edges = []
         right_edges = []
 
-        # for edge in (left_gate_edge, left_connected_edge):
-        #     if edge != None:
-        #         left_edges.append(edge)
-
-        # for edge in (right_gate_edge, right_connected_edge):
-        #     if edge != None:
-        #         right_edges.append(edge)
+        for edge in (*left_gate_edge, left_connected_edge):
+            if edge != None:
+                left_edges.append(edge)
 
-        left_edges.append(new_node.get_all_dangling()[0])
-        left_edges.append(new_node.get_all_dangling()[2])
-        left_edges.append(left_connected_edge)
-
-        right_edges.append(new_node.get_all_dangling()[1])
-        right_edges.append(new_node.get_all_dangling()[3])
-        right_edges.append(right_connected_edge)
+        for edge in (*right_gate_edge, right_connected_edge):
+            if edge != None:
+                right_edges.append(edge)
 
         u, s, vdag, _ = tn.split_node_full_svd(
             new_node, left_edges=left_edges, right_edges=right_edges
@@ -182,11 +201,7 @@ def add_two_qubit_gate(self, gate, operating_qubits):
             |   |
             c   d
 
-            0  1
-            |  |
-            gate'
-            |  |
-            2  3
+
         """
         # transpose
         gateT = tn.Node(np.conj(gate.tensor))
@@ -199,31 +214,38 @@ def add_two_qubit_gate(self, gate, operating_qubits):
         temp_nodesA = tn.connect(mpo_indexA, gate.get_edge(2))
         temp_nodesB = tn.connect(mpo_indexB, gate.get_edge(3))
 
-        left_gate_edge = gate.get_edge(0)
-        right_gate_edge = gate.get_edge(1)
+        left_gate_edge = [gate.get_edge(0)]
+        right_gate_edge = [gate.get_edge(1)]
+
+        left_gate_edge.append(mpo_indexC)
+        right_gate_edge.append(mpo_indexD)
         # left_gate_edgeT = gateT.get_edge(2)
         # right_gate_edgeT = gateT.get_edge(3)
 
         new_node = tn.contract_between(
             self._nodes[operating_qubits[0]], self._nodes[operating_qubits[1]]
         )
+
+        x = self.get_mpo_node(operating_qubits[0], True).get_all_dangling()[1]
+        y = self.get_mpo_node(operating_qubits[1], True).get_all_dangling()[1]
+
         node_gate_edge = tn.flatten_edges_between(new_node, gate)
         new_node = tn.contract(node_gate_edge)
 
-        self.two_qubit_svd(new_node, operating_qubits)
+        self.two_qubit_svd(new_node, operating_qubits, left_gate_edge, right_gate_edge)
 
         # Try connecting edges of gate 2, 3
         # Check transposition of gate, try multiindexing gate
-        temp_nodesC = tn.connect(mpo_indexC, gateT.get_edge(0))
-        temp_nodesD = tn.connect(mpo_indexD, gateT.get_edge(1))
+        # temp_nodesC = tn.connect(mpo_indexC, gateT.get_edge(0))
+        # temp_nodesD = tn.connect(mpo_indexD, gateT.get_edge(1))
 
-        new_node = tn.contract_between(
-            self._nodes[operating_qubits[0]], self._nodes[operating_qubits[1]]
-        )
-        node_gate_edge = tn.flatten_edges_between(new_node, gateT)
-        new_node = tn.contract(node_gate_edge)
+        # new_node = tn.contract_between(
+        #     self._nodes[operating_qubits[0]], self._nodes[operating_qubits[1]]
+        # )
+        # node_gate_edge = tn.flatten_edges_between(new_node, gateT)
+        # new_node = tn.contract(node_gate_edge)
 
-        self.two_qubit_svd(new_node, operating_qubits)
+        # self.two_qubit_svd(new_node, operating_qubits)
 
     def mpo_mps_inner_prod(self, mps):
         mpo = self.get_mpo_nodes(False)
diff --git a/scratchpad/qtensor_MPS/test.py b/scratchpad/qtensor_MPS/test.py
index 6abfbfd7..0043048e 100644
--- a/scratchpad/qtensor_MPS/test.py
+++ b/scratchpad/qtensor_MPS/test.py
@@ -178,72 +178,124 @@ def test_mpo_two_qubit():
     # assert (mps1.get_wavefunction() == mps2.get_wavefunction()).all()
 
 
-def test_mpo_construction_from_single_qubit_gate_function():
-    # IX
-    I = np.eye(2)
-    X = xmatrix
+def test_mpo_single_qubit_gate():
+    mps1 = MPS("q", 2, 2)
+    mps1.apply_single_qubit_gate(xgate(), 0)
+    print("MPS with gate", mps1.get_wavefunction())
 
+    mpo = MPOLayer("q", 2, 2)
+    mpo.add_single_qubit_gate(xgate(), 0)
     mps2 = MPS("q", 2, 2)
-    mps2.apply_single_qubit_gate(xgate(), 0)
-    print("MPS with gate", mps2.get_wavefunction())
+    mps2.apply_mpo_layer(mpo)
+    print("MPS with MPO", mps2.get_wavefunction())
+
+    condition = np.allclose(
+        np.array(mps1.get_wavefunction()),
+        np.array(mps2.get_wavefunction()),
+        rtol=1e-05,
+        atol=1e-08,
+    )
+    message = "Arrays are not equal within tolerance"
+
+    assert condition, message
+
+
+def test_mpo_two_qubit_gate():
+    mps1 = MPS("q", 2, 2)
+    mps1.apply_single_qubit_gate(xgate(), 1)
+    mps1.apply_two_qubit_gate(cnot(), [0, 1])
 
-    gate_func = np.kron(X, I)
     mpo = MPOLayer("q", 2, 2)
-    mpo.construct_mpo(gate_func)
-    mps = MPS("q", 2, 2)
+    mpo.add_single_qubit_gate(xgate(), 1)
+    mpo.add_two_qubit_gate(cnot(), [0, 1])
+    mps2 = MPS("q", 2, 2)
 
-    mps.apply_mpo_layer(mpo)
-    print("MPS with MPO", mps.get_wavefunction())
+    mps2.apply_mpo_layer(mpo)
 
     condition = np.allclose(
-        np.array(mps.get_wavefunction()),
+        np.array(mps1.get_wavefunction()),
         np.array(mps2.get_wavefunction()),
         rtol=1e-05,
         atol=1e-08,
     )
-    message = "Arrays are not equal within tolerance"
+    message = "MPO two qubit gate test failed"
 
     assert condition, message
 
 
-def test_mpo_construction_from_two_qubit_gate_function():
-    I = np.eye(2)
-    CNOT = np.array(
-        [
-            [1.0, 0.0, 0.0, 0.0],
-            [0.0, 1.0, 0.0, 0.0],
-            [0.0, 0.0, 0.0, 1.0],
-            [0.0, 0.0, 1.0, 0.0],
-        ]
-    )
+def test_mpo_construction_from_pauli_string():
+    pauli_string = "IXZXXXZZZZI"
+    n = len(pauli_string)
 
-    mps2 = MPS("q", 2, 2)
-    mps2.apply_single_qubit_gate(xgate(), 1)
-    print("INITIAL MPS", mps2.get_wavefunction())
+    mpo = MPOLayer("q", n, 2)
+    mpo.construct_mpo(pauli_string)
 
-    gate_func = CNOT
-    mpo = MPOLayer("q", 2, 2)
-    mpo.construct_mpo(gate_func)
-    mps = MPS("q", 2, 2)
-    mps.apply_single_qubit_gate(xgate(), 1)
-    mps.apply_mpo_layer(mpo)
-    print("FINAL MPS", mps.get_wavefunction())
+    mps1 = MPS("q", n, 2)
+    for i, ps in enumerate(pauli_string):
+        if ps == "X":
+            mps1.apply_single_qubit_gate(xgate(), i)
+        if ps == "Z":
+            mps1.apply_single_qubit_gate(zgate(), i)
+
+    mps2 = MPS("q", n, 2)
+    mps2.apply_mpo_layer(mpo)
 
     condition = np.allclose(
-        np.array(mps.get_wavefunction()),
+        np.array(mps1.get_wavefunction()),
         np.array(mps2.get_wavefunction()),
         rtol=1e-05,
         atol=1e-08,
     )
-    message = "Arrays are not equal within tolerance"
+    message = (
+        "(gates) on MPS and (MPO pauli string on) MPS are not equal within tolerance"
+    )
 
     assert condition, message
 
 
+# def test_mpo_construction_from_two_qubit_gate_function():
+#     I = np.eye(2)
+#     CNOT = np.array(
+#         [
+#             [1.0, 0.0, 0.0, 0.0],
+#             [0.0, 1.0, 0.0, 0.0],
+#             [0.0, 0.0, 0.0, 1.0],
+#             [0.0, 0.0, 1.0, 0.0],
+#         ]
+#     )
+
+#     mps2 = MPS("q", 2, 2)
+#     mps2.apply_single_qubit_gate(xgate(), 0)
+#     mps2.apply_two_qubit_gate(cnot(), [0, 1])
+#     print("INITIAL MPS", mps2.get_wavefunction())
+
+#     gate_func = CNOT
+#     mpo = MPOLayer("q", 2, 2)
+#     mpo.construct_mpo(gate_func)
+#     mps = MPS("q", 2, 2)
+#     mps.apply_single_qubit_gate(xgate(), 0)
+#     mps.apply_mpo_layer(mpo)
+#     print("FINAL MPS", mps.get_wavefunction())
+
+#     condition = np.allclose(
+#         np.array(mps.get_wavefunction()),
+#         np.array(mps2.get_wavefunction()),
+#         rtol=1e-05,
+#         atol=1e-08,
+#     )
+#     message = "Arrays are not equal within tolerance"
+
+#     # assert condition, message
+
+
 # TODO:
 # Verify for two qubit gates. - (SVD and then apply it) (apply and then SVD)
 # Evolve half of circuit in mps and half in mpo.
 # Evaluate Bell state, epectation value of X at a particular index
 # Evaluate X, Y, Z observable
 # Check for one/two/three qubits (cover edge cases) for both mps/mpo
-# test_mpo_construction_from_single_qubit_gate_function()
+# test_mpo_construction_from_pauli_string()
+# test_mpo_two_qubit_gate()
+# Initilisation given as pauli string and give mpo
+# Test for CNOT - Pauli then apply cnot + decompose cnot
+# expectation of an MPO - <psi| MPO | psi > psi = MPS (for this |psi | G1G2 X G2'G1' | psi>)